Gravity modeling in aerospace applications

Aerospace Science and Technology 13 (2009) 301–315

Contents lists available at ScienceDirect

Aerospace Science and Technology www.elsevier.com/locate/aescte

Gravity modeling in aerospace applications George M. Siouris 1 Consultant: Avionics and Weapon Systems

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 4 December 2007 Received in revised form 15 October 2008 Accepted 22 May 2009 Available online 17 June 2009 Keywords: Gravity model Inertial navigation system Markov process Geopotential Zonal harmonics

The science of inertial navigation has evolved to the point that the traditional gravity model is a principal error source in advanced, precise systems. Specifically, the unmodeled vertical deflections of the earth’s gravitational field are a major contributor to CEP (circular error probable) divergence in precise terrestrial inertial navigation systems (INS). Over the years, several studies have been undertaken to the development of advanced techniques for accurate, real-time compensation of gravity disturbance vectors. More complex on-board gravity models which compute vertical deflection components will reduce the CEP divergence rate, but imperfect modeling due to on-board processing limitations will still cause residual vertical deflection errors. In order to eliminate or reduce gravity-induced errors in the INS requires measurement of gravity disturbance values and in-flight compensation to the inertial navigator. It is assumed in this paper that gravity disturbance values have been measured prior to the airborne mission and various techniques for compensation are to be considered. As part of a screening process in this study, several gravity compensation techniques (both deterministic and stochastic models) were investigated. The screening process involved identification of gravity models and algorithms, and developments of selection criteria for subsequent screening of the candidates. © 2009 Elsevier Masson SAS. All rights reserved.

1. Introduction The objective of this modeling study is to identify and trade-off alternative gravity disturbance compensation techniques in terms of system accuracy. A linear optimal estimation algorithm is proposed for use in estimating and compensating an aircraft inertial navigator or inertially guided missile for gravity disturbance errors. The algorithm is examined both at the gravity disturbance residual level and also at the inertial navigation position error level as a function of gravity data, grid spacing, and measurement error. One of the first requirements of formulating an approach to gravity compensation technique development and evaluation is to carefully distinguish between the terms “model” and “algorithm”. A model is a quantitative description of the gravity field. This description may be deterministic or statistical. Given a description or perhaps more than one description, one can then develop a calculation (i.e., an algorithm) based on the model. Many algorithms could be obtained from a given model and, in fact, several very similar algorithms can be applied to different models. However, some models naturally lend themselves to a rather natural algorithm type. For instance, a direct description of deflections of the vertical could be a three-dimensional grid map with two values associated with each grid point (i.e., the north and east deflections).

E-mail address: [email protected]. Formerly, Adjunct Professor, Dept. of Electrical and Computer Engineering, Wright-Patterson AFB, OH, USA. 1

1270-9638/$ – see front matter doi:10.1016/j.ast.2009.05.005

© 2009 Elsevier Masson

SAS. All rights reserved.

Such a model naturally suggests an interpolation algorithm to span the space between grid points; but many types of interpolation algorithms are available. From the above discussion, it is quite clear that there are many models and many algorithms; many models can be used with a given algorithm, and many algorithms can be used with a given model. Gravity error, which is the difference between actual gravity and reference gravity derived from an ellipsoidal potential function, represents a major source of error for high precision airborne inertial navigation systems (INS). That is, there are areas around the earth where gravity anomalies are so significant that they become the dominant error source in these systems. The gravity model used in most present operational systems is based on a gravity field perpendicular to a reference ellipsoid, which approximates the mean-sea level equipotential surface called the geoid. Since the surface of the geoid deviates from the reference ellipsoid, the gravity error also deviates (in magnitude and direction). The deviation of the actual gravity magnitude from the reference value is known as the gravity anomaly, while the deviation of the gravity vector direction from the reference ellipsoid normal is known as the vertical deflection. These parameters are related to the undulation of the geoid (the difference in height between the geoid and the reference ellipsoid). A large class of inertially guided missiles employ a stable platform to maintain a fixed orientation in inertial space by utilizing gyroscopes to sense deviations in orientations. The inertial guidance system signals are obtained by combining the first and second integrals of the output of accelerometers fixed to such a platform.

302

G.M. Siouris / Aerospace Science and Technology 13 (2009) 301–315

Nomenclature ellipsoidal equatorial radius = 6,378.1372 km semi-minor axis of the reference ellipsoid (or polar radius) = 6,356.752314 km  e eccentricity = 1 − (b2 /a2 ) = 0.08209443799496 f flattening = 1 − (b/a) = 1/298.257223563 G universal constant of gravitation = 6.673 × 10−11 m3 /kg s2 h geodetic altitude (or height above the reference ellipsoid) M mass of earth N undulation of the geoid  height r geocentric radius = xe2 + y e2 + ze2 V normal gravitational potential (r , θ, λ) spherical coordinates (x, y , z) Cartesian ECEF coordinate frame

a b

It should be noted that the gravitational effects of the sun, moon and earth’s atmosphere are estimated and found to be negligible for a ballistic missile. Moreover, in reality, the missile cannot be considered as a point object. An observer on the earth measures the acceleration of the missile by observing a particular point on the missile (e.g., the Doppler antenna, tail flame, etc.). Inertial navigation systems in the 1 nm/hr class do not require modeling of the vertical deflection components of the gravity field, since it is not a dominant error source. However, with the advent of more precise inertial navigation systems in the  0.1 nm/hr class makes the vertical deflection of the gravity field a limiting factor in navigation accuracy performance during flight as well as for weapon delivery. For example, a cruise missile trade-off study concluded that it would be desirable for the cruise missile carrier aircraft to use gravity compensation in order to meet the weapon’s delivery accuracy goals. Furthermore, this study involved the development of advanced techniques for accurate real time compensation of gravity disturbance vectors. Having identified the various gravity disturbance models and estimation algorithms, the next step was to identify candidate gravity compensation techniques from that set. Note that the mere combination of a model and algorithm does not necessarily constitute a gravity compensation technique. In identifying those combinations that do qualify, the first priority is to clearly define exactly what constitutes a gravity compensation technique. Since this study is restricted to gravity compensation for INS, whether aided or unaided, it is reasonable to expect that the question of compatibility with an INS have something to do with the definition of a compensation technique. In its broadest sense, an applicable gravity compensation technique for INS estimates the components of the gravity disturbance vector (δ g ) as a function of position (P) indicated by the INS (see Fig. 1). This, of course, presupposes that the INS software has an acceptable model of the normal gravity field and that the purpose of the gravity compensation technique is to account only for the gravity disturbance field. The gravity compensation technique definition as shown in Fig. 1 should be considered as a high level definition and is primarily mentioned here to exemplify its input (i.e., position from the INS) and output (i.e., estimate of gravity disturbance) relationship. However, in actuality, a gravity compensation technique is more involved and can have several parts to it. For example, a database of gravity anomalies may be used to compute the gravity disturbance

2 Unless otherwise specified, all numerical values are from the WGS84, which represents the current best global geodetic reference system of the earth.

deviation (or deflection) of the normal free air gravity (or theoretical gravity on the surface of the ellipsoid) λ geocentric longitude (measured east of the Greenwich meridian) μ = G M earth’s gravitational constant = 3.986004418 × 1014 m3 /s2 ϕ geodetic latitude ϕc geocentric latitude Ω earth rotation rate = 7.292115 × 10−5 rad/s ECEF Earth-Centered, Earth-Fixed INS Inertial Navigation System WGS84 World Geodetic System 1984 nm nautical mile

δ

γ

Fig. 1. Input/output definition of gravity compensation technique.

Fig. 2. Definition of gravity compensation technique.

components in a flight corridor. These results, computed at the intersections (nodes) of a three-dimensional lattice would represent a direct description of gravity disturbance vector and be stored in the flight computer. In flight, the indicated INS position and then some form of an interpolation algorithm would be used to obtain an estimate of the gravity disturbance vector at the desired INS position. Thus, we see that the definition of the complete gravity compensation technique consists of a pre-flight database such as surface gravity anomalies [19,23]. The gravity compensation technique can be defined by the following: 1) pre-flight database, 2) pre-flight computation, 3) inflight database, and 4) in-flight-computation. The complete definition of gravity compensation is illustrated in Fig. 2. Note that the output of the gravity compensation technique is shown as an input to the INS. This illustrates the fact that the gravity compensation technique and the INS form a closed loop system.

G.M. Siouris / Aerospace Science and Technology 13 (2009) 301–315

Statistical characterization of the gravity field is motivated by its short-wavelength features. The long wavelength features of the earth’s gravity field are known accurately from many years of tracking satellites and space probes. Spherical harmonic expansions are a convenient representation of these long-wavelength features, and the geopotential coefficients have been accurately determined out to degree and order 360; this defines the gravity field with wavelength longer than about 3000 km. However, this “known” part of the gravity field represents only a small percentage of the energy in the gravity field at the surface. For example, NASA’s EGM96 (Earth Gravitational Model 1996) recommends using degree and order of 70 for satellite simulation and navigation [11,16]. Consider now another gravity compensation technique where the gravity disturbance field is represented by some series, such as a high order spherical harmonic series or point masses whose contributions are summed to obtain the field at any desired point. In either case, the database may involve one or more data types, which can be used in some form of pre-flight least-squares algorithm to estimate the coefficients of the series or the point masses [13]. These coefficients or point masses are then stored in the flight computer which are then used in the series representation during flight to compute the gravity disturbance vector at any desired position. 2. Fundamental definitions and relations This section reviews the fundamentals of the physics of gravimetric uncertainties and the earth’s gravitational field. No attempt will be made at completeness or rigor; for details, the reader is referred to the text of Heiskanen and Moritz [5]. The discussion following will begin with a somewhat informal presentation of geodetic/gravimetric terminology and the definition of gravimetric uncertainties, followed by a brief review of geopotential theory and selected mathematical relationships between the various gravimetric uncertainties. To begin with, the distinction must be made between the gravitation vector G and the gravity vector g. Briefly, these definitions are as follows: Gravitation. Gravitation is the force exerted on a body due to mass attraction, and it is described by the well-known inverse square law. Gravity. Gravity is the sum of gravitation and the centrifugal force due to the earth’s rotation. In the sequel, the terms “gravitation” and “gravity” will be used in their formal sense, unless otherwise noted. The gravity field of the earth is commonly described in terms of equipotential surfaces; that is, the force of gravity is perpendicular to these surfaces. The equipotential surface at mean sea level is called the geoid. Because of the mass-density variations within the earth and topographic features such as mountains and ocean basins, the surface geoid is irregular. Consequently, since the geoid is an irregular surface, it is neither suitable as a basis for a coordinate system nor as a basis for gravity compensations. Therefore, an analytic (i.e., smooth) reference shape is chosen that closely approximates the geoid. The reference shape universally chosen is an ellipsoid. Because of its use as a position reference, the ellipsoid (also known as spheroid) is of basic importance in many problems of navigation over the surface of the earth [16]. Thus, it seems appropriate at the outset to summarize a number of the properties of the spheroid in this section. Although some of the items covered may seem to be of little practical importance, it is hoped that

303

these parts of the discussion will at least serve to clarify some of the characteristics of this fundamental surface. 2.1. The reference ellipsoid Mathematically speaking, an ellipsoid is a geometrical surface or solid that is symmetrical about its three axes, and the plane sections of which are circles or ellipses. When its axes coincides with the coordinate axes its equation is

x2 a2

y2

+

b2

+

z2 c2

=1

(1)

where the points (±a, 0, 0), (0, ±b, 0), and (0, 0, ±c ) are intercepts with the x, y , z axes respectively. The surface is symmetric to all three coordinate planes. Also, we note that

|x|  a,

| y |  b,

| z|  c

(2)

Replacing z by k in Eq. (1), and writing it in the form

x2 a2 (1 − k2 /c 2 )

+

y2 b2 (1 − k2 /c 2 )

=1

(3)

we can see that the cross-section cut out by z = k is an ellipse. It is interesting to observe that a sphere may be regarded as a special case of an ellipsoid for which a2 = b2 = c 2 . Also, an ellipsoid having two of the quantities a2 , b2 , and c 2 equal is called an ellipsoid of revolution for it may be obtained by revolving an ellipse about one of its axes. Fig. 3 shows in detail the elements of the reference ellipsoid. Two important parameters of the ellipse are the flattening (also known as ellipticity), f , and the eccentricity, e. Mathematically, they are expressed as

f = (a − b)/a 2



2

2



(4) 2

e = a − b /a = 2 f − f

2

(5)

Finally, we note that the center of the reference ellipsoid is the center of the mass of the earth. The reference meridian plane is parallel to the zero meridian adopted by the BIH (Bureau International de l’Heure) on the basis of astronomic longitudes of the BIH. 2.2. Definitions For the benefit of the reader, we will now present some definitions that are central in understanding this paper. Astronomic latitude — The angle, measured north or south (positive or negative) from the equipotential plane to the direction of a plumb-line at the point, i.e., the normal to the equipotential level surface, which in general is slightly distorted from the simple spheroid surface due to local gravitational anomalies. Or, more simply, the angle between the perpendicular to the surface of the geoid and the plane of the equator (also defined as the elevation of Polaris above the horizon). Astronomic longitude — The angle between the plane of the meridian at Greenwich3 (Prime Meridian) and the astronomic meridian of the point. Deviation of the normal — The deviation of the normal, δ , is commonly defined as the angle between the geocentric and geographic verticals. Thus, from Fig. 4, δ can be expressed in the form

  δ = ϕ − ϕc = f 1 − (h/ R 0 ) sin(2ϕ )

(6)

where 3 Longitude is measured east and west of the Prime Meridian, which is defined to pass through the sight of the Royal Observatory in Greenwich, England.

304

G.M. Siouris / Aerospace Science and Technology 13 (2009) 301–315

Fig. 3. The reference ellipsoid.

Fig. 4. a) Geodetic parameters and deviation of the normal at altitude. b) Relationship between elevation H above the geoid, ellipsoid height h, and the geoid height (undulation) N above the ellipsoid [7,26].

ϕ = geographic latitude ϕc = geocentric latitude f = ellipticity = 1 − (b/a) h = altitude above the reference ellipsoid R 0 = geocentric position Figure of the earth — Figure of the earth refers to its qualities of shape and size. The overall shape of the earth (disregarding topographical features) is determined by the resultant of gravitational force and the centrifugal force associated with the diurnal rotational speed. In order to define the earth geometrically on terms sufficiently exact for navigation theory and practice, it is necessary to introduce the concepts of the geoid and the reference ellipsoids. Geoid — The geoid is the figure of the earth considered as a mean sea level surface extended continuously through the continents. It is a gravity equipotential surface; anything placed on it does not tend to move toward the equator or the poles, but only downward, i.e., normal to the surface. It is not strictly a mathematically simple surface because of the earth’s inhomogeneous mass distribution. It is slightly bulged, for example, over a local mass of exceptionally dense rock. Moreover, it is not strictly a fixed surface because of tides in the oceans and in the earth mass itself. The geoid is a true equipotential surface and a plumb bob would

hang perpendicular to the surface of the geoid at every point. From Fig. 3 we note that a line normal to the elliptical curve at any point can be constructed by bisecting the angle between lines from the foci to the point (i.e., θ1 = θ2 ). The direction of this line is termed the geodetic vertical. A plane normal to the line is the level or horizon plane. Geocentric latitude — The angle between the equatorial plane and the radius vector from the geocenter to a point on the surface of the ellipsoid. The transformation from geocentric latitude, ϕc , is related to the geodetic latitude, ϕ , can be computed via the relation





tan ϕc = 1 − e 2 tan ϕ or





ϕc = tan−1 1 − e2 tan ϕ

(7)



(8)

Furthermore, it can be shown from Eq. (6) that

tan(ϕ − ϕc ) =

[e 2 /(2 − e 2 )] sin 2ϕ 1 + [e 2 /(2 − e 2 )] cos 2ϕ

(9)

so that by using the Fourier expansion

tan−1



x sin θ 1 − x cos θ

 =

∞ n x sin(nθ) n =1

n

(10)

G.M. Siouris / Aerospace Science and Technology 13 (2009) 301–315

we have

ϕc = ϕ −

2

e sin 2ϕ 2 − e2

−

e

2

2 − e2

2

sin(4ϕ ) 2

− ···

(11)

The difference between ϕ and ϕc reaches a maximum value of 11.6 minutes of arc at about 45◦ latitude. There, the difference between geodetic and geocentric latitude corresponds to a meridional displacement of about 20.5 km or 11.6 nm a quantity which must be taken into account even in crude forms of navigation. Geodetic latitude — Geodetic latitude, ϕ , also known as geographic latitude, is the angle between the equatorial plane and the normal to the surface of the ellipsoid. This is the latitude used almost universally in charts, maps, and lists of geographical positions for navigation purposes.

ϕ = tan−1



tan ϕc



(12)

(1 − f )2

Geodesic — A geodesic is the shortest curve between two points on a curved surface lying wholly on the surface. On the surface of a sphere, a geodesic is an arc of a great circle. On a plane, a geodesic is a straight line. Mathematically, the length of any space curve is

S=



dx2 + dy 2 + dz2

(13)

Therefore, in order to find the shortest lines (i.e., geodesics) on any surface, the integral must be minimized subject to the condition F (x, y , z) = 0. Gravity anomaly — Difference between the observed gravity at some point on the earth and the theoretical value of gravity at the same point. Gravity anomalies are caused by local or regional mass variations. Gravity anomaly (free air) — The difference between the observed and theoretical gravity which has been computed for latitude and corrected for elevation of the station above the geoid, by application of the normal rate of change of elevation, as in free air. Note that the elevation correction is for the height above the geoid, not the ellipsoid. Thus, the expression for the resulting freeair reduction in the neighborhood of the reference geoid is given by

γ = γ s + 3.086 × 10−6 h

(14) s

where h is the height above the surface, and γ is the gravity on the surface of an ellipse. International gravity formula — The 1980 Geodetic Reference System gives the international gravity formula on the surface of an ellipsoid by the expression [11,14]





g s = − ge 1 + 0.00523024 sin2 ϕ − 0.0000058 sin2 2ϕ m/s2 (15) where g e = 9.780327 is the surface gravity at the equator and ϕ is the geodetic latitude. The variation from the equator to the poles is 0.27%. Rhumb line (or loxodrome) — A line on the earth’s surface which intersects all meridians at the same angle. Any two places on the surface of the earth may be connected by such a line. (For derivation of the mathematical expression of a loxodrome, see [20].) Undulation (or geoidal undulation) — The difference between the reference ellipsoid and the geoid.

305

in applications (e.g., ballistic missile guidance, high precision INS, etc.) a more accurate theoretical closed-form expression for the gravity on the surface of the earth is needed. Therefore, the theoretical gravity, γ , on the surface of the ellipsoid described above is given by the equation of Somigliana [11] in the form

1 + k sin2 ϕ

γ = γe

(16)

1 − e 2 sin2 ϕ

where

k=b

γp −1 aγe

(17)

In Eqs. (16) and (17), a and b are the semi-major and semiminor axes of the ellipsoid, respectively; γe and γ p are the theoretical gravity at the equator and poles, respectively; and ϕ is the geodetic latitude. Or more simply, we can define gravity flattening in the form [5]

fg =

γp −1 γe

(18)

The DMA (Defense Mapping Agency) in St. Louis, Missouri, developed an equation for calculating the magnitude of gravity on the surface of the WGS84 reference ellipsoid. The value of gravity, referred to as normal gravity γn , is calculated as follows [16]:



γn = 9.7803267715 1 + 0.001931851353 sin2 ϕ



−1/2  m/s2 × 1 − 0.0066943800229 sin2 ϕ

(19)

where ϕ = geodetic latitude. Finally, and as mentioned above, high accuracy guidance system (e.g., aircraft, cruise missiles, air-to-air missiles, and intercontinental ballistic missiles, etc.), demand the best possible gravity model the designer can provide. For more information on gravity requirements, the reader is referred to [12] and Chapters 6 and 7 of [21]. 2.4. Coordinate systems In the last few years, and with the availability of a full GPS constellation, the GPS/INS integration for aircraft and surface vehicle navigation has found widespread use. It is well known that the position errors in an inertial navigation system tend to increase with time in an unbounded manner. One way of halting this growth and bounding the error is to update the inertial system periodically with external position fixes. Such position fix data is commonly derived from GPS signals, ground tracking radars, laser range finders, etc. The effect of the position fixes is to reset the position error in the inertial system to the same level of accuracy inherent in the position fixing technique selected by the user. The advantages of combining the long-term accuracy of the GPS with the short-term stability of an autonomous INS, in an integrated mode, are well known. Several AF programs that demand high accuracy mentioned in the introduction, and high quality navigation, such as the F-16 navigation system, the Combat Talon II (MC-130H aircraft), etc., use integrated GPS/INS systems. Inertial sensor (i.e., accelerometers and gyroscopes) errors, gravity uncertainties, baro-inertial altimeter errors, GPS receiver errors, etc., are commonly modeled as random constants, random walks, and firstorder Markov processes. At this point, it is appropriate to develop and/or define certain expressions that are useful in navigating at or near the surface of the earth. Specifically, the geodetic, inertial, and local-level coordinate systems will be briefly discussed.

2.3. Gravity on the surface of the ellipsoid In many navigation applications, approximate formulae have been used to calculate the gravity effects of the earth. However,

(a) Geodetic coordinates The geodetic ellipsoidal coordinates are convenient for terrestrial navigation because the latitude is a critical variable, fixed

306

G.M. Siouris / Aerospace Science and Technology 13 (2009) 301–315

site locations (e.g., navigation waypoints, runways, ground stations, etc.), and gravity is directed essentially along the ellipsoidal normal. Moreover, because of the particular advantages of each coordinate system, aircraft and spacecraft inertial navigation algorithms frequently require transformations between geodetic ellipsoidal and geocentric Cartesian coordinates. It should be pointed out, however, that earth-centered-inertial (ECI) coordinates are more convenient for orbital navigation. In geodesy, however, the ellipsoidal or geodetic coordinates, ϕ (latitude), λ (longitude), and h (altitude) are commonly used. We are therefore interested in relating these ellipsoidal coordinates to the rectangular coordinates. Fig. 4 shows the meridional section through the rotation ellipsoid. The x-axis is the BIH zero meridian plane, the z-axis coincides with the reference ellipsoid b-axis, and the y-axis points 90◦ east of the positive x-axis. From Fig. 4, the conversion from geodetic ellipsoidal coordinates (ϕ , λ, h) into Cartesian coordinates in an ECEF coordinate frame (x, y , z) can be realized by the following algorithm [7,24, 25]:

x = ( R m + h) cos ϕ cos λ 2

Ω = kΩ, where k is a unit vector in z-direction r = position relative to the center of the earth In a nonrotating inertial coordinate frame, gravitation g constitutes the significant gravitational influence. For instance, satellite motion is best described in terms of gravitation. In an earth-fixed coordinate frame (or one related to it, such as a local-level), gravity g  constitutes the significant gravitational influence. Aircraft motion, for example, is the best described in terms of gravity [4,20]. (b) Inertial coordinates Gravitation g, is comprised of a basic inverse-square, central force field plus higher-order terms representing the nonspherical mass distribution of the earth. These higher-order terms represent expansion of the earth’s gravitation potential in spherical harmonics, the associated coefficients being deduced from their observed effects on satellite orbits. For inertial navigation applications, the geocentric coordinates G c = {−G ϕ , 0, −G r } must be transformed into the inertial frame as follows:

G i = C ci G c

y = ( R m + h) cos ϕ sin λ



Ω = earth’s angular velocity (Ω = 7.292115 × 10−5 rad/sec);



z = R m + h − e R m sin ϕ

(20)

where

⎡

where

R m = meridian radius (also known as north–south radius of curvature)

   3/2 = a 1 − e 2 / 1 − e 2 sin2 ϕ

(24)

C ci

−sin ϕc cos λ −sin λ −cos ϕc cos λ

= ⎣ −sin ϕc sin λ

cos λ

⎤

−cos ϕc sin λ ⎦

cos ϕc 0 ϕc = geocentric latitude

−sin ϕc

λ = celestial longitude

e 2 = eccentricity squared = 0.00669437999

(c) Local-level coordinates In local-level coordinates, gravity is computed quite simply as

h = geodetic height R n = prime vertical radius of curvature (also known as

g l = [0 0 − g ] T

(25)

the east–west radius of curvature)

= a/ 1 − e 2 sin2 ϕ

In inertial navigation system analysis, using a spherical earth model, the radius of curvature R s is given by

R s = ( R m + R n )/2

at ϕ = 45◦

(21)

In order to eliminate the computation of the square root in the east–west radius, R m , the following approximation can be used:



Rm ≈ a 1 +

e 2 sin

2



ϕ

2(1 − 0.75e 2 sin2

ϕ)

(22)

If the latitude varies from 0◦ to 90◦ , then the error is less than 0.2 ft (0.06096 m). As we shall discuss in more detail in the next section, gravitation represents the mass attraction existing between the earth and another object. For the purpose of inertial navigation, gravity represents the gravitational acceleration apparent in a coordinate frame rotating with the earth, that is, the sum of actual gravitation plus centripetal acceleration due to earth rotation. The relation between gravitation g and gravity g  is simply [2]

g  = g − Ω × (Ω × r ) where 

(23)

For more details on coordinate transformation, the reader is referred to [1,19]. Finally, we note that quaternions are widely used in strap-down inertial navigation systems due to their singularity-free characterization of orientation of the inertial sensor axes relative to the local geodetic frame. For more details on coordinate transformations, the reader is referred to [4,20]. 2.5. Newton’s law of universal gravitation If the earth were perfectly spherical and of uniform density, then the earth’s gravity would behave as if the earth were a point mass. Let an object of mass m be located at position vector r in an ECI coordinate system. Then, according to Newton’s universal law of gravitation, the force F G of attraction between two bodies m and M at a distance r between their center-of-mass will be [4]

F G = ma = −

G Mm r3

r

where

a = acceleration of object m = mass of object M = mass of the earth

g = gravity acceleration

G = universal gravitational constant

g = gravitational acceleration

r = separation between their center of mass = |r |

(26)

G.M. Siouris / Aerospace Science and Technology 13 (2009) 301–315

The negative sign on the right of Eq. (26) results from the fact that gravitational forces are always attractive. Eq. (26) can also be written as

d2 r

=−

dt 2

μ r

r 3

(27)

from the positive x-axis. For the reference ellipsoid model, symmetry exists about the polar axis ze . Assuming a spherical harmonic expansion of the world’s gravitational potential, we have [4,20]

V (r , ϕc , λ) =

Note that Eq. (27) is the well-known expression for the socalled “two-body problem”. Furthermore, if the function V measures the true gravitational potential of the earth at an arbitrary point in space, then Eq. (27) can be written as follows:

α=

d2 r dt 2

∇

μ

(28)

⎡ ⎤



x

=−

r

μ⎣ ⎦ μ y =− r

r3

(29)

r3

z

2.6. The earth’s gravitational field We begin this section by defining the gravitation field, g. The gravitation field g, is a vector field which may be expressed as a scalar function, called the normal gravitational potential, V . Furthermore, the gravitational potential is a function of position, r, relative to the attractive body. Thus,

g = ∇ V (r )

(30)

where the operator ∇ is given by ∇ = In Cartesian coordinates,



∂V ∂x

∂V ∂y

∂V ∂z

∂ ∂x i

+

∂ ∂y j

+

∂ ∂ z k.

1+

r

n n ∞ a n =2 m =0

r

P n,m (sin ϕc )



 × C n,m cos(mλ) + S n,m sin(mλ) 

(34)



(31)

∇ × g = 0 and ∇ • g = 0

(32)

The cross product in Eq. (32) demonstrates the symmetry of the gradient tensor, while the dot product in Eq. (32) is Laplace’s equation, requiring that the sum of the diagonal components of the gradient tensor be zero. The normal gravitational potential V is the solution of the Laplacian 2

2

2

∂ V ∂ V ∂ V + + =0 2 2 ∂x ∂y ∂ z2

V = gravitational potential a = semi-major axis of the reference ellipsoid r = geocentric radius =

xe2 + y e2 + ze2

n = degree of harmonic term m = order of the harmonic term C n,m , S n,m = harmonic coefficients

degree m with an argument of sin(ϕc ) Eq. (34) consists of sectoral harmonic (i.e., n = m) and Tesseral harmonics (n = m). The former deals with mass distributions depending on longitude, while the latter depends on mass distributions in both latitude and longitude. The harmonic terms, that is, m = 0, are called zonal harmonics, and are longitude independent. For the reference ellipsoid model, symmetry exists about the polar axis, ze . Moreover, the potential no longer depends on the earthreferenced longitude, λ. Therefore, ignoring order higher than zero (i.e., m = 0), the earth’s gravitational potential at a distance r from the center of mass is defined by the equation [10,19]

V (r , ϕc ) =

Therefore, for a conservative gravitation force,

∇2 V =



P n,m (sin ϕc ) = associated Legendre polynomial of order n and

Finally, note that Newton’s law computes only gravitation, not gravity (see next section).

gT =

μ

where

= ∇V

where V = G M /r = μ/r. Moreover, as it will be demonstrated in the next subsection, the earth’s gravitational potential is modeled by a spherical harmonic series. Therefore, for the two-body motion, V = μ/r. Now, Eqs. (27) and (28) are equivalent. Thus,



307

(33)

with the boundary condition V = V 0 = constant for points on the ellipsoidal surface S 0 (see also Eqs. (1) and (3)): [(x2 + y 2 )/a2 ] + (z2 /b2 ) = 1. 2.7. The gravitational potential The WGS-84 gravitational potential (for a two-body motion) at point P , is evaluated at the system location specified by the geocentric position vector, r, having the general spherical coordinates (r , θ, λ), arising from the gravitational effects of the distributed mass of the earth. In turn, this coordinate system can be transformed into the rectangular (x, y , z) system. Now, let θ be the angle of colatitude, and λ the longitude measured anticlockwise

μ r



1+

∞ n a n =2

r



C n P n (sin ϕc )

(35)

where the terms in Eq. (35) have been defined above. Furthermore, we note that in Eq. (35) the even harmonics are symmetric about the polar axis giving rise to the oblate terms, while the odd harmonics are anti-symmetric giving rise to the so-called pear-shaped term. In this case, the Legendre polynomials (also known as surface zonal harmonics), P n (ζ ) (n = 0, 1, 2, . . .) can be defined by the Rodrigues formula [3,18]

P n (ζ ) =

1 dn (ζ 2 − 1)n 2n n!

dζ n

,

P n (−ζ ) = (−1)n P n (ζ )

(36)

where ζ denotes sin ϕc . Furthermore, the Legendre polynomials satisfy the following identities:

P n (1) = 1,

P n (−1) = (−1)n ,

and

     P n (ζ )  1 |ζ |  1

Each harmonic term in the potential is due to a deviation of the potential from a uniform sphere. In the present analysis and for the purpose of inertial navigation, the second-, third-, and fourth-order terms will be considered. The gravitational acceleration is normally obtained by taking the gradient of the gravitational potential. Thus, taking the gradient of Eq. (35) in Cartesian coordinates, and writing out Eq. (35) for the potential through the C 4 term, results in [1,4,16]:

x¨ = −μ

+

x r3



1+

3C 2 a2 (5z2 − r 2 ) 2r 4

−

5C 4 a4 (3r 4 − 42z2 r 2 + 63z4 ) 8r 8

5C 3 a3 z(3r 2 − 7z2 )



2r 6

308

G.M. Siouris / Aerospace Science and Technology 13 (2009) 301–315

y¨ = x¨

y x

z¨ = −μ

+

z

 1−

r3

3C 2 a2 (3r 2 − 5z2 ) 2r 4

−

5C 4 a (15r − 70z r + 63z ) 4

4

2 2

4

C 3 a3 (30z2 r 2 − 35z4 − 3r 4 ) 2zr 6



(37)

8r 8

where x, y , z = earth-fixed coordinate system. The coefficients C n are determined experimentally from observation of artificial satellite orbital observations. These values are as follows [16]:

C 2 = −1.08262668355 × 10−3

C 3 = 2.53265648533 × 10−6

C 4 = 1.61962159136 × 10−6 It should be noted that in Refs. [4,16] the zonal harmonic coefficients are given by the symbol J n instead of C n,0 , and are published with opposite sign to the C n,0 coefficients. Eq. (37) is the desired analytical expression for the earth’s gravitational potential, and it has been derived under the following assumptions: a) the potential is evaluated at a point that is external to the mass of the earth, b) the mass distribution of the earth is symmetric about its polar axis, and c) the center of the mass coincides with its geometric center and with the origin (or center) of the earth coordinate frame. Ref. [5] shows that the second degree zonal harmonic C 2,0 is a function of the world’s mass, equatorial radius, and moment of inertia. However, C 2,0 cannot be measured directly and can only be empirically derived. Thus, modeling the geopotential using only the first two even zonal harmonics, we have

V (r , ϕc ) =

μ r

+



2

1+

4 a

a

C 2,0

r

2 4

C 4,0

r

(3 sin2 ϕc − 1)

(35 sin ϕc − 30 sin2 ϕc + 3)



8

(38)

In general, making use of Eq. (36) and expanding to the first five terms (note that only the P 1 through P 4 terms are used in computing the geopotential), results in

P 0 (sin ϕc ) = P 0 (ζ ) = 1 P 2 (ζ ) = P 4 (ζ ) =

3ζ 2 − 1

P 1 (sin ϕc ) = sin ϕc

P 3 (ζ ) =

2 35ζ 4 − 30ζ 2 + 3 8

5ζ 3 − 3ζ 2 P 5 (ζ ) =

63ζ 5 − 70ζ 3 + 15ζ 8

where the first two terms (i.e., P 0 and P 1 ) are by definition. From Eq. (36) we see that, since the nth derivative of a polynomial of nth degree is a constant, we have



P nn (x) = C 1 − x2

m/2

= C sinn θ,

and

P nm (x) = 0 for m > n (39)

The nth polynomial which is Rodrigues’ formula, solves Legendre’s differential equation [3]





1 − x2 y  − 2xy  + p ( p + 1) y = 0

whose solution, when p is a natural number, is the pth Legendre polynomial. Finally, taking the gradient of Eq. (38), we have

γr =



2   ∂V μ 3 a =− 2 1+ C 2,0 3 sin2 ϕc − 1 ∂r 2 r r 

4   5 a + C 4,0 35 sin4 ϕc − 30 sin2 ϕc + 3 8

r

Fig. 5. Geopotential and spheropotential surfaces.

γϕ =

1 ∂V r ∂ ϕc

2

= cos ϕc sin ϕc



× 3C 2,0 +

1 2

2 a r

μ a r2

r



C 4,0 35 sin2 ϕc − 15





(40)

It should be noted that the zonal harmonic coefficients are published as J n and are of opposite sign to the C n,0 coefficients [16]. In Eq. (38) the C n,0 coefficients are as follows:

C 2,0 = −1.08262982131 × 10−3 C 4,0 = 2.37091120053 × 10−6 The total “actual” potential W of the earth at a boundary or external point is the sum of the gravitational potential (i.e., mass attraction), V , due to mass attraction, and the potential due to rotation (centrifugal force, non-field force), and is given by the expression [7,11]

W =V +

Ω 2 r 2 cos2 ϕc

(41)

2

where Ω is the earth’s rate of rotation, r is the magnitude of a radius vector from the center of the earth to the point in question, and ϕc is geocentric latitude (i.e., latitude measured with respect to a line through the center of the earth). The total potential, W , is considered to be a “normal” gravity potential, U , plus small irregular variations resulting in an anomalous potential or disturbing potential, T , responsible for the fine structure of the earth’s gravity field: thus

W =U +T

(42)

The total gravity vector, g, the normal gravity vector, gravity disturbance vector δ are given as [14,15]:

g = grad W

γ = grad U

δ = grad T 4

γ , and the (43)

Compare now the geoid W (x, y , z) = W 0 with the reference ellipsoid U (x, y , z) = U 0 of the same potential U 0 = W 0 . A point P of the geoid is projected onto the point Q of the ellipsoid by means of the ellipsoidal normal, as shown in Fig. 5. The distance P Q between geoid and ellipsoid is called the geoidal height, or geoidal undulation, and is denoted by N. The undulation N is related to the anomalous potential T by the well-known Bruns formula [5,8,9,26]

N = T /γ

(44)

4 The reader should note that the δ given in Eq. (43) is not the same as the δ in Eq. (6).

G.M. Siouris / Aerospace Science and Technology 13 (2009) 301–315

309

where the negative sign is chosen so as to be consistent with Eq. (47). The normal component of the gravity disturbance vector δ g is given as

δg = −

∂T (x, y , z) ∂z

(52)

The gravity anomaly, g, is given in terms of the anomalous potential, T , by the fundamental equation of physical geodesy as [5]

g = cT (x, y , z) −

∂T (x, y , z) ∂z

(53)

where

1 ∂γ γ ∂z

Fig. 6. Deflection of the vertical as illustrated by means of a unit sphere with center at O .

c=

where γ is the magnitude of the normal gravity. Consider now the gravity vector g at P and the normal gravity vector γ at Q . The gravity anomaly vector, g, is defined as their difference:

can be considered to be a constant for all practical purposes. The disturbing potential, T , satisfies the Laplace equation above the surface of the earth, that is,

g = g p − γ Q

∇2 T =

(45)

where

  ∂T ∂  2  ∇ T =0 = c∇ 2 T − ∇ 2 ( g ) = ∇ 2 cT − ∂z ∂z

γ Q = theoretical (normal) acceleration of gravity   δ = cos−1 (γ Q • g p )/|γ Q || g p |

A vector is characterized by magnitude and direction. The difference in magnitude is the gravity anomaly

g = gp − γQ

(46)

and difference in direction is the deflection of the vertical (see also Section 2.2). The deflection of the vertical is defined as the angular difference between the true vertical (perpendicular to the geoid) and reference vertical or geodetic normal (perpendicular to the reference ellipsoid). This angle is commonly broken up into two components, that is, a meridian (north–south) component, ξ , and a longitudinal (east–west) component, η , as shown in Fig. 6. These components are given by [20,23]

ξ = ϕa − ϕ

η = (λa − λ g ) cos ϕ

(47)

where ϕa and λa are the astronomic latitude and longitude, respectively, and ϕ and λ g are the geodetic latitude and longitude, respectively. Basically, the geoid, geodetic latitude, and longitude define position with respect to the reference ellipsoid. It is also possible to compare the vectors g and γ at the same point P . Then we get the gravity disturbance vector

δ = gp − γ p

(48)

accordingly, and the difference in magnitude is the normal component of the gravity disturbance vector

δ g = gp − γ p

(49)

The difference in direction, that is, the deflection of the vertical, is the same as before, since the directions of γ p and γ Q practically coincide. For a local coordinate system (x, y , z) with x pointed towards the east, y towards north, and z towards the vertical up, the deflections of the vertical ξ and η are connected to the disturbance potential T as [2,22,23]

ξ =−

η=−

(55)

Now, using Eqs. (52), (53), (54), and (55) yields

g p = actual (or observed) acceleration of gravity



∂2T ∂2T ∂2T + + ∂ x2 ∂ y2 ∂ z2

(54)

1 ∂T

γ ∂y 1 ∂T

γ ∂x

(x, y , z)

(50)

(x, y , z)

(51)

(56)

where c was assumed to be constant. Thus, gravity anomaly itself satisfies Laplace’s equation in rectangular coordinates. 3. Survey of gravity models This section describes the various models of potential gravity compensation methods that were examined to identify candidate techniques for use in INS gravity compensation. Moreover, such models can be used in the design of a Kalman filter for representation of the effect of the gravity disturbance vector. Specifically, the gravity model criteria covered aspects of maturity, screenability, computability, optimality, flexibility, versatility, sensitivity, stability, and integrity. From the literature search and classification process, a set of twelve gravity models were identified as “deterministic” and fifteen as “statistical.” The list is by no means exhaustive. There are many models that will not be discussed in this paper. For more details, the reader is referred to [8,17,22]. 3.1. Deterministic gravity model Gravity models in this category are represented by a list of numerical values, a set of deterministic functions and/or a set of precalculated coefficients. It should be pointed out that although statistical methods may be used to arrive at the model from a real gravity database, the end result is a non-statistical, that is, deterministic model. The twelve deterministic gravity models are: 1) Direct Description, 2) Spherical Harmonics, 3) Ellipsoidal Harmonics, 4) Binary Sample Functions, 5) B-Function Representation, 6) Multiquadratic Functions, 7) Chebychev Polynomials, 8) FiniteElement, 9) Molodenskii Coefficients, 10) Surface Coating, 11) Point Mass, and 12) Bicubic Spline. The following paragraphs discuss briefly four of the models, that is, the direct description, spherical harmonics, finite four element, point models, and statistical gravity models. Where possible, the mathematical representation of the model is given. However, in some cases the model is so mathematically involved that its presentation will not be commented in this paper. In other cases the literature did not divulge sufficient detail to succinctly describe the model.

310

G.M. Siouris / Aerospace Science and Technology 13 (2009) 301–315

ξ( λ, Φ, 1) = (1 − Φ)ξ( λ, 0, 1) + Φξ( λ, 1, 1) (62) f) Finally, we estimate ξ at the indicated position by linearly interpolating between the points ( λ, Φ, 1) and ( λ, Φ, 1)

ξ( λ, Φ, h) = (1 − h)ξ( λ, Φ, 0) + hξ( λ, Φ, 1) Fig. 7. Aircraft location within parallelepiped.

1) Direct Description This model, which is one of the most widely used, stores gravity disturbance data in terms of three-dimensional space lattices. These data may be the disturbance potential, undulation, deflections, anomalies or gravity disturbance components. Thus, each data point may reference up to three disturbance values, for example, three components of the gravity disturbance vector. Based on pre-flight data, the size of the space lattice will have to not only encompass the anticipated trajectory, but also include sufficient boundary to accommodate possible navigation errors. During flight, the navigator’s indicated position would be used to identify the elemental parallelepiped. An interpolation scheme would then be used to obtain the gravity disturbance data. The model under discussion utilizes a three-dimensional space lattice of north and east gravity disturbance components. The disturbances at the INS indicated present position are found by three-dimensional linear interpolation between the nearest lattice points. The disturbances are applied directly to the level axis accelerometer output as corrections. Specifically, suppose after suitable determination of present position, the aircraft is identified to be within a particular parallelepiped. Referring to Fig. 7, the corners are numbered where the data is available and the aircraft location is defined by incremental longitude ( λ), latitude ( Φ ), and altitude ( h) within the parallelepiped. That is, λ, Φ , and h represent functional distances. The linear interpolation proceeds as follows: a) Consider the plane formed by corners 1, 2, 3, and 4 and label the data at these points ξ1 , ξ2 , ξ3 , and ξ4 respectively (see Fig. 7). b) Using linear interpolation between points 1 and 2, we estimate the value ξ at the point indicated by ( λ, 0, 0) as

ξ( λ, 0, 0) = (1 − λ)ξ1 + λξ2

ξ( λ, 1, 0) = (1 − λ)ξ3 + λξ4

A similar process is used to estimate the second component at the indicated present position. 2) Spherical Harmonics In Section 2, the earth’s gravitational potential was derived using a spherical harmonic expansion. In this section we will discuss in some more detail the role the spherical harmonics play in connection with the gravitational potential. We begin by noting that the Laplacian potential function V = 0 leads to another class of function. Specifically, we seek such solutions which are homogeneous and of degree n in the coordinate system (x, y , z). These solutions, expressed in three-dimensional polar coordinates, are of the form

V = r n F n (θ, φ)

(58)

d) Again, using linear interpolation between the points (λ, 0, 0) and ( λ, 1, 0), we estimate the value of ξ at the point ( λ, Φ, 0) as

ξ( λ, Φ, 0) = (1 − Φ)ξ( λ, 0, 0) + Φξ( λ, 1, 0) (59) e) Considering the plane formed by the corners 5, 6, 7, and 8, we obtain (essentially repeating steps b, c, and d)

ξ( λ, 0, 1) = (1 − λ)ξ5 + λξ6

(60)

ξ( λ, 1, 1) = (1 − λ)ξ7 + λξ8

(61)

(64)

Assuming independence of φ , a simpler equation results by setting cos θ = x and F n = y:

d dx





1 − x2

 dy dx



+ n(n + 1) y = 0

(65)

Now, since the anomalous potential T = W − U is a harmonic function, it can be expanded into a series of spherical harmonics [11,15]:

T (r , θ, λ) =

∞ n+1 R n =0

r

T n (θ, λ)

(66)

where r , θ, λ are the spherical coordinates related to the rectangular coordinates by

x = r sin θ cos λ y = r sin θ sin λ z = r cos θ and

(57)

c) Similarly, we estimate the value of ξ at the point ( λ, 1, 0) as

(63)

T n (θ, λ) =

n 



A nm R nm (θ, λ) + B nm S nm (θ, λ)

(67)

m =0

These are fully normalized harmonics defined by:

R nm (θ, λ) = P nm (cos θ) cos mλ S nm (θ, λ) = P nm (cos θ) sin mλ where

P no = 2−n

√

2n + 1

 (−1)k k =0

(2n − 2k)! × t n−2k for m = 0 k!(n − k)!(n − 2k)! and

(68)

G.M. Siouris / Aerospace Science and Technology 13 (2009) 301–315

P nm = 2

−n

×



1−t

 2 m/2

 (−1)k k =0

 2(2n + 1)

(n − m)! (n + m)!

(2n − 2k)! k!(n − k)!(n − m − 2k)!

which makes the approximation unbiased. Each function u i , j ,k is expressed as a series of basis functions. Dropping the indices i , j , k for simplicity, assume

t n−m−2k

for m = 0

u=

In the above expressions,  is the greatest integer  (n − m)/2. A nm and B nm are coefficients obtained from anomalous potential data. Furthermore, it should be recognized that

P nm = 2

×



1−t

 2 m/2

 (−1)k k =0



2(2n + 1)

(n − m)! (n + m)!

(70)

and applied to gravity compensation. In addition, the infinite series must be terminated at some finite value. When the coefficients A nm and B nm are statistically obtained from gravimetric data, this model is called a “Degree Variance” model. Degree variance gravimetric models have been extensively used by Moritz [14,15]. 3) Finite Element Model The Finite Element model is a scheme, which may be used to obtain either a local or global representation of the geopotential. The history of its development appears to be that of a three-dimensional curve fitting method which was originally created for application to cartography and later expanded to geodesy [7]. In this model, the first task is to select a volume, either global or local, for representation. An example of a global volume is a spherical shell covering the surface of the Earth and extending from the surface out to 1.2 earth radii. An example of a local volume is a region with a base on the surface of the earth of 10◦ × 10◦ , and an altitude of 300 km. The next task is to partition this volume into cells or volume elements of finite dimensions, hence the name, finite element method. Within each volume element or cell, the geopotential is represented by a weighted sum of polynomials. (Note that the number of elements used to partition the total volume varies from 100 to 1500.) The chosen scheme imposes stringent continuity conditions on the potential and its derivatives, both within a given cell and from one cell to the next. The shape of the cell is generally chosen to be a rectangular parallelepiped, so each cell has eight vertices. Now, denoting the disturbing potential by u, then each cell is represented by

u =

ωi, j,k u i, j,k

(71)

Each function u i , j ,k (i , j , k = 0, 1) is an approximation of u, which is quite accurate at one of the cell vertices. Each function ωi, j,k (i , j , k = 0, 1) is a weighting function which blends the various approximations together to obtain a representation inside the cell. These weighting functions are chosen to satisfy the constant

i =0 j =0 k =0

ωi, j,k = 1

where the coefficients c α β γ are determined from measurements of the value of u in the specific cell in question using standard least squares method (it is assumed that more measurements are available than the number of coefficients to be determined), and the basis function f α β γ has a set of normalized local coordinates (x1 , x2 , x3 ) for its arguments and separable basis functions are adopted, such that

(74)

f α β γ (x1 , x2 , x3 ) = f α (x1 ) f β (x2 ) f γ (x3 )

(72)

(75)

Three different possibilities for the individual one-dimensional polynomials have been examined. Taking f α (x1 ), to be specific, these possibilities are:  f α (x1 ) = xα 1

f α (x1 ) = T α (x1 )

(76)

where T here represents α Chebychev polynomial, and f α (x1 ) is a specially constructed set of orthogonal polynomials in the local coordinate x1 . For the choice of the weighting function ωi , j ,k , a symmetrical Hermite polynomial is suggested. Clearly, this technique is rather complicated, involving as it does series within series, and multiple applications of least squares solutions of linear simultaneous algebraic equations. Also, the total number of coefficients required for the overall representation can ran into the hundreds of thousands. 4) Point Mass Model In Section 2, we have seen that the earth’s gravity acts as a point mass in accordance with Newton’s law of universal gravitation, and given by Eq. (26). The point mass approach is to model the external disturbance potential field as the response to a finite set of mass points rather than as the response to continuous mass distribution over the entire volume. In contrast to surface layer models, there are no numerical integrations over area elements, only finite summations. It also has the advantage of being able to increase resolution (i.e., bandwidth) by adding more point masses in local regions. In missile guidance, for example, an accurate wide-bandwidth gravity model in the launch region is crucial for precise targeting. Getting the required bandwidth in a spherical harmonic expression requires increasing the number of terms in all regions. The expressions used to compute the anomalous gravity potential T (r ) and the anomalous gravity vector, δ(r ), are



i

T (r ) =

i =0 j =0 k =0

1 1 1

(73)

and

is an explicit formula, convenient for use in digital computer programming, to determine any Legendre function of degree n and order m. For inertial navigation systems, the gravity disturbance vector (s) can be obtained from

1 1 1

c α β γ f α β γ (x1 , x2 , x3 )

γ =n−α −β

(2n − 2k)! t n−m−2k k!(n − k)!(n − m − 2k)!

δ = grad T

n n N −α n=0 α =0 β=0

(69)

−n

311

Mi

r − r i r − r i δ(r ) = − Mi r − r i 3

(77) (78)

i

where M i are the point masses, r i are their respective position vectors, r is the point outside the geoid at which T and δ are being computed, and r − r i represents the distance between r and r i . A statistical version of the point mass model can be formulated wherein the point masses and their positions are treated as random variables and statistical estimation algorithms are used to estimate them.

312

G.M. Siouris / Aerospace Science and Technology 13 (2009) 301–315

cross-correlation function between the process U (x1 , y 1 , z1 ) and V (x1 , y 1 , z1 ) is defined as [17]

3.2. Statistical gravity model Statistical characterization of gravity models has received considerable attention in the last few years. The fifteen statistical gravity models most commonly identified in the literature are: 1) Gaussian, 2) Inverse Distance Square, 3) Reciprocal Distance, 4) Attenuated White Noise, 5) Logarithmic, 6) Algebraic, 7) Exponential, 8) Exponential Cosine, 9) Bessel, 10) Second order Markov, 11) Third-order Markov, 12) Expanded Third-order Markov, 13) Anisotropic Third-order Markov, 14) Two-dimensional Markov, and 15) Karhunen–Loéve Random Field. There are three major benefits of statistically modeling the anomalous gravity field: 1. More than 80% of the energy of the anomalous field is contained in wavelengths that are very short relative to the radius of the earth. Therefore, characterization is physically appealing. 2. Modern gravimetric data are inherently multisensor. Statistical estimation methods optimally combine mixed data bases, taking into account the spectral properties of each signal and noise. Most of these statistical estimation techniques require use of a statistical model of the gravity field. 3. Statistical gravity models facilitate powerful conceptual and computational tools, such as Kalman filter and covariance analysis. As stated in Section 1 (Introduction), gravity modeling is used extensively in modeling of inertial error dynamics for Kalman filter integration. Specifically, the gravity modeling characterization is used mainly in terms of devising analytical models for use in the truth models that bias the INS characteristics. In view of the unknown and short-wavelength character of the gravity field, it is natural to use statistical random process theory to characterize the gravity field. However, the theoretical basis for statistical models of the anomalous gravity field is somewhat tenuous. It is difficult to define statistical ensembles and probability distributions when there is only one earth and gravity at any particular point is a fixed, deterministic quantity (excluding transient effects such as tides). A probability space of “possible earths” may appear physically unnatural, but it is logically correct. Note that the arguments against statistical geodesy are not present when dealing with statistical modeling of the local gravity field. For the purposes of this investigation, local modeling is more appropriate than global modeling since, as stated above, more than 80% of the anomalous field is captured in wavelengths that are much shorter than the radius of the earth. Local statistical models of the anomalous gravity field take the form of autocovariance (or autocorrelation) function (ACF), power spectral densities (PSD) or a linear system of differential equations driven by white noise. Most of the statistical gravity models are specified by autocovariance functions, which are usually represented by fewer parameters than are functional models. Consider now two spatial points given by (x1 , y 1 , z1 ) and (x2 , y 2 , z2 ) in a local rectangular coordinate system (x, y , z). The autocovariance of U at these two spatial points is given by [6,17]

  ΦU U (x1 , y 1 , z1 ; x2 , y 2 , z2 ) = E U (x1 , y 1 , z1 )U (x2 , y 2 , z2 )

(79)

where E is the statistical expectation operator and U can be any gravimetric quantity (e.g., anomalous potential T , undulation N, north–south deflection ξ , east–west deflection η , normal component of the gravity disturbance vector δ g, or gravity anomaly g). If the gravity model is nonhomogeneous and nonisotropic, the ACF given by Eq. (79) will be a function of all six spatial position variables. If more than one random process is considered, the

  ΦU V (x1 , y 1 , z1 ; x2 , y 2 , z2 ) = E U (x1 , y 1 , z1 ) V (x2 , y 2 , z2 ) (80)

Notice that in general, ΦU U = ΦU V . Estimation of the cross correlation between north and east deflection in the horizontal plane shows that the cross correlation are fairly small (i.e., < 0.3) [7]. Therefore, little would be gained by using both the cross and autocorrelation functions as opposed to just using the autocorrelation functions. If the model is homogeneous, the ACF will be a function of the shift distance x, y , z only, that is [17]

ΦU U (x1 , y 1 , z1 ; x2 , y 2 , z2 ) = ΦU U (x, y , z)

(81)

where

x = x2 − x1 y = y2 − y1 z = z2 − z1 If the model is homogeneous and isotropic, the ACF is a function only of the distance between (x1 , y 1 , z1 ) and (x2 , y 2 , z2 ), or

ΦU U (x1 , y 1 , z1 ; x2 , y 2 , z2 ) = ΦU U (r )

(82)

where

r=

x2 + y 2 + z2

(83)

Since all the gravimetric quantities are linear functions of the anomalous potential and since the anomalous potential is harmonic, the elevations z1 and z2 enter into the ACF only through their sum z1 + z2 so we have

ΦU U (x1 , y 1 , z1 ; x2 , y 2 , z2 ) = ΦU U (ρ , z1 , z2 ) where

(84)

ρ is the (x– y ) shift distance given by

ρ = x2 + y 2

(85)

The form of the ACF given by Eq. (84) is considered to be homogeneous and isotropic in the z = 0 plane and harmonic for z  0. Simple harmonic extension of an ACF into outer space z > 0 implies continuity of its spectrum. In the x– y plane where z1 = z2 = 0, the corresponding spectral representation is given by

G (ω) =

1 2π

∞ ΦU U (ρ ) J 0 (ωρ )ρ dρ

(86)

−∞

where J 0 is the Bessel function of zero order and G (ω) is the spectrum of the function ΦU U (ρ ). Inversely,

∞ ΦU U (ρ ) = 2π

G (ω) J 0 (ωρ )ω dω

(87)

−∞

Eqs. (86) and (87) define a Hankel transformation. It is related to the well-known Fourier transformation. The two-dimensional Fourier transform reduces to the Hankel transform when isotropic functions (which depends on ρ but not on the azimuth φ ) are used [17]. If the gravity model is homogeneous but not isotropic in the x– y plane, the two-dimensional Fourier transform can be used to obtain the power spectral density (PSD). If the gravity field is isotropic, then the one- and two-dimensional ACFs are identical but the PSDs are different. This is due to the difference in operations between one-dimensional Fourier and Hankel transforms.

G.M. Siouris / Aerospace Science and Technology 13 (2009) 301–315

Simple harmonic extension into outer space z > 0 implies the existence of a continuous spectrum, and positive definiteness of ACF implies that the spectrum is positive, that is,

G (ω)  0 for all ω  0

(88)

Some of these ACF or PSD models can be expressed in statespace form. However, it is not always possible to find a state space representation since not all PSD are rotational functions and hence cannot be represented as outputs of a shaping filter (linear system driven by white noise) [17]. Even if the PSD is a rational function, state space models require straight line, constant altitude and constant speed trajectories. The time domain correlation functions are obtained by making the replacement

ρ = vt

(89)

where v is the speed of the vehicle and the characteristic distance d is usually transformed to a characteristic time constant β where

β = v /d

(90)

The gravity model characteristic distance, d, varies as a function of altitude and grid spacing. The autocorrelation function given by Eq. (79) or the cross-correlation function given by Eq. (80) are further characterized as homogeneous or isotropic based on the following definitions:

were: the direct description/linear interpolation, attenuated white noise/collocation, Karhunen–Loéve random field estimation, and the first- and second-order Markov models. Because of its tractability and ease of implementation, the first- and second-order Markov models will be discussed here. As stated in the introduction, the effects of gravitational uncertainties in limiting accuracy performance for precision INS have created a need for modeling these uncertainties. A linear optimal estimation algorithm is proposed for use in estimating and compensating the aircraft inertial navigation for gravity disturbance errors. The algorithm has been examined both at the gravity disturbance residual level and also at the inertial navigation position error level as a function of gravity data, grid spacing, and measurement error [15]. The application of modern covariance analysis and Kalman filtering tools for analyzing INS performance makes the linear state space Markov process a convenient tool for analyzing gravity errors. This section outlines statistical analysis of gravity models currently in use. Specifically, we will investigate the firstand second-order Markov gravity models.6 The first (and simplest) models used to evaluate the effects of gravitational anomalies and vertical deflections are based on a first order spatial autocorrelation function of the form [3,6,17]

φ g g (r ) = σ 2 g e −r / D c 2

Homogeneity. The statistical model of a gravity quantity is a homogeneous (or stationary) random process if, and only if, the ACF depends only on the distance vector between any two spatial points given by (x1 , y 1 , z1 ) and (x2 , y 2 , z2 ) and not on the spatial points themselves over the entire region. In other words, a homogeneous random field is invariant with respect to translation, but not rotation.

313

2 −r / D c

φεε (r ) = g σε e φηη (r ) = g 2 ση2 e −r / D c

(91) (92) (93)

where

σε , ση = rms values of the east and north component of the vertical deflection, respectively

Isotropic. The statistical model of a gravity quantity is isotropic if, and only if, its ACF depends only on the scalar distance between any two points (x1 , y 1 , z1 ) and (x2 , y 2 , z2 ) and not on the direction of the distance vector. In other words, an isotropic random field is invariant with respect to rotation, but not translation. Isotropy necessarily applies only to multidimensional random fields and not to one-dimensional stochastic processes, since one-dimensional processes are by definition invariant to rotation. It should be noted that gravity anomalies are neither stationary nor isotropic. Some areas are rough, others smooth, and direction dependent features corresponding to mountain range, continental shelf, fault, etc., are present in most areas. If the gravity model is homogeneous but not isotropic in the x– y plane, the twodimensional Fourier transform can be used to obtain the power spectral density (PSD). If the gravity field is isotropic, then the one- and two-dimensional ACFs are identical but the PSDs are different. This is due to the difference in operations between onedimensional Fourier and Hankel transforms.5 4. Candidate gravity models During the course of this investigation, 35 candidate gravity models have been identified and screened for compatibility, stability, ease of implementation of the models in the airborne navigation computer, and optimality. In the analysis and screening criteria of the candidate models, the candidates selected

5

The Karhunen–Loéve random model is a nonstationary nonisotropic characterization of the disturbance potential in a local region. Such a model is “spatially dynamic” as opposed to a static model specified by autocorrelation functions or power spectral densities. The model is also “bilateral” such that any given spatial point’s dependency upon all its neighboring point is preserved and no unnecessary limitation of causality is enforced on the system.

σ g = rms value of the gravitational anomaly D c = correlation distance g = normal gravity value normal to the reference ellipsoid r = distance parameter The statistics of this process is represented by the rms value of the errors (σ g , σε , ση ) and a correlation distance (D c ). These correlation functions can be represented by a linear firstorder differential equation driven by white noise, as illustrated in Fig. 8. In order to convert the spatial correlation distance to an equivalent correlation time (as required for covariance analysis), the correlation distance7 is divided by the vehicle horizontal velocity magnitude. Typical rms values of gravity errors and correlation distance for uncompensated gravity errors are:

σ g = 0.25 μg σε = ση = 5 arc sec D c = 20 nm These values are derived from a survey along the 35th parallel of the U.S. For analysis purpose, it is useful to approximate the first-order Markov process as an equivalent white noise driving the velocity error states [17] (see Appendix B for a brief discussion of the first-order differential equation Markov process). We will now

6 With regards to the third-order Markov models, Jordan [6,9] examined all the conditions of self-consistency of gravity models and proposed the third-order Markov undulation model that satisfies all the conditions of self-consistency. 7 Correlation distance is defined as the distance where the correlation coefficient becomes 1/e. However, in the horizontal direction, at the distance D, the correlation does not have the value 1/e, and therefore, D is designated the characteristic distance d.

314

G.M. Siouris / Aerospace Science and Technology 13 (2009) 301–315

Fig. 8. First-order Markov gravity error model.

discuss briefly the second-order Markov model. As stated earlier, gravity anomalies are dependent upon the mass distribution of the earth, which tends to be locally constant, so that the derivative of the ACF at zero shift should approach zero, corresponding to unity correlation. The second-order Markov gravity model has an autocorrelation (ACF) given by [6]



φ(ρ ) = σ 2 1 +

ρ d



e −ρ /d

(94)

where d is the characteristic distance (see Eqs. (89) and (90)) and ρ is the (x– y) shift distance given by Eq. (85). Some researchers suggest using the second-order Markov as a model for the gravity “anomaly” whereas others prefer using the gravity “undulation” model. We will select the gravity undulation model. Therefore, the second-order Markov undulation can be expressed by the following relations [6,8,9,17]

  ρ −ρ /d φ N N (ρ ) = σ N2 1 + e d   ρ −ρ /d φ g g (ρ ) = σ 2 g 1 − e 2d   ρ 2 2 φξ ξ (x, y ) = σξ 1 − cos ϕ e −ρ /d d   ρ 2 2 φηη (x, y ) = ση 1 − sin ϕ e −ρ /d d

(95)

  φ N N (t ) = σ N2 1 + β|t | e −β|t |   β|t | −β|t | φ g g (t ) = σ 2 g 1 − e 2   β|t | −β|t | φξ ξ (t ) = φηη (t ) = σξ2 1 − e 2

d2 u dt 2

+ 2β

du dt

− β2u = w

(101)

(102)

where u is either g or N, and w is a white noise process with variance 4β 2 σu2 δ(t ). The initial condition u (0) is a zero-mean random variable with variance σu2 , the initial condition du (0)/dt is a zero-mean random variable with variance β 2 σu2 , and u (0) and du (0)/dt are uncorrelated. The two-dimensional power spectral density (PSD) of any second-order Markov process is obtained by Fourier transforming the ACF; thus,

ΦU U (ω1 , ω2 ) = (97)

where ϕ is the heading angle (i.e., tan ϕ = x/ y), and ξ and η are the north–south and east–west coordinates, respectively. Note that even though the gravity anomaly and/or undulation statistics are assumed to be isotropic, the corresponding deflection statistics are nonisotropic. If a space–time transformation is made, the secondorder Markov undulation model takes the following form:

(100)

where β is the characteristic time constant given by Eq. (90). In either case (i.e., anomaly or undulation), the second-order Markov ACF can be associated with a time-dependent differential equation as follows:

(96)

(98)

(99)

−6π β 3 σu2 (β 2 + ω12 + ω22 )5/2

(103)

where ω1 and ω2 are frequency variables. The undulation model is favored because it consistently derives the anomaly statistics from the undulation statistics, whereas the undulation statistics are missing from the anomaly model. 5. Conclusions A variety of different techniques for in-flight gravity compensation of an INS have been examined, and techniques for evaluation of gravity residual errors (i.e., errors after compensation) and

G.M. Siouris / Aerospace Science and Technology 13 (2009) 301–315

methodologies for evaluating the effects of compensated gravity disturbances on INS performance have been presented. From the evaluation of the different gravity models, it is concluded that: 1. If the accurate gravity deflection data is available, navigation error due to gravity disturbances can be made small. 2. For compensation of gravity errors to be small compared to other INS errors neither the data storage nor the computation rate requirement is excessive based on current airborne computer technology. That is, a gravity compensation scheme requires a “Model” for the gravity disturbance vector and an “Algorithm” to derive from that model, the gravity disturbance vector at the computed position of the INS that is then used for compensating the accelerometer measurements of specific force. 3. It is important that the data stored and used for compensating the INS not contain any mid- or low-frequency gravity errors. Gravity data should not have any errors in the Schuler frequency (note that the Schuler frequency which appears in inertial platform √ leveling, is given by ωs = g /r where r is the radius of the earth, and g is the gravitational acceleration) region or below which, for vehicle velocities of 300 kn8 or more, corresponds to gravity frequencies of 1 cycle/400 miles or less. Constraints need to be placed on both the magnitude and spectral content of the error in the gravity data. Furthermore, this may be adequate in some applications (e.g., high speed vehicles), where a substantial amount of the unmodeled gravitational energy resides outside the Schuler bandwidth of the INS. The ultimate objective in gravity modeling efforts will be to demonstrate the effectiveness of gravity modeling in improving the accuracy of the high accuracy INS in flight tests. The flight tests used for comparison of accuracy with or without gravity compensation must be over the area where gravity map data is available. Finally, post processing of flight test data should be used to minimize the actual number of flight tests required. Acknowledgement The author would like to acknowledge the many helpful suggestions of the following individuals: Dr. Chaoyong Li, a research associate at the Dept. Electrical Engineering and Computer Science, University of Central Florida, FL, and Major Lynn R. Herche, USAFR, of NOAA. Dr. Li provided the author with many valuable suggestions for improving the quality of this paper, while Major Herche, a former student of mine at the AFIT, Dept. Electrical and Computer Engineering, WPAFB, OH, meticulously plotted numerous statistical cases and validated the results. Last but not least, the author gratefully acknowledges the anonymous reviewers for their perceptive comments and their time and effort spent with the manuscript. Appendix A In addition to the WGS84, the following geopotential models are commonly used in many aerospace applications: 1. EGM96 (Earth Gravitational Model 1996) — Developed by NASA, Ohio State University (OSU), and the National Imagery and Mapping Agency (NIMA). 2. GeoForschungs-Zentrum Potsdam (GFZ, Germany). 3. Goddard Earth Model (GEM) — Developed by the NASA Goddard Space Flight Center. 4. Joint Gravitational Model (JGM) — Developed by NASA, the University of Texas (UT), the OSU, and the Centre National d’Etudes Spatiales (CNES). 5. OSU-91A. 8

kn = knot (1 kn = 1 nm/hr = 1.8532 km/hr).

315

6. Texas Earth Gravity (TEG) — Developed by the UT. Appendix B The first-order Markov model is driven by a zero-mean, white Gaussian noise of strength q, and is modeled as follows:

dx/dt = −(1/τ )x(t ) + w (t )





E x(t )2 = (1/2)qt where

w (t ) = zero-mean Gaussian noise

τ = time constant (also referred to as correlation time) 

σ = E x(t )2



This model is used to represent exponentially time-correlated noises (i.e., band-limited errors). References [1] R. Bate, D. Mueller, E. White, Fundamentals of Astrodynamics, Dover, New York, 1971. [2] G. Bomford, Geodesy, 2nd ed., Oxford University Press, London, 1962. [3] B. Carnahan, H.A. Luther, J.O. Wilkes, Applied Numerical Methods, John Wiley & Sons, Inc, New York, 1969. [4] A.B. Chatfield, Fundamentals of High Accuracy Inertial Navigation, vol. 174, AIAA Progress in Astronautics and Aeronautics, Reston, VA, 1997. [5] W.A. Heiskanen, H. Moritz, Physical Geodesy, W.H. Freeman Co., 1967. [6] W.G. Heller, S.K. Jordan, Attenuated white noise statistical gravity model, Journal of Geophysical Research 84 (Nb9) (1979) 4680–4688. [7] D.Y. Hsu, An accurate and efficient approximation to the normal gravity, in: IEEE Position Location and Navigation Symposium, 20–23 April, 1998, pp. 38– 44. [8] S.K. Jordan, Statistical models for vertical deflection from gravity-anomaly models, Journal of Geophysical Research 77 (5) (1972) 971. [9] S.K. Jordan, Self-consistent statistical models for gravity anomaly – vertical deflections, and undulation of geoid, Journal of Geophysical Research 77 (20) (1972) 3660. [10] W.M. Kaula, Theory of Satellite Geodesy, Blaisdell Publishing Company, Waltham, MA, 1966. [11] X. Li, H.J. Goetze, Ellipsoid, geoid, gravity, geodesy, and geophysics, Geophysics 66 (6) (2001) 1660–1668. [12] C. Li, W. Jing, Fuzzy PID controller for 2D differential geometric guidance and control problem, IET Control Theory and Applications 1 (3) (2007) 564–571. [13] C. Li, W. Jing, C. Gao, Adaptive backstepping-based flight control system using integral filters, Aerospace Science and Technology 13 (2–3) (2009) 105–113. [14] H. Moritz, Geodetic reference system 1980 (vol. 66, pg. 2, 1992), Journal of Geodesy 74 (1) (2000) 128–133. [15] H. Moritz, Covariance functions in least-squares collocation, Dept. of Geodetic Science, Ohio State University, Report No. 240, June 1976. [16] National Imagery and Mapping Agency, Dept. of Defense World Geodetic System, NIMA, TR8350.2, 3rd ed., Amendment 1, January 2000. [17] A. Papoulis, Probability, Random Variables, and Stochastic Processes, 2nd ed., McGraw-Hill, New York, 1984. [18] C.E. Pearson, Handbook of Applied Mathematics, Van Nostrand Reinhold Co., New York, 1974. [19] D.T. Sandwell, W.H.F. Smith, Marine gravity anomaly from Geosat and ERS 1 satellite altimetry, Journal of Geophysical Research – Solid Earth 102 (B5) (1997) 10039–10054. [20] G.M. Siouris, Aerospace Avionics Systems: A Modern Synthesis, Academic Press, Inc., San Diego, CA, 1993. [21] G.M. Siouris, Missile Guidance and Control Systems, Springer-Verlag, New York, 2004. [22] P. Vanicek, E. Krakiwsky, Geodesy: The Concepts, 2nd ed., Elsevier, Amsterdam, 1986. [23] G.P. Woollard, The new gravity system – changes in international gravity base values and anomaly values, Geophysics 44 (8) (1979) 1352–1366. [24] J.J. Zhu, Exact conversion of earth-centered, earth-fixed coordinates to geodetic coordinates, Journal of Guidance Control and Dynamics 16 (2) (1993) 389–391. [25] J.J. Zhu, Conversion of earth-centered earth-fixed coordinates to geodetic coordinates, IEEE Transactions on Aerospace and Electronic Systems 30 (3) (1994) 957–962. [26] K.P. Schwarz, M.G. Sideris, Heights and GPS, in: GPS WORLD, February 1993, pp. 50–56.