Coarse Initial Orbit Determination for a Geostationary Satellite Using Single-Epoch GPS Measurements

A practical algorithm is proposed for determining the orbit of a geostationary orbit (GEO) satellite using single-epoch measurements from a Global Positioning System (GPS) receiver under the sparse visibility of the GPS satellites. The algorithm uses three components of a state vector to determine the satellite’s state, even when it is impossible to apply the classical single-point solutions (SPS). Through consideration of the characteristics of the GEO orbital elements and GPS measurements, the components of the state vector are reduced to three. However, the algorithm remains sufficiently accurate for a GEO satellite. The developed algorithm was tested on simulated measurements from two or three GPS satellites, and the calculated maximum position error was found to be less than approximately 40 km or even several kilometers within the geometric range, even when the classical SPS solution was unattainable. In addition, extended Kalman filter (EKF) tests of a GEO satellite with the estimated initial state were performed to validate the algorithm. In the EKF, a reliable dynamic model was adapted to reduce the probability of divergence that can be caused by large errors in the initial state.


Introduction
Many ground-tracking networks and facilities are required to track the position of geostationary orbit (GEO) satellites if a ground-tracking system is used [1,2]. The Global Positioning System (GPS) receiver can provide high-accuracy position data at a low cost; thus, it is reasonable to use GPS receivers in GEO satellites. However, position accuracy, which is calculated using a single-point solution (SPS) algorithm with snapshot measurements, is low compared with that of low Earth orbit (LEO) satellites or ground users; in certain cases, no result is produced. These disadvantages occur because the GPS signal power is lower than that of the LEO satellite, and the geometrical configuration of the GPS satellites relative to the GEO satellites is unfavorable. Therefore, an orbit determination (OD) filter is needed to overcome these unfavorable circumstances [3].
There are two main types of orbit determination (OD) filters: post-processing and real-time processing techniques. A real-time OD filter is needed to operate a geostationary orbit (GEO) satellite instantly and properly; one such filter is the extended Kalman filter (EKF), which is well known for its accuracy and efficiency. The EKF convergence time is determined by the initial state and conditions [4]. However, most studies do not consider how to determine the initial state and its conditions [5]. More than four GPS satellites are observable in low Earth orbit (LEO); thus, single-point solutions (SPS) can be applied at any time and the result can be used as the initial states of the EKF with a high level of accuracy. Typically, fewer than four GPS satellites are observable at the GEO; therefore, SPS are not always applicable. In these situations, alternative solutions must be employed, such as the short arc batch technique, which is a post-processing technique that requires measurements over a long period of time; thus, this technique does not provide navigation data instantly [6,7].
In this paper, we developed a coarse initial orbit determination algorithm to improve the accuracy of the initial EKF states under the sparse visibility of GPS satellites. The application of the proposed algorithm is illustrated in Figure 1. We used the characteristics of GEOs to develop our algorithm: GEO is an almost circular orbit, and its inclination angle is nearly zero. We set the minimum number of state components to calculate the state of the satellite using snapshot measurements under the sparse visibility of GPS satellites. This algorithm can determine the GEO satellite's state vector even when fewer than four GPS satellites are visible, and it is very practical because it does not require long-term measurements. The remainder of this paper is organized as follows. In Section 2, the details of the coarse initial orbit determination algorithm are explained. In Section 3, the general EKF-based orbit determination scheme is discussed. In Section 4, simulations of the proposed algorithm and EKF are performed to test the accuracy and the availability of the proposed algorithm. Finally, the results of the developed algorithm and EKF are discussed.

Coarse Initial Orbit Determination Algorithm
Four unknown variables appear in the classical SPS algorithm using the GPS signal: position components ( , , x y z ) and receiver clock bias error ( r δ ). Thus, measurements from at least four GPS satellites are required. The user position can typically be calculated using the GPS at a LEO satellite at any time because more than four GPS satellites are observed at the LEO satellite's altitude. The error of the calculated single point position of a LEO is less than several dozens of meters and can be used as the initial state of the EKF. However, it is impossible to calculate the position using the classical SPS algorithm when fewer than four GPS satellites are visible at the GEO. More than four GPS satellites are infrequently visible at the GEO satellite; thus, the point position is not always determined. Thus, we must wait until more than four GPS satellites are visible to obtain navigation data.
We developed a coarse initial orbit determination algorithm to calculate a point solution using measurements obtained from the receiver with two or three GPS satellites. The algorithm uses a minimum number of state variables, which were selected by considering the characteristics of the GEO. The classical orbital elements of the ideal GEO are shown in Figure 2. The right ascension of the ascending node (RAAN) of the ideal GEO cannot be defined because the inclination angle is zero. The eccentricity is also zero, and thus, the argument of perigee cannot be defined either. Thus, we can define the state of the GEO satellite using only its true longitude value. The true longitude is a useful term when defining circular and equatorial orbits, and its equation is given by If we assume that the satellite rotates in the ideal GEO, the states of a satellite can be approximately expressed by one element: the true longitude. We included two additional variables (the clock bias and clock bias rate of the receiver) in the state vector of the GEO satellite because these error components are included in the C/A code pseudorange and pseudorange rate. Thus, the state vector of the GEO satellite is defined as where c is the speed of light, r δ is the receiver's clock bias, and r δ  is the receiver's clock bias rate.
The state vector in Equation (2) can be converted into Cartesian coordinates and given by where G is the gravitational constant and M ⊕ is the mass of the Earth. The value of a was set to a constant, and thus, the differential of a is zero. Then, the measurement vector can be defined as where M is the number of visible GPS satellites (two or three in this paper); t r  and t v  are the position and velocity vectors, respectively, of the GPS satellite; r  and v  are the position and velocity vectors, respectively, of the GEO satellite; t δ and t δ  are the GPS satellite clock bias and clock bias rate, respectively; r δ and r δ  are the receiver clock bias and clock bias rate, respectively; ê is the unit direction vector; ρ and ρ  are the pseudorange and pseudorange rate, respectively; and n is the measurement noise from the assumed Gaussian distribution. The ionospheric error is not included in Equations (6) and (7) because we assume that the ionospheric error can be removed by the dual-frequency GPS receiver. We The equations of the measurement matrix Kep H and cart H are given by If we define the state vector as Equation (2), we can calculate the position of the GEO satellite using single-epoch measurements of two or three GPS satellites. In the definition of the state vector in Equation (2), the inclination angle and eccentricity are not included in the variables; however, these variables exist in the real GEO orbit and could increase the error of the state vector calculated iteratively. The geometric state error estimated by the proposed algorithm can be defined as 3 3 3 cos( ) where Ω is the ascending node, i is the inclination, ω is the argument of perigee, ν is the true anomaly, r is the radius, est λ is the estimated true longitude, and est cδ is the estimated receiver clock bias.
In Equation (20), ( ) 2 z is a constant bias term that cannot be reduced or removed by the proposed algorithm; therefore, the error could be increased if the z-component has great value. However, we can still determine the state vector as accurately as possible given the sparse visibility of the GPS satellites. After calculating the state vector using Equations (9)-(18), the vector can be used as the initial state value for the EKF with a practical level of accuracy.

EKF Scheme
The EKF is well known for its accuracy and speed; thus, it is used in non-linear system applications, such as real-time OD. We will not explain the entire algorithm in detail, as it is not the focus of this paper. Thus, the equations of the EKF algorithm scheme are summarized in Table 1 [5,8]. Table 1. Extended Kalman filter (EKF) processing scheme using GPS measurements.

Time update
Predicted state estimation: Linear approximation: Predicted covariance matrix:

Measurement update
Predicted Measurement: Measurement matrix linearization: Kalman gain computation: Updating posteriori covariance matrix: We used two-body gravity and other perturbations to predict the state of the satellite. The equation of the general accelerations acting on a satellite is [5]: where r  is the position vector of the satellite in the Earth-centered inertial (ECI) coordinate system, μ is the gravitational constant of the Earth and perturbed a  is the other perturbed acceleration.
We included gravitational attraction by a nonspherical central body, third-body effects (Sun and Moon) and solar-radiation pressure in the perturbed accelerations [5,9]. The gravitational attraction by the nonspherical central body is expressed as: where R ⊕ is the radius of the Earth, p is the associated Legendre polynomials, sat φ is the latitude of the satellite, sat λ is the longitude of the satellite, and C and S are the gravitational coefficients.
The equation of the third-body effect is expressed as: where ⊕ μ is the gravitational constant of the Earth, 3 μ is the gravitational constant of the third body, The acceleration by solar radiation pressure is included and is one of the significant accelerations acting on GEO satellites. The solar radiation pressure is expressed as: where SR P is the solar pressure, R C is the reflectivity, A  is the effective area of the satellite and m is the mass of the satellite. We introduced some errors into the dynamics model of the EKF filter compared with the measurement dynamics model to better approximate a real situation, and the errors are listed in Table 2. The degree and order of the geopotential were lowered from 20 to 10, and the satellite area for the solar pressure model was adjusted to create a 10% error [2].
The errors of the initial state were assumed to be approximately 10 km and 0.01 km in position and velocity, respectively. With these values, we set the initial covariance as follows: where ( ) The noise of the range and that of the range rate are Gaussian distributions with standard deviations of 0.01 km and 0.0001 km/s, respectively, in the simulation. We set R as follows for the case of two or three GPS satellites, respectively: where ( )  The elements of the process noise covariance must be approximately 10 −16 , with the same level of error as the dynamics model of the filter; however, we tuned these values with consideration for the initial state error and convergence time [10,11]. We determined Q as follows: where ( )

Simulation Procedure
We chose the Communication, Ocean and Meteorological Satellite (COMS) launched by the Korea Aerospace Research Institute (KARI) as our GEO satellite for simulation. The COMS is located at 128.2° east, and COMS missions are related to Ka-band communication services, meteorological monitoring, and ocean observation [2]. We simulated the satellite's orbit to validate and test the developed algorithm. First, we generated the position and velocity data of COMS for 24 h using a numerical orbit propagation algorithm. Then, the pseudorange and pseudorange rate, which were obtained from the receiver in the GEO satellite, were generated with simulated GPS signals. The developed algorithm was tested using the generated measurements under the condition that two or three GPS satellites were observable. Then, the calculated state vector was set to the initial state of the EKF, and the EKF was processed. The overall simulation procedure is summarized in Figure 3.

Generation of Measurements
The first simulation step was to propagate the GEO satellite's orbit. We propagated the GEO using Cowell's method, which propagates the position and velocity of the satellite by integrating the accelerations caused by perturbations at each time step [5]. We included the geopotential, solar pressure, and third-body gravity (the Sun and Moon). We chose the EGM-96 model as the geopotential, and the degree and order were both set to 20. The Runge-Kutta 68 algorithm was chosen as the numerical integrator, and the integral step was set to 10 s. The initial orbital elements for propagation are listed in Table 2. The GPS satellites' position and velocity were generated every 10 s for 24 h using the Almanac data. After generating the positions and velocities of the GEO and GPS satellites, we checked whether the GPS satellites were observable at each epoch. A GPS satellite was only visible when it was not blocked by the Earth and its signal power was sufficiently strong to be processed by the GPS receiver.
After determining the visibility of the GPS satellites, the C/A code pseudorange and pseudorange rate were calculated. The C/A code pseudorange is given by [12] where ρ is the C/A code pseudorange in L1, RX is the signal reception time, TX is the signal emission time, r σ is the receiver clock bias, t σ is the GPS satellite clock bias, n is noise, τ is the time delay, 0 t is the reference time, 0 a and 1 a are the polynomial coefficients of the GPS satellite clock bias, and 0 b and 1 b are the polynomial coefficients of the receiver clock bias.
The geometric distance from the GPS satellite to the receiver was calculated using ( ) . The GPS signal travels through space at the speed of light, and thus, the receiver does not instantly receive the signal emitted from the GPS satellite. Therefore, the signal reception time is later than the signal emission time. The equation of the elapsed time from the emission to reception is given by τ appears on both sides of Equation (30). Thus, an iteration technique was used to calculate the proper τ . First, ( ) t r RX  was used instead of ( ) t r RX − τ  to calculate the temporary τ on the left side of the equation, Then, the temporary τ was used on the right side of Equation (30) to update the temporary τ , and the iterations continued until τ converged [13,14].
The equation for the pseudorange rate is similar to that for the pseudorange and is given by where r δ  and t δ  are the clock bias rates of the receiver and GPS satellite, respectively. These variables are calculated from the derivatives of r σ and .

Simulation Results
The algorithm was tested using data from the four points (A, B, C and D) selected from the 24 h simulated GEO orbit. We selected four points at 6 h intervals to rigorously validate the algorithm. The simulation times of the selected points A, B, C and D were UTC 00:00:00, 06:00:00, 12:00:00 and 18:00:00, respectively, on 1 January 2006. The position and velocity vector of each point is given in Table 3.  The observable GPS satellites over 24 h are depicted in Figure 4; however, the simulations were performed under controlled conditions in which only two or three GPS satellites were visible. The initial error in the reference position was set to 1000 km at each of the four points, and the reference receiver clock bias was set to 100 km. The noise terms in Equations (28) and (31) were selected from Gaussian distributions with standard deviations of 0.01 km and 0.0001 km/s, respectively.
The positions of the COMS and its visible GPS satellites at each point are shown in Figures 5-12. The red spot represents the COMS, and the yellow spots represent GPS satellites. At each point, we produced scenarios such that two or three GPS satellites were visible at the COMS. Thus, we intentionally chose GPS satellites among those visible to control the number of visible satellites if more than three GPS satellites were visible. For consistency, we simply removed one satellite from the scenario where three GPS satellites were visible such that only two GPS satellites were visible. As shown in the figures, the visible GPS satellites are located behind the Earth and are located closely together.         Tables 4 and 5. We also tested the influence of the position of the third GPS satellite at point D; more than three GPS satellites are observable at point D. We ran simulations with two fixed GPS satellites and a third satellite placed at several different positions. The results are summarized in Table 6. The residual refers to the difference between the calculated and true values. The results indicate that the residual of the calculated position is less than 40 km in range, and three of the four points have residuals of less than several kilometers in range when using three visible GPS satellites. The residuals increase as the z-component of orbital state increases. The maximum error occurs when the z-component of the orbital state is greatest, and the errors are small when the z-components of the orbital state are small. This relationship occurs because the developed algorithm does not include the z-component in the state vector, and thus, the z-component in the real orbit influences the x-and y-components in the state vector. The z-component value increases with the inclination and the relationship between the inclination and error at point D are presented in Figure 13. Based on Figure 13, we can conclude that the accuracy level of the proposed algorithm is high when the GEO satellite's inclination is small.   No significant difference occurs when using measurements from two or three visible GPS satellites, except for point B, where the residual decreases when using the measurements from three GPS satellites. Furthermore, the geometric relationship among the GPS satellites and the GEO also affected the accuracy of the algorithm, as demonstrated by the simulation results for point D. The EKF was tested using the initial state vector calculated by the developed algorithm at each point. The time update was processed using the Runge-Kutta method. The simulation results of the EKF are shown in Figures 14-18. The time for filter convergence varied across the simulation points and the error of the initial state; however, the filter converged within 120 min with an accuracy of 100 m in all simulations. The errors expressed in the RIC frame after the filter is stabilized are shown in Figures 19 and 20, and the error norm is bounded at 30 m when using three GPS satellites. These results were quite acceptable because the EKF filter for a GEO converges very slowly due to the orbit's characteristics. For example, the convergence time of the EKF filter for a GEO is approximately one or two hours under the sparse visibility of GPS satellites [15]. The convergence rate of the filter depends on the geometric location of the GPS satellites and the changing visibility of the GPS satellites, and thus, the convergence rate varies even though the accuracies of the initial conditions are not significantly different.

Conclusions
The main goal of the algorithm is to calculate the state of the GEO satellite using a single-epoch measurement under the sparse visibility of the GPS satellites, which is usually impossible, even when applying the classical SPS algorithm. The proposed algorithm can calculate the position, velocity and receiver's clock bias using only a single-epoch measurement of two or three GPS satellites without data from external sources. Therefore, the calculated result can be used as the initial state as soon as the measurements are generated by the receiver. The resulting maximum range error is less than 40 km, and when using the result as the initial state of the EKF, which uses a very accurate dynamic model, the filter converges to an error of 100 m within 120 min in the worst case. The maximum error norm of the EKF after filter stabilization is bounded at 30 m when three GPS satellites are observable. Given this acceptable result and the benefits of the algorithm, we expect that our algorithm is sufficiently accurate, robust, efficient, and practical for determining the initial orbit of GEO satellites.