Abstract

 

The global positioning system (GPS) is the first and only full functional satellite navigation system that provides autonomous geo-spatial positioning with global coverage.

 

For computing the user position on the Earth, it is essential to determine the accurate GPS satellites orbit. One applicable way is to calculate the satellite orbit using the broadcasted ephemeris, the merit of which is that one could determine the orbit information at any given time. However this broadcasted ephemeris could only provide a limited accuracy of about 1.6 meters. Another approach is to obtain the satellite orbit information directly from the International GNSS Service (IGS) which is also called the precise orbits or precise ephemeris, the precision of which can reach the millimeter level. However the precise orbits are reported at an interval of 15 minutes thus requires to be interpolated into the time of positioning.

 

This thesis introduces these two approaches as how the GPS satellite orbits are determined or interpolated. It also analyzes the advantages and disadvantages of using the broadcasted and precise ephemeris. Finally, the performance of GPS positioning is presented and compared by using the satellite orbit information obtained in each method.

 


1. Introduction

 

The global positioning system (GPS) is a satellite based navigation system, developed by US Department of Defense since 1960s. It is initially in full operation in December 1993, and currently the only fully functional GNSS system in the world.

 

In calculating the correct user position, the user needs to know the orbit of the GPS satellites, and the satellite clock corrections. Generally, there are two approaches to obtain the satellite orbit and the clock correction.

 

The broadcast ephemerides are produced by the GPS Control Segment that uses the Kalman filter to estimate satellite position, velocity, solar radiation, pressure coefficient, clock bias, clock drift and clock drift rate [1]. The estimated parameters are used to propagate the satellite orbit and clock correction in the future. Then, the propagated values are converted into a set of keplerian elements and uploaded into the navigation message. Therefore, by decoding the navigation message in the GPS receiver, the user could obtain the ephemeris and calculate satellite orbit and clock correction at any given time. Besides, for post-processing applications, users could also obtain broadcast ephemeris from the Crustal Dynamics Data Information System (CDDIS), and they are stored in RINEX format. The broadcast ephemeredes are uploaded daily and are valid for two hours, the accuracies of which are 160 cm for the orbit and 7 ns for the satellite clock [2].

 

An alternative and more accurate approach to obtain the satellite orbit is through the IGS website managed by the Jet Propulsion Laboratory. The data is stored in the SP3 format. The accuracy of the satellite orbit and clock correction is much higher than the broadcast orbits, which are 5 cm and 0.1 ns respectively. The precise orbits are provided every 15 minutes, thus interpolation is required when calculating user position at a given time.

 

This paper will analyze the two approaches of obtaining the satellite orbit. For the broadcast orbits, detailed algorithm of calculating satellite position, velocity and clock correction from the ephemeris will be given based on GPS ICD [3] and reference 4. As for the precise orbits, an interpolation method is provided based on reference 5. The satellite position and clock correction obtained in two approaches will then be compared and eventually sent to the least square program to calculate their respective user position.

 

The RINEX and SP3 files used in this test are downloaded from the IGS website. The tracking site is called NIST which is located at Boulder, CO. The time of the files starts from 2008/11/07 00:00:00 UTC, with 30 seconds interval and the duration is 24 hours.

 


2. Broadcast Ephemeris

 

The Broadcast ephemeris is a set of modified Keplerian elements which are produced and managed by the GPS control segment. As introduced in the section 1, the Broadcast ephemeris can be either obtained by decoding the navigation message in the GPS receiver or directly downloaded in the IGS website. The former is usually for the real time positioning in the receiver, and the latter is for the differential GPS positioning or GPS post-processing.

 

Usually, one only needs to calculate satellite position to meet requirements of the least square GPS positioning. However, for some advanced application such as receiver clock drift estimation, obtaining the satellite velocity is also significant. GPS ICD [3] listed the equations of calculating satellites position from the Keplerian elements, and Benjamin W. Remondi [4] derived the satellites velocity from the derivative of certain Keplerian elements. Those Keplerian elements are listed in the tables below.

 

 

 

Parameters

Definition and Values

 

 

 

 

 

 

 

Satellite Orbit

Mean anomaly at a reference time

Eccentricity of the satellite orbit

Mean Motion Difference from Computed Value

Square Root of the Semi-Major Axis of the satellite orbit

Longitude of Ascending Node of Orbit Plane at Weekly Epoch

Inclination Angle at Reference Time

Argument of Perigee

Rate of Right Ascension

Rate of Inclination Angle

Amplitude of the Sine Harmonic Correction Term to the Argument of Latitude

Amplitude of the Cosine Harmonic Correction Term to the Argument of Latitude

Amplitude of the Sine Harmonic Correction Term to the Orbit Radius

Amplitude of the Cosine Harmonic Correction Term to the Orbit Radius

Amplitude of the Sine Harmonic Correction Term to the Angle of Inclination

Amplitude of the Cosine Harmonic Correction Term to the Angle of Inclination

Reference Time Ephemeris

Time of the week

 

Table 1. Definition of Keplerian Elements for Calculating Satellite Orbit

 

 

 

 

Parameters

Definition and Values

 

 

 

Clock Correct

Satellite group delay differential

Satellite clock correction

Clock Reference time

order parameter

order parameter

Speed of Light

 

 

Constants

WGS-84 value of the Earth's universal

WGS-84 value of the Earth's rotation rate

PI.

Table 2. Definition of Keplerian Elements for Calculating Clock Correction

 

The detailed algorithm of calculation and testing results are shown in the following subsections.

 

2.1. Calculating the Satellite Orbit and Clock Correction

 

All symbols of parameters and constants are listed in the Table 1 and Table 2 or based on their derivatives.

 

Firstly, the Semi-major axis is found as

The mean motion is computed as

The time of transmissionis modified from TOW from navigation data and the correction of is performed as

The mean anomaly is found as

The change rate of mean anomaly is calculated as

Eccentric anomaly can be found through iterative calculation as

The change rate of E is calculated as

The relativistic correction term is given as

Where F is a constant with the value of

 

The clock correction is them performed as

The GPS time at time of transmittion is corrected again as

The true anomaly is calculated as

The change rate of the true anomaly is calculated as

The argument of latitude is found as

The change rate of the latitude is

The following correction terms are generated fromandas

 

The argument of latitude, radius and inclination and their change rate are corrected as

 

The positions in orbital plane is

The change rate of which is

The corrected longitude of ascending node is

The change rate of which is

The position and velocity of satellites are the calculated as

 

2.2. Plot of Satellite Orbit and Clock Correction

 

The RINEX files are downloaded from IGS site located at Boulder, CO. Though the files are 24 hours long, certain GPS satellites are only visible for several hours. The broadcast ephemeris for a specific satellite is also valid for several hours.

 

For the RINEX files downloaded for this report, GPS satellite PRN 17 is chosen to be presented. The valid time duration of this satellite is from 2008/11/07 15:35:00 UTC to 2008/11/07 21:53:00 UTC, which corresponds to GPS time of 488100 seconds to 510780 seconds at GPS week of 1504.

 

The later analysis and comparisons of the GPS satellite orbits and clock correction are all based on PRN 17 in the above time duration. By using the equations listed in the section 2.1, the satellite orbit and clock corrections of PRN 17 are depicted below.

Figure 1. Plot of GPS satellite orbit calculated from the broadcast ephemeris

 

A significant jump of about 0.05 ms is seen in the plot of the satellite clock correction, the same discrepancy also happens in the plot of satellite position and velocity. The reason of this discrepancy is that an update of broadcast ephemeris happens at the jumping time and the calculations of the satellite orbit and clock correction before and after the jumping time use different ephemeris data.

 

3. Precise Ephemeris

 

Different from the broadcast orbits, the precise orbits are estimated from Jet Propulsion Laboratory with direct 3-D satellites position at ECEF coordinate and their clock corrections. Therefore, no keplerian computation is required in obtaining the satellite orbits. The precise satellite orbits are stored in the standard SP3 files, and directly read from such files, the satellite orbits and clock correction of PRN 17 for the 24 hours time span are depicted below.

Figure 2. Plot of GPS satellite orbit read from the SP3 file

 

However the limit of precise orbits is that the reported satellite position and clock correction are at an interval of 15 minutes. Thus if users want to have the satellite position at a random given time, they have to interpolate the orbits through the time.

 

3.1 Methods of Interpolation

 

One could easily the N order polynomial interpolation which is shown below as [5]

 

 

Where C represents the X, Y, Z coordinate values, T is the GPS time and A0 through An are coefficients of the polynomial.

 

Low order of interpolation using the form above would be sufficient for the satellite clock corrections since they are not as dynamic as position data as shown in the figure 2. However it would bring a significant degradation for the satellite position information. Increased order of polynomial interpolation could ease the degradation, but due to the rising complexity of calculation, it will begin to fail in realizing the interpolation in the computer if the order of the polynomial increases to a certain level.

Considering the characteristic of the satellite 3-D position over the 24 hours time span shown in the figure 2, the X and Y have a cyclic period of 24 hours and Z has a period of 12 hours. Mark Schenewerk [5] suggests a trigonometric form of data interpolation on the satellite 3-D position data, shown as

 

 

In the equation shown above, C, A and T have the same definition as the first equation, W is a characteristic frequency that equals totimes the ratio of a solar to sidereal days length.

 

Based on the tests made by Mark Schenewerk [5], 9 terms of trigonometric interpolation would be adequate to achieve a satisfactory accuracy. Matlab codes are written in implementing the 9 terms trigonometric interpolation at 30 seconds interval and results are plotted below for PRN 17 from the GPS time of 488100 seconds to 510780 seconds at the GPS week of 1504, as same as section 2.2.

 

Figure 3. Plot of GPS satellite orbit interpolated from the precise ephemeris

 

4. Comparison

 

4.1 Satellite Orbit and Clock Correction

 

The satellite position and clock correction are similar in scale as shown in figure 1 and 3. However, as described in previous sections, the accuracy of the satellite position calculated from broadcast ephemeris is about 160 cm and 5cm for the precise ephemeris. Therefore the satellite position in the figure 3 is more accuracy than that in the figure 1, so is the satellite clock corrections. The differences are compared in the figure 4 shown below.

Figure 4. Comparison of GPS satellite orbit and clock correction between the precise and broadcast orbits

 

The discontinuity is seen in the figure 4 between the GPS time of 489600 seconds and 489630 seconds. As explained in the section 2.2, the reason of this discontinuity is due to the update of broadcast ephemeris and the calculation of the satellite orbit and clock correction use different ephemeris data before and after GPS time 489630 seconds, the latter one is more accurate.

 

Figure 5 depicts the above difference again from the GPS time 489630 seconds, in which the discontinuity still exists with smaller magnitude. This is discontinuity is also due to the regular update of broadcast ephemeris every two hours, it is also shown that the precise orbit does not suffer from this issue and turns to be continuous in the full time duration.

Figure 5. Comparison of GPS satellite orbit and clock correction between the precise and broadcast orbits modified.

 

4.2 User Position

 

In this section, a comparison is made to evaluate the performance of satellite orbit and clock correction calculation in affecting the final GPS user position solution.

 

The satellites 3-D positions are obtained either from broadcast orbits or precise orbits as introduced in previous sections. Satellites observation measurements such as Pseudorange are obtained from the observation RINEX files. By knowing those information from more than four GPS satellites, the user position could be calculated using the least square solution.

 

Details of the least square solution is not the focus of this report, but can be referred to Misra P & Enge P [7]. Matlab codes of user position least square solution are given in the attachment of this report, based on least square solution, the user positions computed from both broadcast and precise satellite orbits are shown in the figure below. The time range is from 2008/11/07 00:00:00 UTC to 2008/11/07 06:14:30 UTC, or GPS time of 432000 seconds to 454470 seconds at GPS week of 1504. The total time duration is 22500 seconds at 30 seconds time interval.

 

Figure 6. Estimated ECEF coordinates of user position based on broadcast and precise Orbits

 

Figure 6 shows that the 3-D user positions calculated from broadcast and precise orbits are very close with each other, and it is hard to judge which one has better accuracy.

 

Taking the approximate XYZ position in the header of the observation file as the reference user position, the relative user position to the reference can be found and plotted as

Figure 7. Relative user position to the reference position

 

The statistics of the user position in the figure 6 are listed in the table below.

 

 

Broadcast orbit

Precise Orbit

Mean

Std

Mean

Std

X

-4.55 m

1.88 m

-4.32 m

1.55 m

Y

-16.09 m

3.16 m

-16.06 m

2.79 m

Z

13.56 m

2.62 m

13.24 m

2.37 m

 

Table 3.Statistics of relative user positions

 

From the table 3, one may notice that the accuracy of estimated user position using precise orbits is only slightly better than using the broadcast orbits. Both user positions are very noisy as shown in the figure 6 and 7. The reason is that the least square solution uses code pseudorange measurements rather than carrier phase measurements to calculate user position. Besides, the program combines code pseudorange at L1 and L2 frequency to obtain an Ionospheric free measurement that further increases the noise. The carrier phase based least square solution could achieve an accuracy of centimeter level, but the implementation is also more complicated than the using the code pseudorange. Details about how to obtain the Ionospheric free measurements and carrier phase based least square solution can be found at Misra P & Enge P [7].

 

One may also find that the means of user positions are several to tens meters biased from the reference coordinates, the reason is that tropospheric errors are not removed but remain in the corrected measurement, bringing a certain amount of bias in the final positioning solution. Besides, multipath, noise and other un-modeled error also contribute to this bias.

 

It could be seen that when using the code measurements to calculate user position, though the satellite orbits and clock corrections are more accurate comparing with those using broadcast orbit, the former brings little improvements to the final user position, due to the noises of using the code measurements. The precise orbit will help when computing the precise user position using carrier phase measurement, the accuracy of which is in centimeters level. It is suggested to use the precise orbits rather than the broadcast orbits when dealing with some even more precise GPS applications such as seismic research and analysis.

5. Conclusion

 

This report introduces the methods of obtaining GPS satellites orbits and clock corrections from broadcast and precise orbits.

 

The steps and algorithms of calculation through broadcast ephemeris are given, by using which satellites orbits and clock corrections could be obtained at any given time. However, the limit of using the broadcast orbits is that the calculated orbits and clock corrections are not very precise for some applications. Besides, the broadcast ephemeris requires frequent update every two hours, or the precision would degrade dramatically. The receiver will need to make a choice of using a relatively more accurate ephemeris data at updating time, an inappropriate choice will result in computed orbits at large errors.

 

For using the precise orbits which maintains a high precision, the issue is to interpolate the reported satellite orbits and clock corrections from 15 minutes interval to any given time, Mark Schenewerk [5] suggested a 9 orders trigonometric form of data interpolation which proved to be satisfactory.

 

Comparisons of satellite orbits, clock corrections and user position, calculated from both broadcast orbits and precise orbits were made in this report, results showed that though precise orbits has higher precision, it does little improvements when calculating user position using code measurements. But the precise orbits will be beneficial in the GPS applications that require high precision.

 

Actually, in the real GPS receiver design, the single frequency receiver and real time kinematics receiver will decode the broadcast ephemeris from the satellites navigation messages. The precise orbits, storing in the SP3 files are usually used in post-processing applications.

 

Reference

 

1.       David L.M. Warren, John F. Raquet (2003), Broadcast vs. precise GPS, ephemerides: a historical perspective, GPS Solution 2003 7:151-156

2.       International GNSS Service, http://igscb.jpl.nasa.gov/components/prods.html

3.       ICD-GPS-200, October, 1993.

4.       Benjamin W. Remondi (2004), Calculate Satellite Velocity Using The Broadcast Ephemeris, GPS Solutions 2004 8:181-183.

5.       Mark Schenewerk (2003), A brief review of basic GPS orbit interpolation strategies, GPS Solution 2003 6: 265-267

6.       http://earth-info.nga.mil/GandG/sathtml/sp3format.html

7.       Misra P & Enge P (2001), Global Positioning System: Signals, Measurements, and Performance. Ganga-Jamuna Press ISBN 0-9709544-0-9

 

 

 

 

Appendix

 

Some of the Matlab codes attached in this report are modified from the Matlab toolbox for course ASEN 5090 and ASEN 6090. The codes directly used from the tool are indicated in the codes description below. Codes without tag are original works.

 

RINEX and SP3 files are downloaded from the IGS website.

 

1.     Files Description

 

(1). NISU.OBS

RINEX observation file

 

(2). NISU.NAV

RINEX navigation file

 

(3). NISU.SP3

Precise orbits file

 

1.     Codes Description

 

(1). calSatUserPos.m

This is the main code file to parse the broadcast and precise ephemeris, calculate or interpolate satellites position and store them in the proper data structure. Besides, it also does the user positioning.

 

(2). cal2gpstime.m (toolbox)

This function converts calendar date/time to GPS week/time.

 

(3). comp_pos.m

Least square solution for user position

 

(4). check_rinex_line_length.m (toolbox)

This function takes a string and checks if it < 80 characters long. If so the string is "padded" with zeros until it is 80 characters long. Useful when reading RINEX files.

 

(5). get_broadcast_orbits.m

Calculate satellites orbits based on broadcasted ephemeris

 

(6). get_clk_corr.m

Get clock corrections from SP3 file

 

(7). get_precise_orbits.m

Get satellites orbits from the precise ephemeris/SP3 file

 

(8). GPSweek.m (toolbox)

This function finds GPS week and GPS second of the week based on the input calendar date and time.

 

(9). Parsef.m (toolbox)

parsef parse string value using FORTRAN formatting codes

 

(10). plotBroadcast.m

Plot the satellites orbits and clock corrections from the broadcast orbits

 

(11). plotCompare.m

Plot the comparison of satellites orbits obtained from the precise orbits and broadcast orbits

 

(12). plotPrecise.m

Plot the satellites orbits and clock corrections from the precise orbits

 

(13). plotSP3Orbits.m

Directly plot the contents in the SP3 file (no interpolations)

 

(14). plotUserPos.m

Plot of user position solved by least square solution

 

(15). precise_orbit_interp.m

The 9th order interpolation of precise orbits based on the algorithm introduced in the reference 5.

 

(16). read_rinex_nav.m (toolbox)

Read the RINEX navigation file

 

(17). read_rinex_obs.m (toolbox)

Read the RINEX observation file

 

(18). read_sp3.m (toolbox)

Read the SP3 file

 

(19). set_constants (toolbox)

This script contains many useful constants for GPS and related work.

 

(20). wgsxyz2lla.m (toolbox)

Convert the user position from ECEF coordinates to WGS84 latitude, longitude and height.

 

(21). createObs.m

Create the data structure for storing the satellite orbits for both broadcast and precise orbits.

 

2.     Start the program

 

[satOrbits,userPos,alt1,alt2]=calSatUserPos;