Abstract
This paper represents orbit propagation and determination of low Earth orbit (LEO) satellites. Satellite global positioning system (GPS) configured receiver provides position and velocity measures by navigating filter to get the coordinates of the orbit propagation (OP). The main 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 filtering method 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 estimation method, it is measured at the estimated measurements on more and more applications. This paper will be the first study on UKF microsatellites 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.
1. 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.
2. 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.
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 and 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.
2.1. 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 and covariance matrix from discrete observations. The set of observations includes the measurements 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 whose orbit has been determined from past observations. Moreover, the covariance matrix is propagated. A numerical integration of the dynamic model is applied for orbit prediction.
The OP solution 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 for algorithm initialization/reinitialization;(ii)the time moment to new OP solution calculation;(iii)the set of observations for new estimation calculation.
The following input data are obtained from the previous OP solutions calculation:(i)the last OP solution ;(ii)the time of last calculation of estimation ;(iii)the covariance matrix .
The maximal time of continuous propagation , maximal integration time step , minimal count of available Navigation SV , and default covariance matrix are used for algorithm control.
2.2. 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 , , and 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, is the Sun and Moon gravitation forces, is the acceleration due to geopotential, is a perturbing force due to aerodynamic drag, and is the angular velocity of Earth rotation, hire and below 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 Earth-fixed frame (ECEF).
According to [4] magnitude ratio of atmospheric drag and solar radiation for average size spherical objects with 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).
3. 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 (), is the Earth’s gravitational constant, hire and below , is the distance from the Earth’s center of mass , is the Semimajor axis of the WGS 84 Ellipsoid, hire and below , and are the degree and order, respectively; is the geocentric latitude; is the geocentric longitude; and are the normalized gravitational coefficients defined in EGM-96 model [5], is the normalized associated Legendre function, is the associated Legendre function, and is the Legendre polynomial. Consider
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 and calculation:
3.2. -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 : astronomical unit, hire and below ; : 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: : is Julian centuries from J2000.0; : is radius vector of Moon mass center in Ecliptic Earth Centered Inertial Frame of J2000.0 which is defined as: : is the distance between Earth and Moon mass centers defined as: : is the Moon geocentric longitude. It can be defined as: : is the mean longitude of the Moon; : is the Moon geocentric latitude as defined in the following equation: , , , , : are fundamental arguments of Moon motion theory, they are defined in [6].
All equations of this item are given according to [6, 7].
3.3. Atmospheric Drag
A near-Earth space vehicle of arbitrary shape moving with some velocity in an atmosphere will experience drag force. The drag force acceleration can be modeled as [4] where is the atmospheric density, 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 cross-sectional area of the main body perpendicular to .
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 .
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 is the reference height, ; is the atmospheric density on reference height, ; 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.
4. 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 (see Figure 1 for details); is the predicted ECEF velocity vector of the object; is the predicted ECEF radius vector of the object; is the predicted receiver clock bias; is the predicted receiver clock drift; 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 is a state vector that includes the spacecraft position and velocity vectors and the receiver clock bias and drift.
The propagation of from the previous state for covariance matrix propagation is generated by the following reduced equation: where is the integration step, limited by (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:
4.1. Orbit Determination Algorithm Description
The orbit determination algorithm applies an UKF to estimate the state vector which comprises 8 components:(i)object velocity ,(ii)object position ,(iii)receiver clock bias ,(iv)receiver clock drift ,(v).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: is a system condition vector in the moment, is unscented system model, is a dynamic mixed signal in the moment, is a measuring dynamic vector in the moment, is a unscented system measuring model, and is dynamic measuring mixed signal in the moment.
The measurement vector is denoted by , and the process noise and the measurement noise are assumed to be zero-mean white noise. The process noise covariance matrix and the measurement noise covariance are given by and , respectively.
The system error covariance matrix is as follows:
The measuring error covariance matrix is as follows:
Here, and are independent and unrelated:
4.2. Unscented Kalman Filter Processing
(1) Producing Standard Points and Calculation Value. Consider
The parameter is a scaling parameter defined as
The positive constants , are used as parameter of the method, and a priori and a posteriori estimates of the state are denoted by and .
(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:
(4) Calculating Kalman Gain Value. Consider
(5) Updating Estimated Value to Measurement Value. Consider
(6) Updating Condition Error Covariance Matrix. Consider
4.3. Unscented Kalman Filter Flow Chart
The flow chart of the UKF is as shown in Figure 2.
The time update phase of the UKF includes the propagation of state vector from the previous object state and the state covariance matrix , which is defined in detail in Section 2.1.
The subsequent measurement update adjusts the state vector components and state covariance matrix to best fit the GPS pseudorange and pseudorange rate measurement data.
4.4. 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; are the positional states of the user object at reception time; , , and are the positional states of the th GPS satellite according to item; is the receiver clock offset; and 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 , are the th GPS satellite velocity vector and radius vector according to item; is the user object velocity; is the coordinates of the user object at reception time; is the receiver clock drift; and 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 , , and represent ; , , and represent , , and ; represent , , ; , , and represent , , and ; and represents .
If both pseudorange and Doppler measurement are used, the sensitivity matrix will be composed of and matrices in the following way: where and are matrixes size of defined as The measurement residuals or innovations sequence is where and are the measured pseudoranges and pseudorange rates which are computed according to section; and 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 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.
5. GPS Measurement Models
The basic measurement types that will be employed in the current project are GPS pseudorange and Doppler in frequency. The equation of the pseudorange in frequency is given by where is the pseudorange in , 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, 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 frequency because of 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 frequency is given by where is the th GPS satellite velocity at the observation instant , is the receiver velocity, is the line-of-sight to the th GPS satellite at , and MHz 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 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.
6. OP Algorithm Settings
6.1. Integration Settings
The maximal time of continuous propagation seconds (it is specified in 0).
The maximal integration time step 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.
6.2. Dynamic Model Settings
The maximal half interval of multibody acceleration fixing 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 , where is time of acceleration fixing.
6.3. Estimation Settings
The minimal count of available Navigation SV (it is defined by test results).
The discrete state-noise covariance matrix , the discrete measurement noise covariance , and the initial covariance matrix of can have different components values and structure for special receiver application. According to the requirements 0 in current protect them can be defined, for example, as the follows:where is the pseudorange measurement dispersion and is the pseudorange rate measurement dispersion.
7. 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.
7.1. Simulated Conditions
Microsatellite altitude 500 km, longitude 108° and latitude 35°, sampling time 4 sec, on-modulator magnitude = 2, satellite attitude motion trajectory is shown in Figure 3. The UKF of direct observation equation is used in the simulation.
(a)
(b)
Figure 3(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. 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.
Disclosure
This paper represents original, unpublished material that is not under editorial consideration elsewhere, in whole or in part, other than in abstract form.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.