# Orbit Propagation and Determination of Low Earth Orbit Satellites

Authors: Ho-Nien Shou
Journal: International Journal of Antennas and Propagation
 (2014)

---

## 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.


---

## Body

## 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. Section2 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 P-codes at both the L1 and L2 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 Figure1.Figure 1
Orbit propagation algorithm diagram.The orbit determination algorithm is based on Unscented Kalman Filter (UKF) and estimates the object state vectorx^k+1 and covariance matrix P^k+1 from discrete observations. The set of observations includes the measurements {rawdata} 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 vehiclex~k+1 whose orbit has been determined from past observations. Moreover, the covariance matrix P~k+1 is propagated. A numerical integration of the dynamic model is applied for orbit prediction.The OP solutiontk+1,xk+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 vectortinit,xinit for algorithm initialization/reinitialization;(ii)
the time momenttk+1 to new OP solution xk+1 calculation;(iii)
the set of observationsNk+1SV,{rawdata} for new estimation x^k+1 calculation.The following input data are obtained from the previous OP solutions calculation:(i)
the last OP solutiontk,xk;(ii)
the timetest of last calculation of estimation x^k+1;(iii)
the covariance matrixPk.The maximal time of continuous propagationTmax⁡prp, maximal integration time step hmax⁡, minimal count of available Navigation SV Nmin⁡SV, and default covariance matrix Pdef 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(1)f(x,t)=(v˙xv˙yv˙zx˙y˙z˙b˙d˙)=(aGEOx+aNbodx+Fxm+Ωe2x+2Ωey˙aGEOy+aNbody+Fym+Ωe2y-2Ωex˙aGEOz+aNbodz+Fzmvxvyvzd0),
where vx, vy, and vz are the ECEF velocity components of object, x, y, and z are the ECEF radius vector components of object, b is the receiver clock bias, d is the receiver clock drift, aN-body is the Sun and Moon gravitation forces, aGEO is the acceleration due to geopotential, Fdrag={Fx,Fy,Fz} is a perturbing force due to aerodynamic drag, and Ωe is the angular velocity of Earth rotation, hire and below Ωe=7.2921151467e-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 Earth-fixed frame (ECEF).According to [4] magnitude ratio of atmospheric drag and solar radiation for average size spherical objects with cx=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 (Table1).Table 1
SV altitude, km
400
500
600
700
800
Solar radiation to atmospheric drag ratio
0.018
0.08
0.27
0.8
2.1


## 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 Figure1.Figure 1
Orbit propagation algorithm diagram.The orbit determination algorithm is based on Unscented Kalman Filter (UKF) and estimates the object state vectorx^k+1 and covariance matrix P^k+1 from discrete observations. The set of observations includes the measurements {rawdata} 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 vehiclex~k+1 whose orbit has been determined from past observations. Moreover, the covariance matrix P~k+1 is propagated. A numerical integration of the dynamic model is applied for orbit prediction.The OP solutiontk+1,xk+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 vectortinit,xinit for algorithm initialization/reinitialization;(ii)
the time momenttk+1 to new OP solution xk+1 calculation;(iii)
the set of observationsNk+1SV,{rawdata} for new estimation x^k+1 calculation.The following input data are obtained from the previous OP solutions calculation:(i)
the last OP solutiontk,xk;(ii)
the timetest of last calculation of estimation x^k+1;(iii)
the covariance matrixPk.The maximal time of continuous propagationTmax⁡prp, maximal integration time step hmax⁡, minimal count of available Navigation SV Nmin⁡SV, and default covariance matrix Pdef 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(1)f(x,t)=(v˙xv˙yv˙zx˙y˙z˙b˙d˙)=(aGEOx+aNbodx+Fxm+Ωe2x+2Ωey˙aGEOy+aNbody+Fym+Ωe2y-2Ωex˙aGEOz+aNbodz+Fzmvxvyvzd0),
where vx, vy, and vz are the ECEF velocity components of object, x, y, and z are the ECEF radius vector components of object, b is the receiver clock bias, d is the receiver clock drift, aN-body is the Sun and Moon gravitation forces, aGEO is the acceleration due to geopotential, Fdrag={Fx,Fy,Fz} is a perturbing force due to aerodynamic drag, and Ωe is the angular velocity of Earth rotation, hire and below Ωe=7.2921151467e-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 Earth-fixed frame (ECEF).According to [4] magnitude ratio of atmospheric drag and solar radiation for average size spherical objects with cx=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 (Table1).Table 1
SV altitude, km
400
500
600
700
800
Solar radiation to atmospheric drag ratio
0.018
0.08
0.27
0.8
2.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 functionU is defined as [5]
(2)U=GMr[1+∑n=2N∑m=0n(Rr)nP-nm(sinϕ)nnnnn×(C-nmcos⁡mλ+S-nmsinmλ)∑n=2N∑m=0n(Rr)nP-nm(sinϕ)],
where U is the gravitational potential function (m2/s2), GM is the Earth’s gravitational constant, hire and below  GM=3986004.418e8m3/s2, r is the distance from the Earth’s center of mass (m), R is the Semimajor axis of the WGS 84 Ellipsoid, hire and below R=3986004.418e8m, n and m are the degree and order, respectively; ϕ is the geocentric latitude; λ is the geocentric longitude; C-nm and S-nm are the normalized gravitational coefficients defined in EGM-96 model [5], P-nm(sinϕ) is the normalized associated Legendre function, Pnm(sinϕ) is the associated Legendre function, and Pn(sinϕ) is the Legendre polynomial. Consider
(3)P-nm(sinϕ)=[(n-m)!(2n+1)k(n+m)!]1/2Pnm(sinϕ),Pnm(sinϕ)=(cos⁡ϕ)mdmd(sinϕ)m[Pn(sinϕ)],Pn(sinϕ)=12nn!dnd(sinϕ)n(sin2ϕ-1)n.In (2) and (3), k=1 for m=0 and k=2 for m≠0.The series is theoretically valid forr≥R.The acceleration due to geopotential can be defined as(4)aGEOx=dU(r,ϕ,λ)dx=∂U(r,ϕ,λ)∂rdrdx+∂U(r,ϕ,λ)∂ϕdϕdx+∂U(r,ϕ,λ)∂λdλdx,aGEOy=dU(r,ϕ,λ)dy=∂U(r,ϕ,λ)∂rdrdy+∂U(r,ϕ,λ)∂ϕdϕdy+∂U(r,ϕ,λ)∂λdλdy,aGEOz=dU(r,ϕ,λ)dz=∂U(r,ϕ,λ)∂rdrdz+∂U(r,ϕ,λ)∂ϕdϕdz+∂U(r,ϕ,λ)∂λdλdz,
where
(5)∂U(r,ϕ,λ)∂r=-GMR[1+∑n=2∞(n+1)(Rr)n+2nnnnnnnln×∑m=0nP-nm(sinϕ)(C-nmcos⁡mλ+S-nmsinmλ)],∂U(r,ϕ,λ)∂ϕ=GMR∑n=2∞(Rr)n×∑m=0ncos⁡ϕd(P-nm(sinϕ))d(sinϕ)(C-nmcos⁡mλ+S-nmsinmλ),∂U(r,ϕ,λ)∂λ=GMR∑n=2∞(Rr)n×∑m=0nmP-nm(sinϕ)(-C-nmsinmλ+S-nmcos⁡mλ).
Variables r, ϕ, and λ are related to object ECEF radius vector components x, y, and z by
(6)r=x2+y2+z2;sinφ=zr;sinλ=yρ;ρ=x2+y2.
According to (6),
(7)drdx=xr;drdy=yr;drdz=zr;dϕdx=-zxr2ρ;dϕdy=-zyr2ρ;dϕdz=ρr2;dλdx=-yρ2;dλdy=xρ2;dλdz=0.
Projections of Earth gravitation force in ECEF are
(8)aGEOx=∂U(r,ϕ,λ)∂rxr-∂U(r,ϕ,λ)∂ϕxzρr2-∂U(r,ϕ,λ)∂λyρ2,aGEOy=∂U(r,ϕ,λ)∂ryr-∂U(r,ϕ,λ)∂ϕyzρr2+∂U(r,ϕ,λ)∂λxρ2,aGEOz=∂U(r,ϕ,λ)∂rzr+∂U(r,ϕ,λ)∂ϕρr2.
The following recurrence equations can be applied to cos⁡mλ and sinmλ calculation:
(9)cos⁡((m+1)λ)=cos⁡mλcos⁡λ-sinmλsinλ,sin((m+1)λ)=cos⁡mλsinλ+sinmλcos⁡λ.


### 3.2.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, theN-body accelerations can be expressed as [4]
(10)aN-body=∑iGMi[ri|ri|3-ri-r|ri-r|3],
where G is the universal gravitational constant, Mi mass of the ith perturbing body (Sun or Moon), ri position vector of the ith perturbing body in ECEF, ri-r position vector of the ith perturbing body with respect to the object mass center in ECEF, and i planet index, where i=S for the Sun and i=M for the Moon.The values of the Sun and Moon position vectorsri can be obtained from the following equations:
(11)rS=ΩCrSe,rM=ΩCrMe,
where Ω is a transfer matrix from current Equatorial Earth Centered Inertial Frame to ECEF defined as
(12)Ω=(cos⁡ΩetsinΩet0-sinΩetcos⁡Ωet0001),
with Ωe the angular velocity of Earth and t time in seconds from the beginning of current sidereal day defined as
(13)t=0.997269566329084·GSMT
(Greenwich Sidereal Mean Time (GSMT) (see [6] for details)); C 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:
(14)C=(1000cos⁡ε-sinε0sinεcos⁡ε),
with ε the mean obliquity of the ecliptic as defined in [6]; rSe is radius vector of Sun mass center in Ecliptic Earth Centered Inertial Frame of J2000.0 defined as
(15)rSe=ρS(cos⁡λSsinλS0)ρS mean distance between Earth and Sun mass centers
(16)ρS=AU(1.00014-0.01673cos⁡l′-0.00014*cos⁡2l′)AU: astronomical unit, hire and below AU=1.49597870e11m; l′: is the mean anomaly of the Sun, see [6, 7] for detail; λS: is the ecliptic longitude of the Sun defined as:
(17)λS=LS+1.9171·sinl′+0.02·sin2l′+0.0003·sin3l′LS: the mean longitude of the Sun as defined below:
(18)LS=180+(100.46457166+35999.37244981·T)T: is Julian centuries from J2000.0; rMe: is radius vector of Moon mass center in Ecliptic Earth Centered Inertial Frame of J2000.0 which is defined as:
(19)rMe=ρM(cos⁡λMcos⁡βMsinλMcos⁡βMsinβM)ρM=1/Q: is the distance between Earth and Moon mass centers defined as:
(20)Q=1+0.0545cos⁡l+0.0100*cos⁡(l-2D)+0.0082cos⁡2D+0.0030cos⁡2*l+0.0009cos⁡(l+2D)+0.0006cos⁡(l′-2D)+0.0004cos⁡(l+l′-2D)+0.0003cos⁡(l-l′)λM: is the Moon geocentric longitude. It can be defined as:
(21)λM=L+6.289sinl-1.274sin(l-2D)+0.658sin2D+0.214sin2l-0.186sinl′-0.114sin2F-0.059sin(2l-2D)-0.057sin(l+l′-2D)+0.053sin(l+2D)-0.046sin(l′-2D)+0.041sin(l-l′)-0.035sinD-0.031sin(l+l′)-0.015sin(2F-2D)+0.011sin(l-4D)L=Ω+F: is the mean longitude of the Moon; βM: is the Moon geocentric latitude as defined in the following equation:
(22)βM=5.128sinF+0.281sin(l+F)-0.278sin(F-l)-0.173sin(F-2D)+0.055sin(F+2D-l)-0.046sin(l+F-2D)+0.033sin(F+2D)+0.017sin(2l+F)l, l′, F, D, Ω: 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 velocityv in an atmosphere will experience drag force. The drag force acceleration can be modeled as [4]
(23)Fdragm=-12ρ(CdAm)vr|vr|,
where ρ is the atmospheric density, vr is the object velocity vector relative to the atmosphere, m is the mass of the object, Cd is the drag coefficient for the object, and A is the cross-sectional area of the main body perpendicular to vr.The parameterσx=CdA/m 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 σx=0.01m2/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(24)ρ=ρ1e(h-h1)/H;H=70000+0.075(h-h1);h=r-R,
where h1is the reference height, h1=500000m; ρ1 is the atmospheric density on reference height, ρ1=2.e-13kg/m3; h is the object height; H is the height scale of the uniform atmosphere. r is the distance from the Earth’s center of mass; R is the semimajor axis of the WGS 84 Ellipsoid.


## 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 functionU is defined as [5]
(2)U=GMr[1+∑n=2N∑m=0n(Rr)nP-nm(sinϕ)nnnnn×(C-nmcos⁡mλ+S-nmsinmλ)∑n=2N∑m=0n(Rr)nP-nm(sinϕ)],
where U is the gravitational potential function (m2/s2), GM is the Earth’s gravitational constant, hire and below  GM=3986004.418e8m3/s2, r is the distance from the Earth’s center of mass (m), R is the Semimajor axis of the WGS 84 Ellipsoid, hire and below R=3986004.418e8m, n and m are the degree and order, respectively; ϕ is the geocentric latitude; λ is the geocentric longitude; C-nm and S-nm are the normalized gravitational coefficients defined in EGM-96 model [5], P-nm(sinϕ) is the normalized associated Legendre function, Pnm(sinϕ) is the associated Legendre function, and Pn(sinϕ) is the Legendre polynomial. Consider
(3)P-nm(sinϕ)=[(n-m)!(2n+1)k(n+m)!]1/2Pnm(sinϕ),Pnm(sinϕ)=(cos⁡ϕ)mdmd(sinϕ)m[Pn(sinϕ)],Pn(sinϕ)=12nn!dnd(sinϕ)n(sin2ϕ-1)n.In (2) and (3), k=1 for m=0 and k=2 for m≠0.The series is theoretically valid forr≥R.The acceleration due to geopotential can be defined as(4)aGEOx=dU(r,ϕ,λ)dx=∂U(r,ϕ,λ)∂rdrdx+∂U(r,ϕ,λ)∂ϕdϕdx+∂U(r,ϕ,λ)∂λdλdx,aGEOy=dU(r,ϕ,λ)dy=∂U(r,ϕ,λ)∂rdrdy+∂U(r,ϕ,λ)∂ϕdϕdy+∂U(r,ϕ,λ)∂λdλdy,aGEOz=dU(r,ϕ,λ)dz=∂U(r,ϕ,λ)∂rdrdz+∂U(r,ϕ,λ)∂ϕdϕdz+∂U(r,ϕ,λ)∂λdλdz,
where
(5)∂U(r,ϕ,λ)∂r=-GMR[1+∑n=2∞(n+1)(Rr)n+2nnnnnnnln×∑m=0nP-nm(sinϕ)(C-nmcos⁡mλ+S-nmsinmλ)],∂U(r,ϕ,λ)∂ϕ=GMR∑n=2∞(Rr)n×∑m=0ncos⁡ϕd(P-nm(sinϕ))d(sinϕ)(C-nmcos⁡mλ+S-nmsinmλ),∂U(r,ϕ,λ)∂λ=GMR∑n=2∞(Rr)n×∑m=0nmP-nm(sinϕ)(-C-nmsinmλ+S-nmcos⁡mλ).
Variables r, ϕ, and λ are related to object ECEF radius vector components x, y, and z by
(6)r=x2+y2+z2;sinφ=zr;sinλ=yρ;ρ=x2+y2.
According to (6),
(7)drdx=xr;drdy=yr;drdz=zr;dϕdx=-zxr2ρ;dϕdy=-zyr2ρ;dϕdz=ρr2;dλdx=-yρ2;dλdy=xρ2;dλdz=0.
Projections of Earth gravitation force in ECEF are
(8)aGEOx=∂U(r,ϕ,λ)∂rxr-∂U(r,ϕ,λ)∂ϕxzρr2-∂U(r,ϕ,λ)∂λyρ2,aGEOy=∂U(r,ϕ,λ)∂ryr-∂U(r,ϕ,λ)∂ϕyzρr2+∂U(r,ϕ,λ)∂λxρ2,aGEOz=∂U(r,ϕ,λ)∂rzr+∂U(r,ϕ,λ)∂ϕρr2.
The following recurrence equations can be applied to cos⁡mλ and sinmλ calculation:
(9)cos⁡((m+1)λ)=cos⁡mλcos⁡λ-sinmλsinλ,sin((m+1)λ)=cos⁡mλsinλ+sinmλcos⁡λ.


## 3.2.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, theN-body accelerations can be expressed as [4]
(10)aN-body=∑iGMi[ri|ri|3-ri-r|ri-r|3],
where G is the universal gravitational constant, Mi mass of the ith perturbing body (Sun or Moon), ri position vector of the ith perturbing body in ECEF, ri-r position vector of the ith perturbing body with respect to the object mass center in ECEF, and i planet index, where i=S for the Sun and i=M for the Moon.The values of the Sun and Moon position vectorsri can be obtained from the following equations:
(11)rS=ΩCrSe,rM=ΩCrMe,
where Ω is a transfer matrix from current Equatorial Earth Centered Inertial Frame to ECEF defined as
(12)Ω=(cos⁡ΩetsinΩet0-sinΩetcos⁡Ωet0001),
with Ωe the angular velocity of Earth and t time in seconds from the beginning of current sidereal day defined as
(13)t=0.997269566329084·GSMT
(Greenwich Sidereal Mean Time (GSMT) (see [6] for details)); C 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:
(14)C=(1000cos⁡ε-sinε0sinεcos⁡ε),
with ε the mean obliquity of the ecliptic as defined in [6]; rSe is radius vector of Sun mass center in Ecliptic Earth Centered Inertial Frame of J2000.0 defined as
(15)rSe=ρS(cos⁡λSsinλS0)ρS mean distance between Earth and Sun mass centers
(16)ρS=AU(1.00014-0.01673cos⁡l′-0.00014*cos⁡2l′)AU: astronomical unit, hire and below AU=1.49597870e11m; l′: is the mean anomaly of the Sun, see [6, 7] for detail; λS: is the ecliptic longitude of the Sun defined as:
(17)λS=LS+1.9171·sinl′+0.02·sin2l′+0.0003·sin3l′LS: the mean longitude of the Sun as defined below:
(18)LS=180+(100.46457166+35999.37244981·T)T: is Julian centuries from J2000.0; rMe: is radius vector of Moon mass center in Ecliptic Earth Centered Inertial Frame of J2000.0 which is defined as:
(19)rMe=ρM(cos⁡λMcos⁡βMsinλMcos⁡βMsinβM)ρM=1/Q: is the distance between Earth and Moon mass centers defined as:
(20)Q=1+0.0545cos⁡l+0.0100*cos⁡(l-2D)+0.0082cos⁡2D+0.0030cos⁡2*l+0.0009cos⁡(l+2D)+0.0006cos⁡(l′-2D)+0.0004cos⁡(l+l′-2D)+0.0003cos⁡(l-l′)λM: is the Moon geocentric longitude. It can be defined as:
(21)λM=L+6.289sinl-1.274sin(l-2D)+0.658sin2D+0.214sin2l-0.186sinl′-0.114sin2F-0.059sin(2l-2D)-0.057sin(l+l′-2D)+0.053sin(l+2D)-0.046sin(l′-2D)+0.041sin(l-l′)-0.035sinD-0.031sin(l+l′)-0.015sin(2F-2D)+0.011sin(l-4D)L=Ω+F: is the mean longitude of the Moon; βM: is the Moon geocentric latitude as defined in the following equation:
(22)βM=5.128sinF+0.281sin(l+F)-0.278sin(F-l)-0.173sin(F-2D)+0.055sin(F+2D-l)-0.046sin(l+F-2D)+0.033sin(F+2D)+0.017sin(2l+F)l, l′, F, D, Ω: 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 velocityv in an atmosphere will experience drag force. The drag force acceleration can be modeled as [4]
(23)Fdragm=-12ρ(CdAm)vr|vr|,
where ρ is the atmospheric density, vr is the object velocity vector relative to the atmosphere, m is the mass of the object, Cd is the drag coefficient for the object, and A is the cross-sectional area of the main body perpendicular to vr.The parameterσx=CdA/m 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 σx=0.01m2/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(24)ρ=ρ1e(h-h1)/H;H=70000+0.075(h-h1);h=r-R,
where h1is the reference height, h1=500000m; ρ1 is the atmospheric density on reference height, ρ1=2.e-13kg/m3; h is the object height; H is the height scale of the uniform atmosphere. r is the distance from the Earth’s center of mass; R 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 Section2.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
(25)v~(tk+h)=∫tktk+hv˙(t)dt+v(tk),r~(tk+h)=∫tktk+hv(t)dt+r(tk),b~(tk+h)=∫tktk+hd(t)dt+b(tk),d~(tk+h)=d(tk),
where h is the integration step limited by hmax⁡ (see Figure 1 for details); v~(tk+h)={vx,vy,vz} is the predicted ECEF velocity vector of the object; r~(tk+h)={x,y,z} is the predicted ECEF radius vector of the object; b~(tk+h) is the predicted receiver clock bias; d~(tk+h) is the predicted receiver clock drift; v(tk),r(tk),b(tk),d(tk) define last object state.The object state derivatives vector is defined using dynamic motion model which is described in Section2.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(26)x˙=f(x,t),
where x=(vx,vy,vz,x,y,z,b,d)T is a state vector that includes the spacecraft position and velocity vectors and the receiver clock bias and drift.The propagation ofx~k+1 from the previous state for covariance matrix propagation is generated by the following reduced equation:
(27)x~k+1=fm(xk)=xk+∫tktk+hfr(x,t)dt=xk+hfr(x,t),
where h=tk+1-tk is the integration step, limited by hmax⁡ (see Figure 1 for detail), xk is the state vector from the previous step, and fr(x,t) is the reduced dynamic model of notion witch is defined as follows:
(28)fr(x,t)=(v˙xv˙yv˙zx˙y˙z˙b˙d˙)=(-GMr3x+Ωe2x+2Ωey˙-GMr3y+Ωe2y-2Ωex˙-GMr3zvxvyvzd0).


### 4.1. Orbit Determination Algorithm Description

The orbit determination algorithm applies an UKF to estimate the state vector which comprises 8 components:(i)
object velocityv^={v^x,v^y,v^z},(ii)
object positionr^={x^,y^,z^},(iii)
receiver clock biasb^,(iv)
receiver clock driftd^,(v)
x=[v^Tr^Tb^d^].The diagram of orbit determination algorithm is described as follows.The following process and measurement models can be established:(29a)xk+1=f(xk)+wk,(29b)zk=h(xk)+vk.The variables in the above equation will be described:xk is a system condition vector in the k moment, f(·) is unscented system model, wk is a dynamic mixed signal in the kmoment, zk is a measuring dynamic vector in the k moment, h(·) is a unscented system measuring model, and vk is dynamic measuring mixed signal in the k moment.The measurement vector is denoted byzk, and the process noise wk and the measurement noise vk are assumed to be zero-mean white noise. The process noise covariance matrix and the measurement noise covariance are given by Qj and Rj, respectively.The system error covariance matrixQj is as follows:
(30)E{wiwjT}={Qj,i=j,0,i≠j.The measuring error covariance matrixRj is as follows:
(31)E{vivjT}={Rj,i=j,0,i≠j.Here,wk and vk are independent and unrelated:
(32)E{wivjT}=0,∀i,j.


### 4.2. Unscented Kalman Filter Processing

(1) Producing Standard Points and Calculation Value. Consider(33)χ0=x^k-,χi=x^k-+((n+ε)Pk)iT,i=1,2,…,n,χi+n=x^k--((n+ε)Pk)iT,i=1,2,…,n,Wim=ε(n+ε),W0c=W0m+(1-ε12+ε2),Wim=Wic=12(n+ε).The parameter is a scaling parameter defined as(34)ε2=ε12(n+ε3)-n.The positive constantsεi, i=1,2,3 are used as parameter of the method, and a priori and a posteriori estimates of the state are denoted by x^k- and x^k.(2) Time Updating. Condition predicted value is
(35)(ςk-)i=f((χk)i).Condition predicted average is(36)x^k-=∑i=02nWim(ςk-)i.Covariance matrix is(37)Pk-=∑i=02nWic[(ςk-)i-x^k-][(ςk-)i-x^k-]T+Qk.(3) Observation Updating. Observation measurement predicted value is
(38)(zk-)i=h((ςk-)i).Observation measurement predicted average is(39)z^k-=∑i=02nWim(zk-)i.P
x
y and Pyy are updated as follows:
(40)Pxy=∑i=02nWic[(ςk-)i-x^k-][(ςk-)i-z^k-]T,Pyy=∑i=02nWic[(zk-)i-x^k-][(zk-)i-z^k-]T+Rk.(4) Calculating Kalman Gain Value. Consider
(41)Kk=PxyPyy-1.(5) Updating Estimated Value to Measurement Value. Consider
(42)x^k=x^k-+Kk(zk-z^k-).(6) Updating Condition Error Covariance Matrix. Consider
(43)Pk+1=Pk--KkPyyKkT.


### 4.3. Unscented Kalman Filter Flow Chart

The flow chart of the UKF is as shown in Figure2.Figure 2
The algorithm flow chart of UKF.The time update phase of the UKF includes the propagation of state vectorxk+1- from the previous object state and the state covariance matrix Pk+1, which is defined in detail in Section 2.1.The subsequent measurement update adjusts the state vectorxk+1- components and state covariance matrix Pk+1 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(44)ypi(x~k+1,tk+1)=ρ~k+1i+b~k+1+βk+1i,
where
(45)ρ~k+1i=(xGPSi-x~k+1)2+(yGPSi-y~k+1)2+(zGPSi-z~k+1)2
is the geometric range; x~k+1,y~k+1,z~k+1 are the positional states of the user object at reception time; xGPS, yGPS, and zGPS are the positional states of the ith GPS satellite according to item; b~k+1 is the receiver clock offset; and βk+1i 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 vectorx:
(46)Hpi(tk+1)=[000-(xGPSi-x~k)ρ~ki-(yGPSi-y~k)ρ~ki-(zGPSi-z~k)ρ~ki10].
The model for a GPS Doppler measurement is given by
(47)yvi(x~k+1i,tk+1)=(vk+1i-v~k+1)·rk+1i-r~k+1∥rk+1i-r~k+1∥+d~k+1-εk+1i,
where vk+1i, rk+1i are the ith GPS satellite velocity vector and radius vector according to item; v~k+1 is the user object velocity; r~k+1 is the coordinates of the user object at reception time; d~k+1 is the receiver clock drift; and εk+1i 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 vectorx:
(48)Hvi(tk+1)=[∂Di∂vx∂Di∂vy∂Di∂vz∂Di∂x∂Di∂y∂Di∂z01],
where
(49)∂D∂vx=-(xGPS-x)ρ,∂D∂vy=-(yGPS-y)ρ,∂D∂vz=-(zGPS-z)ρ,∂D∂x=-(vxGPS-vx)(1ρ-(xGPS-x)2ρ3)+(vyGPS-vy)(xGPS-x)(yGPS-y)ρ3+(vzGPS-vz)(xGPS-x)(zGPS-z)ρ3,∂D∂y=(vxGPS-vx)(xGPS-x)(yGPS-y)ρ3-(vyGPS-vy)(1ρ-(yGPS-y)2ρ3)+(vzGPS-vz)(yGPS-y)(zGPS-z)ρ3,∂D∂z=(vxGPS-vx)(xGPS-x)(zGPS-z)ρ3+(vyGPS-vy)(yGPS-y)(zGPS-z)ρ3-(vzGPS-vz)(1ρ-(zGPS-z)2ρ3),
where xGPS, yGPS, and zGPS represent rk+1i={xGPSi,yGPSi,zGPSi}; vxGPS, vyGPS, and vzGPS represent vk+1i={(vx)GPSi, (vy)GPSi, and (vz)GPSi}; x,y,andz represent r~k+1={x~k+1, y~k+1, z~k+1}; vx, vy, and vz represent v~k+1i={(v~x)k+1i, (v~y)k+1i, and (v~z)k+1i}; and D represents Di.If both pseudorange and Doppler measurement are used, the sensitivity matrix will be composed ofHp and Hv matrices in the following way:
(50)Hk+1=[Hp(tk+1)Hv(tk+1)],
where Hp and Hv are matrixes size of [Nk+1SV×8] defined as
(51)Hp=[Hp1⋮HpNSV];Hv=[Hv1⋮HvNSV].
The measurement residuals or innovations sequence is
(52)Δyk+1=[yp(xk+1,tk+1)yv(xk+1,tk+1)]-[y~p(x~k+1,tk+1)y~v(x~k+1,tk+1)],
where yp(xk+1,tk+1) and yv(xk+1,tk+1) are the measured pseudoranges and pseudorange rates which are computed according to section; y~p(x~k+1,tk+1) and y~v(x~k+1,tk+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:(53)Kk+1=P~k+1Hk+1T(Hk+1P~k+1Hk+1T+Rk+1)-1,x^k+1=x~k+1+Kk+1Δyk+1,P^k+1=(I-Kk+1Hk+1)P~k+1(I-Kk+1Hk+1)T+Kk+1Rk+1Kk+1T,
where Rk+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.


## 4.1. Orbit Determination Algorithm Description

The orbit determination algorithm applies an UKF to estimate the state vector which comprises 8 components:(i)
object velocityv^={v^x,v^y,v^z},(ii)
object positionr^={x^,y^,z^},(iii)
receiver clock biasb^,(iv)
receiver clock driftd^,(v)
x=[v^Tr^Tb^d^].The diagram of orbit determination algorithm is described as follows.The following process and measurement models can be established:(29a)xk+1=f(xk)+wk,(29b)zk=h(xk)+vk.The variables in the above equation will be described:xk is a system condition vector in the k moment, f(·) is unscented system model, wk is a dynamic mixed signal in the kmoment, zk is a measuring dynamic vector in the k moment, h(·) is a unscented system measuring model, and vk is dynamic measuring mixed signal in the k moment.The measurement vector is denoted byzk, and the process noise wk and the measurement noise vk are assumed to be zero-mean white noise. The process noise covariance matrix and the measurement noise covariance are given by Qj and Rj, respectively.The system error covariance matrixQj is as follows:
(30)E{wiwjT}={Qj,i=j,0,i≠j.The measuring error covariance matrixRj is as follows:
(31)E{vivjT}={Rj,i=j,0,i≠j.Here,wk and vk are independent and unrelated:
(32)E{wivjT}=0,∀i,j.


## 4.2. Unscented Kalman Filter Processing

(1) Producing Standard Points and Calculation Value. Consider(33)χ0=x^k-,χi=x^k-+((n+ε)Pk)iT,i=1,2,…,n,χi+n=x^k--((n+ε)Pk)iT,i=1,2,…,n,Wim=ε(n+ε),W0c=W0m+(1-ε12+ε2),Wim=Wic=12(n+ε).The parameter is a scaling parameter defined as(34)ε2=ε12(n+ε3)-n.The positive constantsεi, i=1,2,3 are used as parameter of the method, and a priori and a posteriori estimates of the state are denoted by x^k- and x^k.(2) Time Updating. Condition predicted value is
(35)(ςk-)i=f((χk)i).Condition predicted average is(36)x^k-=∑i=02nWim(ςk-)i.Covariance matrix is(37)Pk-=∑i=02nWic[(ςk-)i-x^k-][(ςk-)i-x^k-]T+Qk.(3) Observation Updating. Observation measurement predicted value is
(38)(zk-)i=h((ςk-)i).Observation measurement predicted average is(39)z^k-=∑i=02nWim(zk-)i.P
x
y and Pyy are updated as follows:
(40)Pxy=∑i=02nWic[(ςk-)i-x^k-][(ςk-)i-z^k-]T,Pyy=∑i=02nWic[(zk-)i-x^k-][(zk-)i-z^k-]T+Rk.(4) Calculating Kalman Gain Value. Consider
(41)Kk=PxyPyy-1.(5) Updating Estimated Value to Measurement Value. Consider
(42)x^k=x^k-+Kk(zk-z^k-).(6) Updating Condition Error Covariance Matrix. Consider
(43)Pk+1=Pk--KkPyyKkT.


## 4.3. Unscented Kalman Filter Flow Chart

The flow chart of the UKF is as shown in Figure2.Figure 2
The algorithm flow chart of UKF.The time update phase of the UKF includes the propagation of state vectorxk+1- from the previous object state and the state covariance matrix Pk+1, which is defined in detail in Section 2.1.The subsequent measurement update adjusts the state vectorxk+1- components and state covariance matrix Pk+1 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(44)ypi(x~k+1,tk+1)=ρ~k+1i+b~k+1+βk+1i,
where
(45)ρ~k+1i=(xGPSi-x~k+1)2+(yGPSi-y~k+1)2+(zGPSi-z~k+1)2
is the geometric range; x~k+1,y~k+1,z~k+1 are the positional states of the user object at reception time; xGPS, yGPS, and zGPS are the positional states of the ith GPS satellite according to item; b~k+1 is the receiver clock offset; and βk+1i 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 vectorx:
(46)Hpi(tk+1)=[000-(xGPSi-x~k)ρ~ki-(yGPSi-y~k)ρ~ki-(zGPSi-z~k)ρ~ki10].
The model for a GPS Doppler measurement is given by
(47)yvi(x~k+1i,tk+1)=(vk+1i-v~k+1)·rk+1i-r~k+1∥rk+1i-r~k+1∥+d~k+1-εk+1i,
where vk+1i, rk+1i are the ith GPS satellite velocity vector and radius vector according to item; v~k+1 is the user object velocity; r~k+1 is the coordinates of the user object at reception time; d~k+1 is the receiver clock drift; and εk+1i 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 vectorx:
(48)Hvi(tk+1)=[∂Di∂vx∂Di∂vy∂Di∂vz∂Di∂x∂Di∂y∂Di∂z01],
where
(49)∂D∂vx=-(xGPS-x)ρ,∂D∂vy=-(yGPS-y)ρ,∂D∂vz=-(zGPS-z)ρ,∂D∂x=-(vxGPS-vx)(1ρ-(xGPS-x)2ρ3)+(vyGPS-vy)(xGPS-x)(yGPS-y)ρ3+(vzGPS-vz)(xGPS-x)(zGPS-z)ρ3,∂D∂y=(vxGPS-vx)(xGPS-x)(yGPS-y)ρ3-(vyGPS-vy)(1ρ-(yGPS-y)2ρ3)+(vzGPS-vz)(yGPS-y)(zGPS-z)ρ3,∂D∂z=(vxGPS-vx)(xGPS-x)(zGPS-z)ρ3+(vyGPS-vy)(yGPS-y)(zGPS-z)ρ3-(vzGPS-vz)(1ρ-(zGPS-z)2ρ3),
where xGPS, yGPS, and zGPS represent rk+1i={xGPSi,yGPSi,zGPSi}; vxGPS, vyGPS, and vzGPS represent vk+1i={(vx)GPSi, (vy)GPSi, and (vz)GPSi}; x,y,andz represent r~k+1={x~k+1, y~k+1, z~k+1}; vx, vy, and vz represent v~k+1i={(v~x)k+1i, (v~y)k+1i, and (v~z)k+1i}; and D represents Di.If both pseudorange and Doppler measurement are used, the sensitivity matrix will be composed ofHp and Hv matrices in the following way:
(50)Hk+1=[Hp(tk+1)Hv(tk+1)],
where Hp and Hv are matrixes size of [Nk+1SV×8] defined as
(51)Hp=[Hp1⋮HpNSV];Hv=[Hv1⋮HvNSV].
The measurement residuals or innovations sequence is
(52)Δyk+1=[yp(xk+1,tk+1)yv(xk+1,tk+1)]-[y~p(x~k+1,tk+1)y~v(x~k+1,tk+1)],
where yp(xk+1,tk+1) and yv(xk+1,tk+1) are the measured pseudoranges and pseudorange rates which are computed according to section; y~p(x~k+1,tk+1) and y~v(x~k+1,tk+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:(53)Kk+1=P~k+1Hk+1T(Hk+1P~k+1Hk+1T+Rk+1)-1,x^k+1=x~k+1+Kk+1Δyk+1,P^k+1=(I-Kk+1Hk+1)P~k+1(I-Kk+1Hk+1)T+Kk+1Rk+1Kk+1T,
where Rk+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.


## 5. GPS Measurement Models

The basic measurement types that will be employed in the current project are GPS pseudorange and Doppler inL1 frequency. The equation of the pseudorange in L1 frequency is given by
(54)Pki=ρki+Iki+c[dtGPS(tk)-dtU(tk)]-εk,
where Pki is the pseudorange in L1, ρki is the geometric range to the ith satellite at the observation instant tk is given by
(55)ρki=(xGPSi-x(tk))2+(yGPSi-y(tk))2+(zGPSi-z(tk))2,Iki is the ionospheric delay, c is the vacuum speed of light, dtGPS(tk) is the GPS satellite clock offset, dtU(tk) is the receiver clock offset, tk is the observation instant in GPS time, and εk 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 nominalL1 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 L1 frequency is given by
(56)Dki=-(vki-vuc·LOSki)L1,
where vki is the ith GPS satellite velocity at the observation instant tk, vu is the receiver velocity, LOSki is the line-of-sight to the ith GPS satellite at tk, and L1=1575.42 MHz is the transmitted frequency.The Doppler can be converted to a pseudorange rate observation given by(57)ρ˙ki=(vki-vu)·rki-ru∥rki-ru∥+fk+εki,
where f is the receiver clock drift in m/s and εki is the error in observation in m/s.Another possible GPS measurement model is a linear combination of GPSL1 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 propagationTmax⁡prp=2400 seconds (it is specified in 0).The maximal integration time stephmax⁡=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.Table 2
h, sec
Altitude, km
500
1200
d
R, m
d
V, m/sec
d
R, m
d
V, m/sec
1
6.00
×
10
-
7
8.00
×
10
-
8
4.00
×
10
-
7
4.00
×
10
-
8
10
0.009
9.00
×
10
-
6
5.00
×
10
-
3
4.50
×
10
-
6
20
0.16
1.60
×
10
-
4
0.09
8.50
×
10
-
5
30
0.9
9.00
×
10
-
4
0.5
4.50
×
10
-
4
50
9
9.00
×
10
-
3
5
4.50
×
10
-
3
100
150
0.15
80
0.075


### 6.2. Dynamic Model Settings

The maximal half interval of multibody acceleration fixingτmax⁡mb=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.Table 3
Altitude, km
Half interval of fixing
10
30
60
300
d
R, m
d
V, m/sec
d
R, m
d
V, m/sec
d
R, m
d
V, m/sec
d
R, m
d
V, m/sec
500
0.169
0.000185
0.551
0.000602
1.168
0.00128
5.806
0.00635
800
0.189
0.000193
0.615
0.000631
1.304
0.00134
6.489
0.00666
1200
0.216
0.000204
0.706
0.000668
1.497
0.00142
7.454
0.00705The fixed multibody acceleration components are available on time intervaltfix±τmax⁡mb, where tfix is time of acceleration fixing.


### 6.3. Estimation Settings

The minimal count of available Navigation SVNmin⁡SV=2 (it is defined by test results).The discrete state-noise covariance matrixQ, the discrete measurement noise covariance R, and the initial covariance matrix of P0 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:(58)P0=[0.25000000000.25000000000.25000000000.25×102000000000.25×102000000000.25×1020000000010600000000103],Q=[2×10-7000000002×10-7000000002×10-7000000006×10-2000000006×10-2000000006×10-2000000005×10-1000000005×10-3],where Rpr=55.9mis the pseudorange measurement dispersion and Rd=0.037m2/sec2 is the pseudorange rate measurement dispersion.


## 6.1. Integration Settings

The maximal time of continuous propagationTmax⁡prp=2400 seconds (it is specified in 0).The maximal integration time stephmax⁡=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.Table 2
h, sec
Altitude, km
500
1200
d
R, m
d
V, m/sec
d
R, m
d
V, m/sec
1
6.00
×
10
-
7
8.00
×
10
-
8
4.00
×
10
-
7
4.00
×
10
-
8
10
0.009
9.00
×
10
-
6
5.00
×
10
-
3
4.50
×
10
-
6
20
0.16
1.60
×
10
-
4
0.09
8.50
×
10
-
5
30
0.9
9.00
×
10
-
4
0.5
4.50
×
10
-
4
50
9
9.00
×
10
-
3
5
4.50
×
10
-
3
100
150
0.15
80
0.075


## 6.2. Dynamic Model Settings

The maximal half interval of multibody acceleration fixingτmax⁡mb=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.Table 3
Altitude, km
Half interval of fixing
10
30
60
300
d
R, m
d
V, m/sec
d
R, m
d
V, m/sec
d
R, m
d
V, m/sec
d
R, m
d
V, m/sec
500
0.169
0.000185
0.551
0.000602
1.168
0.00128
5.806
0.00635
800
0.189
0.000193
0.615
0.000631
1.304
0.00134
6.489
0.00666
1200
0.216
0.000204
0.706
0.000668
1.497
0.00142
7.454
0.00705The fixed multibody acceleration components are available on time intervaltfix±τmax⁡mb, where tfix is time of acceleration fixing.


## 6.3. Estimation Settings

The minimal count of available Navigation SVNmin⁡SV=2 (it is defined by test results).The discrete state-noise covariance matrixQ, the discrete measurement noise covariance R, and the initial covariance matrix of P0 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:(58)P0=[0.25000000000.25000000000.25000000000.25×102000000000.25×102000000000.25×1020000000010600000000103],Q=[2×10-7000000002×10-7000000002×10-7000000006×10-2000000006×10-2000000006×10-2000000005×10-1000000005×10-3],where Rpr=55.9mis the pseudorange measurement dispersion and Rd=0.037m2/sec2 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 Figure2, 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 Figure3. The UKF of direct observation equation is used in the simulation.(a)
(b)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. 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.


## 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 Figure3. The UKF of direct observation equation is used in the simulation.(a)
(b)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. 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.



---

*Source: 903026-2014-05-08.xml*