Orbit determination in satellite geodesy

Available

Pergamon

SCIENCE

www.elsevier.com/locate/asr

ORBIT

online at www.sciencedirect.com DIRECT.

doi: lO.l016/SO273-1177(03)00171-6

DETERMINATION G. Beutlerl,

1Astronomical

Institute,

IN SATELLITE

T. Schildknechtl,

U. Hugentobler’,

GEODESY

W. Gurtner 1

University of Berne, Sidlerstrasse 5, 3012 Berne, Switzerland ABSTRACT

For centuries orbit determination in Celestial Mechanics was a synonym for the determination of six socalled Keplerian elements of the orbit of a minor planet or a comet based on a short series of (three or more) astrometric places observed from one or more observatories on the Earth’s surface. With the advent of the space age the problem changed considerably in several respects: (1) orbits have to be determined for a new class of celestial objects, namely for artificial Earth satellites; (2) new observation types, in particular topocentric distances and radial velocities, are available for the establishment of highly accurate satellite orbits; (3) even for comparatively short arcs (up to a few revolutions) the orbit model that has to be used is much more complicated than for comparable problems in the planetary system: in addition to the gravitational perturbations due to Moon and planets higher-order terms in the Earth’s gravity field have to be taken into account as well as non-gravitational effects like atmospheric drag and/or radiation pressure; (4) the parameter space is often of higher than the sixth dimension, because not only the six osculating elements referring to the initial epoch of an arc, but dynamical parameters defining the (a priori imperfectly known) force field have to be determined, as well. It may even be necessary to account for stochastic velocity changes. Orbit determination is not a well-known task in satellite geodesy. This is mainly due to the fact that orbit determination is often imbedded in a much more general parameter estimation problem, where other parameter types (referred to station positions, Earth rotation, atmosphere, etc.) have to be determined, as well. Three examples of “pure” orbit determination problems will be discussed subsequently: l

The first problem intends to optimize the observation process of one Satellite Laser Ranging (SLR) observatory. It is a filter problem, where the orbit is improved in real time with the goal to narrow down the so-called range-gate, defining the time interval when the echo of the LASER pulse is expected.

l

Secondly we highlight orbit determination procedures (in particular advanced orbit parametrization techniques) related to the determination of the orbits of GPS satellites and of Low Earth Orbiters (LEOs) equipped with GPS receivers.

We conclude by discussing the problem of determining the orbits of passive artificial space debris using high-precision astrometric CCD-observations of these object. 0 2003 COSPAR. Published by Elsevier Science Ltd. All rights reserved. l

satellites or of

INTRODUCTION The Satellite Orbit and its Parametrization The classical orbital parameters, called orbital elements, are illustrated in Figure 1. In satellite geodesy these elements are referred to the (mean) equatorial (X, Y)-plane of a standard epoch (usually 52000.0). The X-direction is the (mean) direction of the vernal equinox at the reference epoch. The right ascension R of the ascending node and the inclination i w.r.t. the equatorial plane define the orbital plane in the equatorial coordinate system. The perigee distance from the ascending node, called the argument of perigee w, defines the orientation of the conic section (the ellipse) in the orbital plane. The semi-major axis a and the eccentricity e define size and shape of the ellipse, the perigee passing time To initializes the orbital motion. SpaceRes. Vol. 31, No. 8, pp. 1853-1868.2003 Q 2003 COSPAR. Published by Elsevier Science Ltd. All rights Printed in Great Britain

Adv.

0273-1177103 $30.00+ 0.00

reserved

1854

G. Beutler

Fig. 1. The six classical

orbital

et al.

elements.

These elements, or a subset thereof, are the “final” orbit parameters to be determined in any orbit determination process. The attribute “final” expresses the fact that intermediary parameters may be better suited to characterize the problem. But eventually, the parameters of Figure 1 (or a modified version thereof), have to be determined. If the Earth is assumed to be the only attracting body, if its mass distribution is approximated as spherically symmetric, and if all non-gravitational forces are neglected, the six classical elements uniquely describe the satellite’s orbit, which is then a solution of the one-body version of the equations of motion: +=.-GMT

(1)

r3 ’

where G is the constant of gravitation, M the Earth’s mass, r the geocentric position vector of the satellite (and r its length), and i: its second time derivative. If the length of the satellite arc is short (a small fraction of the revolution period), the orbit model (1) may be sufficient to represent the observations within the arc with sufficient accuracy. If the arc length is growing, the equations of motion become more elaborate. Symbolically, we write them as follows:

i:=-GM

;+bf(t;r,i,pl,pz

,...,

pm),

(2)

where Sf) is the perturbing acceleration and pi, i = 1,2,. . . ,m are additional unknown parameters, socalled dynamical parameters of the orbit; we may think of them as parameters characterizing radiation pressure or atmospheric drag. Usually, one may assume that the absolute value of Sf (. . .) is much (typically by a factor of the order of 103) smaller than the absolute value y of the main term (which justifies their neglect for src lengths much shorter than one revolution period). If the satellite motion is governed by equations of motion of type (2), we still use the classical orbital elements of Figure 1 as the unknowns characterizing the initial values of the parameter estimation process, but we have to interpret them as osculating orbital eZements referring to a particular epoch to (usually selected as the initial epoch of the arc). The osculating elements referring to the epoch to are derived from the position- and velocity-vectors ~0 G r(tc) and +e A i(to) using the formulae of the two-body problem. The orbit to be determined may then be written as the solution of an initial value problem of the following type (as long as the orbit model is assumed to be purely deterministic)

i:

-GM

r(t0)

1

Wt0)&0j,.

+(to)

=

i.(dto),

q +

fif(t;r,+:,pl,~z,.

. . ,pm>

(3)

. . Jo(t0)) e(to),

. . . ,To(to))

.

Orbit

Fig. 2. The Zimmerwald

Observatory

(left)

Determination

in Satellite

and the Zimmerwald

Geodesy

l-m

Telescope

(right)

Observations in Satellite Geodesy Figure 2 shows the instruments of classical astrometry and of modern satellite geodesy. Figure 2 (left) gives an overview of the Zimmerwald observatory, Switzerland’s contribution to classical astrometry and to the modern global space geodetic networks. The dome to the left (west) contains the Schmidt telescope (on the same mount as the observatory’s Cassegrain telescope). More than 100 minor planets and five comets were discovered with this instrument in the second half of the twentieth century, among them comet Wild-2, the target of the US-American Stardust mission. The dome to the right is devoted to satellite laser ranging (SLR) and to CCD-astrometry (CCD=charge coupled device). The SLR-instrument is composed of an astronomical telescope of 1 m aperture and of a Titanium-Sapphire Laser with a repetition rate of 10 Hz, a pulse width of 100 ps (corresponding to a length of 3 cm of the individual Laser pulse), and two wavelengths of Xi = 423 nm and X2 = 846 nm. The CCD-observations are performed with the same telescope using a cooled 2k x 2k pixel CCD. The 10-m mast in Figure 2 (left) hosts the antenna of the observatory’s dual-band GPS-receiver. Zimmerwald is one of more than 300 global stations of the International GPS Service (IGS) network (IGS homepage http://igscb.jpl.nasa.gov) and one of about 30 stations of the International Laser Ranging Service (ILRS) (homepage http://ilrs.gsfc.nasa.gov). The new Zimmerwald telescope (Figure 2 (right)) was designed as a multipurpose instrument. The astrometric places of the artificial satellites to be analyzed subsequently were made with this telescope. The setup allows it, in principle, to observe directions and distances to.one and the same satellite within one and the same satellite pass with state-of-the-art accuracies. Figure 3 (left) shows that the astrometric place is the geometric direction of the observer at observation time t to the satellite at time t - At A t - $, where p is the distance between the observatory at time t and the satellite at time t - f. The observed functions (Y (right ascension) and b (declination) (i.e., the observables) are defined as follows: P

+

e

A

o!

= =

6

1 T(t q&$)-R(t)

f) P

-

R(t)

1

(4)

arctan 2 arcsine, ,

where e,, eY and e, are the components of the unit vector e pointing from the observer at observation time to the satellite at time t - f in the equatorial system. Atmospheric refraction effects were neglected in the above equations. Figure 3 (middle) explains the SLR observable Atslr, which is the light travel time Atslr G 5 + e from

1856

G. Beutler

et al.

Fig. 3. Topocentric astrometric place e(t) at observation time t of a satellite S at a distance p from an observer as R, and the geocentric position vector r(t - At) w.r.t. the geocenter (where At A z) (left); SLR observation light travel time At A 9 from observer at signal emission time tl to satellite at time (tl + Q/2 to observer at signal receptidn time t2 (center); GPS-observable as difference of clock readings t, of receiver clock at signal reception time and reading tS of satellite clock at signal emission time (right).

the observatory at pulse emission time tl to the satellite and back to the observatory at light reception time t2. With only minor neglects the SLR observable At, may be written as c A&r

= 2dhn)

G2

I d&n) - Wm) I ,

(5)

where t,,, G 9. Atmospheric refraction was not considered in the above treatment. (Subsequently the subscript sir will be omitted if no confusion is possible). The GPS-observable is illustrated by Figure 3 (right). It is the difference of the readings of two different clocks, namely the reading of the receiver clock at signal reception time t, and the reading ts of the satellite clock at signal emission time. The difference of the two clock readings, multiplied by the velocity of light, is called pseudorange. It may be written as the sum of the range between observer (at time tr) and the satellite (at time t,) and the two clock corrections (both multiplied by the speed of light c) w.r.t. a perfect time frame, in practice realized by the so-called system time: pA

p(tp) + c (At, - Ats) .

(6)

The atmospheric propagation delays were neglected in the above equation. The GPS-observable (6) and the SLR-observable (5) would be identical, if the clock errors could be neglected. Equation (6) governs the so-called GPS-code observation. Code observations are accurate to about 20 cm to 3 m depending on the method and the type of signal used. The so-called GPS-phase observation is much more accurate (of the order of a few mm), but an additional term NX (where N is an unknown and X the carrier wavelength), the so-called ambiguity term, has to be taken into account in addition to the two clock parameters: 4G

p(tr) + c (At, - At’) + NX .

(7)

For one triplet (satellite, station, carrier) the ambiguity term NX retains its value as long as the phase-lock can be maintained. N must be estimated as a real-valued number; because it is in reality an integer number, it may be fixed to the correct integer value under certain conditions (see Teunissen and Kleusberg (1998)).

Orbit

Determination

in Satellite

Geodesy

1857

Orbit Improvement and First Orbit Determination Orbit determination is a non-linear parameter estimation problem. If approximate values uK, eK, iK, flK , uK, Z’eK, and p:, i = 1,2 7”‘7 m are available for the (osculating) orbital elements and for the dynamical parameters, the problem may be linearized. In this case we speak of an orbit improvement process. Introducing the general notation

where (Y is the right ascension, 6 the declination, At,,, the light-travel time, p the code pseudorange, $ the corresponding phase pseudorange, and phi the phase observable, for the observed functions and

the linearized observation equations may be written in the form ‘gg(tk)

(qi

-4:)

-

(0:

-0:)

=Vk

, k=l,&:..:n,

2

whereoK_o(tk;q~,qZK,...,qgK++,)isth e ob served function at observation time tl, (which, according to the above definitions, has slightly different meanings for different observation techniques), OL is the observation itself, and Vk is the residual. The linear system of observation equations (10) is solved by the method of least squares, which minimizes a (possibly weighted) sum of the residual squares vk. 2 Obviously the number n of observation equations must be equal to or exceed the number of unknowns in order to guarantee that the system has a unique solution. This is a necessary, but not a sufficient condition. The observation equations approximate the non-linear observed functions oi (41, q2 , . . .) by Taylor series, which are truncated after the terms of first order. Due to these neglects the orbit improvement process (in principle) must be an iterative process, where the resulting parameters in the iteration step K are used as the (known) a priori values in the next step. The iteration process may be stopped, if the increments qi - q; are “infinitesimal”. In most cases only few iteration steps are required in practice. For more details and other termination criteria we refer to Beutler (2003). The partial derivatives of the observed functions w.r.t. the unknown parameters may, in essence, be written as the scalar product of the position gradient of the observed function and the partial derivative of the geocentric position vector r(t) w.r.t. the unknown parameters: Z(t)

= z

VoK(t).

g(t)

. 2

Due to the light travel times, small velocity dependent terms should be applied, as well (for details we refer to Beutler (2003)). The position gradient of the observed function (first vector of the product on the right-hand side of eqn. (11) depends on the observation type - its computation is a purely algebraic problem (see, e.g., Beutler (2003)). The second vector in this equation is obtained by taking the partial derivative of eqns. (3) w.r.t. the parameter qi considered. The initial value problem for the determination of the vector function F(t) may be written as:

(12)

The matrix Ao is the Jacobian of the acceleration vector w.r.t. the Cartesian components of the position vector of the satellite, Al the Jacobian of the same vector w.r.t. the Cartesian components of the velocity vector (which is the zero matrix if no velocity dependent accelerations occur). Obviously, the partial derivative of Sf w.r.t. qi is zero, if qi is one of the osculating elements, whereas the partial derivatives of the initial position and velocity vectors are zero w.r.t. the dynamical parameters. The differential equation system in eqns. (12) is called the system of variational equations. It is a linear system

1858

G. Beutler et al. Logeos-1:

Pulse Trovkl Times

Logeos

:. 20

I, 30

G.

1: Observed

,..

.. ’ I ,., 40

- Predicted

. ,

I 50

I.,.

L 61

Fig. 4. Observed (left), observed (ms) minus predicted (in ns, right) light-pulse travel times to Lageos 1, observed at Zimmetwald Observatory on July 7, 2002, 07 ’ 19m - 08h OOm (left); time on horizontal pxis in min relative to

07h

the exact form of which is given in Beutler (2003). Despite the fact that the system is linear, the variational equations usually have to be solved numerically. If the equations of motion are, however, those of the one-body problem (l), the partial derivatives of the orbit may be calculated analytically. This is a rather important fact, because in most cases it is sufficient to approximate the partial derivatives of the orbit by those of the one-body problem, even if the orbit itself is governed by the initial value problem (3). If no approximate values are available, we speak of first orbit determination, which is (at least) one order of magnitude more complicated than the orbit improvement problem. We will discuss this problem in the section on orbit determination based on astrometric places.

REAL-TIME

ORBIT

IMPROVEMENT

USING

SLR OBSERVATIONS

OF ONE STATION

Figure 4 (left) shows the measurement signature of a typical LASER satellite pass (of Lageos 1 in the particular case) over a space-geodetic observatory. When Lageos 1 was first observed in this pass, the light travel times were about 53 ms, corresponding to a distance of about 7950 km. Afterwards, the light travel times continually decreased to reach a minimum at the epoch of the closest approach (the minimum At M 39 ms corresponds to a distance of about 5850 km). After closest approach, the light travel times grow again. The entire pass lasted for about 40 minutes. A total of about 11200 Laser pulses were sent out, about 2800 candidate echoes were detected, a real-time analysis during the pass accepted 1882 measured light travel times as candidate echoes. A screening procedure (briefly described below) eventually accepted about 1817 measurements as real echoes. This performance is typical for day-light passes. The night-time performance (as judged from the ratio of accepted echoes and the total number of pulses) is much better. The gaps between the recorded data in Figures 4 are not due to instrument failures. They are caused by the fact that two other laser satellites were tracked during the time period of the Lageos 1 pass, too. Figure 4 (right) illustrates the real-time parameter estimation process performed during each satellite pass: The predicted orbit of the satellite and the (rather well) known geocentric coordinates of the observatory allow it to compute the difference between the predicted ranges and the actual measurements. Predictions are, however, never 100% accurate. In the case of a satellite orbit the uncertainty is mainly along track, i.e., all orbital elements may be assumed to be perfectly known, except one, the perigee passing time TO (actually the mean anomaly at to is used as an orbit parameter in this application). The predictions allow it to define a so-called range gate. Only echoes within a (observer-defined) time interval centered at the predicted ranges (originally f70 ns in Figure 4 (right)) are considered as candidate echoes. With a relatively high degree of probability it is then possible to decide in real time whether a registered echo within this range gate is real

Orbit

Determination Lageos-1:

Fig. 5. Screened

residuals

(ns) of Lageos

in Satellite

Residuals

1 observations

after

Geodesy

1859

Screening

at Zimmetwald

Observatory

on July 7, 2002, 07h19m

-

08hOOm. by looking for “identical” values in a list of recently established values “observed minus predicted”. If three or more coincidences are encountered (by a majority voting technique), the accepted observations may be used to determine an improved value for To. This is realized by a conventional orbit improvement process with only one unknown. From this time onwards, the predictions are much more reliable, allowing it to narrow down the range gate. Obviously, this was done in the example illustrated by Figure 4 (right). After each longer break, the range gate is reset to the original value. Unnecessary to say that the determination of TOwas based on the variational equation for the element TO. Real-time decisions have to be based on a limited amount of information. This is why a more correct analysis, based on a full orbit improvement process with all six elements (usually no other parameters have to be introduced for the short arc of one pass), is performed after the satellite pass before the data are sent to the ILRS data centers. It is in general not possible to accurately determine all six osculating elements using the range observations of one observatory only, even if the individual observations are of “infinite” accuracy, because the orbital plane is very poorly defined (this plane may in essence rotate about a straight cone with the axis pointing from the geocenter to the observatory). In practice this means that the orbital elements i and R have to be slightly constrained in the analysis. The orbit improvement is done iteratively, where in principle the flag of every observation (indicating whether the observation shall be used for the next step) may be reconsidered in each step. This process, the result of which is documented in Figure 5, is fully automated. The rms-error of this particular pass was 0.15 ns, corresponding to a mean error of about 2 cm in the measured ranges. ORBITS OF LOW EARTH ORBITERS USING THE GPS Kinematic and Deterministic Orbits All code observations (6) recorded simultaneously at time t by a space-borne GPS receiver on a Low Earth Orbiter (LEO) may be used to determine the current LEG position with an accuracy of a few meters. This estimate is completely independent of any orbit model. If the GPS satellite clock term At3 in eqn. (6) is taken over from analyses performed by the IGS Analysis Centers, four unknowns, namely the three Cartesian coordinates and the receiver clock term have to be determined. A minimum of four GPS satellites have to be observed simultaneously for this purpose. Due to the presence of the ambiguity term NX it is not possible to determine the satellite position pertaining to one observation time using only the phase observation equations (7) of this particular epoch only. It is, however, possible to combine all code observations (6) and all phase observations (7) of a certain time period (e.g., of one day) to determine all satellite positions, receiver clock corrections, and ambiguity terms within this period.

1860

G. Beutler

Fig. 6. Determination

of LEO position

differences

et al.

using the phase difference

observations

of subsequent

epochs.

Assuming that the LEO GPS receiver’s measurement frequency is 0.1 Hz we obtain a parameter estimation problem with 4.8640 unknown coordinates and clock parameters plus a certain (rather large) number of ambiguity parameters to be determined simultaneously for a time interval of one day. Problems of this dimension are not easily dealt with, but they promise a highly accurate result, namely a purely kinematic LEO trajectory with an accuracy of only a few cm per coordinate for each satellite position (the receiver’s clock corrections and the (real-valued estimates of the) ambiguity terms as byproducts. These satellite positions may be further analyzed by interpreting them as pseudo-observations in an orbit determination process using a model of type (3). U sually a large number of dynamical orbit parameters (gravity field parameters and parameters characterizing the non-gravitational forces) have to be determined in this latter process. The result of this latter orbit improvement problem is a pure deterministic orbit. Alternatively it is, of course, also possible to combine the two orbit determination problems (the determination of a kinematic trajectory and the deterministic problem using the satellite positions as pseudoobservations) into one process and to establish directly a deterministic trajectory by writing the observation equations (6) and (7) as a function of the orbit parameters (9) and by solving directly for these parameters. The result is a very bulky orbit determination process. For pure research purposes, where highest accuracy is aimed at, such procedures must be considered. For attempts in this direction we refer, e.g., to Svehla and Rot hacher (2002). There are many applications where the accuracy requirements are less stringent, but where efficiency and robustness are central issues. It may be sufficient in many cases to determine the LEO orbit with an accuracy of a few dm, but (almost) in real time. In these cases it would be nice to make use of the phase observations without invoking an orbit determination process with thousands of unknowns to be determined simultaneously. This is indeed possible: one may set up a parameter estimation process based on the phase observation equations (7) of two subsequent epochs ti and ti+l with the goal to determine the position F i g ure 6 illustrates the observation situation for the two epochs tl and difference Ari+i,i A r(ti+i) - r(ti). tz.

The principle is simple: instead of analyzing the ph+se observations of only one epoch one analyzes the phase difference observations of one and the same satellite for two subsequent epochs. These phase differences are ambiguity free. From eqns. (7) one obtains: A&+l,i

G &+I - 4i = p(ti+l)

-

p(h)

+

c

(At,,i+i,i

- Ats~i+‘~i) ,

(13)

where At,,i+i,i A At,i+i - At,i is the error of the receiver clock difference between the two epochs. At first sight the equations (13) are not very useful: whereas it is true that the ambiguity term was eliminated, and that we only have one unknown clock term (assuming that the GPS satellite clock difference

Orbit

Determination

in Satellite

Geodesy

1861

term is taken over as known from IGS analyses) we have now six instead of only three unknowns related to position, namely the Cartesian coordinates of the vectors I and r(ti+r). One may easily show, however, that the observable is orders of magnitude more sensitive to the position difference Ar(ti+,,i) 2 r(ti+i)-I than to either of the two position vectors of the differences. This sensitivity difference is indeed such that the dependence of eqn. (13) on the position may be taken as known from the analysis of the code observations (or from an a priori orbit) and to retain only the position difference as unknown in eqn. (13). This approach was first presented by Bock, (2000) and it is briefly described in Beutler (2003). Having determined the position differences in the way outlined above and the positions in the conventional way, one may use the positions (accurate to about l-2 m) and the position differences (accurate to a few cm) as pseudo-observations in a deterministic (or a generalized) orbit determination process to be discussed in the next section. It is possible to reconstruct the LEO trajectory with an accuracy of one to a few dm using this approach. Empirical Accelerations and Pseudo-Stochastic Pulses Emnirical Accelerations When analyzing long arcs (5 to 20 days) of GPS satellites at the Center for Orbit Determination (CODE) Analysis Center of the IGS it became evident that it is necessary to introduce empirical acceleration terms in order to represent the GPS satellite positions determined routinely by the IGS Analysis Centers at the subdm level. The models are described in detail in (Beutler et al., 1994). The perturbing acceleration, which is attributed to radiation pressure, may be written as (after having taken into account the deterministic model to the extent possible): %-ad

=

1

~~~~~+ X(t) 0

ex + Y(t) ey + Z(t) ez

; ;

in sunlight in Earth shadow

(14)

where ano& is the a priori model as documented in (Fliegel et al., 1992). ex, ey. and ez are the unit vectors in three orthogonal directions. According to (Beutler et al., 1994) the three components are defined as X(t) Y(t) Z(t)

A A G

X0 + X, cos(u + 4~) Yo +Yp cos(u+q&) 20 + 2, cos(u + 4~)

= = =

X0 + X, cost + X, sin21 YO+Yc cosu+Ys sinu 20 + 2, cosu + 2, sinu ,

(15)

where ‘u is the satellite’s argument of latitude. The unit vectors ex, ey and ez are defined as the unit vector from the satellite to the Sun (Z), the unit vector along the space vehicle’s solar panel axes (Y), and the third vector X completing the orthogonal system. The empirical model (14), (15) p roved to be an essential improvement over the predecessor which just allowed for a scale factor of the a priori constituent a&,& and a constant y-bias. Based on this empirical model Springer et al. (1999) came up with an improved radiation pressure model for GPS satellites. Other applications may ask for other decompositions. As an alternative we also use the decomposition into the radial, along-track, and out-of plane directions - a decomposition leading to an approximately orthogonal system in the case of low-eccentricity orbits, The second option of decomposition may be preferable if atmospheric drag instead of radiation pressure is the principal source for the unmodelled accelerations. The empirical model (15) may be easily generalized by adding periodic terms with arguments Ic IL, k=2,3,... (the model (15) may be viewed as a harmonic series with the first two terms corresponding to k = 0,l). Such a generalization is considered for implementation in the near future. Pseudo-Stochastic Pulses Usually, the attitude of active scientific satellites is maintained by the operators of the satellite. It may also be necessary to perform small orbital manoeuvres from time to time to optimize the orbital characteristics of a satellite. When determining the orbit of a LEO one has moreover the problem that the forces acting on the satellite may vary rapidly in space and time. It is, e.g., close to impossible to account for atmospheric drag with highest accuracy using deterministic models. One radical method of curing modelling problems is to break up the original arc into many shorter arcs. The modelling deficiencies are absorbed by the (frequent) set up of initial state vectors. This well-known

1862

G. Beutler

et al.

method usually is called the short-arc method. One should be aware of the fact, however, that the procedure in essence multiplies the number of parameters by the number of arcs - what considerably weakens the solutions. It would therefore be preferable to allow from time to time for instantaneous velocity changes bvi in predefined directions, where these velocity changes are considered (and determined) as orbit parameters in conventional least squares procedures. Subsequently, an instantaneous velocity change SIJ~at time t in a predetermined direction ei will be called a pseudo-stochastic pulse. Depending on the application it may be necessary to set up many of these pulses. It is possible to introduce up to three pulses at one and the same epoch t (in linearly independent directions). Usually, pseudo-stochastic pulses are set up either in the radial, the along-track, or the out-of-plane directions or in a combination thereof. The velocity changes may or must be constrained to “reasonably small” values by introducing artificial observations of the velocity changes with a weight proportional to IJ-~(&u~), where o(6vi) is the user-specified rms-scatter of the pulses. Modern orbit determination procedures for LEOs therefore require the capability of a stochastic component in the orbit model. One might in principle replace the deterministic differential equation systems describing orbital motion by stochastic di$erential equation systems, containing on top of all deterministic effects so-called stochastic accelerations which are characterized by the (known) mean values (usually the zero vector) and the associated (known) variance-covariance matrix. Stochastic modelling is, e.g., enabled, if the classical least-squares approach is replaced by digital filter methods, in particular Kalman- and KalmanBucy filters. In the framework of this more general theory every deterministic parameter may be replaced by a stochastic process. Not only the measurement noise, but also the process noise (represented by the stochastic differential equations in the case of orbit parameters) is considered in this case. Alternatively, the approach of pseudo-stochastic modelling may be used, which was introduced by (Beutler et al., 1994). The procedure may be characterized as follows: l

Each orbital arc is continuous.

l

Each arc is represented piecewise by ordinary differential equation systems, the deterministic of motion.

l

At predetermined epochs (e.g., every five minutes) the orbit is allowed to change the velocity instantaneously in (up to three) pre-determined directions. These velocity changes are called pseudo-stochastic pulses.

l

The pseudo-stochastic

l

Each pseudo-stochastic

pulses 6vi are parameters of the classical least-squares adjustment

equations

process.

pulse is characterized by an a priori variance.

There is a certain degree of arbitrariness in such an approach (exactly as in the case of the Kalman-Bucy filtering): the number of epochs and the associated variances may be arbitrarily selected. The procedure is on the other hand very flexible: depending on the actual definition of the number and the properties of the pulses the orbit may be either purely deterministic, purely kinematic, or something in-between, usually referred to as reduced dynamics orbit. The size of the velocity changes are steered by artificial observation equations of the kind 6v; = 0 )

(16)

with prescribed weight

. 00” wi=o2(6ui)*

(17)

The scalar velocity change Svi thus is constrained as a random variable with expectation value zero and variance a2(Swi). 00 is the mean error of the weight unit (observation with a mean error of o = 1). If (~(bvi) is big, the weight wi is small, which allows 6vi to assume rather big values. If a(bvi) is small, only minor velocity changes are allowed, roughly within the range f3 ~(Svi).

Orbit 6

Determination

R-

0

ml

404

600

.soo mill

1oDo

ml0

in Satellite

Geodesy

1863

0.08

R$ ...*...

_..*

1400

1600

0

200

400

600

800 Illin

looo

1200

1400

1600

Fig. 7. Residuals of positions (left) and position differences (right) (in radial (R), along-track (S) and out-of-plane (W) directions) in an orbit determination process of one day of CHAMP data (data sampling frequency 0.1 Hz, six osculating elements, nine empirical accelerations, and 429 pseudo-stochastic pulses estimated).

The partial derivatives associated with the pseudo-stochastic pulses may be represented as linear combinations of the six variational equations associated with the original six osculating elements (Beutler et al., 1994, 2003). Even if many parameters have to be set up it is possible to solve for these parameters in a very efficient way. Figure 7 shows the residuals of the estimated positions (left) and position differences (right) of one day of 10-s observations of the GPS receiver onboard the CHAMP (homepage http://op.gfz-potsdam.de/champ/index-CHAMP.html) spacecraft. The root mean square error of the kinematically determined satellite positions is about 1.3 m, that of the position differences of the order of 1 cm. The orbit quality could be compared with other solutions based on a more elaborate and mathematically correct analysis. According to Beutler (2003) it is about 20 cm rms. The orbit determination procedure is very efficient and robust. FIRST ORBIT DETERMINATION USING ASTROMETRJC PLACES The difference between orbit improvement and first orbit determination is similar - qualitatively - to the difference between craftsmanship and fine arts. In first orbit determination we do not try to reduce the non-linear orbit determination problem to a linear one but make instead the attempt to solve the non-linear problem directly. In principle it would be possible to define a first orbit determination problem for all observation types (and all possible combinations thereof). Here we confine ourselves to the determination of a first orbit from a short series of astrometric places, where the attribute “short” characterizes the time interval containing all observations (4): We assume that this time interval is short compared to the revolution period of the observed object. The number of astrometric places must exceed two, if we want to determine the six osculating elements. The method for first orbit determination outlined below is illustrated using the observations from a search campaign for space debris (Schildknecht et al., 2001,2002). Figure 8 illustrates two observation scenarios for the two methods for objects in the geostationary belt. The CCD-pictures were taken with the Zimmerwald l-m telescope (see Figure 2 (right)) on November 2, 2001. In Figure 8 (left) the satellites are moving w.r.t. the star background, because the telescope tracking is compensating for the diurnal motion of the stars, in Figure 8 (right) the telescope is in the “staring mode”, i.e., it points into one and the same Earth-fixed direction; this mode implies that geostationary objects are mapped as “points” and the stars as dashes. For space debris search campaigns in the geostationary belt the second mode is preferable (better signal-to-noise ratio for satellites). For more information concerning the astrometric observation technique developed for moving objects or for the description of search campaigns for space debris we refer to (Schildknecht et al.,

G. Beutler

1864

et al.

Fig. 8. Astrometric CCD observation of seven geostationary Astra satellites at 19.2” eastern longitude with the Zimmerwald observatory (left: sidereal tracking, right: instrument Earth-fixed). Note the satellite motion w.r.t. the stars within the 14 minutes between the two pictures

2001, 2002).

Assuming short observation time spans we may usually approximate the orbit as a particular solution of the differential equations (1) of the one-body problem. For the method presented below, this is not a requirement - we might as well assume that the orbit is governed by a differential equation system of type (2). We will always assume, however, that no dynamical parameters have to be determined. Under these assumptions we have to solve the six-dimensional non-linear parameter estimation problem of determining the six osculating orbital elements shown in Figure 1 from n 2 3 astrometric places characterized by eqn. (4). In this “naive” formulation the parameter estimation problem is (probably) unsolvable - or the computational effort to solve it would be tremendous. The non-linear parameter estimation problem is solved in the following steps: 1. The dimension of the problem may be immediately reduced from six to two by formulating the problem not as an initial value problem (as indicated by eqns. (3)), but as a boundary value problem. i: f(ta - %)

= =

-GM 5 + df(t; T, +) R(t,) + pa e,

T(tb

=

R(tb)

-

$1

+

Pb eb

(18)

,

where the boundary epochs t, and tb must be selected as two observation epochs (not necessarily the first and the last one). The indices a and b thus characterize two observation numbers. e,, eb are the observed unit vectors, and pa, & are the (originally unknown) distances. 2. The osculating elements as unknowns are replaced by the following six parameters knowns to be determined first): {Pl,PZ>.

. . ,pS}

=

{pa,

/‘b, %,

abr

da, 66)

.

(as auxiliary un-

(19)

3. For any given set of parameters (19) the boundary value problem (18) is solved by numerical integration. In the program package included in Beutler (2003) a collocation method is used for that purpose.

Orbit

Determination

in Satellite

Geodesy

1865

After the numerical solution of the boundary value problem (18) the orbital arc is available as a Taylor series

r(t) = -& $ @(t - ta)i )

(W

i=O

where the development origin t, could be chosen arbitrarily. values are 10 5 q 5 14.

q is the order of the Taylor series, typical

4. Formula (20) allows us to compute the position- and velocity-vectors, r(t) and t(t), for any epoch t, 2 t 5 tb. As there is a one-to-one correspondence between the osculating elements referring to epoch t and the position- and velocity-vectors of the same epoch, we are also in a position to calculate the osculating elements corresponding to parameters (19) - after having solved the boundary value (18) problem numerically. 5. Excellent approximate values are available for the latter four of the above six unknowns (19): After all, the right ascensions and the declinations were observed directly. We may therefore use the approximations {Pl,PZ,

(P3,.

. ’ ,pS)}

=

{Pdb,

@‘,,

&

$7

(21)

$1)

in order to establish a first crude orbit. 6. The two unknowns left are the two topocentric distances pa and pb. They have the following properties: l

pa and pb have to be positive and

l

they cannot differ by too much - at least if the time interval 1 tb - t, ) is short (“infinitesimal”).

7. We might define a two-dimensional search, where pa and 1 Pb - pa 1 are varied within reasonable (e.g., user-defined) limits. For each pair of values the mean error of the observed angle is computed as

2n The factor cos2 6: in front of the residuals in right ascension is introduced to make small variations in (Y and S comparable on the celestial sphere). The minimum (or possibly the minima) of function B(p,, Pb) are found by a two-dimensional search in the table of the function values l3 (pa, pb) This iterative search is “fool proof’ broad enough) way.

as long as the search ranges are defined in an appropriate

(i.e.,

8. In view of the fact that I tb - t, 1 is assumed to be a small quantity and that therefore / pb - pa / must be small as well, we may replace the two-dimensional search (22) by varying only ra systematically and by determining for each selected value of pa that value for the parameter POwhich minimizes the sum (22). For that purpose a one-dimensional and solved. It is initialized by

orbit improvement

problem with pb as only unknown

is set up

(23) i.e., the observations al,, o$, and bh, 6; replace the observed functions, orbit improvement process.

and pb 2 pa to initialize the

1866

G. Beutler

et al.

1

J looo4

0

ml00

3oooo

4oooo

.5oooO

6mo

7oooo

within

lh 26m (solid line) or 5 places

km

Fig. 9. Function log(B(p,)) within 4m (dotted line).

using 23 astrometric

places of a GTO-object

The resulting sum of residuals squares after the solution of the orbit improvement as

process then reads

2 [cos2 s: (a;- cYi)2 + (6;- Q2] &a)

A

%‘a,k’b(h))

Sz

d i=l

(24

2n - 1

Figure 9 shows graphs of the above defined function for two orbit determinations for the same object. The solid line stands for the case where all astrometric places were used, the dashed line for the case where only five places were used. 9. The minimum (or the minima) of B(p,) is (are) easily found by a one-dimensional search in the resulting table for function B(p,). Observe that B(p,) is the mean error of one observed angular coordinate a posteriori (i.e., of (Ycos2 6’ or 6) of the orbit improvement process. 10. The values pa and pb corresponding to the minimum’of the function B(p,) together with the observed unit vectors ea,b define the first reference orbit for a full orbit improvement process, where all six parameters (19) are determined. Let us mention that the six partial derivatives of the position vector r(t) w.r.t. the orbit parameters are calculated analytically using the formulae of the one-body problem (see Beutler (2003)). After completion of the final orbit determination

step the osculating elements are computed.

The orbit determination procedure outlined above resembles in some (essential) aspects the Gaussian orbit determination procedure. The formulation of the orbit determination problem as a boundary value problem of type (18) and not as an initial value problem of type (3) is due to Gauss, and, admittedly, this is the key to the solution. It should be pointed out, however, that from there onwards our solution strategy differs considerably from that originally proposed by Gauss. For more information we refer to Beutler (2003). Let us apply the theory outlined above to an object found in the context of a search for space debris in geostationary transfer orbits (GTO) using the ESA 1 m telescope in Teneriffe. The search strategy and the observatory are described in Schildknecht (2001). 23 astrometric places were measured in the night of September 8 and 9, 2002 between 23h and Oh30m for an object suspected to be in a geostationary transfer orbit. Five observations were made within five minutes, then other areas in the sky were searched. The same object was re-observed after a break of about 45 minutes nine times within four minutes, then, after another break of about thirty minutes nine times within eight minutes.

Orbit

Table 1. Orbit

I = E I = NODE PER = TPER =

elements

23819403.496614 0.7146729612 6.6971388 92.3293039 -248.1323982 -14M11.4262976

for the previously

+/+I+I+I+I+,-

2604.810940 O.OWO696G42 0.0006263 0.0032306 0.0024963 3.0647606

Determination

unknown

in Satellite

object

Geodesy

1867

geo-004

II m.G) (OEG) m.0) SEC

GTO objecl geo-004 (WI-2002OSOSS4-16)

-.0

40

20

30

40

60

60

70

*o

90

100

Epoch [minutes]

Fig. 10. Residuals

of orbit

determination

for object

geo-0004

Table 1 gives the final results of the orbit determination for the object considered as an example in this article, Figure 10 the corresponding residuals. The orbit determination was based on all 23 astrometric places. The conventional elements of Figure 1 including their mean errors are provided as a result. The orbit model included the Earth’s oblateness term and lunisolar gravitation. The estimated mean error of the one observation is 0.27” which is a remarkable accuracy for astrometric places of moving objects and sufficient for search purposes. The orbit quality would be sufficient to find the object in subsequent nights.

SUMMARY

AND

CONCLUSIONS

We reviewed the case of pure orbit determination encountered today in satellite geodesy. We identified three areas, namely the real-time optimization of SLR observations, the determination of LEO orbits using data from spaceborne GPS receivers, and the determination of satellite and/or space debris orbits from astrometric places as prominent applications. In the introductory section the equations of motion and the orbit parametrization were defined. Then, three essential observation techniques, namely astrometric places, SLR light travel times, and GPS (code and phase) observations, were introduced. We made the distinction between orbit determination and orbit improvement. In all applications orbit improvement plays a key role. In this problem type we assume that approximate values for the parameters are available. These values allow us to linearize the nonlinear functional dependence of the observed functions from the orbit parameters. According to eqn. (11) this relationship is characterized by the partial derivative of the observed functions w.r.t. orbit parameters, which may be written as scalar products of two vectors. One of these vectors is the partial derivative of the reference trajectory r(t) w.r.t. orbit parameters. This vector is a solution of the so-called variational equations (12). Approximate analytical solutions (based on the one-body motion) are known and may be used in all cases considered here.

1868

G. Beutler

et al.

We distinguished three types of orbit parameters, namely the osculating elements, socalled dynamical parameters, and so-called pseudo-stochastic pulses. The osculating elements (or a subset thereof) as illustrated in Figure 1 have to be set up in all cases. Dynamical parameters define the force (acceleration model) on the right-hand sides of the equations of motion. For many applications, e.g., for gravity field determination, these parameters are of primary interest. We also gave an example for empirical parameters as they are commonly used to characterize radiation pressure for navigation satellites. When considering the orbits of low Earth orbiters it may be required to introduce in addition so-called pseudostochastic pulses, i.e., instantaneous velocity changes in pre-determined directions in order to absorb non-modelled effects (e.g., due to an imperfectly known gravity field or atmospheric drag). The partial derivatives w.r.t. these parameters are linear combinations of the six partials associated with the osculating elements of an arc. It is therefore possible to implement very efficient procedures to solve for such parameters - even if hundreds of them are introduced. Exactly one case of first orbit determination was encountered in this article, namely the determination of an orbit from a short series of astrometric places. The problem is, e.g., encountered in optical search campaigns for space debris. A very robust and efficient method, essentially based on Gauss’s idea to transform the orbit determination problem into a boundary value problem, was outlined and illustrated with (only) one example. Unnecessary to say that the same method may be applied to minor planets, comets, etc. It is in particular planned to adapt the method to near Earth asteroids (NEAs) in the near future. REFERENCES Beutler, G., Methods of Celestial Mechanics: Basic Principles and Applications, Springer-Verlag Heidelberg, 2003. Beutler, G., E. Brockmann, W. Gurtner, et al., Extended orbit modeling technique at the CODE processing center of the International GPS Service for Geodynamics (IGS): Theory and initial results, Manvscripta Geodaetica, 19, 367-386, 1994. Bock, H., G. Beutler, S. Schaer, et al., Processing aspects related to permanent GPS arrays, Earth Planets Space, 52, pp. 657-662, 2000. Fliegel, H. F., T. E. Gallini, and E. R. Swift, Global Positioning System radiation force model for geodetic application, Journal of Geophysical Research, 97, No. Bl, 559-568, 1992. Schildknecht, T., M. Ploner, and U. Hugentobler, The search for debris in GEO, Adv. Space Research, 28, No. 9, pp.1291-1299, 2001. Schildknecht, T., R. Musci, M. Ploner, et al., Optical observation of space debris in GE0 and in highlyeccentric orbits, COSPAR Scientific Assembly, Ott 10 - 19, 2006, Houston, Texas, USA, submitted to Adv. Space Research, ZOO&. Springer, T., G. Beutler, and M. Rothacher, A new solar radiation pressure model for GPS satellites, GPS Solutions, 2, No. 3, 50-62, 1999. Springer, T., Modeling and validating orbits and clocks using the Global Positioning System, Geodiitischgeophysikalische Arbeiten in der Schweiz, 60, Schweizerische Geodatische Kommission, Institut fiir Geod&sie und Photogrammetrie, Eidg. Technische Hochschule Zurich, Zurich, 2000. Svehla, D., and M. Rothacher, Kinematic orbit determination of LEOs based on zero or double-difference algorithms using simulated and real SST data, in Vistas for Geodesy for the New Millenium, IAG Symposia Vol. 125, edited by J. Adam and L.P. Schwarz, pp. 322-328, 2002. Teunissen, P.J.G. and A. Kleusberg (Eds.), GPS for Geodesy, 2nd Edition, Springer- Verlag Heidelberg, 1998. E-mail address of G. Beutler: [email protected] Manuscript received 3 December 2002; revised 7 March 2003; accepted: 14 March 2003