Keplerian orbit

For the moment, it is assumed that both the position and the velocity vector of the satellite have been derived from observations. Now, the question arises of how to use these data for the derivation of the Keplerian parameters.

The position and velocity vector given at the same epoch t define an initial value problem which has been treated in Sect. 4.2.1. Recall that the two given vectors contain altogether six components which enable the calculation of the six Keplerian parameters.

Position vectors are preferred for orbit determination since they are more accurate than velocity vectors. Let us assume that two position vectors gS(tl) and fl(t2) at epochs tl and t2 are available. These data correspond to boundary values in the solution of the basic differential equation of second order, cf. Eq. (4.1). An approximate method for the derivation of the Kep-lerian parameters makes use of initial values defined for an averaged epoch t

= !

^(tl

+

t2):

(}s(t)

=

^gS(t2)

+

gS(tt)

- 2

iJS(t)

=

^{(}S(t2) -} (}s(tt) .

- t2 -

tt

( 4.49)

The rigorous solution starts with the computation of the geocentric distances

rl

=

r(tt)

=

lI(}s(t1)1I

r2

=

r(t2)

=

IIgS(t2)1I. ^{( 4.50)}

4.3 Orbit determination 55 The unit vector ~, orthogonal to the orbital plane, is obtained from a vector product by

gS(tt) ^X gS(t2)

~ =

^{IIgS(tl) X gS(t2)1I (4.51)}

and produces the longitude l and the inclination angle i, cf. Eqs. (4.11) and (4.13). As demonstrated earlier, the argument of the latitude ^U = w

+

v is defined as the angle between the satellite's position and the ascending node vector

Ii

= (cosl, sinl, O)T. Hence, the relation

i

=

^{1,2 (4.52)}

holds from which the Uj with U2

>

Ul can be deduced uniquely. Now, there are two equations, cf. Eq. (4.7),

a(1- e^{2 )}

Tj = ₁

+

^{ecos(uj - w)} , ^{i = 1,2 ( 4.53)}

where the parameters a, e, ware unknown. The system can be solved for a and e after assigning a preliminary va1m~ such as the nominal one to w, the argument of the perigee. With the assumed wand the ^Uj the true anomalies

Vi and subsequently the mean anomalif~s Mi are obtained. Therefore, the mean angular velocity n can be calculat·:!d twice by the formulas

(iL M2 - Ml

n

= V;;3 =

^{t2 - tl ' (4.54)}

cf. Eqs. (4.2) and (4.3). The equivalence is achieved by varying w. This iterative procedure is typical for boundary value problems. Finally, the epoch of perigee passage To follows from the relation

To

=

ti - - ' . M·

n (4.55)

For a numerical solution of the boundary value problem assume two position vectors g(tl) ^and g(t2), both represented in the earth-fixed equatorial system

X4

and with At

=

t2 - tl

=

^{1 hour:}

g(tl) = (11465, 3818, -20 923l [km]

g(t2)

= (

5220,16754, -18421l [km].

The application of the Eqs. (4.50) through (4.55), apart from rounding errors, results in the following set of parameters for the associated Kepler ellipse:

a

=

^{26000km, e}

=

^{0.1, w}

=

^{-140°, i}

=

^{60°, l}

=

^{110°, and} To

=

tl-1~3183.

56 4. Satellite orbits H there are redundant range observations, the parameters of an instan-taneous Kepler ellipse can be improved. The position vector

gff

associated to the reference ellipse can be computed, and each observed range gives rise to an equation

I ^·s gff - fR

^S

e

=

^eo

+

de

= I

~

- gRIl + lieS _ ell'

dg . so -R

(4.56) The vector deS can be expressed as a function of the Keplerian parameters, d. Eq. (4.29). Thus, Eq. (4.56) actually contains the differential increments for the six orbital parameters.

An orbit improvement is often performed in the course of GPS data processing when in addition to terrestrial position vectors

gR

the solution also allows for the determination of the increments dPoi. Of course the procedure becomes unstable or even fails for small networks. Sometimes only three degrees of freedom are assigned to the orbit by a shift vector.

This procedure is called orbit relaxation.

4.3.2 Perturbed orbit

Analytical solution. As known from previous sections, the perturbed mo-tion is characterized by temporal variamo-tions of the orbital parameters. The analytical expressions for these variations are given by the Eqs. (4.37) or (4.39).

In order to be applicable for Lagrange's equations, the disturbing poten-tial must be expressed as a function of the Keplerian parameters. For the effects of the nonsphericity of the earth the Legendre functions in Eq. (4.40) must be transformed which first has been performed by Kaula (1966). The resulting relation is

00 n n

R =

E E

An{a) EFnmp{i)

n=2m=O FO

00 (4.57)

. E

Gnpq{e)Snmpq{w,n,Mj00,Jnm,Knm)'

q=-oo

Recall that n denotes the degree and m the order of the spherical har-monics in the disturbing potential development. Each of the functions

An, Fnmp, Gnpq contains only one parameter of a Kepler ellipse. However, the function Snmpq is composed of several parameters and can be expressed by

Snmpq

=

Jnm cos

t/J +

^{Knm sin}

^t/J

Snmpq = -Knm cos

t/J +

Jnm ^sin

t/J

for n - m even

for n - m odd. (4.58)

4.3 Orbit determination

Table 4.4. Perturbations due to the earth's gravity field Parameter Secular Long-period Short-period

a no no yes

e no yes yes

t no yes yes

n

^{yes yes yes}

w yes yes yes

M yes yes yes

The relation of the frequency

tb

associated to the argument 'IjJ is

tb

^{= (n - 2p)w}

⁺

^{(n - 2p}

⁺

^q)

^{if +} men -

^{00 )}

which is a measure for the spectrum of the perturbations.

57

( 4.59) The conditions (n-2p) = (n-2p+q) = m = 0 lead to

tb

⁼ 0 and thus to secular variations. Because of m = 0, they are caused by zonal harmonics.

If ( n - 2p) " 0, then the variations depend on wand are, therefore, generally long-periodic. Finally, the conditions (n - 2p

+

^q) " 0 and/or m " 0 result in short-period variations. The integer value (n - 2p+ q) gives the frequency in cycles per revolution and m the frequency in cycles per day.

A rough overview of the frequency spectrum in the Keplerian parameters due to the gravity field of the earth is given in Table 4.4. Summarizing, one can state that the even degree zonal coefficients produce primarily secular variations and the odd degree zonal coefficients give rise to long-periodic perturbations. The tesseral coefficients are responsible for short-periodic terms. From Table 4.4 one can read that short-period variations occur in each parameter. With the exception of the semimajor axis, the parameters are also affected by long-periodic perturbations. However, secular effects are only contained in

n,

w, M. The analytical expression for the secular variations in these parameters due to the oblateness term J₂ is given as an example:

. 3 ² cosi

n =

^{-2naE a2(1- e2)2}

h

3 ² 5 cos²i - 1

w =

4"

n aE a2 (1 _ e2)2

h

^(4.60)

. 3 ² 3 cos²i - 1

m =

4"

n aE a2 J(l _ e2)3

h.

58 4. Satellite orbits

Fig. 4.3. Secular perturbations caused by the oblateness term J2

The first equation describes the regression of the node in the equatorial plane, the second equation expresses the rotation of the perigee, and the third equation contributes to the variation of the mean anomaly by

M =

n

+ m.

In the case of GPS satellites, the numerical values

n

^{~ -0.°03}

per day,

w

^{~ 0~01} per day, and

m

~ 0 are obtained. The result for

m

is verified immediately since for the nominal inclination i = 55° the term 3 cos²i - 1 becomes approximately zero. A graphical representation of the secular perturbations is given in Fig. 4.3. Special attention must be paid to resonance effects which occur when the period of revolution corresponds to a harmonic in the gravity potential. This is why the GPS satellites are raised into slightly higher orbits to preclude an orbital period very close to half a sidereal day.

The tidal potential also has a harmonic representation. Analogously, the tidal perturbations can be analytically modeled. This was done first by Kozai (1959) and, analogously to the earth's potential effect, has led to analytical expressions for the secular variations of the node's right ascension

n

and of the perigee's argument w. The reader is referred to the publication of Kozai (1959) for formulas specified.

Numerical solution. If the disturbing acceleration cannot be expressed in analytical form, one has to apply numerical methods for the solution. There-fore, in principle, with initial values such as the position and velocity vectors

4.3 Orbit determination 59 g(to) and g(to) at a reference epoch to, a numerical integration of Eq. (4.30) could be performed. This simple concept can be improved by the intro-duction of a Kepler ellipse as a reference. By this means only the smaller difference between the total and the central acceleration must be integrated.

The integration results in an increment lig which, when added to the position vector computed for the reference ellipse, gives the actual position vector.

The second-order differential equation is usually transformed to a system of two differential equations of the first order for the numerical integration.

This system is given by

t

get) = g(to)

+ J

^{g(to) dt}

to t

= g(to)

+ J

^{[dg(to) -}

^e3~to)

g(to)] dt

to

(4.61)

t

get) = g(to)

+ J

^{g(to) dt.}

to

The numerical integration of this coupled system can be performed by apply-ing the standard Runge-Kutta algorithm, cf. for example Kreyszig (1968), p.89. This method is shortly outlined. Let y(x) be a function defined in the interval ^Xl ~ X ~ X2 and denote by y' = dy/dx its first derivative with respect to the argument x. The general solution of the ordinary differential equation of the first order

I dy '( )

y

= - =

_dx ^{y y,x ( 4.62)}

follows from integration and the particular solution is found after assigning the given numerical initial value YI = Y(XI) to the integration constant.

For the application of numerical integration, first the integration interval is subdivided into n equal and sufficiently small .6.x = (X2 - xI)/n where n is an arbitrary integer> O. Then, the difference between successive functional values is obtained by the weighted mean

.6.y = y( X

+

.6. x ) - y( X )

=

~

[.6.y(1)

+

2 (.6.y(2)

+

.6.y(3))

+

.6.y(4)] ( 4.63)

60 where

~y(1)

=

^y'(y, X) ~X

~y(1)

~y(2)

=

^y'(y

+ __ ,

2

~y(2)

~y(3)

=

^y'(y

+ __ ,

2

~y(4) = y'(y

+

^~y(3),

x

+

~X 2)~x x+2)~x ~x

x

+

^{~x )~x.}

4. Satellite orbits

Hence, starting with the initial value YI for the argument ^x}, the function can be calculated for the successive argument Xl

+

^~x and so on.

Numerical methods can also be applied to the integration of the disturb-ing equations of Lagrange or Gauss. These equations have the advantage of being differentials of the first order and, therefore, must only be integrated once over time.

Let us conclude this paragraph with an example for numerical integra-tion. Assume the ordinary first-order differential equation y'

=

y - x

+

1 with the initial value YI

=

1 for ^Xl

=

^O. The differential equation should be solved for the argument ^X2

=

1 with increments ~x

=

0.5. Starting with the initial values for the first interval, successively the values ~y(1)

=

^1.000,

~y(2)

=

1.125, ~y(3)

=

1.156, and ~y(4)

=

1.328 are obtained. The corre-sponding weighted mean thus is ^~y

=

1.148. Replacing in the second interval the initial values by ^X

=

^Xl

+

^~x

=

0.5 and y

=

YI

+

^~y

=

2.148 and pro-ceeding in an analogous manner as before yields ^~y = 1.569 and thus the final result Y2

=

^y(l)

=

^Y

+

^~y

=

2.717. Note that the function y

=

^eX

+

^x,

satisfying the differential equation, gives a true value y(l) = 2.718. With an increment ~x

=

0.1 the numerical integration would provide an accuracy in the order of 10-^{6 •}

Trong tài liệu B. Hofmann-Wellenhof, H. Lichtenegger, and 1. Collins (Trang 70-76)