Orbit Propagation and Determination of Low Earth Orbit Satellites

This paper represents orbit propagation and determination of low Earth orbit (LEO) satellites. Satellite global positioning system (GPS) configured receiver provides position and velocitymeasures by navigating filter to get the coordinates of the orbit propagation (OP).Themain contradictions in real-time orbit which is determined by the problem are orbit positioning accuracy and the amount of calculating two indicators.This paper is dedicated to solving the problem of tradeoffs. To plan to use a nonlinear filteringmethod for immediate orbit tasks requires more precise satellite orbit state parameters in a short time. Although the traditional extended Kalman filter (EKF) method is widely used, its linear approximation of the drawbacks in dealing with nonlinear problems was especially evident, without compromising Kalman filter (unscented Kalman Filter, UKF). As a new nonlinear estimationmethod, it is measured at the estimatedmeasurements onmore andmore applications.This paper will be the first study onUKFmicrosatellites in LEO orbit in real time, trying to explore the real-time precision orbit determination techniques. Through the preliminary simulation results, they show that, based on orbit mission requirements and conditions using UKF, they can satisfy the positioning accuracy and compute two indicators.


Introduction
The satellite orbit determination (OD) estimates discrete observation of the position and velocity of an orbiting object.The set of observations includes the measurements from the space based GPS receiver (GPSR) that is located on the object itself.Satellite orbit propagation (OP) estimates the future state of motion of an object, whose orbit has been determined from past observations.The satellite's motion is described by a set of approximate equations of motion.The degree of approximation depends on the intended use of orbital information.Observations are subject to systematic and random uncertainties; therefore, orbit determination and propagation are probabilistic.
The satellite is influenced by a variety of external forces, including terrestrial gravity, atmospheric drag, multibody gravitation, solar radiation pressure, tides, and spacecraft thrusters.Selection of forces for modeling depends on the accuracy and precision required from the OD process and the amount of available data.The complex modeling of these forces results in a highly nonlinear set of dynamical equations.Many physical and computational uncertainties limit the accuracy and precision of the object state that may be determined.Similarly, the observational data are inherently nonlinear with respect to the state of motion of the object and some influences might not have been included in models of observation of the state of motion.
The remainder of this paper is organized as follows.Section 2 describes the methodology of GPSR based orbit determination; Section 3 is a brief introduction of the disturbance mathematical model; Section 4 legends the orbit determination algorithm description; Section 5 describes the GPS measurement models; description of OP algorithm settings is included in Section 6; and Section 7 offers the conclusion.

Methodology of GPSR Based Orbit Determination
Three basic strategies are presently used to determine precise LEO orbits with GPS.They are the dynamic, kinematic or nondynamic, and hybrid or reduced-dynamic strategies.

International Journal of Antennas and Propagation
The dynamic orbit determination approach [1] requires precise models of the forces acting on user object.This technique has been applied to many successful space vehicle missions and has become the mainstream of precision OD (POD) approach.Dynamic model errors are the limiting factor for this technique, such as the geopotential model errors and atmospheric drag model errors, depending on the dynamic environment of the user space vehicle.With the continuous, global, and high precision GPS tracking data, dynamic model parameters, such as geopotential parameters, can be tuned effectively to reduce the effects of dynamic model error in the context of dynamic approach.The dense tracking data also allows for the frequent estimation of empirical parameters to absorb the effects of unmodeled or mismodeled dynamic errors.
The kinematic or geometric approach does not require the description of the dynamics except for possible interpolation between solution points for the user object, and the orbit solution is referenced to the phase center of the on-board GPS antenna instead of the space vehicle's center of mass.Yunck and Wu [2] proposed a geometric method that uses the continuous record of object position changes obtained from the GPS carrier phase to smooth the position measurements made with pseudorange.This approach assumes the accessibility of -codes at both the  1 and  2 frequencies.Byun [3] developed a kinematic orbit determination algorithm using double-and triple-differenced GPS carrier phase measurements.Kinematic solutions are more sensitive to geometrical factors, such as the direction of the GPS satellites and the GPS orbit accuracy, and they require the resolution of phase ambiguities.
The previous two strategies each have counterbalancing disadvantages: various mismodelling errors in dynamic OD and GPS measurement noise in kinematic OD.A hybrid dynamic and kinematic OD strategy would down weight the errors caused by each strategy but still utilize the strengths of each.One such strategy has been devised and is referred to as reduced dynamic orbit determination.The reduced-dynamic approach [1] uses both geometric and dynamic information and weighs their relative strength by solving local geometric position corrections using a process noise model to absorb dynamic model errors.

Orbit Propagation Algorithm Description.
The orbit propagation algorithm can be divided into two main tasks: orbit determination and orbit prediction (propagation).The general diagram of orbit propagation algorithm is described in Figure 1.
The orbit determination algorithm is based on Unscented Kalman Filter (UKF) and estimates the object state vector x+1 and covariance matrix P+1 from discrete observations.The set of observations includes the measurements {raw data} from the space based GPS receiver that is located on the space vehicle itself.The orbit determination algorithm includes the orbit prediction task as time update stage of UFK.
Orbit prediction algorithm calculates the future state of motion of a vehicle x+1 whose orbit has been determined from past observations.Moreover, the covariance matrix P+1 is propagated.A numerical integration of the dynamic model is applied for orbit prediction.
The OP solution  +1 ,  +1 is output data of orbit propagation algorithm.The OP solution and covariance matrix can be obtained as from prediction task as from determination task.The following external data are required for OP solution calculation: (i) init time and state vector  init ,  init for algorithm initialization/reinitialization; (ii) the time moment  +1 to new OP solution  +1 calculation; (iii) the set of observations  SV +1 , {raw data} for new estimation x+1 calculation.
The following input data are obtained from the previous OP solutions calculation: (i) the last OP solution   ,   ; (ii) the time  est of last calculation of estimation x+1 ; (iii) the covariance matrix   .
The maximal time of continuous propagation  prp max , maximal integration time step ℎ max , minimal count of available Navigation SV  SV min , and default covariance matrix  def are used for algorithm control.

Dynamic Model.
A dynamic model of the object motion essentially adds a priori knowledge from the equations of the orbital motion to the kinematic position knowledge as obtained from the raw GPS measurements.In our case, the dynamic model incorporates the complex Earth gravity field (EGM 96) truncated to order and degree 18.Furthermore, the Sun and Moon gravitation and atmospheric drag are accounted.
The differential dynamic equation of motion is given by where V  , V  , and V  are the ECEF velocity components of object, , , and  are the ECEF radius vector components of object,  is the receiver clock bias,  is the receiver clock drift,  -body is the Sun and Moon gravitation forces,  GEO is the acceleration due to geopotential, F drag = {  ,   ,   } is a perturbing force due to aerodynamic drag, and Ω  is the angular velocity of Earth rotation, hire and below Ω  = 7.2921151467 − 5 rad/sec.The user coordinates are in the rotating Earth-fixed frame (ECEF).Although this adds some complexity, especially due to the Coriolis and centrifugal acceleration in the dynamic model, no reference system transformations are required in the main program since input (initial position and velocity) and OP output are consistently referring to the Earth-fixed frame.In this way, reference system transformations may completely be encapsulated in the dynamic model.Moreover, some dynamic algorithms, which compute the acceleration due to the Earth's gravity field and the atmospheric drag, may be formulated simpler in an Earth-fixed than in an inertial frame.
The integration is carried out by using the simple fourth order Runge-Kutta algorithm without any mechanism of step adjustment or error control.The fourth order Runge-Kutta is considered an adequate numerical integrator due to its simplicity, fair accuracy, and low computational burden.Numerical integration is performed in the rotating Earthfixed frame (ECEF).
According to [4] magnitude ratio of atmospheric drag and solar radiation for average size spherical objects with   = 2.4 moving along the circular LEO, solar radiation can be neglected because its effect on total object acceleration is much smaller than acceleration due to perturbing geopotential, the third body forces from the Sun and the Moon, and atmospheric drag.
We can see that for altitudes less than 600 km solar radiation pressure is significantly smaller than atmospheric drag.Furthermore, atmospheric drag decreases with altitude and it becomes negligible for altitudes higher than 700 km (Table 1).

Disturbance Mathematical Model
3.1.Earth Gravity.The gravitational potential function for the solid-body mass distribution of the Earth is generally expressed in terms of a spherical harmonic expansion, referred to as the geopotential in the Earth-fixed reference frame (ECEF).The gravitational potential function  is defined as [5] where  is the gravitational potential function (m In ( 2) and ( 3),  = 1 for  = 0 and  = 2 for  ̸ = 0.The series is theoretically valid for  ≥ .The acceleration due to geopotential can be defined as where Variables , , and  are related to object ECEF radius vector components , , and  by According to (6), Projections of Earth gravitation force in ECEF are The following recurrence equations can be applied to cos  and sin  calculation:

N-Body Perturbation.
The gravitational perturbations of the Sun, Moon, and other planets can be modeled with sufficient accuracy using point mass approximations.In the geocentric inertial coordinate system, the -body accelerations can be expressed as [4] where  is the universal gravitational constant,   mass of the th perturbing body (Sun or Moon),   position vector of the th perturbing body in ECEF,   −  position vector of the th perturbing body with respect to the object mass center in ECEF, and  planet index, where  =  for the Sun and  =  for the Moon.The values of the Sun and Moon position vectors   can be obtained from the following equations: where Ω is a transfer matrix from current Equatorial Earth Centered Inertial Frame to ECEF defined as with Ω  the angular velocity of Earth and  time in seconds from the beginning of current sidereal day defined as (Greenwich Sidereal Mean Time (GSMT) (see [6] for details));  is a transfer matrix from Ecliptic Earth Centered Inertial Frame of J2000.0 to current Equatorial Earth Centered Inertial Frame defined in the following equation: with  the mean obliquity of the ecliptic as defined in [6];    is radius vector of Sun mass center in Ecliptic Earth Centered Inertial Frame of J2000.0 defined as mean distance between Earth and Sun mass centers AU: astronomical unit, hire and below AU = 1.4959787011 m;   : is the mean anomaly of the Sun, see [6,7] for detail;   : is the ecliptic longitude of the Sun defined as: : the mean longitude of the Sun as defined below: ,   , , , Ω: are fundamental arguments of Moon motion theory, they are defined in [6].All equations of this item are given according to [6,7].

Atmospheric Drag.
A near-Earth space vehicle of arbitrary shape moving with some velocity V in an atmosphere will experience drag force.The drag force acceleration can be modeled as [4] where  is the atmospheric density, V  is the object velocity vector relative to the atmosphere,  is the mass of the object,   is the drag coefficient for the object, and  is the crosssectional area of the main body perpendicular to V  .The parameter   =   / is sometimes referred to as the ballistic coefficient.It is varied during an orbital motion due to an object attitude and an object mass center evolution and others factors.The middle (typically) value of a ballistic coefficient is used because this factors are unknown for OP algorithm.Realistically, chosen values   = 0.01 m 2 /kg.
Different empirical dynamical atmospheric models can be used for computing the atmospheric density.These include the Jacchia 71 [8], Jacchia 77 [9], the drag temperature model (DTM) [10], DTM-2000 [11], MSIS-90 [12], and NRLMSISE-00 [13].The density computed by using any of these models could be in error anywhere from 10% to over 200% depending on solar activity.A deal settings are used for aforementioned atmospheric models computation, for example, geomagnetic activity index, and daily and average solar flux index.They are fluctuated during orbital flight and must be monitored.Sizeable density errors can be acquired otherwise.Furthermore, all abovementioned models require appreciable computation resources.
According to aforesaid in the current project local atmosphere density model [4] is employed.It is a rough density model relative to dynamical models, but this model is very easy for computation and requires no settings monitoring.
Equations for density calculation are where ℎ 1 is the reference height, ℎ 1 = 500000 m;  1 is the atmospheric density on reference height,  1 = 2.−13 kg/m 3 ; ℎ is the object height;  is the height scale of the uniform atmosphere. is the distance from the Earth's center of mass;  is the semimajor axis of the WGS 84 Ellipsoid.

Prediction Algorithm Description
Orbit prediction algorithm calculates the future state of motion of a vehicle whose orbit has been determined from past observations.Moreover, the covariance matrix is propagated.
To construct the future object trajectory, the orbit prediction algorithm uses the dynamic equation of motion given in Section 2.1.This fundamental equation of mechanics provides the dynamic constraint governing the orbit solution.The true acceleration at any instant depends on the space vehicle position and velocity at that instant and on many other parameters that characterize the forces at work.The predicted trajectory is then generated by integration of where ℎ is the integration step limited by ℎ max (see Figure 1 for details); Ṽ(  + ℎ) = {V  , V  , V  } is the predicted ECEF velocity vector of the object; r(  + ℎ) = {, , } is the predicted ECEF radius vector of the object; b(  + ℎ) is the predicted International Journal of Antennas and Propagation  receiver clock bias; d(  + ℎ) is the predicted receiver clock drift; V(  ), (  ), (  ), (  ) define last object state.The object state derivatives vector is defined using dynamic motion model which is described in Section 2.1.As the numerical integrator, we will use a simple fourth order Runge-Kutta algorithm without any mechanism of step adjustment or error control.Numerical integration is performed in an ECEF reference frame.
The covariance matrix propagation is defined below.
The differential dynamic equations of motion are given by ẋ =  (, ) , where  = (V  , V  , V  , , , , , )  is a state vector that includes the spacecraft position and velocity vectors and the receiver clock bias and drift.
The propagation of x+1 from the previous state for covariance matrix propagation is generated by the following reduced equation: where ℎ =  +1 −  is the integration step, limited by ℎ max (see Figure 1 for detail),   is the state vector from the previous step, and   (, ) is the reduced dynamic model of notion witch is defined as follows: ) .
(28) The diagram of orbit determination algorithm is described as follows.
The following process and measurement models can be established: The variables in the above equation will be described: x  is a system condition vector in the  moment, (⋅) is unscented system model, w  is a dynamic mixed signal in the moment, z  is a measuring dynamic vector in the  moment, ℎ(⋅) is a unscented system measuring model, and k  is dynamic measuring mixed signal in the  moment.
The measurement vector is denoted by z  , and the process noise w  and the measurement noise k  are assumed to be zero-mean white noise.The process noise covariance matrix and the measurement noise covariance are given by Q  and R  , respectively.
The system error covariance matrix Q  is as follows: The measuring error covariance matrix R  is as follows: Here, w  and k  are independent and unrelated: The parameter is a scaling parameter defined as The positive constants   ,  = 1, 2, 3 are used as parameter of the method, and a priori and a posteriori estimates of the state are denoted by x−  and x .(2) Time Updating.Condition predicted value is Condition predicted average is Covariance matrix is (3) Observation Updating.Observation measurement predicted value is Observation measurement predicted average is and   are updated as follows:  2.
The time update phase of the UKF includes the propagation of state vector x − +1 from the previous object state and the state covariance matrix  +1 , which is defined in detail in Section 2.1.
The subsequent measurement update adjusts the state vector x − +1 components and state covariance matrix  +1 to best fit the GPS pseudorange and pseudorange rate measurement data.

The Measurement Update
Phase.The measurement residual and sensitivity matrix are found by forming the computed observation equation.
The model for a GPS pseudorange measurement is given by where is the geometric range; x+1 , ỹ+1 , z+1 are the positional states of the user object at reception time;  GPS ,  GPS , and  GPS are the positional states of the th GPS satellite according to item; b+1 is the receiver clock offset; and   +1 accumulates all unmodeled errors.
Using the abovementioned nonlinear measurement equation, the sensitivity matrix is given by the Jacobian matrix of partial derivatives of nonlinear measurement vector with respect to the state vector : The model for a GPS Doppler measurement is given by where V  +1 ,   +1 are the th GPS satellite velocity vector and radius vector according to item; Ṽ+1 is the user object velocity; r+1 is the coordinates of the user object at reception time; d+1 is the receiver clock drift; and   +1 accumulates all unmodeled errors of the Doppler observation.
Using the abovementioned nonlinear measurement equation, the sensitivity matrix is given by the Jacobian matrix of partial derivatives of nonlinear measurement vector with respect to the state vector : where where  GPS ,  GPS , and GPS , (V  )  GPS , and (V  )  GPS }; , , and  represent r+1 = {x +1 , ỹ+1 , z+1 }; V  , V  , and , and (Ṽ  )  +1 }; and  represents   .If both pseudorange and Doppler measurement are used, the sensitivity matrix will be composed of   and  V matrices in the following way: where   and  V are matrixes size of [ SV +1 × 8] defined as The measurement residuals or innovations sequence is where   ( +1 ,  +1 ) and  V ( +1 ,  +1 ) are the measured pseudoranges and pseudorange rates which are computed according to section; ỹ ( x+1 ,  +1 ) and ỹV ( x+1 ,  +1 ) are the predicted pseudoranges and pseudorange rates which are computed according to section.The measurement update phase uses the Kalman equations to incorporate the information given by the measurements themselves and obtains improved estimates of the state and of the covariance as follows: where  +1 is the discrete measurement noise covariance which is basically a measurement weight matrix.The QR-decomposition algorithm is applied to inverse matrix calculation.The general idea of this algorithm is described in detail.

GPS Measurement Models
The basic measurement types that will be employed in the current project are GPS pseudorange and Doppler in  1 frequency.The equation of the pseudorange in  1 frequency is given by where    is the pseudorange in  1 ,    is the geometric range to the th satellite at the observation instant   is given by is the ionospheric delay,  is the vacuum speed of light,  GPS (  ) is the GPS satellite clock offset,   (  ) is the receiver clock offset,   is the observation instant in GPS time, and   is a remnant error supposed random gaussian.
The numerically controlled oscillator (NCO), which controls the carrier tracking loop, provides an indication of the observed frequency shift of the received signal.This observed frequency differs from the nominal  1 frequency because of where V   is the th GPS satellite velocity at the observation instant   , V  is the receiver velocity, LOS   is the line-of-sight to the th GPS satellite at   , and  1 = 1575.42MHz is the transmitted frequency.
The Doppler can be converted to a pseudorange rate observation given by where  is the receiver clock drift in m/s and    is the error in observation in m/s.
Another possible GPS measurement model is a linear combination of GPS  1 C/A code and carrier phase.Since both data types are affected by systematic ionospheric errors with the same magnitude but opposite signs, their arithmetic mean is free of ionospheric errors.This approach, as proposed by Yunck in 1996 [1], removes the dominant systematic error source for raw GPS data, which may amount to 10-20 m [14] at low elevations.As a matter of fact, the resulting so-called Group and Phase Ionospheric Calibration (GRAPHIC) data provide a low-noise biased range with an accuracy of half the C/A code noise.A drawback of using the GRAPHIC data type originates from the employed carrier phase data which introduce range biases for each of the twelve receiver channels.As consequence, twelve range biases have to be adjusted as part of the estimation process which significantly complicates the orbit determination algorithms.Finally, it has to be noted that GPS broadcast ephemeris errors with a mean standard deviation of about 3 m (3D position) and 1 m (User Equivalent Range Error, UERE) are still present in real-time applications [15], if no counter measures, such as the upload of precise ephemeris, are taken.

OP Algorithm Settings
6.1.Integration Settings.The maximal time of continuous propagation  prp max = 2400 seconds (it is specified in 0).The maximal integration time step ℎ max = 30 seconds.It was defined according to Table 2 which describes maximal Runge-Kutta method errors, respectively, to integration step.The period of dynamic model integration is one turn.

Dynamic Model Settings.
The maximal half interval of multibody acceleration fixing  mb max = 30 sec.It was defined according to Table 3 which describes integration errors, respectively, to the half interval.The period of dynamic model integration is one turn.
The fixed multibody acceleration components are available on time interval  fix ±  mb max , where  fix is time of acceleration fixing.

Estimation Settings. The minimal count of available
Navigation SV  SV min = 2 (it is defined by test results).The discrete state-noise covariance matrix , the discrete measurement noise covariance , and the initial covariance matrix of  0 can have different components values and structure for special receiver application.According to the requirements, 0 in current protects them and can be defined, for example, as the follows:  where  pr = 55.9 m is the pseudorange measurement dispersion and   = 0.037 m 2 /sec 2 is the pseudorange rate measurement dispersion.

Simulation and Analysis
Using the Kalman algorithm to estimate orbit propagation and determination in this section has been proposed to simulate and validate the derivation of the formula.Simulation results are shown in Figure 2, and the initial conditions are selected at the beginning of the track after leaving a balance within 200 seconds after convergence.The simulation results as expected are as follows.Deliberately made a real track star with an initial value is not the same as kalman filtering, 200 seconds after the kalman algorithm converges within 200 sec.Figure 3(b) shows time response of microsatellite velocity measured, estimated and difference between measured and estimated.The same as in Figure 3(a), a willfully made real track star with an initial value is not the same kalman filtering; at 200 seconds, using the kalman algorithm, it converges within 200 sec.
Simulation results show that, when the microsatellite attitude and orbit maintain balance, satellite orbital position and velocity are in the estimated and the measured value via GPS satellite computer considerable amount.Initial value changes in a relatively small error, the maximum error of 10 ∘ .If the number of GPS by the change, the number of GPS consists of four changes into three or less, the error is relatively large when the attitude changes, the maximum error is instantaneous 32 ∘ .In the academic theory and engineering practice, a systematic analysis is generally considered to be consistent with the conclusion that it is feasible.

Figure 2 :
Figure 2: The algorithm flow chart of UKF.

7. 1 .
Simulated Conditions.Microsatellite altitude 500 km, longitude 108 ∘ and latitude 35 ∘ , sampling time 4 sec, onmodulator magnitude = 2, satellite attitude motion trajectory is shown in Figure3.The UKF of direct observation equation is used in the simulation.

Figure 3 (
Figure3(a), time response of microsatellite measured, estimated, and difference between measured and estimated.Deliberately made a real track star with an initial value is not the same as kalman filtering, 200 seconds after the kalman algorithm converges within 200 sec.Figure3(b)shows time response of microsatellite velocity measured, estimated and difference between measured and estimated.The same as in Figure3(a), a willfully made real track star with an initial value is not the same kalman filtering; at 200 seconds, using the kalman algorithm, it converges within 200 sec.

Table 2
Doppler shifts produced by the GPS satellite and user motion, as well as the frequency error or drift of the satellite and user clocks.The equation of the Doppler shift in  1 frequency is given by