3  Earth models

Author

David Gianazza

This chapter provides a few notions of geodesy that are useful to the computation of aircraft trajectories. For a more complete documentation, the reader may refer to the first chapters of the book of Michel Capderou  on satellites, and also to the book of Dominic J. Diston .

There are several possible models of the Earth’s surface, among which the geoid, the ellipsoid of revolution and the sphere. The geoid is defined as the equipotential surface of the gravity field, conforming to the shape defined by the actual mean sea level. This equipotential surface is not easy to compute and is often approximated by a simpler model: the ellipsoïd of revolution, where the Earth is considered as a sphere flattened at the poles. In an ellipsoid of revolution, each section in a meridian plane is an ellipse of parameters a for the major axis (lying in the equatorial plane), and b for the minor axis (between the South and North poles). The parameters a and b of the ellipse are constant, whatever the meridian. For some applications or computations, an even simple model can be used: a spherical Earth model.

Let us briefly describe the spherical model and the ellipsoid of revolution in the rest of this chapter, starting with the sphere for simplicity’s sake.

3.1 The Spherical Earth Model

3.1.1 The Earth-Centered, Earth-Fixed (ECEF) Coordinate System

The Earth is modeled here as a sphere of radius RT. Let us define a reference frame fixed to the Earth and for which the axis ze passes through the poles and is oriented from the South pole to the North pole. The xe axis is chosen in the equatorial plane and passes through the center of the Earth and through the Greenwich meridian (an arbitrarily chosen meridian). The ye completes this system and is chosen so as to form a direct orthonormal coordinate system centered on O, the Earth’s center.

In the following, we shall denote (ıe,ȷe,ke) the orthonormal vectors of the ECEF reference frame Oxeyeze.

Figure 3.1: Spherical Earth model

In this coordinate system, the position of a point P is given by its latitude, longitude, and distance from the center of Earth, or more simply its altitude h above the surface of the globe, as illustrated in . The latitude, denoted μ in the following, is the angle between the equatorial plane xeOye and OP. The longitude, denoted λ, is the angle between the Greenwich meridian plane xeOze and the meridian plane containing OP.

3.1.2 The North-East-Down (NED) reference frame

For given point P, located on or at proximity of the Earth’s surface, let us define another system of axes, called the NED system, or the local horizontal reference frame centered on point P. In this NED frame, the xh axis is in the local horizontal plane and passes through P, pointing to the North. The yh axis is also in the horizotal plane and passes through P, but it points to the East. Finally, the vertical axis zh passes through P and points downward, toward the center of the Earth. The axes system Oxhyhzh is represented on .

An orthonormal basis of vectors (ıh,ȷh,kh) is associated to this system.

Figure 3.2: NED reference frame, on a spherical earth

3.1.3 Coordinates of a Point on the Sphere

The Cartesian coordinates – in the ECEF system – of a point P at altitude h are given by :

(3.1){OP}ECEF|xe=(RT+h)cosμcosλye=(RT+h)cosμsinλze=(RT+h)sinμ

In the NED system where the unit vector kh points downward from P, the position vector can be expressed very simply by , and the coordinates by .

(3.2)OP=(RT+h)kh

(3.3){OP}NED|xh=0yh=0zh=(RT+h)

3.1.4 Velocity of a Point with Respect to the ECEF Frame

Let Oxmymze be the orthonormed direct system such that xmOze is the meridian plane containing P, the position of the mobile agent. This referential is obtained simply by rotating the ECEF axes Oxeyeze of an angle λ around the axis of the poles ze. Let (ım,ȷm,ke) be the orthonormal basis associated with the reference frame fixed too the meridian plane passing through P.

The velocity with respect to the ECEF reference frame, considered as fixed, is given in , taking into account the angular speed ΩOxmymze/ECEF=λ˙ke of the Oxmymze system around the axis Oze:

(3.4)dOPdt/ECEF=dOPdt/Oxmymze+ΩOxmymze/ECEFOP

The derivative dOPdt/Oxmymze is obtained simply by deriving assuming Oxmymze is fixed. With this assumption, we place ourselves in the meridian plane and only take into account the rotational motion of kh in this plane when expressing its derivative :

(3.5)dOPdt/Oxmymze=(RT+h)dμdtdkhdμh˙kh=μ˙(RT+h)ıhh˙kh

By combining and , the velocity with respect to the ECEF reference frame is then expressed (see ) in a simple form in the basis of vectors (ıh,ȷh,kh) associated with the local horizontal frame NED.

(3.6)dOPdt/ECEF=μ˙(RT+h)ıhh˙kh+(λ˙ke)((RT+h)kh)=μ˙(RT+h)ıh+λ˙(RT+h)cosμȷhh˙kh

Adopting the slightly unwieldy but explicit notation convention where “/ECEF” means “with respect to the Earth’s reference frame” and {.}NED means “in the NED coordinate system”, the coordinates of the velocity with respect to the Earth are given by .

(3.7){dOPdt/ECEF}NED=(μ˙(RT+h)λ˙(RT+h)cosμh˙)

3.2 The Ellipsoid of Revolution

A more accurate approximation than the sphere of the Earth’s surface is the ellipsoid of revolution. In the ellipsoid model, each section in a meridian plane is an ellipse of major axis a and minor axis b. The parameters a and b of the ellipse are the same in all meridian planes. The approximate surface of the Earth is thus determined by a single ellipse, which is rotated around its minor axis, which is coincident with the polar axis.

In the following, we will denote e the eccentricity of the ellipse, which verifies the following equation :

(3.8)e2=1b2a2

The equation of the ellipse is recalled in , where z is the minor axis (the polar axis), and x is the major axis of the ellipse:

(3.9)x2a2+z2b2=1

3.2.1 Several Definitions of Latitude

Figure 3.3: Geodetic latitude (μ) and geodetic altitude (h) on an ellipsoid of revolution

Historically, latitude was calculated by measuring angles between the local horizontal plane and the direction of reference stars (the North Star for example). The horizontal plane is determined as being perpendicular to the direction of the plumb line. This direction is therefore the normal to the equipotential surface represented by the local geoid. This latitude measured “on the ground” is called the “astronomical latitude”.

The geodetic latitude is defined as the angle between the equatorial plane and the line perpendicular, at the considered point N, to the surface of the ellipsoid of revolution (see ). Note that this perpendicular line does not pass through the center O, but intersects the polar axis at a point I. The distance IN is sometimes called the “great normal” and is denoted N.

The geocentric latitude is the angle between the equatorial plane and the axis going from the center of the Earth O to the point N considered. Finally, the parametric latitude or reduced latitude is used for calculation purposes and corresponds to the geocentric latitude of a point N (see ) on a circle with center O and radius a. This point N is obtained by following the parallel to the axis Oz passing through the point N. The geocentric and parametric latitudes are illustrated in .

In the rest of this chapter, μ will denote the geodetic latitude, χ the geocentric latitude, and u the parametric latitude.

Figure 3.4: Geocentric (χ), geodetic (μ), and parametric (u) latitudes

3.2.2 Some Useful Characteristics of the Ellipse

3.2.2.1 Relations between geodetic, geocentric, and parametric latitudes

Let us remind that the ellipse is an affine transformation, of axis Ox of direction Oz and of ratio ba, of the circle of center O and radius a. In other words, the points of the ellipse are obtained from the circle of center O and radius a by applying a multiplicative factor ba to the coordinate z.

Therefore, in , we observe that HNOH=baHNOH. We immediately deduce the relation between the geocentric latitude χ and the parameterized latitude u.

(3.10)tanχ=batanu

Moreover, by differentiating , we obtain the slope of the tangent to the ellipse at point N, expressed in :

(3.11)dzdx=b2a2xz

The direction of the normal is then given by :

(3.12)tanμ=a2b2zx

From the definition of geocentric latitude, we can also see that:

(3.13)tanχ=zx

From and we deduce the relation between the geocentric latitude χ and the geodetic latitude μ :

(3.14)tanχ=b2a2tanμ

By combining and , we can easily deduce the relation between the parametric latitude u and the geodetic latitude μ:

(3.15)tanu=batanμ

3.2.2.2 Coordinates of a point on the ellipse

The point N of is on a circle of radius a. Its coordinates are simply expressed as a function of the parameterized latitude, as follows:

{ON}Oxz|x=acosuz=asinu

The point N on the ellipse is obtained by an affine transformation of the point N, of ratio ba. The z-coordinate of N is therefore simply ba(asinu). The coordinates of N are given by :

(3.16){ON}Oxz|x=acosuz=bsinu

Furthermore, considering the geodetic latitude μ and the great normal N=IN, we observe that x=Ncosμ. Taking into account this expression and , the coordinate z is rewritten as a function of the geodetic latitude as follows:

z=bsinu=batanu(acosu)=b2a2tanμ(Ncosμ)=Nb2a2sinμ

By introducing the eccentricity , we finally obtain the coordinates of N as a function of the geodetic latitude μ, expressed in the following equation :

(3.17){ON}Oxz|x=Ncosμz=N(1e2)sinμ

3.2.2.3 Expression of the great normal N

The expression of the great normal N is simply found from and the equation of the ellipse :

x2a2+z2b2=1N2cos2μa2+N2(1e2)2sin2μb2=1N2[1sin2μ+a2b2(1e2)2sin2μ]=a2N2[1sin2μ+(1e2)sin2μ]=a2N2(1e2sin2μ)=a2

The great normal can be expressed as follows: (3.18)N=a1e2sin2μ

3.2.3 The ECEF and Other Reference Frames for the Ellipsoid of Revolution

Figure 3.5: ECEF and mobile reference frames, on the ellipsoid of revolution

In the case of the ellipsoid of revolution, we use an ECEF reference frame similar to the spherical model. Its system of axes is noted here Oxeyeze, with a base of orthonormal vectors (ıe,ȷe,ke). For calculation purposes, we will also use an orthonormal reference frame (ım,ȷm,ke) obtained by a rotation of angle λ around the polar axis ze. The vectors ım and ke are thus in the meridian plane passing through N (and P). We will denote Oxmymze the associated axis system.

The local horizontal reference frame NED, with orthonormal base (ıh,ȷh,kh), is similar to the one of the spherical model, except that the vertical direction determined by vector kh does not pass through the center of the Earth but is defined as the normal to the surface of the ellipsoid. The vector ıh points to the North, and ȷh to the East. These different reference points and axis systems are described in .

3.2.4 Coordinates of a Point on the Ellipsoid

The position of a point P located at an altitude h above the surface of the ellipsoid of revolution can be broken down as follows (see ), where the point N (nadir) is the point on the surface of the ellipsoid located at the vertical of the point P: OP=ON+NP

3.2.4.1 Coordinate of the Nadir Point N

(3.19){ON}Oxmymze|xm=Ncosμym=0ze=N(1e2)sinμ

(3.20){ON}ECEF|xe=Ncosμcosλye=Ncosμsinλze=N(1e2)sinμ

In these equations, N is the “great normal” at point N, given by .

3.2.4.2 Coordinates of point P at altitude h

(3.21){OP}Oxmymze|xm=(N+h)cosμym=0ze=[N(1e2)+h]sinμ

(3.22){OP}ECEF|xe=(N+h)cosμcosλye=(N+h)cosμsinλze=[N(1e2)+h]sinμ

3.2.5 Velocity of a Point with Respect to the ECEF Frame

dONdt/Oxmymze=μ˙(dNdμcosμNsinμ)ım+(1e2)μ˙(dNdμsinμ+Ncosμ)ke

withdNdμ=ae2sinμcosμ(1e2sin2μ)32

dONdt/Oxmymze=μ˙[ae2sinμcos2μ(1e2sin2μ)32a(1e2sin2μ)12sinμ]ım+μ˙(1e2)[ae2sin2μcosμ(1e2sin2μ)32+a(1e2sin2μ)12cosμ]ke=μ˙a(1e2sin2μ)32sinμ[e2cos2μ(1e2sin2μ)]ım+μ˙a(1e2)(1e2sin2μ)32cosμ[e2sin2μ+1e2sin2μ]ke=μ˙a(1e2)(1e2sin2μ)32[sinμım+cosμke]=μ˙a(1e2)(1e2sin2μ)32ıh

(3.23)dONdt/Oxmymze=μ˙RMıh

where RM is the meridian radius of curvature (see ):

(3.24)RM=a(1e2)(1e2sin2μ)32

Figure 3.6: Meridian radius of curvature RM, in blue, and radius of curvature RP in the local plane of latitude μ, in red

In addition, we have:

dNPdt/Oxmymze=h˙khhdkhdt=h˙khhμ˙dkhdμ=hμ˙ıhh˙kh

We deduce that:

dOPdt/Oxmymze=μ˙(RM+h)ıhh˙kh

Taking into account the rotation ΩOxmymze/ECEF=λ˙ke of the reference frame Oxmymze around the axis Oze, we obtain: dOPdt/ECEF=dOPdt/Oxmymze+ΩOxmymze/ECEFOP=μ˙(RM+h)ıhh˙kh+λ˙ke(ON+NP)=μ˙(RM+h)ıh+λ˙(N+h)cosμȷhh˙kh

By introducing RP=Ncosμ, the radius of the parallel passing through the point N (see ), we finally obtain the expression

(3.25){dOPdt/ECEF}NED=(μ˙(RM+h)λ˙(N+h)cosμ)h˙)=(μ˙(RM+h)λ˙(RP+hcosμ)h˙)

In this equation, the meridian radius of curvature RM is given by , and the radius RP of the parallel circle is given by the following equation :

(3.26)RP=Ncosμ=acosμ1e2sin2μ

List of Acronyms

ECEF
Earth-Centered, Earth-Fixed
NED
North-East-Down