## ZZi ZZ

Where p (X - X )uX + (Y - Y )uY + (Z - Z )uZ, the position vector of the satellite with respect to the instrument. The relative velocity is P (X - Xi )ux + (Y- Y )uy + (Z - Z )uz. Note that the position and velocity of a ground-based instrument expressed in the nonrotating (X, Y, Z) system will be dependent on the rotation of the Earth. Equation (3.2.6) can readily be interpreted as the component of the relative velocity in the direction de ned by the relative position vector, p. In other...

## State Noise Compensation Algorithm

In addition to the effects of the nonlinearities, the effects of errors in the dynamical model can lead to divergence in the estimate. See, for example, the discussion in Schlee etal, (1967) . As pointed out previously, for a sufficiently large number of observations the elements of the covariance matrix Pk will asymptotically approach zero and the estimation algorithm will be insensitive to any further observations. This condition can lead to fiter divergence. One approach to preventing this...

## The Role of the Observation

As an example of the orbit determination procedure, consider the situation in which the state is observed at some time epoch, j. Then if Xj, Yj, Xj and Yj are given at tj, Eq. (1.2.2) can be used to form four equations in terms of four unknowns. This system of equations can be used to determine the unknown components of the initial state. For example, from Eq. (1.2.2) it follows that Then, the initial state can be determined as follows Providing that the matrix inverse in Eq. (1.2.5) exists,...

## F2 State Noise Compensation

The State Noise Compensation (SNC) algorithm (see Section 4.9) provides a means to improve estimation performance through partial compensation for the unknown sinusoidal acceleration. SNC allows for the possibility that the state dynamics are influenced by a random acceleration. A simple SNC model uses a two-state fiter but assumes that particle dynamics are perturbed by an acceleration that is characterized as simple white noise where u(t) is a stationary Gaussian process with a mean of zero...

## Qk Vfc vT

A Givens algorithm to compute an upper (or lower) triangular matrix, Wk, is given in Section 5.8, Eqs. (5.8.6) through (5.8.13). If there is no process noise, set Vk 0 in this algorithm. Methods for maintaining the measurement update in triangular form are discussed next. Triangular Square Root Measurement Update The algorithm for performing the measurement update on the state error co-variance matrix square root, W, can be expressed as If W is lower triangular and if A is lower triangular,...

## F1 Introduction

Consider a target particle that moves in one dimension along the x axis in the positive direction.* Nominally, the particle's velocity is a constant 10 m sec. This constant velocity is perturbed by a sinusoidal acceleration in the x direction, which is unknown and is described by The perturbing acceleration, perturbed velocity, and position perturbation (perturbed position nominal position) are shown in Fig. F. 1.1. A measurement sensor is located at the known position x 10 m. This sensor takes...

## Amsterdam Boston Heidelberg London New York Oxford Paris San Diego San Francisco Singapore Sydney Tokyo

Senior Publishing Editor Project Manager Editorial Coordinator Marketing Manager Composition Printer Frank Cynar Simon Crump Jennifer Hele Linda Beattie Author 200 Wheeler Road, Burlington, MA 01803, USA 525 B Street, Suite 1900, San Diego, California 92101-4495, USA 84 Theobald's Road, London WC1X 8RR, UK This book is printed on acid-free paper. (oo) Copyright 2004, Elsevier Inc. All rights reserved. No part of this publication may be reproduced or transmitted in any form or by any means,...

## N trtr n

N x - x, ej y - Hjx. (5.4.34) Consequently, ej j that is, the elements of the error vector e mx1 contain a contribution from errors in the a priori information as well as the observation residuals. The RMS of the observation residuals, j, is given by If the procedure is initialized with n observations in place of a priori information, x and P, Eq. (5.4.33) becomes and again ej < j because the fist n observations serve the same function as a priori values of X and P. Here we have assumed that...

## Example 2241

The ECI position and velocity of a Space Shuttle has been determined to be the following Y 3984001.40 m Z 2955.81 m X - 3931.046491 m sec Y 5498.676921 m sec With 3.9860044 x 1014m3 s2 and the preceding equations, it follows that the corresponding orbit elements are 5616.2198 sec 6767394.07 m ra 6890552.40 m f 44.608202

## Variation of Parameters

As illustrated in the previous section, the influence of a perturbing force on an orbit can be illustrated through time variations of the elements. The temporal variations illustrated for the Moon suggest the development of differential equations for the elements, analogous to the process known as variation ofparameters in the solution of ordinary differential equations. The differential equations of motion for the perturbed satellite problem can be expressed as where r is the position vector...

## TT dU 1 d U 1 dUn010

Although Eq. (2.3.18) contains a singularity at the poles, which poses a problem for exactly polar orbits, alternate formulations exist that remove this singularity (e.g., Pines, 1973). If the body-ix ed system is chosen to represent r , then r will include the Coriolis, centripetal, tangential, and relative acceleration terms. Alternatively, expressing r in the nonrotating expression is simpler. Since the gradient of U given by Eq. (2.3.18) provides force components in spherical coordinates, a...

## Observations

In the late 1950s, a ground-based tracking system known as Minitrack was installed to track the Vanguard satellites. This system, which was based on radio interferometry, provided a single set of angle observations when the satellite passed above an instrumented station's horizon. These angle measurements, collected from a global network of stations over several years, were used to show that the Earth had a pear-shape. They also provided the means to investigate the nature of forces acting on...

## State Transition Matrix

The matrix multiplying the initial state vector in Eq. (1.2.12) is referred to as the state transition matrix, (i,io). For the state deviation vector, x, defined in Eq. (1.2.17), this matrix is given by The state transition matrix maps deviations in the state vector from one time to another. In this case, deviations in the state are mapped from to to t. Note that we have arbitrarily set t0 0. For this linear example, the mapping is exact, but in the general orbit determination case, the state...

## Range Rate

Most range-rate systems in current use are based on a single propagation path, either uplink or downlink. The following discussion treats the problem from two viewpoints an instrument transmitting a short duration pulse at a known interval and a beacon that transmits a signal with a known frequency. Assume a satellite-borne instrument transmits a pulse at a specifed and fix ed interval, JtT. Hence, the pulses are transmitted at a sequence of times (e.g., tTi,tT2, etc). The transmitted pulses...

## GfM211

Where G is the Universal Constant of Gravitation and d is the distance between the two point masses or particles, Mi and M2. The Laws of Newton contain some important concepts. The concept of a point mass is embodied in these laws, which presumes that the body has no physical dimensions and all its mass is concentrated at a point. The concept of an inertial reference frame is crucial to Newton's Second Law that is, the existence of a coordinate system that has zero acceleration. The concept of...

## References

Kockarts, and G. Thuiller, A thermospheric model based on satellite drag data, Annales de G ophysique, Vol. 34, No. 1, pp. 9-24, 1978. Battin, R., An Introduction to the Mathematics and Methods of Astrodynamics, American Institute of Aeronautics and Astronautics, Reston, VA, 1999. Beutler, G., E. Brockmann, W. Gurtner, U. Hugentobler, L. Mervart, M. Rothacher, and A. Verdun, 'Extended orbit modeling techniques at the CODE Processing Center of the...

## The Gauss Markov Process

A fist-or der Gauss-Markov process is often used for dynamic model compensation in orbit determination problems to account for unmodeled or inaccurately modeled accelerations acting on a spacecraft. A Gauss-Markov process is one that obeys a Gaussian probability law and displays the Markov property. The Markov property means that the probability density function at tn given its past history at tn_2, --- is equal to its probability density function at tn given its value at A Gauss-Markov process...

## Example 225

The position and velocity of the Space Shuttle at a later time can be determined from the state given in Example 2.2.4.1. For example, at 30 minutes after the epoch used in Example 2.2.4.1, the eccentric anomaly is 159.628138 . It follows that the position and velocity at 30 minutes is An example of an ephemeris of positions, or table of positions, for the interval 30 min to 34 min is Figure 2.2.6 Shuttle ground track the latitude-longitude ephemeris of Example 2.2.4.1. The latitude-longitude...

## Summary

Following the introduction of the basic concepts of nonlinear parameter estimation in Chapter 1, Chapter 2 provides a detailed discussion of the characteristics of the orbit for Earth-orbiting satellites. This is the primary dynamical system of interest in this treatment. Chapter 3 describes the various observations used in the Earth-orbiting satellite orbit determination problem. Chapter 4 describes the properties of the Weighted Least Squares and Minimum Variance Estimation algorithms as they...

## Review of Matrix Concepts

The following matrix notation, definitions, and theorems are used extensively in this book. Much of this material is based on Graybill (1961). A matrix A will have elements denoted by aj, where i refers to the row and j to the column. At will denote the transpose of A. A-1 will denote the inverse of A. A will denote the determinant of A. The dimension of a matrix is the number of its rows by the number of its columns. An n x m matrix A will have n rows and m columns. If m 1, the matrix will be...

## Bibliography

Kinoshita, 'Note on the relation between the equinox and Guinot's non-rotating origin, Celest. Mech., Vol. 29, No. 4, pp. 335-360, 1983. 2 Ashby, N., 'Relativity and the Global Positioning System, Physics Today, Vol. 55, No. 5, pp. 41-17, May 2002. 3 Barlier, F., C. Berger, J. Falin, G. Kockarts, and G. Thuiller, A thermo-spheric model based on satellite drag data, Annales de G ophysique, Vol. 34, No. 1, pp. 9-24, 1978. 4 Battin, R., An Introduction to the Mathematics and...

## Exercises

Date time 31 March 1993 02 00 00.000 (receiver time) x -2107527.21m y 6247884.75 z -4010524.01 PRN21 x 10800116.93m y 23914912.70 z 1934886.67 where these positions are expressed in the same reference frame as T P. Determine the T P GPS receiver clock offset using both satellites. Ignore the GPS transmitter clock correction. State and justify any assumptions. 3. Repeat the example in Section 3.3.1 but use orbit altitudes of The following data apply to Exercises 4-13 The following ephemerides...

## Transformation between ECF and ECI

If the equations of motion given by Eqs. (2.3.46) and (2.3.47) are expressed in an ECI system, such as the J2000 system, the gravitational force given by Eq. (2.3.18) must be transformed from the Earth-iX ed system used to describe the potential into the J2000 system. For the simpliied case of simple ECF rotation about the body-iX ed z axis, the necessary transformation was given in Eq. (2.3.20). In reality, with consideration given to precession, nutation, polar motion, and UT1, the...

## State Transition Matrix and Error Propagation

The initial state for a particular problem, such as the example state vectors in the Section 2.2.4 examples, can be expected to be in error at some level. If the error in position, for example, is Ar(t0), this error will produce an error Ar(t) at a later time t, as well as an error in velocity, Ar(t). The propagation of the error at t0 to an error at t can be approximately expressed by where the matrix of partial derivatives is referred to as the state transition matrix, also known as the...

## T br

This equation verifes Kepler's Third law. From the conic section equation, Eq. (2.2.14), if the angle f is known, the distance r can be determined. However, for most applications, the time of a particular event is the known quantity, rather than the true anomaly, f. A simple relation between time and true anomaly does not exist. The transformation between time and true anomaly uses an alternate angle, E, the eccentric anomaly, and is expressed by Kepler's Equation, which relates time and E as...

## Consider Covariance Analysis

The design and subsequent performance of the statistical estimation algorithms used for orbit determination, parameter identification, and navigation applications are dependent on the accuracy with which the dynamic system and the measurements used to observe the motion can be modeled. In particular, the design of orbit determination algorithms usually begins with the defhi tion of the important error sources and a statistical description of these error sources. The effect of erroneous...

## Linearization Procedure

To illustrate the linearization procedure, return to the flat Earth example discussed previously. Assuming that there are errors in each of the initial position and velocity values and in g, the state vector is XT X, Y, X, Y,g . It is a straightforward procedure to replace X by X * + X, Y by Y* + Y, , and g by g* + g in Eq. (1.2.2) to obtain the equation for perturbation, or deviation, from the reference trajectory. The ()* values represent assumed or specified values. If the equations for the...

## Bj Rj xj

By multiplying Eq. (5.10.27) by Rj, we obtain where nj Rj j. Further, nj satisfes the condition E(nj) 0 and it is easily demonstrated that E(njnJ) I. Now assume that we have the new observation, yj, where and j satisfes the conditions given in Eq. (5.10.18). To obtain a best estimate of Xj, the least squares performance index to be minimized is J(Xj) RjXj - bj 2 + (HjXj - yj)2 . Following the procedure discussed in the previous section, J(xj) can be written as Multiplying by an orthogonal...

## Doppler Systems

A classic example of this technique has been the Navy Navigation Satellite System (NNSS) or Transit, which consisted of a constellation of several satellites in polar orbit. The NNSS began operation in the 1960s with the purpose that the constellation was to support navigation of U.S. Navy submarines. Operations in support of Navy applications ceased in the 1990s with the operational use of GPS. Each satellite carried a beacon that broadcasted signals with two frequencies for ionosphere...

## Maximum Likelihood And Bayesian Estimation

The method of Maximum Likelihood Estimation for determining the best estimate of a variable is due to Fisher (1912). The Maximum Likelihood Estimate (MLE) of a parameter O-gi ven observations yi, y2, yk and the joint probability density function is defiled to be that value of O that maximizes the probability density function (Walpole and Myers, 1989). However, if O is a random variable and we have knowledge of its probability density function, the MLE of O is defiled to be the value of 0, which...

## Minimum Variance Estimate with A Priori Information

If an estimate and the associated covariance matrix are obtained at a time tj, and an additional observation or observation sequence is obtained at a time tk, the estimate and the observation can be combined in a straightforward manner to obtain the new estimate Xk. The estimate, Xj, and associated covariance, Pj, are propagated forward to tk using Eqs. (4.4.19) and (4.4.22) and are given by Xk (tk,tj)Xj, Pk (tk,tj)Pj T(tk,tj). (4.4.23) The problem to be considered can be stated as follows...

## Em 84381448 468150t 000059t2 0001813t3 H32

The expression for the nutation in longitude and in obliquity is given by McCarthy (1996), for example. The IAU 1980 theory of nutation is readily available in a convenient Chebychev polynomial approximation with the planetary ephemerides (Standish, et al., 1997). Corrections to the adopted series, based on observations, are provided by the IERS in Bulletin B. -CemSA > CemCetCA + SemSet CemSetCA - SemCet -SemSA SemCetCA - CemSet , SemSetCA + CemCet

## Two Way Range slr

The basic concept of a two-way range was described in Section 3.3.1. Examples of this measurement type include ground-based satellite laser ranging (SLR) systems and satellite-borne altimeters. SLR was applied for the first time by Figure 3.5.3 The illustrated SLR-2000 is designed for autonomous operation, whereas previous designs have required operators. The SLR concept is based on the measurement of time-of-flght between laser pulse transmission and receipt of echo photons. The telescope...

## TYz Tcf WSN PH11

The transformation matrices W, S', N, and P account for the following effects W the offset of the Earth's angular velocity vector with respect to the z axis of the ECF (see Section 2.4.2) S' the rotation of the ECF about the angular velocity vector (see Sections 2.3.3 and 2.4.2) N the nutation of the ECF with respect to the ECI (see Section 2.4.1) P the precession of the ECF with respect to the ECI (see Section 2.4.1). The ECI and ECF systems have fundamental characteristics described in...

## Sin q cos q

4.16.1 Transformation of the Covariance Matrix between Coordinate Systems Sometimes it is desirable to transform the state vector and the estimation error covariance into alternate coordinate systems. For example, it may be of interest to view these quantities in a radial, transverse, and normal (RTN) system. Here the transverse direction is normal to the radius vector and the normal direction lies along the instantaneous angular momentum vector. RTN forms a right-hand triad. The general...

## W drt drt drt

MMo dr(i) l(i'io) + + 93(t,io)- < l(to,to) < 2 (to, to) 0 3 x3 J 3x3 0 3xml - We could solve Eq. (G.2.5), a 3 x n system of second-order differential equations, instead of the 6x n fist-orde r system given by Eq. (G.2.2). Recall that the partial derivatives are evaluated on the reference state and that the solution of the m x n system represented by < > 3(t, to) 0 is trivial, In solving Eq. (G.2.5) we could write it as a system of n x n first-order equations,

## Iio eAii

(d) Calculate (to,t) by direct inversion. (e) Let -1(i,io) 6(i,io) and show that by integration and comparison with the results of d. (f) Calculate (t2 ,t 1) by finding the product (t2 ,to ) (to ,t 1). (g) Compare this result with the result obtained by integrating the equation (t, t1) A (i, t1), with the condition (i _, i) I. where xT (x x x) and x (xo xco Xo). (9) The equations of motion for a satellite moving in the vicinity of a body with a homogeneous mass distribution can be expressed as...

## Contents

1 Orbit Determination Concepts 1 1.2 Uniform Gravity Field 1.2.1 Formulation of the 1.2.2 The Equation of the 1.2.3 The Role of the 1.2.4 Linearization 1.2.5 State Transition 1.3 Background and 1.6 Exercises 2.1 Historical 2.2 Problem of Two Bodies General 2.2.1 Motion in the 2.2.2 Motion in Space 2.2.3 Basic Coordinate 2.2.4 Orbit Elements and 2.2.5 Position Velocity 2.2.6 State Transition Matrix and Error Propagation 40 2.3 Perturbed Motion 44 2.3.1 Classical Example Lunar Problem 45 2.3.2...

## Gravitational Perturbations Mass Distribution

Consider the gravitational attraction between two point masses, M1 and M2, separated by a distance r. The gravitational potential between these two bodies, U, can be expressed as The gravitational force on M2 resulting from the interaction with Mi can be derived from the potential as the gradient of U that is, where r is the position vector of M2 with respect to M1 given by Figure 2.3.4 Defiiit ion of position vectors and differential mass for a body with arbitrary mass distribution. Figure...

## H6 References

Kinoshita, 'Note on the relation between the equinox and Guinot's non-rotating origin, Celest. Mech., Vol. 29, No. 4, pp. 335-360, 1983. Kaplan, G. H. (ed.), The IAUResolutions on astronomical constants, time scales, and the fundamental reference frame, USNO Circular No. 163, U.S. Naval Observatory, 1981. Lambeck, K., Geophysical Geodesy, Clarendon Press, Oxford, 1988. Lieske, J. H., T. Lederle, W. Fricke, and B. Morando, 'Expressions for the precession quantities based upon...

## Earth Fixed and Topocentric Systems

Since the Earth is not a rigid body, the concept of an Earth-fix ed reference frame is not precisely defiled. As previously noted, the Earth deforms in response to luni-solar gravity, which produces a redistribution of mass. The deformations associated with this redistribution, or tides, produce not only changes in gravity but changes in the coordinates of points on the Earth's surface. In addition, the Earth's crust consists of several moving tectonic plates. As a consequence, a relative...

## H y J ydx

Is the marginal density function of Y. Hence, the marginal density function of a random variable is obtained from the joint density function by integrating over the unwanted variable. A.ll INDEPENDENCE OF RANDOM VARIABLES We have previously defined the independence of two events A and B by p (A, B) p (A) p (B). For the case of random variables X and Y, we say that they are independent if we can factor their joint density function into where g (x) and h (y) are the marginal density functions of...

## Hyxgx h y

The last equations for g (x y) is the most elementary form of Bayes theorem. It is a useful starting point in the development of statistically based estimation criteria. If we defile X to be the state vector and Y to be the observation vector, then g (x y) a posteriori density function h (y x) a priori density function. From a Bayesian viewpoint, we wish to develop a filter to propagate as a function of time the probability density function of the desired quantities conditioned on knowledge of...

## E[ei 0 E[eeT R

A priori estimates for x and c are given by x and c, where x x + n, c c + fl . The errors, n and have the following statistical properties E n E 0 E nnT Px p cc It is convenient to express this information in a more compact form, such as that of a data equation. From Eqs. (6.3.4) and (6.3.7) It follows that the observation equations can be expressed as y Hzz + g e (0, R) We wish to determine the weighted least squares estimate of z, obtained by choosing the value of z, which minimizes the...

## E[nfc 0 E[VkV

From these conditions it follows that k has mean E k k E(HifeXfe + efc - HkXk) P E ( 3fe - k)( k - Jj E 3fe T E (yfe - HfcXfc)(yfe - -HfeXfe)T E (efc- ffkr fcXek- Hiknk)T P Rk + HHkPk Hf (4.7.34) Hence, for a large prediction residual variance-covariance, the Kalman gain will be small, and the observation will have little influence on the estimate of the state. Also, large values of the prediction residual relative to the prediction residual standard deviation may be an indication of bad...

## The Sequential Computational Algorithm

The algorithm for computing the estimate sequentially is summarized as Given Xk-1, Pk-1, XQ-1, and Rk, and the observation Yk, at (at the initial time t0, these would be XQ, X0, and P0). X* F(X*,i), X*(ife-i Xk-i (4717) < & (i,ifc-i) A(i) (i,ifc-i), (ifc_i,ifc_i) J. Xfc (ife,ife-i)Xfe-i Pfe (ife,ife-i)Pfe-i T(tfc, tfe-i)- yfe Yfe-G(X*fe,ife) HHfc dG(g'tfc) - (4) Compute the measurement update Kfe Pfe JiT HfcPfc HT + Rfe -1 (4.7.18) Xfc Xfc + Kfe yfe - HfeXfc Pk J - KfeJife Pfe- (5) Replace...

## Precession and Nutation

Consider the motion of a spherical Moon orbiting the nonspherical Earth. Based on the discussion in Section 2.3, the gradient of the gravitational potential given by Eq. (2.3.12), multiplied by the mass of the Moon, yields the gravitational force experienced by the Moon. By Newton's Third Law, the Earth will experience an equal and opposite force however, this force will not, in general, be directed through the center of mass of the Earth. As a consequence, the Earth will experience a...

## Differenced Measurements

In some cases, the measurements made by a particular technique are differenced in special ways. The differencing of the measurements removes or dimin- Figure 3.5.11 The TOPEX altimeter data collected during a long pass over the Pacific TOPEX Poseidon Repeat Cycle 303 are shown here. ishes one or more error sources. The most common use of such differencing is with GPS pseudorange and carrier phase measurements and with altimeter measurements.

## Gravitational Perturbations Oblateness and Other Effects

Consider an ellipsoid of revolution with constant mass density. It is apparent that Cim and Slm will be zero for all nonzero m because there is no longitudinal dependency of the mass distribution. Furthermore, it can readily be shown that because of symmetry with respect to the equator there will be no odd-degree zonal harmonics. Hence, the gravitational potential of such a body will be represented by only the even-degree zonal harmonics. Furthermore, for a body like the Earth, more than 95 of...

## D3 Lunar Solar And Planetary Masses

Additional parameters required for the description of satellite motion include the gravitational parameters for the Sun, Moon, and planets. The values that are used with the planetary ephemerides, such as DE-405 Standish, et al., 1997 , are given in Table D.2. Additional information can be found in McCarthy 1996 , Seidelmann 1992 , and Standish et al. 1997 . All mass parameters have been determined from observations consult the references for uncertainties. 5 0 0.68658987986543E 07 G 0...

## Y

14 Given a simple pendulum with range observations p , from a fk ed point, as shown in the igure. a Write the equations of motion and form the observation-state matrix H and the state propagation matrix A . Assume the state vector is b Assume small oscillations that is, sinO O, cosO 1. Solve the equations of motion and derive expressions for the state transition matrix. c How does the assumption that O is small differ from a linearized formulation of this problem d Write a computer program to...

## Orbit Elements and Position Velocity

Assume that an artificial satellite orbits a planet whose mass, Mi, is known i.e., p GMi is known . Assume further that a technique and or procedure has been utilized to determine the position and velocity of the satellite at a time, to, in a nonrotating coordinate system. The position and velocity at t0 are ro XoUX YoUY Zouz Vo Xoux YoUY ZoUz . 2.2.32 It follows that the inclination i and node location can be determined from Eqs. 2.2.28 , 2.2.29 , 2.2.30 , and 2.2.33 . The energy per unit mass...

## Example 2261

Assuming an initial error in Example 2.2.4.1 given by the Goodyear 1965 formulation of the state transition matrix can be used to map this error to 30 min, which produces Note that the error at t0 propagated to 30 min can also be obtained by adding the above initial error to the state, which results in X 5492001.34 m Y 3984003.40 m Z 2958.81m and no change in velocity. Following the procedure in Section 2.2.5, the predicted state at 1800 sec is Y 2729258.37 m Z 2973906.50 m Y 6300.787892 m s Z...

## Basic Coordinate Systems

The X, Y, Z coordinate system used in the preceding sections has no preferred orientation for the simple description of the two-body motion. Nevertheless, a specific, well-defined orientation of the X, Y, Z axes is required in practice. Consider the motion of the Earth about the Sun. Since the mass of the Sun is 328,000 times the mass of the Earth, it is quite natural to describe the motion of this system as the motion of the Earth relative to the Sun. But from the description of the two-body...