Investigation on possible gravity field extraction from satellite orbits.
Arsov, Kirco; Virtanen, Jenni
Finnish Geodetic Institute, FINLAND

Earth satellites orbits are to a great extent controlled by the Earth's gravity field. It is the biggest force contributor in the equations of motion. Any irregularities in the Earth's gravity field are directly mapped in the orbital shape. We may say then that satellite orbit irregularities are the fingerprints of the irregularities in the Earth's gravity field. It is then very much clear that if one is able to find appropriate mathematical and physical relation between these orbit irregularities and gravity field parameters, the gravity field might be determined indirectly from these orbital perturbations on a global scale. In doing that, the spatial data distribution, data quality and parameters chosen control the gravity field extent we recover in terms of wavelength, for example if as parameters we choose the solid spherical harmonics. Theoretically, an infinite number of spherical harmonics coefficients will be needed to describe the full structure of the Earth's gravity field. However, in reality we deal with finite data sets which implies only limited information of the gravity field to be recovered. For example, the GOCE mission enables us to recover 100 km half-wavelength of the gravity field leading to ~ 200 - 250 degree spherical harmonics expansion (including gradiometry and SST data). Once having these parameters, its derivatives might be computed, for example gravity anomalies or geoid undulations etc .

So, the basic principle we follow is deriving physical quantities (spherical harmonics coefficients representation of the Earth's gravity field) from linear data (orbit discrepancies with respect to reference orbit). There are also other methods ex. energy integral approach, but here we stay on the traditional orbit/variational equations integration. This means that by proper integration of the satellite orbit by using best available dynamical models or empirical dynamical satellite data together with the variational equations, we are able to compute physical quantities (spherical harmonics coefficients) from these linear orbit discrepancies. This paper(presentation) presents continuation of author's work on ESA funded "From Eotvos to mgal+" project. We want to see here how far we may push the border in determining the Eart's gravity field from satellite orbital data (SST alone). We therefore implement some enhancements in the theory derived in the course of the abovementioned project following the today's IT improvements and developments. We recover the normalized spherical harmonics coefficients based on the orbit perturbations theory that will be elaborated and explained. For this we follow this methodology:

Nowadays there are satellites flying in pretty low orbits. These dedicated gravity satellite missions allow us to determine relatively high degree spherical harmonics coefficients on a global basis. Having in mind that the higher degree unnormalized spherical harmonics coefficients are relatively small numbers, as well as their correlation with other force parameters controlling the orbit and the height of the satellite, it might lead to certain numerical instabilities if one is using the unnormalized harmonics and goes to higher degrees and orders. Especially if one tries to estimate other force parameters together with spherical harmonics. Also, the already available orbital analysis software packages are developed over few decades, and understanding of the code, models used etc might be somehow difficult. Not in all of these packages normalized spherical harmonics are used as parameters. This is especially true for high flying satellites where the gravity field is used only up to some tens of degree. Only in the past decade some gravity dedicated missions were launched, such as CHAMP, GRACE and GOCE for example. But the orbital dynamics software packages have been developed over many decades as mentioned above. On the other hand force models used in orbital dynamics are in constant improvement and the intention of this study is also to implement the most up to date force models in the equations of motion. This is very important not only for the upcoming new dedicated gravity field missions but also for other satellite missions that do not have its main goal to be gravity field determination.

Based on all said above, the theoretical foundations of the numerical integration of variational equations together with the orbit will be presented and elaborated. Also, the formulation of the observational equations with accuracy study of this method will be presented together with some useful overall suggestions. Where necessary, an explanation on possible optimization will be given, since the method under certain scenarios could be very slow.

Few test runs will be computed and presented furthermore, dealing with the numerical accuracy, stability of the system, possible errors affecting the overall solution, linearization errors etc. All computations presented here are done with our own software ARCSST. There are two versions of ARCSST. One is written in FORTRAN90, and in all real number computations uses 16 significant digits 8 bytes double precision and works only on 32 bit systems. The other version is written in C++ and uses for all real numbers either 16 or 32 significant digits (8 or 16 bytes) based on parameter file settings. It is also implemented to work on both 32 and 64 bit operating system. Both ARCSST versions are freely available from the author to interested readers for scientific use.

So, we will assess first the numerical integrator accuracy of ARCSST and will present that it is for 365 days integration within few milimeters accurate with respect to analytical keplerian orbit. For 24h long arcs we will be good within few nanometers. This is important for the upcoming interferometric satellite missions that require relatively good integrator accuracy.
We will also show some computations with some generated orbits with 32 and 16 significant digits and elaborate how iteratively we may push the limits of parameters estimation as compared with todays max SST degree 140. We will present results based on pure iterative SST up to degree and order 200. Also, some studies on optimal arc length data gaps and normal equations matrix stability will be performed and explained.