2.1. Satellite State Vector in ECI Using Keplerian Orbit Elements

Keplerian orbit elements are determined using the equatorial and orbital plane. Figure 1 shows the Keplerian orbit elements in the ECI and the perifocal coordinate frame. The elements are the following:

• a is the semi major axis;

• $\vec e$ is the eccentricity vector;

• i is the inclination angle;

• Ω is the right ascension of ascending node (RAAN);

• ω is the argument of perigee;

• ν is the true anomaly;

• E is the eccentricity anomaly; and

• U is the argument of latitude.

Figure 1. Classical orbit elements in the ECI and the perifocal coordinate frame.

There are two main parts to determining the satellite position and velocity. The first part is the calculation of the position and velocity in the orbital plane (perifocal coordinate frame), and the second part is to transform the position and velocity from the orbital plane to the ECI coordinate frame. The position determination of the satellite in the orbital plane can be expressed as Equation (1) (Montenbruck and Gill, Reference Montenbruck and Gill2005; Misra and Enge, Reference Misra and Enge2006).

(1)$$\vec r_{orb} = \left[ \matrix{x_{orb} \cr y_{orb} \cr 0 } \right] = \left[ {\matrix{ {r\cos \nu} \cr {r\sin \nu} \cr 0 \cr}} \right] = \left[ {\matrix{ {a\cos E - ae} \cr {a\sqrt {1 - e^2} \sin E} \cr 0 \cr}} \right]$$

The velocity of the satellite in the orbital plane can be derived as:

(2)$$ \dot {\vec r}_{orb} = \displaystyle{{a \cdot n} \over {1 - e\cos E}}\left[ {\matrix{ { - \sin E} \cr {\sqrt {1 - e^2} \cos E} \cr 0 \cr}} \right]$$

where n is the mean motion of the satellite and is expressed as:

(3)$$ n = \sqrt {\displaystyle{{GM} \over {a^3}}} $$

G is the gravitational constant and M is mass of the Earth.

The acceleration of the satellite in the orbital plane can be expressed as:

(4)$$ \ddot {\vec r} = - \displaystyle{{GM} \over {r^3}} \cdot {\vec r}$$

The state vector of the satellite in the orbital coordinate can be transformed to the ECI coordinate frame using the rotation matrix. The transform matrix from the orbital plane to the ECI coordinate frame is:

(5)$$ \vec r_{ECI} = R_3 ( - \Omega )R_1 ( - i)R_3 ( - \omega )\vec r_{orb} $$

The R₁, R₂ and, R₃ rotation matrices are:

(6)$$ \eqalign{R_1 (\theta ) & = \left[ {\matrix{ 1 & 0 & 0 \cr 0 & {\cos \theta } & {\sin \theta } \cr 0 & { - \sin \theta } & {\cos \theta } \cr } } \right],R_2 (\theta ) = \left[ {\matrix{ {\cos \theta } & 0 & { - \sin \theta } \cr 0 & 1 & 0 \cr {\sin \theta } & 0 & {\cos \theta } \cr } } \right], \cr R_3 (\theta ) & = \left[ {\matrix{ {\cos \theta } & {\sin \theta } & 0 \cr { - \sin \theta } & {\cos \theta } & 0 \cr 0 & 0 & 1 \cr } } \right]} $$

where θ is the arbitrary angle that you want to rotate. The position of the satellite in the ECEF coordinate frame can be transformed using Equation (7).

(7)$$ \vec r_{ECEF} = \left[ W \right]\left[ R \right]\left[ N \right]\left[ P \right]\left[ B \right]\vec r_{ECI} $$

W, R, N, P and B are the polar motion, the rotation, the nutation, the precession and the bias transformational matrices (Vallado and McClain, Reference Vallado and McClain2001; IERS, 2010; Bradley, Vallado et al., Reference Bradley and Vallado2011). A detailed description of the ECEF to the ECI transformation algorithm will be discussed in Section 4.

2.2. GPS Satellite State Vector in the GPS Orbital Plane

The orbital coordinate system described in IS-GSP-200F is different from the natural (or classical) orbital coordinate system. In the natural orbital coordinate system, the direction of the x-axis is defined by the direction of the perigee, but the direction of the x-axis in the IS-GPS-200F orbital coordinate system is defined by the direction of ascending node (Zhang, Zhang et al., Reference Zhang and Zhang2006). In IS-GPS-200F, the position of the GPS satellite is calculated in the GPS orbital plane and is then transformed into the ECEF coordinate frame by applying two rotation matrices; however, the position in the classical orbital plane is first transformed into the ECI and is then transformed into the ECEF if needed.

The reason that the IS-GPS-200F system produces a position in the ECEF coordinates is because, in most cases, ground-based users want the position only in the ECEF coordinate frame. The ECEF-to-ECI coordinate transform algorithm requires Earth Orientation Parameters (EOP) and additional complicated numerical calculations. Therefore, this orbital plane to ECEF direct transformation method reduces the computational burdens to the ground-based user who does not need positional data in the ECI coordinate frame.

However, this method is a weak point for the space-based user who wants the positional data in both the ECI and the ECEF frames. The user that wants the position in the ECI frame using the LNAV message must have additional EOP data, which are provided by other agencies. This requirement means that position in the ECI frame cannot be calculated using the LNAV message alone and EOP data must be obtained from other sources through an external data link. To improve this burdensome process, the EOP data will be included in the CNAV message type 32.

There are two main procedures in the position calculation algorithm described in IS-GSP-200F. The first step is calculating the position in the CNAV orbital plane, and the second step is rotating the position from the GPS (or CNAV) orbital plane to the ECEF frame using equation (8). $\vec r_{CNAVorb} $ is the position vector of the GPS satellite in the CNAV orbital plane, and $\vec r_{ECEF} $ is the position vector in the ECEF coordinate frame. R _l(−i _k) and R ₃(−Ω_k) are the rotation matrices that rotate the position vector by the CNAV angle of inclination and the CNAV longitude of ascending node. It is important to note that the CNAV longitude of the ascending node is completely different from the classical orbit element. The inclination angle of the CNAV is slightly different from the inclination of the classical orbit elements in the ECI coordinate frame.

(8)$$\vec r_{ECEF} = R_3 ( - \Omega _k )R_1 ( - i_k )\vec r_{CNAVorb} = R_i^e \vec r_{CNAVorb} $$

(9)$$R_i^e = \left[ {\matrix{ {\cos \Omega _k} & { - \sin \Omega _k \cos i_k} & {\sin \Omega _k \sin i_k} \cr {\sin \Omega _k} & {\cos \Omega _k \cos i_k} & { - \cos \Omega _k \sin i_k} \cr 0 & {\sin i_k} & {\cos i_k} \cr}} \right]$$

The velocity determination equation is not described in IS-GPS-200F, but the velocity can be calculated by differentiating the position determination equation analytically. The simple format of the analytical velocity determination equation is represented as:

(10)$$ \dot {\vec r}_{ECEF} = \dot R_i^e \vec r_{CNAVorb} + R_i^e \dot {\vec r}_{CNAVorb}$$

The acceleration determination algorithm can be calculated by differentiating the velocity equation. The simple format of the analytical acceleration equation is represented as:

(11)$$ \ddot {\vec r}_{ECEF} = \ddot R_i^e {\vec r}_{CNAVorb} + 2 \dot R_i^e \dot {\vec r}_{CNAVorb} + R_i^e \ddot {\vec r}_{CNAVorb} $$

2.3. GPS SATELLITE POSITION IN ECEF

The position determination algorithm described in IS-GPS-200F is divided into the position calculation part and the rotation from the CNAV orbital coordinate frame to the ECEF coordinate frame. The position determination equation in the CNAV orbital coordinate frame is:

(12)$$\vec r_{CNAVorb} = \left[ {\matrix{ {X_{CNAVorb}} \cr {Y_{CNAVorb}} \cr 0 \cr}} \right] = \left[ {\matrix{ {r_k \cos u_k} \cr {r_k \sin u_k} \cr 0 \cr}} \right]$$

where r _k is the radius of the satellite in the CNAV orbital plane, which is calculated using A _k, e_n and E_k. The value of e _n is a constant value that is broadcast in the CNAV ephemeris message. The equation to calculate the radius:

(13)$$r_k = A_k \left( {1 - e_n \cos E_k} \right) + \delta r_k $$

where A _k is determined by:

(14)$$A_k = A_0 + \dot A \cdot t_k {\rm} \Leftarrow {\rm} A_0 = \left( {A_{REF} + \Delta A} \right)$$

A _REF is predefined in IS-GPS-200F, $\dot A$ and ΔA are broadcast parameters in the CNAV ephemeris and are new parameters that were not included in the LNAV. In equation (13), E _k is an eccentricity anomaly at epoch time t _k that is determined by M _k and e _n. The relationship of M _k, E_k and e _n is:

(15)$$M_k = E_k - e_n \sin E_k $$

where M _k is the mean anomaly that changes in time, n _A is the rate of M _k, and the relationship between M _k and n _A is defined in Equation (16), where M ₀, Δn ₀ and $\Delta \dot n_0 $ are broadcast in the CNAV ephemeris.

(16)$$\eqalign{n_0 & = \sqrt {\displaystyle{\mu \over {A_0^3 }}} \cr \Delta n_A & = \Delta n_0 + \displaystyle{1 \over 2}\Delta \dot n_0 \cdot t_k \cr n_A & = n_0 + \Delta n_A {\rm } \cr M_k & = M_0 + n_A \cdot t_k } $$

δr _k in Equation (13) is the argument of latitude correction term for the radius in the CNAV orbital plane. δr _k is determined by the following algorithm:

(17)$$\delta r_k = C_{rs - n} \sin (2\Phi _k ) + C_{rc - n} \cos (2\Phi _k )$$

where C _rs−n and C _rc−n are the harmonic perturbation correction coefficients that are broadcast in CNAV ephemeris. Φ_k is the argument of the latitude and is determined by:

(18)$$\eqalign{\Phi _k & = \nu _k + \omega _n \cr \nu _k & = \tan ^{ - 1} \left( {\displaystyle{{\sqrt {1 - e_n^2 } \sin E_k } \over {\cos E_k - e}}} \right)} $$

where ν _k is the true anomaly, and ω _n is the argument of the perigee in CNAV orbital plane. The value of ν _k changes in time, and ω _n is a broadcast value in the CNAV ephemeris.

In Equation (12), μ _k is the corrected argument of the latitude given by:

(19)$$\eqalign{u_k & = \Phi _k + \delta u_k \cr \delta u_k & = C_{us - n} \sin (2\Phi _k ) + C_{uc - n} \cos (2\Phi _k )} $$

where C _us−n and C _uc−n are the amplitude of the harmonic correction terms for the argument of latitude and are broadcasted in the CNAV ephemeris.

In Equation (9), $R_i^e $ is the rotation matrix from the CNAV orbital plane to the ECEF coordinate frame. The rotation matrix$R_i^e $ is:

(20)$$\eqalign{\Omega _k & = \Omega _0 + ({\rm \dot \Omega} - {\rm \dot \Omega} _e ) \cr {\rm \dot \Omega} & = {\rm \dot \Omega} _{REF} + \Delta {\rm \dot \Omega} } $$

where:

• Ω₀ is the reference of the right ascension;

• Ω is the rate of right ascension;

• Ω_REF is the predefined reference rate of right ascension; and

• ΔΩ is the correction of the rate of right ascension.

The i _k term and its related terms in Equation (9) have had their notation changed for clarity because $\dot i_k $ (the rate of inclination angle) can be misread as i _k (the inclination angle). The modified notations are:

(21)$$\eqalign{I_k & = i_k \cr I_{0 - n} & = i_{0 - n} \cr \dot I_{0 - n} & = i_{0 - n} - DOT \cr \delta I_k & = \delta i_k } $$

The inclination angle I _k is calculated using the rate of inclination angle $\dot I_{0 - n} $ and the inclination correction δI _k that is given in Equation (22). In Equation (22), C _is−n and C _ic−n are the amplitudes of the harmonic correction term to the inclination and are broadcasted in the CNAV ephemeris.

(22)$$\eqalign{I_k & = I_{0 - n} + ( \dot I_{0 - n} ) \cdot t_k + \delta I_k \cr \delta I_k & = C_{is - n} \sin 2\Phi _k + C_{ic - n} \cos 2\Phi _k } $$

2.4. GPS SATELLITE VELOCITY IN ECEF

The GPS satellite velocity $\dot \vec {r}_{CNAVorb} $ can be determined using Equation (23). In Equation (23), $\dot \vec {r}_{CNAVorb} $ and $\dot R_i^e $ are needed to calculate the velocity.

(23)$$ \dot {\vec r}_{CNAVorb} = \left[ {\matrix{ { \dot X_{CNAVorb}} \cr { \dot Y_{CNAVorb}} \cr 0 \cr}} \right] = \left[ {\matrix{ { \dot r_k \cos u_k - r_k \sin u_k \cdot \dot u_k} \cr { \dot r_k \sin u_k + r_k \cos u_k \cdot \dot u_k} \cr 0 \cr}} \right]$$

where $\dot r_k $ is the differential of the radius in the CNAV orbital plane and is calculated by:

(24)$$ \dot r_k = \dot A_k (1 - e_n \cos E_k ) + A_k \cdot e_n \cdot \sin E_k \cdot \dot E_k + [\delta r_k ]\prime $$

In Equation (24), $\dot A_k $, $\dot E_k $ and $\left[ {\delta r_k} \right]\prime $ are needed to calculate derivative of radius $\dot r_k $. $\dot A_k $ can be set to broadcast CNAV ephemeris term $\dot A$ as:

(25)$$ \dot A_k = \dot A$$

$\dot E_k $ can be obtained by differentiating Equation (15) and the differential of the eccentricity anomaly equation is:

(26)$$\eqalign{& \dot E_k = \displaystyle{{n_A} \over {1 - e_n \cos E_k}} \cr & n_A = n_0 + \Delta n_0 + \displaystyle{1 \over 2}\Delta \dot n_0 \cdot t_k} $$

The mean anomaly value in Equation (26) is different from that of the LNAV, as the rate of the mean motion is added in the CNAV.

In Equation (24), the derivative of the radial correction is:

(27)$$[\delta r_k ]\prime = 2\left( {C_{rs - n} \cos (2\Phi _k ) - C_{rc - n} \sin (2\Phi _k )} \right) \dot \nu _k $$

where $\dot \nu _k $ is the derivative of the true anomaly. The derivative of the true anomaly can be defined by:

(28)$$ \dot \nu _k = \displaystyle{{A_k ^2 \sqrt {1 - e_n^2} \cdot n_A} \over {r_k^2}} $$

The derivative of the corrected argument of latitude $\dot u_k $ is:

(29)$$\eqalign{ \dot u_k = & {\rm \dot \Phi} _k + 2\left( {C_{us - n} \cos (2\Phi _k ) - C_{uc - n} \sin (2\Phi _k )} \right){\rm \dot \Phi} _k \cr = & \dot \nu _k + 2\left( {C_{us - n} \cos (2\Phi _k ) - C_{uc - n} \sin (2\Phi _k )} \right) \dot \nu _k} $$

${\rm \dot \Phi} _k $ can be obtained by differentiating Equation (18) under the condition that ω _n is a constant value.

$\dot R_i^e $ in Equation (10), can be obtained by differentiating the rotation matrix $R_i^e $. This is a simple process that only requires the derivative of each component and the result is:

(30)$$ \dot R_i^e = \! \left[ {\matrix{ { - \sin \Omega _k \cdot {\rm \dot \Omega} _k } & { - \cos \Omega _k \cos I_k \cdot {\rm \dot \Omega} _k + \sin \Omega _k \sin I_k\cdot \dot I_k } & {\cos \Omega _k \sin I_k \cdot {\rm \dot \Omega} _k + \sin \Omega _k \cos I_k \cdot \dot I_k } \cr {\cos \Omega _k \cdot {\rm \dot \Omega} _k } & { - \sin \Omega _k \cos I_k \cdot {\rm \dot \Omega} _k - \cos \Omega _k \sin I_k \cdot \dot I_k } & {\sin \Omega _k \sin I_k \cdot {\rm \dot \Omega} _k - \cos \Omega _k \cos I_k \cdot I_k } \cr 0 & {\cos I_k \cdot \dot I_k } & { - \sin I_k \cdot \dot I_k } \cr } } \right]$$

where ${\rm \dot \Omega} _k $ and $\dot I_k $ are:

(31)$${\rm \dot \Omega} _k = {\rm \dot \Omega} _{REF} + \Delta {\rm \dot \Omega} - {\rm \dot \Omega} _e $$

2.5. GPS Satellite Acceleration in ECEF

To obtain the GPS satellite acceleration, two additional variables are needed: the second order differential of the rotation matrix $ \ddot R_i^e $ and the second order differential of the radius $ \ddot \vec r_{ECEF} $ in the CNAV orbital plane. To obtain the second order differential of the rotation matrix $ \ddot R_i^e $, each component of the first order derivative of the rotation matrix $\dot R_i^e $ is simply differentiated. In this differential calculation, $\ddot I_k $ and $\ddot \Omega _k $ are set to zero, because these values have only first differential values in the CNAV ephemeris. The second order of differential of the rotation matrix is:

(32)$$ \ddot R_i^e = \left[ {\matrix{ { - \cos \Omega _k \cdot {\rm \dot \Omega} _k^2} & {\sin \Omega _k \cdot {\rm \dot \Omega} _k^2 \cdot \cos I_k} & { - \sin \Omega _k \cdot {\rm \dot \Omega} _k^2 \cdot \sin I_k} \cr { - \sin \Omega _k \cdot {\rm \dot \Omega} _k^2} & { - \cos \Omega _k \cdot {\rm \dot \Omega} _k^2 \cdot \cos I_k} & {\cos \Omega _k \cdot {\rm \dot \Omega} _k^2 \cdot \sin I_k} \cr 0 & 0 & 0 \cr}} \right]$$

The second order differential of the radius $\ddot \vec {r}_{CNAVorb} $ in the CNAV orbital plane is defined in Equation (4). The acceleration of the GPS satellite is:

(33)$$ \ddot {\vec r}_{ECEF} = \ddot R_i^e {\vec r}_{CNAVorb} + 2 \dot R_i^e \dot {\vec r}_{CNAVorb} + R_i^e \ddot {\vec r}_{CNAVorb} $$

4.1. Virtual CNAV Message

To test the proposed algorithm, the virtual CNAV ephemeris data were constructed using the precise GPS position and velocity data from the NGA and the LNAV ephemeris data from International GNSS Service (IGS). The procedure of generating the virtual CNAV is described in Figure 2. The main difference between LNAV and CNAV is that CNAV has two new additional terms: the rate of mean motion and the rate of semi-major axis motion. These two terms are generated using the NGA precise GPS data and the LNAV ephemeris correction terms. In this paper, we used data from 1 March 2011 from time 04:00:00 to 06:00:00.

Figure 2. Virtual CNAV generation procedure using the LNAV ephemeris and the precise GPS ephemeris

4.2. Correctness Test

To check the correctness of the proposed position, velocity and acceleration determination algorithm, we compared these values to the numerical velocity and acceleration. The analytical velocity $\dot \bar {r}_{ECEF} $ was determined using Equation (10), and the numerical velocity $\dot \bar {r}_{Num\_ECEF} $ was generated using Equation (34). The same procedure was performed for the correctness test on the acceleration. The analytical acceleration $ \ddot \bar {r}_{ECEF} $ and the numerical acceleration $ \ddot \bar {r}_{Num\_ECEF} $ were compared using the virtual CNAV. The analytical acceleration was generated using Equation (11), and the numerical acceleration $\ddot \bar r_{Num\_ECEF} $ was generated using Equation (35).

(35)$$ \dot {\bar r}_{Num\_ECEF} (t) = \displaystyle{{\bar r_{ECEF} (t + \Delta t) - \bar r_{ECEF} (t - \Delta t)} \over {2\Delta t}}$$

(36)$$ \ddot {\bar r}_{Num\_ECEF} (t) = \displaystyle{{ \dot {\bar r}_{ECEF} (t + \Delta t) - \dot {\bar r}_{ECEF} (t - \Delta t)} \over {2\Delta t}}$$

where Δt was set to 1 second.

The comparison of the analytical and the numerical position and velocity is presented in Figure 3 and Figure 4. The difference between the analytical and the numerical velocity is less than 0·00001 m/s for each axis, and the difference in the acceleration is less than 0·00008 m/s². With these results, we conclude that the proposed algorithm has continuity and correctness in terms of GPS satellite position, velocity and acceleration determination.

Figure 3. Velocity difference between the analytical and the numerical methods.

Figure 4. Acceleration difference between the analytical and numerical methods.

4.3. Accuracy of GPS Satellite State Vector Using CNAV

The position, velocity and acceleration of the PRN5 GPS satellite were calculated using the virtual CNAV ephemeris, and these results were compared to the NGA precise position, velocity and acceleration of the PRN5 GPS satellite. The simulation start time was set as 1 March 2011 04:00:00 and the simulation had a duration of two hours and step intervals of 15 minutes. The NGA precise GPS position and velocity data were assumed to be the true data, and the numerically calculated acceleration data using the NGA precise GPS velocity were also assumed to be the true acceleration data.

During the same time and PRN conditions, the GPS satellite position was calculated using the LNAV and the CNAV. The residual of the position is shown in Figure 5. Over the two hours of the simulation, we observed that the calculated position data using the virtual CNAV is more accurate than using the LNAV. The maximum error of the position calculated using the virtual CNAV is 0·34 m and that of the position calculated using LNAV is 0·69 m. By considering the maximum value of the residual at each epoch, the virtual CNAV shows more accurate results than the LNAV.

Figure 5. Position residuals of the LNAV and CNAV ephemerides.

Next, we checked the accuracy of the calculated velocity using the LNAV and the virtual CNAV. The residuals of the velocity, which are the differences between the true and calculated velocity of the satellite are shown in Figure 6. During the two hours of the simulation, the residual of the velocity using the virtual CNAV is less than that of the LNAV except for three epochs. This result is caused by the limited accuracy of the virtual CNAV ephemeris.

Figure 6. Velocity residuals of the LNAV and CNAV ephemerides.

The residuals of the acceleration calculated using the LNAV and the virtual CNAV are shown in Figure 7. The result shows that the acceleration error residuals of the LNAV and the virtual CNAV are almost the same. This result is caused by the fact that the virtual CNAV is not a complete CNAV; for instance, the inclination and the right ascension of ascending node are adapted from the LNAV and the LNAV derivation value as well. These terms are used when determining the acceleration in Equations (29) and (31); as a result, the accelerations of these two values are approximately the same when using the virtual CNAV ephemeris. The residual of the calculated acceleration is less than 0·0001 m/s².

Figure 7. Acceleration residuals of the LNAV and CNAV ephemerides.

Note that the CNAV in each figure is the virtual CNAV that was constructed using the LNAV data and the NGA precise GPS data. The accuracy of the calculated position, velocity and acceleration data is not exactly the expected value that is generated using the broadcast CNAV; however, the accuracy of the virtual CNAV is greater than that of LNAV, and we expect that the accuracy of the broadcast CNAV will be better than that of the virtual CNAV.

4.4. COORDINATE TRANSFORMATION SIMULATION

IAU-2000 Resolution, which is the equinox-based ITRF-to-GCRF transformation, is the standard reference algorithm; therefore, we assumed that the transformed result using the equinox-based transformation is true. If not using EOP parameters, the ECEF to the ECI transformation can be performed by a simple rotation about the z-axis neglecting effects due to wobble, precession and nutation. When applying simple rotation, the rotation angle is set as the Greenwich Apparent Sidereal Time (GAST) (GPSD, 2011). This simple rotation method was compared to the IAU-2000 equinox-based transformation to evaluate the accuracy improvement when using the EOP data. The EOP and UT1 data from the International Earth Rotation and Reference Systems Service (IERS) are used in this simulation (IERS, 2012). The position, velocity and acceleration transformation results are shown in Figures 8 to 10. The maximum position transformation error using the simple rotation is approximately 70 km in the x-axis and several tens of kilometres along each axis. The maximum velocity transformation error is approximately 9 m/s in x-axis, and the maximum acceleration transformation error is 0·015 m/s in x-axis. Therefore, we concluded that when not using EOP data and the IAU-2000 resolution transformation algorithm, the accuracy of the frame transformation result will be degraded.