Orbit determination modeling analysis by GPS including perturbations due to geopotential coefficients of high degree and order , solar radiation pressure and lunisolar attraction

The purpose of this paper was to analyze the modeling of an artificial satellite orbit, using signals of the GPS constellation and least squares algorithms as the method of estimation, with the aim of analyzing the performance of the orbit estimation process. One pursues to verify how differences of modeling can affect the final accuracy of orbit determination. To accomplish that, the following effects were considered: high degree and order for the geopotential coefficients; direct solar radiation pressure; and Sun-Moon attraction. The measurements were used to feed the batch least squares orbit determination process, in order to yield conclusive results about the orbit modeling issue. An application has been done, using GPS data of the TOPEX/Poseidon satellite, whose accurate ephemeris are available on the Internet. It is shown that from a poor but acceptable modeling up to all effects included, the accuracy can vary from 28 to 9 m in the long-period analysis.


INTRODUCTION
The problem of orbit determination consists essentially of estimating parameter values that completely specify the trajectory of a body in space, processing a set of information (measurements) from this body.Such observations can be collected through a tracking network on Earth or through sensors, like the GPS receiver onboard TOPEX/Poseidon (T/P).
The Global Positioning System (GPS) is a powerful and low cost means to allow computation of orbits for artificial Earth's satellites.The T/P satellite is an example of using this system for space positioning.
The orbit determination of artificial satellites is a nonlinear problem in which the disturbing forces are not easily modeled, like geopotential and direct solar radiation pressure.Through an onboard GPS receiver, it is possible to obtain measurements (pseudo-ranges) that can be used to estimate the state of the orbit.
Usually, the iterative improvement of the orbit parameters of a satellite is carried out using the least squares methods.On a simple way, the least squares estimation algorithms are based on the data equations that describe the linear relation between the residual measurements and the estimation parameters.In this work, the algorithm was implemented through orthogonalization using Givens rotations.

LEAST SQUARES METHODS
Parameters estimation aims at estimating things that are constant along the estimation process.It is necessary a set of measurements to mathematically shape the relationship between these measurements and the parameters or state to be estimated.
One of the most used parameter estimator is the least squares algorithm.Basically, the algorithm minimizes the cost function of the residuals squared (Kuga, 2005).The recursive least squares algorithm, when applied to parameters or state estimation, presents two advantages: avoids matrix inversion in the presence of uncorrelated measurement errors and requires smaller matrices, which means less need of memory storage.

Recursive Least Squares using sequential Givens rotations
The Givens rotations are used when it is fundamental to cancel specific elements of a matrix.Alternative formulations were developed, based on the QR factorization methods, to solve this deficiency.Using orthogonal transformations, the equation matrix of data can be transformed on an upper triangular form, to which the least squares solution is obtained by a simple back substitution.The aim of applying orthogonal transformations in matrices and vectors on the least squares problem is to substitute the matrix inversion for a stronger method, with less numerical errors.The Givens rotations (Givens, 1958) are a method to solve recursive least squares through orthogonal transformations (Silva, 2001).
In this procedure, a given matrix H becomes triangular by the multiplication of a series of orthogonal matrices.The full transformation generically can be given by: denotes the sequence of rotations made in order to eliminate the elements of the sub-diagonal on the i-th row of H, and Q is an orthogonal matrix.At each step, the orthogonalization of the H matrix is performed (producing a transformed measurement vector d and r), and the results are stored to the next set of measurements.At the end, the final solution is computed.See details in Silva (2001), and Montenbruck and Suarez (1984).

MODEL OF FORCES -DISTURBING EFFECTS CONSIDERED
The main disturbing forces that affect the orbit of an Earth's artificial satellite are: the non-uniform distribution of Earth's mass; ocean and terrestrial tides and the gravitational attraction of Sun and Moon.There are also the non-gravitational effects, such as: Earth's atmospheric drag; direct and reflected solar radiation pressure; electric drag; emissivity effects; relativistic effects and impacts of meteorites.
The disturbing effects are in general included according to the physical situation presented and to the accuracy that is intended for the orbit determination.

Perturbations due to Earth Gravitational Field
The Earth is not a perfect sphere with homogeneous mass distribution, and cannot be considered as a material point.Such irregularities disturb the orbit of an artificial satellite and the keplerian elements that describe the orbit do not stay constant.The potential function can be given by (Kaula, 1966): (2) where µ is Earth's gravitational constant; R T is Earth's radius; r is the spacecraft radial distance; φ is the geocentric latitude; λ is the longitude on Earth fixed coordinates system; C nm and S nm are the normalized harmonic spherical coefficients, of degree n and order m; P nm are the associated normalized Legendre functions.The constants µ, R T , C nm , and S nm determine a particular gravitational potential.

Perturbations due to Direct Solar Radiation Pressure
The solar radiation pressure is a force of non-gravitational origin that disturbs the motion of an artificial satellite.The way as the perturbations due to solar radiation pressure will affect the keplerian elements depends on the pressure model adopted (if it includes or not shadow, for example).Meanwhile, in the general case, it causes secular and periodic perturbations on the variables (Ω, ω, and M) and periodic perturbations on the metric variables (a, e, and i).
The components of solar radiation pressure force can be expressed in several systems.Throughout these systems, the orbital elements of the satellite can be connected with Sun's position.This procedure was used here, for the direct solar radiation pressure model adopted for the TOPEX/Poseidon satellite (Marshall, Antreasian, Rosborough and Putney, 1991).Since the force due to the emerging radiation flux on the surface of the satellite depends on the angle of incidence, the attitude of the satellite must be also taken into account (Antreasian and Rosborough, 1992).
According to Marshall and Luthcke's model (Marshall and Luthcke, 1994), the total force acting on T/P is: where G is solar radiant flux (W/m 2 ); A is the surface area of each plate (m 2 ); δ is difusive reflectivity, percentage of the total incoming radiation; ρ is specular reflectivity, percentage of the total incoming radiation; n is surface normal vector; s is source incidence vector; θ is the angle between surface normal and solar incidence; and c is the speed of light (m/s).Subscript k varies from 1 to 8, representing each plate, and F is the total direct solar radiation force acting on the satellite.

Perturbations due to Luni-Solar Gravitational Attraction
These perturbations are due to Sun and Moon attraction force and they can be meaningful if the satellite is far from the Earth.As the orbital variations are of the same type, whatever is the Sun or the Moon the attractive body, they should be studied without distinguishing the third body.The luni-solar gravitational attraction mainly acts on Ω and ω, what causes precession of the orbit and the orbital plane.
The general three-body problem model is here simplified by the circular restricted three-body problem, where the orbital motion of a third body, which mass can be neglected, around two other massive bodies is studied.The motion equation that provides the third body acceleration can be expressed as (Prado and Kuga, 2001): , r r r 23 3 2 , and r i i , , , 1 2 3 , is the i-th body distance to the center of mass of the system

DYNAMICAL MODEL
The problem of orbit determination is essentially nonlinear due to the fact that the orbital motion description is based in ordinary differential equations in the form (Kuga, 1989) (Eq.5): and written in the inertial reference frame.In Eq. 5, r corresponds to the position vector; v is the velocity vector; w is the white noise vector with zero mean and covariance Q; a represents the modeled perturbations vector; and b is the vector of the user clock deviation polynomial model coefficients.

MEASUREMENT MODEL
The nonlinear equation that represents the measurements is given by (Eq.6): where y k is the m-observations vector; h x k k is the state x k nonlinear function of dimension m; and v k is the measurement error vector of dimension m.
The pseudorange observations are a measurement between the GPS satellites and the receiver antenna.These measurements will be used in the orbit determination problem via GPS and they can be written as (Eq.7): which corresponds to the real pseudorange ρ i plus the relativistic corrections c dt dT i , the ionospheric and tropospheric deviations D ION and D TROP , and the noise v.

RESULTS
Here, the tests and the analysis from the algorithm developed to compute direct solar radiation pressure are presented.
On the analysis of direct solar radiation pressure is already included TOPEX GPS antenna location that, lately, consists of the influence of the satellite attitude motion in the orbit determination process (Binning, 1996).The algorithm was implemented through FORTRAN 77 language (Pardal, 2007).
To validate and to analyze the proposed method, real data from the T/P satellite were used.Position and velocity to be estimated were compared with T/P precise orbit ephemeris (POE), from the Jet Propulsion Laboratory (JPL) of NASA.
The test conditions considered pseudorange real data, collected by GPS receiver onboard TOPEX, on November 19, 1993.The tests covered the same day, for a short (2 hours) and a long (24 hours) period of orbit determination.
The force model included perturbations due to high order geopotential (50x50), with harmonic coefficients from JGM-2 model, due to direct solar radiation pressure (Vilhena de Moraes, Raimundo and Kuga, 2008), and due to luni-solar attraction.These model improvements were also investigated for other estimation techniques, aiming at the same orbit determination problem (Pardal, Kuga and Vilhena de Moraes, 2010).The measurement model considered ionospheric correction, although the accuracy on position and velocity is not meaningful (Chiaradia, Kuga and Prado, 2000).
The obtained data were evaluated through one parameter: error in position, which represents the difference between the POE/JPL reference and the estimated position components.Such parameter is given by: which is after translated to radial, normal, and transverse components of orbit fixed system.
First, only geopotential was considered for the mentioned periods of orbit determination.After, the direct solar radiation pressure force acting on TOPEX center of mass and the way how such force acts on satellite orbit determination were analyzed.And, finally, perturbations due to Sun-Moon gravitational attraction were added to the model of forces.
The direct solar radiation pressure analysis already includes TOPEX GPS antenna location, which is one of the steps to determine the direct solar radiation pressure effects.
The obtained data were after translated to RTN (radial, transverse, and normal) system.In this system, "R" points to the nadir direction, "N" is perpendicular to orbital plane, and "T" is orthogonal to "R" and "N", and is also the velocity component.Thus, it is possible to analyze what happens with the orbital system components, and with the orbit evolution too.This is better than carry out an analysis from an Earth referential, where is more difficult to interpret the physical situation.
Figure 1 shows the behavior of the error in position, given in meters, along time, given in seconds, considering only geopotential.In Fig. 2, the geopotential and direct solar radiation pressure effects were considered, shown in two different curves.In Fig. 3, the Sun-Moon gravitational attraction in the model of forces was also taken into account.In the legends of Figs. 1 to 3, "R" means radial component; "N", normal component; and "T", transverse component of orbit fixed system.The subscript "geo" means perturbations due to geopotential only, "srp", perturbations due to geopotential and direct solar radiation pressure, and "sm", the effects due to Sun-Moon attraction.Table 1 shows the maximum and minimum amplitudes of the curves from Figs. 1 to 3, for short period (2 hours) and long period (24 hours) of orbit determination.
long periods of orbit determination, when Sun-Moon effects are not considered.Table 1 also confirms the information of the graphics: solar radiation pressure has more meaningful effects in transverse component, and less effects in normal one.It was expected, because solar radiation pressure acts especially on along track component, which is here represented by transverse component, the one associated with velocity of the satellite.
As Table 1 shows, for the short period, the solar radiation pressure decreases up to 43% the radial component value and up to 16% the transverse one.For the long period, solar radiation pressure reduces up to 42% the radial component value and up to 30% the transverse one.In both cases, such perturbation does not act favorably on the normal component.However, the Sun-Moon attraction has more meaningful effects on normal component.For the short period, it decreases up to 91% the normal component value, and for the long period, it reduces up to 94%.
Figure 4 shows residuals of pseudorange evolution along time for the short 2-hour period (at left), and for the long 24-hour period (at right) of orbit determination on November 19, 1993.These results were obtained considering the complete model of forces, including perturbations from geopotential up to Sun-Moon gravitational attraction.The behavior showed next is the same for all the analyzed cases, as it is possible to verify in Table 2, which presents the residuals of pseudorange statistics, mean and standard deviation, for the three models of forces considered in this analysis, and for the short and the long periods of orbit determination.
As can be seen, the statistics of the residuals are similar in all cases.Figure 1, especially for long period graphics, shows a sinusoidal behavior of the errors in RTN coordinates, with a period near once per revolution of the satellite orbit (around 2 hours).Following verification of all know dynamic models, there may exist a residual signature in the orbit as a result of unmodeled accelerations, which come in many forms (CNES, 2007;Soyka and Davis, 2001).In the case of geopotential, the acceleration is due to truncation of geopotential field.The used model is JGM-2 50×50, while the full model is 70×70, which is computationally intensive and may cause numerical problems, due to the order of some terms (10 -127 ), which was not suitable in this work.
The same sinusoidal behavior appears in long period graphics of Fig. 2, in the curves corresponding to geopotential and direct solar radiation pressure.It is also result of unmodeled accelerations, but solar radiation pressure is responsible for another acceleration, caused by limitations in the modeling of solar pressure as a function of the satellite attitude, surface properties and space environment.So, in Fig. 2, there are two unmodeled accelerations in the curves associated within the two disturbing effects.
The previously detected sinusoidal behavior also appears in long period graphics of Fig. 3, in the curves corresponding to the model of forces up to direct solar radiation pressure and to the completed model including Sun-Moon gravitational attraction.It is also result of unmodeled accelerations explained before.

CONCLUSIONS
The principal aim of this paper was to determine the orbit of an artificial satellite, using signals of the GPS constellation and least squares algorithms using sequential Givens rotations as a method of estimation.The analysis period covered a short period (near once T/P period) and a long period orbit determination.Pseudorange measurements were corrected from ionospheric effects, although the accuracy on orbit determination was not expressive.Real time requirements were not present; nevertheless, it was appropriate to keep low computational cost, with accuracy enough to satellite positioning at 10 m level for one day.
The results were compared with real data from TOPEX POE/JPL (Precision Orbit Ephemeris/Jet Propulsion Laboratory), available on the Internet.For short period orbit determination, the magnitude of error in position varied from 5.1 to 3.7 m, and for long period, the magnitude varied from 27.7 to 9.3 m, according to the complexness increase of the model.As the numbers show, the model that includes direct solar radiation pressure decreases at most around 5% the precision in position.Geopotential, direct solar radiation pressure, and Sun-Moon gravitational attraction were taken into consideration and the analysis occurred without selective availability on the measurements of the signals.
Throughout the results, it was found that least squares method through sequential Givens rotations and positioning using GPS showed reliability and accuracy enough for artificial satellites orbit determination.

Figure 1 :
Figure 1: Errors in position, given in RNT coordinates, for 2 hours (on the left) and 24 hours (on the right), considering perturbations due to geopotential (geo).R: radial component; N: normal component; T: transverse component.!

Figure 3 :
Figure 3: Errors in position, given in RNT coordinates, for 2 hours (on the left) and 24 hours (on the right), considering perturbations due to geopotential and direct solar radiation pressure (srp), and Sun-Moon attraction (sm).R: radial component; N: normal component; T: transverse component.

Figure 4 :
Figure 4: Residuals of pseudorange for 2 hours (on the left)and 24 hours (on the right), considering perturbations due to geopotential and direct solar radiation pressure and Sun-Moon gravitational attraction.

Table 1 :
Maximum and minimum values of the errors in position, for the two cases of test, onNovember 19, 1993

Table 2 :
Table 1, it is possible to see that the minor amplitude variation occurs in radial component, and the higher, in normal component, such for short and Residuals of pseudorange statistics for the three studied cases