Chapter 2 Methods for gravity field modeling using GPS observations
2.3 Theory of residual acceleration approach
A residual acceleration method is employed to determine the time variation of the Earth’s gravity field. In this method, the accelerations of LEOs are determined by numerical differentiations of the positions of LEOs. All perturbing forces caused by the static gravity field, Earth's non-sphericity, N-body, solid Earth tide, ocean tide, air drag, solar radiation pressure, Earth radiation and relativity are modeled first. After removing accelerations other than the Earth’s gravity-induced accelerations, linear relations between LEO accelerations and gravity coefficients can be established.
Empirical parameters can be used to model the residual non-gravitational accelerations.
The total gravitational force is the gradient of V. The acceleration vector expressed in the Earth-fixed coordinate system is (GSFC 1989)
⎥⎥ coordinates. The partial derivatives of the non-spherical portion of the Earth’s potential with respect to r,φ , and λ are given by
∑
∞∑
(2-32) through eq. (2-34) into eq. (2-31) and transforming the Earth-fixed coordinates (xb , yb , zb) into the spherical coordinates (r , φ , λ), eq. (2-31) can be rewritten as which are the radial, latitudinal, and longitudinal accelerations, respectively. M is an orthogonal matrix. The transformation of eq. (2-35) can be simplified by neglecting precession, nutation, and polar motion to obtain the acceleration vector aI in the inertial coordinate system (Hwang and Lin 1998):
⎥⎥ latitude are calculated as follows:
−
where GAST is Greenwich apparent sidereal time.
Chapter 3
Force modeling and precise orbit determination for COSMIC
3.1 Introduction
Precise orbits of a satellite position are important for positioning the satellite and for estimation of the Earth’s gravitation. Owing to the development of GPS, the spacecraft equipped with hl-SST receiver can collect position data continuously. The dynamic method (using force models) and the kinematic method (not using force models) (Švehla and Rothacher 2003; Jäggi et al. 2006 and 2007; Hwang et al. 2009) are two popular methods using GPS data applied for POD of LEOs. In addition, the so-called reduced-dynamic orbit determination, which is a compromising method between the dynamic method and the kinematic method, requires simplified force models. All above procedures for POD require GPS satellite ephemeris, Earth rotation information, and LEO GPS observation data as input for processing.
In Section 3.2, we focus on the perturbation force models acting on a COSMIC spacecraft. Without an accelerometer on the COSMIC spacecraft, the non- gravitational accelerations due to atmospheric drag and solar radiation need to be modeled. Section 3.3 and 3.4 describe the kinematic and dynamic orbit determination methods, complete with the principles, computation procedures and some analysis of results. The procedure of observed orbit data compression applied in gravity recovery is presented in the last part of this chapter.
3.2 Orbit dynamics of COSMIC satellite
The perturbing forces (accelerations) can be classified into gravitational forces and surface forces (or non-gravitational forces). The gravitational forces include the
Earth’s non-sphericity, N-body, solid Earth tide, ocean tide, and relativistic effect, and the surface forces include atmospheric drag, solar radiation pressure and the Earth’s radiation pressure. General accelerations, or empirical accelerations, are used to absorb the mis-modeled and un-modeled gravitational and surface forces. The algorithms of N-body, solid Earth tide, ocean tide, and relativistic effects can be found in a standard textbook of orbit dynamics such as Seeber (2003), and they will not be elaborated here. In this section, we focus on the parameters related to Earth’s non-sphericity, atmospheric drag and solar radiation pressure.
3.2.1 Equations of motion and perturbing potential force
In a geocentric inertial rectangular coordinate system, the equations of motion of an artificial Earth satellite such a COSMIC spacecraft can be expressed as (Long et al. 1989; Montenbruck and Gill 2001; Seeber 2003)
Pert
ns a
a r
r =− 3 + +
r
&& GM (3-1)
where
r: vector of satellite coordinates in the inertial frame r&&: acceleration vector
GM: Earth’s gravitational constant
ans: acceleration due to Earth’s non-sphericity
aPert: accelerations due to other perturbing forces
The first term in eq. (3-1) is called the point mass effect of the Earth given by Newton’s law of gravity and is 1000 times larger than any other acceleration. Eq.
(3-1) contains a system of second order differential equations which can be integrated to obtain the satellite positions and velocity at any epoch giving the initial state vector. The direct integration known as Cowell’s method (Balmino 1989) is selected for the simplicity and capacity for incorporating additional perturbations easily (Pavlis et al. 1996). The accelerations in eq. (3-1) are associated with certain parameters and can be adopted from existing values or estimated by satellite tracking data.
The acceleration aearth due to the geopotential is the gradient vector of the geopotential:
The first term in eq. (3-2) corresponds to the potential of a spherically symmetric Earth and the potential can be considered static. In fact, due to mass re-distribution in the Earth system, the geopotential is time-dependent so that the time-varying part should be considered. This is equivalent to dividing the coefficients in eq. (2-5) into a static part and a time-varying part as
)
where t is time. Thus, by recovering temporal gravity we mean estimating the
time-varying coefficients(ΔCnm(t),ΔSnm(t)) with respect to a mean gravity field.
3.2.2 Atmospheric drag and solar radiation effects on COSMIC satellites The acceleration vector of a LEO due to atmospheric drag is given by
) r r ( r r
adrag =− &−&d &−&d m
CDρ Ad 2
1 (3-4)
where CD is drag coefficient depending on the LEO shape and atmospheric composition, ρ is atmospheric density at the LEO position, Ad is effective (cross-sectional) area and m is the mass of the LEO, r r& &, d are the velocity vectors of
the LEO and atmosphere in the inertial frame, and (r& −r&d) is the velocity vector of the LEO relative to atmosphere. For each of the six COSMIC spacecrafts, the mass (with full thrust fuel) has been determined in the chamber test before the launch (April 2006). The remaining thrust fuel during the flight is observed and is used to adjust the time-dependent mass after the launch (Hwang et al. 2006 and 2009). The effective area Ad is the projected area of the area of the satellite in the flight direction onto a plane perpendicular to the direction(r && −rd). A COSMIC spacecraft travels in a manner that the POD+X antenna points to the flight direction (Fig. 3-1). Therefore, the total area in the flight direction is
panel main
T A A
A = +
(3-5)
where Amain is the area of the main body and Apanel is the area of two solar panels, which are computed as
)
where θ is the rotational angle of the solar panel (Fig. 3-1). Here we assume the thickness of the solar panels is negligible. The effective area was then computed by
)
The velocity vector of atmosphere was computed as
⎥⎥
where x and y are geocentric coordinate components of the LEO in the inertial frame, and ωh is the rotational velocity of atmosphere at an altitude of h computed as (King-Hele and Walker 1983)
)
The acceleration vector due to solar radiation pressure is
3
where ν is the eclipse factor, Ps is solar flux at one astronomical unit (au) (4.560 10× −6N/m2), Cr is reflectivity coefficient depending on the characteristics of the LEO, As is the effective area (different from the effective area for atmospheric drag) and r is the position vector of the Sun. The method to compute the effective s area for the solar radiation pressure is the same as that used in the atmospheric drag.
In this case, the effective area lies in a plane perpendicular to the vector (r− ). Crs r
can be expressed as (1+ε), where ε is reflectivity (from 0 to 1), which depends on the material of satellite parts. The eclipse factor depends on the position of LEO; ν=0 when the LEO is in the Earth’s shadow, and ν=1 when the LEO is illuminated by the Sun. The orbit dynamic modeling software we used is able to determine the eclipse factor in the cases of umbra and penumbra based on the ratio of the sunlight received at the LEO location, so that in practice the eclipse factor for a COSMIC satellite varies from zero to one.
Fig. 3-1: The dimensions of the main part and solar panels of a COSMIC LEO (top), velocity vector r& and LEO-to-atmosphere vector (r& −r&d) (Hwang et al. 2008)
3.3 Kinematic orbit determination using Bernese 5.0
The precise kinematic orbits of LEOs in this research were determined by the Bernese Version 5.0 GPS software (Dach et al. 2007). The reduced dynamic and kinematic approaches are available in Bernese 5.0 for POD with GPS observations.
The reduced dynamic approach estimates orbit arc-dependent parameters including the initial state vector (6 Keplerian elements), 9 solar radiation coefficients and three stochastic pulses in the radial, along-track and cross-track directions. The kinematic approach estimates the kinematic parameters of an orbit arc, including epoch coordinate components, receiver clock errors and phase ambiguities. Both the reduced dynamic and kinematic orbit determinations require high precision GPS satellite orbits and clocks. The GPS satellite precise orbits and high-rate clocks can be downloaded on the website provided by the Center for Orbit Determination in Europe (CODE, http://www.aiub.unibe.ch/igs.html). In the kinematic orbit determination with Bernese 5.0, the reduced dynamic orbit serves as a priori orbit for the kinematic orbit.
Fig. 3-2 shows the steps of precise kinematic orbit determination using real GPS data.
The zero-differenced ionosphere-free GPS measurements are usually used for point-wise calculation of the satellite positions in the kinematic approach. The limitations of orbit accuracy associated with kinematic orbits are based on the GPS satellite observation numbers and relative GPS-LEO geometry (Byun and Schutz 2001). Satellite coordinates are estimated together with one GPS receiver clock parameter every epoch. Comparing with SLR observations, the kinematic POD with accuracy of 1–3 cm was demonstrated for the GRACE mission (Švehla and Rothacher 2004). Using an overlapping analysis, the orbit accuracy of COSMIC is about 3 cm, compared to 1 cm in the case of GRACE satellites (Hwang et al. 2009). We find that the quality of GPS data depends on the quality of satellite attitude data. For the case
measurements over the equator and the polar regions are 0.5o and 3o , respectively, and are larger at a lower altitude. When the attitude of a satellite is poorly determined, the uncertainties in the estimated GPS phase ambiguities are relatively large, leading to degraded orbital accuracy. The detail of GPS-determined orbits of GRACE and COSMIC satellites using the kinematic approaches by Bernese have been documented by Švehla and Rothacher (2005), Jäggi et al. (2007) and Hwang et al. (2009).
Fig. 3-2: Steps of precise kinematic orbit determination using GPS data
3.4 Dynamic orbit determination using GEODYN II software
The dynamic orbit determination strategy for LEO based on GPS double- difference tracking data is demonstrated for the first time for TOPEX/Poseidon mission (Bertiger et al. 1994; Schutz et al. 1994). The equations of motion are solved by numerical integration. The dynamic model errors will lead to systematic errors growing with the arc length (Bock 2003).
In this research, we use NASA Goddard GEODYN II software to model the perturbing forces described in Section 3.2.1. GEODYN II is used extensively for satellite orbit determination, geodetic parameter estimation, tracking instrument calibration, satellite orbit prediction, as well as for other applied research in satellite geodesy using virtually all types of satellite tracking data (Pavlis et al. 1996). For direct numerical integration, GEODYN II uses Cowell’s summation method to obtain the position and velocity at epoch and uses the Bayesian’s least-squares method for parameter estimation. The mathematical models of the perturbing forces used in GEODYN II can be found in Long et al. (1989), Pavlis et al. (1996) and McCarthy (1996). GEODYN II has been used for precise orbit determination/prediction and force modeling in various Earth resource satellites such as TOPEX/Poseidon and GRACE. Temporal gravity fields from such satellite tracking data as SLR and GRACE KBR have been derived with GEODYN II (Cox and Chao 2002; Luthcke et al. 2006).
GEODYN II is divided into three major components: the Tracking Data Formatter (TDF), GEODYN IIS and GEODYN IIE. The flow of running GEODYN II can be found in Hwang (2002). The TDF program takes in the one of several tracking data forms. In this research, the tracking data format is PCE data format containing the information of satellite position and velocity as a priori orbit.
In preparation for the execution of GEODYN II, a file containing the ephemeris of the planets and a file containing A1UTC, polar motion, solar flux and magnetic flux must be made ready. A1UTC is the difference between the atomic time (A1) used in GEODYN II and the universal time (UTC), which is available from http://hpiers.obspm.fr/eop-pc/. Solar flux and magnetic flux are obtained from NOAA’s web site ftp://ftp.ngdc.noaa.gov under the directory STP/
GEOMAGNETIC_DATA/INDICE. These raw data are processed to produce binary files suitable for input to GEODYN IIS program.
The GEODYN IIS program is mainly used to read and process the option cards, the input observation data, optional gravity model, station geodetics, area/mass files, ephemeris and table data. The default gravity model is stored in the file ftn12. In this research, we choose GGM03S gravity models up to degree/order 70 derived from GRACE GPS and KBR observations.
The JPL export ephemeris as input file ftn01 is used in GEODYN IIS for nutations, positions, and velocities of the Moon, Sun, and planets. In this research, we use the JPL binary DE-403 ephemeris. GEODYN II interpolates for the ephemeris in the mean of 1950.0 reference system by Chebyshev interpolation. This step gives greater accuracy than interpolating in the reference system (the true of date coordinate system) because the high-frequency perturbations due to nutations are absent. After interpolation, the coordinates are then rotated to the true of date system using the precession and nutation matrices (Seeber 2003). More details on the transformations between different coordinate systems can be found in the textbooks or manuals like Long et al. (1989), Pavlis et al. (1996) and Seeber (2003).
The file ftn05 contains all option cards determining the force and non-force model parameters to be used in the program execution. These option cards are divided
2002). The Global Set consisting of four groups provides all of the common arcs processing information. The first group, Global Set Mandatory Cards, is the mandatory run description on other three cards. The second set of Global Set Option cards is used to define and/or estimate conditions which are common to all the arcs being processed. This group contains the information and estimations of the Earth's gravitational potential and/or new Earth constants, application and/or adjustment of time dependent gravity coefficients, dynamic polar motion, the third body gravitational potential and/or new constants, solid Earth tide model, ocean tide model, atmosphere drag model, solar flux, and tectonic plate motion. The third group is the Position Card Group containing the information of tracking stations. The last group is the Global Set Terminator to end the Global Set.
One or more Arc Set contains information defining its arc in ftn05 file. The Arc Set also has four groups in this order. The first group, Arc Set Mandatory Cards, is to decide the reference coordinate system and time and spacecraft parameters in this arc.
The second group, Arc Option Set cards, is specified to make use of GEODYN II's individual arc capabilities. The third group is the Data Selection/Deletion Subgroup used to edit input observations. The last group is the Arc Set Terminator to end the Arc Set.
The parameters in the Global Set for the force models of COSMIC given in Table 3-1 are defined by input option cards in the file ftn05. For surface forces in the Arc set, we solve for atmospheric drag coefficient, radiation coefficient and 9 empirical coefficients of general acceleration along the radial, along-track and cross-track directions every 1.5 hours (one orbital period) using COSMIC kinematic orbits. As an example, Fig. 3-3 shows the estimated atmospheric drag coefficients and reflectivity coefficients for FM5. These estimated coefficients vary over time, and the
and 1.23/0.30, respectively. The general accelerations for COSMIC (Table 3-1) at an altitude of 520 km are on the order of 10-11ms-2 . GEODYN IIE performs the computation of satellite orbit and geodetic parameter estimations. The output from IIE contains all necessary information for data analysis.
The main purpose of precise dynamic orbit determination is to generate a reference orbit to compute the residual orbit perturbations for gravity field recovery.
The residual orbit perturbation is a function of the perturbing force due to the perturbing geopotential. Like satellite position, satellite acceleration contains the effect due to the geopotential, plus other perturbing forces which must be modeled for gravity field recovery.
Table 3-1: Standards for the orbit dynamics of COSMIC satellites
Model/parameter Standard Conventional inertial
reference frame
J2000
N-body JPL DE-403
Earth gravity model GGM03S
Polar motion IERS standard 2000
Reference ellipsoid ae = 6378136.3 m, f =1/298.257
GM 396800.4415 km3s-2
Ocean tides GOT00.2 (Ray 1999)
Solid Earth tides IERS standard 2000
Atmosphere density
Mass Spectrometer Incoherent Scatter (MSIS) Empirical Drag Model (Hedin, 1991)
Earth radiation pressure
Second-degree zonal spherical harmonic model (Knocke et al. 1988)
Solar radiation pressure one coefficient every 1.5 hours Atmosphere drag one coefficient every 1.5 hours General accelerations 9 parameters every 1.5 hours
Fig. 3-3: Estimated atmospheric drag coefficients (top) and solar reflectivity
3.5 Normal-point reduction
The original sampling interval of COSMIC and GRACE GPS POD carrier-phase and code observables is 1 second. In practice, a 5-s (0.2 Hz) is used in the reduced-dynamic and kinematic orbit computations. To reduce noises and data volume, the 5-s kinematic orbits can be compressed and filtered at a greater item interval by an algorithm similar to that used in the normal-point reduction of satellite laser ranging. In this research, we adopt the Herstmonceux algorithm (Sinclaire 1997) to compress the COSMIC and GRACE precise orbits. Specifically, we use the following steps to generate normal-point kinematic orbits:
(1) Use the reduced-dynamic orbit as the reference orbit to generate differenced orbit. A differenced orbit component is
3 , 2 , 1 , =
−
=x x i
pi ik ir (3-11)
where x and ik x are components of kinematic and reduced-dynamic orbits, ir respectively.
(2) Remove large outliers in the kinematic orbit, which will not be used in the subsequent computations. An outlier is defined as pi ≥20cm.
(3) Within a bin (a window containing many differenced orbits), the differenced orbits are fitted by a polynomial in time using least-squares. The polynomial is called the trend function f(t)
(4) For each orbit component, compute the residuals at the times of observations as
) ( i
i
i p f t
v = − (3-12)
(5) Compute the root-mean-square value RMS of the residuals. Identify outliers using a rejection level of 2.5 times of RMS, and neglect these outliers in step (3) of the next iteration
(6) Repeat steps (3)-(5) until no outlier is found
(7) Divide the accepted residuals into bins starting from 0h UTC.
(8) Compute the mean value v and the mean time of the accepted residuals within m each bin. The number of accepted residuals within bin m is denoted asn . m
(9) For each orbit component, locate the kinematic orbit x and its residual mk v , m
whose observation time t is nearest to the mean time of the accepted residuals in m bin m.
(10) Compute the normal-point kinematic orbit as
m m k m
m x v v
NP = − + (3-13)
(11) Compute the standard error of normal points as (if nm =1, this bin is neglected)
) 1 (
1 2
=
∑
−m m
n j
m n n
v
m
σ (3-14)
The bin size can be adjusted according to the desired spatial resolution of gravity solution and data compression ratio. The degree of the fitted polynomial increases with the bin size. For a one-minute bin, a second-degree polynomial is found to be optimal. Statistically, the standard errors of normal points will be smaller than those of raw orbits. For example, Fig. 3-4 shows the normal-point residuals (differences between reduced-dynamic and kinematic orbits) are smaller than raw residuals in Y-direction of FM5 satellite in DOY 216, 2006.
Fig. 3-4: Raw and normal-point residuals in Y-direction (FM5, DOY216)
Chapter 4
Recovery of temporal gravity field using analytical orbital
Recovery of temporal gravity field using analytical orbital