Naval Research Laboratory

Washington, DC 20375-5320

NRL/FR/5344— 16-10,279

An Overview of Major Terrestrial, Celestial, and Temporal Coordinate Systems for Target Tracking

David Frederic Crouse

Surveillance Technology > Branch Radar Division

August 10, 2016

Approved for public release; distribution is unlimited.

REPORT DOCUMENTATION PAGE

Form Approved OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the data needed, and completing and reviewing this collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this burden to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number. PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY) 10-08-2016

2. REPORT TYPE

Formal Report

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

An Overview of Major Terrestrial, Celestial, and Temporal Coordinate Systems for Target Tracking

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

David Frederic Crouse

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

Naval Research Laboratory 4555 Overlook Avenue, SW Washington, DC 20375-5320

8. PERFORMING ORGANIZATION REPORT NUMBER

NRL/FR/5 344- 1 6- 1 0,279

9. SPONSORING / MONITORING AGENCY NAME(S) AND ADDRESS(ES)

Office of Naval Research One Liberty Center 875 North Randolph Street Suite 1425

Arlington, VA 22203-1995

10. SPONSOR / MONITOR’S ACRONYM(S)

ONR

11. SPONSOR / MONITOR’S REPORT NUMBER(S)

12. DISTRIBUTION / AVAILABILITY STATEMENT

Approved for public release; distribution is unlimited.

13. SUPPLEMENTARY NOTES

14. ABSTRACT

Coordinate systems useful for those designing target tracking algorithms are discussed. Strict internationally recognized definitions of terrestrial and celestial coordinate systems are addressed, as well as a number of other lesser coordinate systems that arise when dealing with data typically used in target tracking algorithms. International organizations defining standard systems are indicated, as well as where to obtain data. Limita¬ tions in the data are demonstrated and methods for determining one’s orientation using gravitation and magnetic measurements, as well as using stars, are presented. The report focuses on implementing the coordinate-system standards using freely available data and software but does not provide specific details on the strict definitions of all of the standards. More detailed information sources are provided in an extensive reference list. So readers can better reproduce the results in this report, summaries of the sources of the data and external code used in creating the plots in the report are provided. All online sources in this report were verified as of April 2015 or later. Even if high-precision global coordinate system definitions are not necessary to meet the basic requirements of a tracking network, they can be necessary in evaluating the performance of tracking algorithms in such a network.

15. SUBJECT TERMS

Coordinate systems Standards organization

Position measurement Orientation estimation

Geomagnetism

Standardization

Time

Tides

Gravity

IERS

WGS 84

16. SECURITY CLASSIFICATION OF:	17. LIMITATION OF ABSTRACT	18. NUMBER OF PAGES
a. REPORT Unclassified Unlimited	b. ABSTRACT Unclassified Unlimited	c. THIS PAGE Unclassified Unlimited	Unclassified Unlimited	180

19a. NAME OF RESPONSIBLE PERSON

David Frederic Crouse

19b. TELEPHONE NUMBER (include area code)

202-404-8106

1

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39.18

This page intentionally left blank

CONTENTS

EXECUTIVE SUMMARY . E-l

1. INTRODUCTION . 1

2. CELESTIAL AND TERRESTRIAL COORDINATE SYSTEMS . 6

3. TEMPORAL COORDINATE SYSTEMS . 19

3.1 Defining Basic Temporal Coordinate Systems . 19

3.2 Time and Earth- Orientation Parameters . 26

3.3 Units of Time . 28

3.4 Time for Celestial Observations . 29

4. COORDINATE SYSTEMS FOR DESCRIBING LOCATIONS ON AND NEAR THE EARTH . 30

4.1 A Number of Common Spherical and Direction-Cosine Systems . 30

4.2 Coordinates on the Reference Ellipsoid . 33

5. VERTICAL HEIGHT AND GRAVITY . 42

5.1 Definitions of Height . 42

5.2 Computing Orthonretric Heights . 44

5.3 Considering Pressure Altitudes . 50

5.4 Gravitational Models and the Local Vertical . 53

6. THE EFFECTS OF TIDES ON COORDINATE SYSTEMS . 57

6. 1 Why Tides Can Matter . 57

6.2 Tides in Terms of Ground and Ocean Motion . 58

6.3 Tides in Terms of the Gravitational Field . 61

7. THE EARTH’S MAGNETIC FIELD . 68

8. DETERMINING THE ORIENTATION OF A SENSOR . 73

8. 1 Orientation Estimation Using Sets of Common Observations . 74

8.2 Orientation Estimation from Magnetic and Gravitational Measurements . 79

8.3 Orientation Estimation Using Star's . 83

9. CONCLUSION . 95

ACKNOWLEDGMENTS . 96

APPENDIX A — Accelerating Reference Frames and the Coriolis Effect . 97

APPENDIX B — Converting from Cartesian to Geodetic Ellipsoidal Coordinates . 101

iii

APPENDIX C — Explaining the Relativistic Time Dilation Plot . 103

APPENDIX D — Converting to and from Ellipsoidal Harmonic Coordinates . 105

APPENDIX E — Principal Axes, Precession, Nutation, and the Pole Tide . 107

APPENDIX F — Evaluating Spherical-Harmonic Potentials and their Gradients . Ill

F. 1 Pines’ Singularity-Free Method . 114

F.2 The Modified Forward Row Algorithm . 117

F.3 Synthesizing the EGM2008 Error Terms . 121

F.4 The Ellipsoidal Approximation and the Ji Model . 123

APPENDIX G — Umeyama’s Algorithm for General Transformation Estimation . 125

APPENDIX H — A Derivation of the Score Function for Star-to-Catalog Assignment . 127

APPENDIX I — Getting a Direction Vector from the SPICE Toolkit . 133

APPENDIX J — Quantifying Errors Due to Aberration . 137

APPENDIX K— Why the Fisher Distribution is Bad for High-Precision Measurements . 141

APPENDIX L — Drawing the Outline of a Celestial Object . 143

APPENDIX M— Acronyms and Abbreviations . 145

REFERENCES . 149

iv

EXECUTIVE SUMMARY

When tracking targets over long ranges, it is essential that all sensors in a network translate their measure¬ ments or tracks into a common global coordinate system. This report addresses how international standards setting organizations and real data arc used to define usable realizations of common terrestrial and celestial coordinate systems. It quantities the magnitudes of errors that arise when establishing a coordinate sys¬ tem due to a number of model limitations and physical effects (other than measurement errors), such as the ground physically moving due to tides. This report brings together concepts across multiple fields to demon¬ strate how one’s orientation can be determined using magnetic, gravitational, and celestial measurements. Links and references where one can find real data and code necessary for simulating and using standard celestial and terrestrial coordinate systems arc provided.

E-l

This page intentionally left blank

AN OVERVIEW OF MAJOR TERRESTRIAL, CELESTIAL, AND TEMPORAL COORDINATE SYSTEMS FOR TARGET TRACKING

1. INTRODUCTION

Three previous tutorials have focused on the use of nonlinear three-dimensional (3D) monostatic and bistatic measurements, in general [55], and in refractive environments [54], as well as the use of arbitrary nonlinear continuous-time dynamic models [57]. However, a prerequisite to performing accurate tracking in a global environment is the establishment of a common coordinate system across sensors. Indeed, compen¬ sation for the effects of atmospheric refraction discussed in [54] might be completely useless if the pointing direction of the sensor cannot be determined to a sufficient precision.

This report looks at strict, internationally recognized definitions of terrestrial and celestial coordinate systems as well as a number of other lesser coordinate systems that arise when dealing with data typically used in target tracking algorithms. Section 2 reviews Cartesian celestial and terrestrial coordinate systems, while Section 3 goes over time standards. Differing time standards might initially seem to be a trivial detail when performing target tracking, except for the fact that the temporal standards arc intertwined with celestial coordinate systems due to general relativity.

When simply tracking aircraft for air defense purposes, one would assume that celestial coordinate systems arc of little importance. However, celestial objects can play a role in establishing the global look directions of each radar in a network to a high precision (networked sensor registration). Additionally, one should be able to identify celestial objects to avoid confusion with true targets. For example, as shown in real weather-radar data in [139], the Sun can increase the clutter background of a radar in a particular direction.1 Knowing the location of the Sun can be useful in advanced target tracking algorithms that try to account for reduced target detection probabilities of targets near interference sources as in [185, Ch. 7.2, 7.3]. However, to identify the location of the Sun relative to an observer, one might be able to convert between celestial and terrestrial coordinates.

Section 4 provides an overview of basic coordinate systems for determining the location of an observer with respect to the surface of the Earth. This covers geodetic, geocentric, and ellipsoidal harmonic coor¬ dinates.2 Geodetic coordinates are how one typically reports locations on a map in terms of latitude and longitude. However, geodetic coordinates are taken with respect to a reference ellipsoid, whereas many think that it is with respect to a sphere, as is the case with geocentric coordinates. Ellipsoidal harmonic components shall play a role in gravitational models later in the report.

Manuscript Approved July 22 2016.

'This is because, as discussed in [98, Ch. 3.1], for wavelengths longer than 0.4pm, the spectrum of the Sun is approximately the spectrum of a black body having temperature 5760 K, meaning that within the radio bands typically used by radar, solar interference increases with increasing frequency.

^Section 4 does not cover different projections that one might use when creating maps, as the different projections tend to be of limited use in target-tracking systems except for the display of data to the user; see [65, 302] for information on basic map projections.

1

2

David Frederic Crouse

On the other hand, a lot of confusion can arise when making more intelligent tracking algorithms that utilize terrain data. This could be, for example, to determine whether a missile will hit a mountain or go over a mountain, or to determine whether a mountain is obscuring the view of a target. Digital terrain elevation data (DTED) is often given in terms of orthometric heights, whereas one typically expects a height reported to match an elevation reported by a Global Positioning System (GPS) receiver. Section 5 describes different vertical datums that one might encounter when designing a target- tracking system. Such datums arc often defined in terms of the Earth’s gravitational field, requiring more physics than one would expect. Section 5 also reviews the definition of pressure altitudes, which one must be familial- with when using data from aircraft transponders that report pressure altitude.

One generally wants to have to account for as few physical effects as possible while still meeting any requirements for the accuracy of a system. Sections 5, 6, and 7 consider the effects of the strict definitions of the Earth’s gravitational field, the effects of tides, and models for the Earth's magnetic field. At first glance, each of these effects can appeal- to be insignificant. However, all of them can be significant in certain applications. For example, as discussed in Section 6, tides can move points on the surface of the Earth by tens of centimeters. Some global navigation satellite system (GNSS) receivers try to “correct” for tides, producing coordinates for stationary objects that do not appeal- to move over time but which will then not represent “true” coordinates in an Earth-centered Earth-fixed (ECEF) coordinate system [264]. 3 Depending on the application and level of precision desired in a sensor network, such corrections may or may not be necessary or desirable. For example, for high-precision tracking of orbital debris, such shifts might matter. The magnetic field plays a role when tracking ships using transponder data as ships often report their heading with respect to magnetic North. However, one might be tempted to also use the magnetic field to align radars in a network. Section 7 demonstrates that magnetic field models are not suitable for high-accuracy global direction estimation.

The establishment of a common coordinate system in a sensor network is an ostensibly simple task: One can use a receiver for signals from a global navigation satellite system (GNSS), such as the United States’ GPS, whose interface is described in [11], to determine the location of all of the sensors to an accuracy better than 5 cm [113]. 4 However, the determining the orientations of the sensors is much more difficult. In the target tracking literature, numerous algorithms for determining the orientation of a sensor in a global coordinate system exist. However, such algorithms generally require an initial estimate of a sensor’s orientation. For example, the algorithms presented in [59, 204, 324] and [372] all use linear approximations to deal with small errors in the range, elevation, and azimuth measurements of multiple sensors in three dimensions. Such approximations are useless without an initial estimate. Moreover, if the local vertical of the radar is itself biased, the corrections for the bias cannot be expressed as constant additive or multiplicative modifications to the measured azimuthal angle.5 The vertical bias problem is worse when considering algorithms such as those in [1, 247, 250, 371] and [373] that tty to establish a common two-dimensional (2D) coordinate system among multiple sensors, as might be the case when using multiple radars that are incapable of making elevation measurements. If the radars are designed such that their local verticals align with the local gravitational verticals, then the planes in which their measurements are taken do not coincide

3For example, the open-source GPS Toolkit (GPStk) [120] includes a class for solid-Earth tides, as documented at http: //www. gpstk.org/ doxygen/ SolidEarthTides_8hpp.html .

4Some of the highest-precision commercially available receivers can be found by searching for "real-time kinematic GPS re¬ ceiver” online.

5 Consider a sensor with the boresight direction along the z-axis that is observing a target with zero elevation and azimuthal angle 6. In the extreme case, if the boresight direction is rotated 90°, the target will have zero azimuthal angle regardless of the original angle 6 , since the projection of the vector into the rotated plane is zero.

Coordinate Systems for Target Tracking

3

and it is not possible to establish a common 2D coordinate system among the radars. Over short distances, an ad hoc 2D coordinate system can be established, but it will never be of high precision. On the other hand, algorithms for aligning a 2D sensor with 3D sensors do exist [128, 130]. However, the nonlinear bias estimation algorithms of [132, 189] and the linearized algorithms of [129,284,321], which explicitly consider biases in the 3D sensor orientation, all require initial orientation estimates.

To fill the gap in the tracking literature for establishing a 3D global coordinate system without any prior information. Section 8 brings together all of the general concepts discussed in this report to present two techniques for determining a sensor’s orientation in a global coordinate system. One method utilizes geomagnetic and gravitational measurements, the other uses measurements of stars;6 both techniques arc based on the same underlying algorithm. Section 8 demonstrates how a wide variety of physical effects must be taken into account for high-precision orientation estimation. The problem of sensor localization is not addressed, because one can simply purchase a commercially available GPS receiver.7

The report is concluded in Section 9, where specific contributions of the report not already present in the literature arc highlighted. Appendix M lists abbreviations and acronyms used in the report.

This report focuses on implementing the coordinate-system standards using freely available data and software and does not provide specific details on the strict definitions of all of the standards. More detailed information can be found in the references in Table 1. So that readers can better reproduce the results in this report. Table 2 summarizes the sources of the data and external code used in creating the plots in this report. All online sources in this report were verified as of April 2015 or later.

Even if high-precision global coordinate system definitions are not necessary to meet the basic require¬ ments of a tracking network, they can be necessary in evaluating the performance of tracking algorithms in such a network. For example, Automatic Dependent Surveillance-Broadcast8 (ADS-B) signals sent from aircraft or Automatic Identification System9 (AIS) messages broadcast y ships can be used as “truth” data against which estimation and tracking algorithms can be compared, because such transponders provide identifying numbers and the Global Positioning System (GPS) locations of planes/ ships, among other in¬ formation. Such data is widely available.10

6Whereas the majority of algorithms for sensor registration for tracking applications focus on the use of common observations of aircraft, the explicit use of observed aircraft is not considered here for the simple reason that such algorithms will not work if very few aircraft are present. For example, after the terrorist attacks in the United States on 1 1 September 2001, the airspace over the United States and Canada was closed for a few days [201]. Additionally, aircraft and satellites are not point targets. In the radar community, it has long been known that the apparent location of an aircraft can vary significantly depending on its shape and reflectivity [143],

7 Commercial GPS solutions are often significantly cheaper than designing a GPS receiver from the ground up - a task that is much more complicated than it seems. To design a GPS receiver, a number of possibly time-consuming, ad hoc fixes are necessary to correct for poor satellite design, which varies from satellite to satellite. For example, as described in [310], the L5 signal on satellite SVN49 is corrupted with a multipath return because of reflections from an auxiliary port to the satellite’s antenna array and its broadcast ephemerides are biased by 150m. Failure to include satellite-specific corrections will result in degraded performance.

8The standards for ADS-B transponders are maintained by the Radio Technical Commission for Aeronautics (RCTA). A list of the documents defining the standard is available at http : //adsb.tc.f aa. gov/ ADS-B. htm The documents can be expensive to those who are not members of RCTA. A detailed description of the data (ADS-B and otherwise) broadcast by Mode S transponders, which tend to be the most common internationally, is available (for a price) from the International Civil Aviation Organization (ICAO) [150],

9 Specifications on AIS transponders are available for free from the International Telecommunication Union (ITU) at [155].

l0Most modern aircraft are equipped with transponders supporting ADS-B. and all aircraft operating in Class A, B, and C airspace in the United States (generally all commercial aircraft) are required by U.S. law to carry such transponders as of 1 lanuary

4

David Frederic Crouse

Table 1: References for the Definitions of Modem Celestial and Terrestrial Coordinate Systems

Standard	Description	Reference
IERS Conventions 2010	Standard for terrestrial and celestial coordinate systems	http://maia.usno.navy.mil/conv2010/ http://tai.bipm.org/iers/conv20 1 0/
IERS Conventions Updates	Errata and changes since the last conventions release	http://maia.usno.navy.mil/conv2010/convupdt.html http://tai.bipm.org/iers/convupdt/listupdt.html
WGS 84 Standard	Coordinate-system standard for the United States Department of Defense	http://earth-info.nga.mil/GandG/publications/ NGA STND 0036 1 0 0 WGS84/NGA.STND.0036 1 . 0.0 WGS84.pdf
Explanatory Supplement to the Astronomical Almanac	Explains the use of the many coordinate systems in detail	Book: [332]
Errata for the Explanatory Supplement		http://aa.usno.navy.mil/publications/docs/exp_supp_errata.pdf
The Astronomical Almanac for the Year 2014	The Navy’s astronomical almanac	Book: [70] Associated Sites: http://asa.usno.navy.mil and http://asa.hmnao.com
Errata 2014	Errata for the Astronomical Almanac	http://asa.usno.navy.mil/Errata/Errata 2014.html
Time from Earth Rotation to Atomic Physics	Describes international time standards	Book: [226]

Note that the Explanatory Supplement to the Astronomical Almanac in [332] refers to the 2014 edition of the almanac, not the frequently cited but out-of-date 1992 edition [285].

2020 [89] (similar regulations exist in Europe [81]). Similarly, the International Maritime Organization’s (IMO) Safety of Life at Sea (SOLAS) Convention [152] mandates that a transponder supporting AIS be carried by all passenger ships, ships of 300 gross tonnage(Gross tonnage is not a measure of weight). Rather, it is a measure of the volume of a ship [151],

Coordinate Systems for Target Tracking

5

Table 2: Sources of the External Data and Code Related to the Plots in this Report

Description	Reference	Uniform Resource Locator (URL)
Gravitational Model Data and Documentation
EGM2008 Tide-Free and Zero-Tide Spherical-Harmonic Coefficients		http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/first release.html
EGM2008 Height Anomaly to Geoid Undulation Coefficients		http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/egm08 wgs84.html
Detailed Report on the EGM96 Model (Predecessor to EGM2008)	[198]	http://bowie.gsfc.nasa.gov/697/staff/lemoine/EGM96 NASA-TP-1998-206861.pdf
Document on the Development of the EGM2008 Model & Errata	[252,253]
Tide Model Data and Documentation
FES2004 Ocean Tide Coefficients		ftp://tai.bipm.org/iers/convupdt/chapter6/tidemodels/fes2004_Cnm-Snm.dat http://maia.usno.navy.mil/conv2010/convupdt/chapter6/tidemodels/ fes2004 Cnm-Snm.dat
Document on the Development of FES2004	[209]	ftp://tai.bipm.org/iers/temp/FES2004 tides Biancale/FES2004.pdf
IERS Subroutines for Tides (Chapter 7), including step2diu and step21on with dehanttideinel		http://maia.usno.navy.mil/conv2010/software.html ftp://tai.bipm.org/iers/convupdt/chapter7/
IERS Conventions 2010 Tide Model	[257, Ch. 6,7]	http://www.iers.org/IERS/EN/Publications/TechnicalNotes/tn36.html
A Solid Tide Potential Model Simpler than in the IERS Conventions	[332, Ch. 5]
Earth Magnetic Field Data, Documentation, and Software
EMM2010 Spherical-Harmonic Magnetic Coefficients		http ://w w w. ngdc .noaa. gov/geomag/EMM/
WMM2010 Spherical-Harmonic Magnetic Coefficients		http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml
WMM Documentation (Relevant to the EMM too)	[223]	http://www.ngdc.noaa.gov/geomag/WMM/data/WMM2010/WMM2010 Report.pdf
Documentation for a Predecessor to the EMM	[219]	http://citeseerx.ist.psu.edu/viewdoc/ download ?doi= 1 0. 1 . 1 .380. 1 3 1 9&rep=rep 1 &type=pdf
Software Created for the Simulations in this Report Performing Spherical-Harmonic Magnetic Field Synthesis (Choose David F. Crouse)		http://www.ngdc.noaa.gov/geomag/WMM/thirdpartycontributions.shtml
Earth Terrain Models and Documentation
Earth2012 Spherical-Harmonic Terrain Coefficients with Water		http://geodesy.curtin.edu.au/research/models/Earth2012/
Article on Models Used Deriving Earth2012	[135]	http://www.cage.curtin.edu.au/ will/201 UB008878.pdf
DTM2006.0 Spherical-Harmonic Terrain Coefficients		http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/first release.html
Document on the Development of the DTM2006.0	[251]	http://earth-info.nga.mil/GandG/wgs84/gravitymod/new_egm/EGM08_papers/ NPavlis&al S8 Revisedl 1 1606.pdf
Reference Station Positions
ITRF2008 Station Positions and Velocities		http://itrf.ign.fr/ITRF solutions/2008/
Open Access Article Describing the ITRF2008	[7]	http://link.springer.com/article/10. 1007/s00 190-01 1-0444-4
Description of the Solution (Software/Technique) Independent Exchange (SINEX) Format		http://www.iers.org/IERS/EN/Organization/AnalysisCoordinator/SinexFormat/ sinex cont.html
Planetary Ephemeris Data, Software, and Documentation
SPICE Toolkit Software (with Documentation)		http://naif.jpl.nasa.gov/naif/toolkit.html
The DE430 Ephemerides for SPICE		http://naif.jpl.nasa.gov/pub/naif/generic kernels/spk/planets/
Documentation for the DE430 Ephemerides	[95]	http://naif.jpl.nasa.gov/pub/naif/generic_kemels/spk/planets/ DE430%20and%20DE431.pdf
Star Catalog Data and Documentation
The Hipparcos Catalog (The New Reduction, 1/311)		http://cdsarc.u-strasbg.fr/viz-bin/Cat ?cat=I/3 1 1
Description of the New Reduction	[346]	http://arxiv.org/abs/astro-ph/0505432
Validation of the Catalog	[347]	http://arxiv.org/abs/0708. 1 752
Extensive Documentation on the Original Hipparcos Catalog	[82]	http://www.rssd.esa.int/index. php?project=HIPPARCOS&page=Overview
Errata for the Original Hipparcos Catalog		http://cdsarc.u-strasbg.fr/viz-bin/getCatFile Redirect/?-plus=-%2B&I/239/errata.htx
Earth-Orientation Parameters, Software, and Documentation
Earth-Orientation Parameter Data According to the IERS Conventions		http://www.iers.org/IERS/EN/DataProducts/EarthOrientationData/eop.html http://www.usno.navy.mil/USNO/earth-orientation http://hpiers.obspm.fr/eop-pc/
The NGA’s Earth-Orientation Parameter Predictions (Not used in the Simulations)		http://earth-info.nga.mil/GandG/sathtml/eopp.html
Explanation of IERS Bulletins A and B Parameters		http://hpiers.obspm.fr/eoppc/bul/bulb/explanatory.html
Documentation on the Combined Solution C04	[30]	http://hpiers.obspm.fr/iers/eop/eopc04/C04.guide.pdf
Code to Add Tides when Interpolating EOPs (interp and Subroutines)		http://hpiers.obspm.fr/eop-pc/index.php?index=models
General Software for Astronomy and Time Conversions
The IAU’s Standards of Fundamental Astronomy Software	[146]	http://www.iausofa.org
Software for Optimal 2D Assignment
An overview of 2D assignment algorithms; the appendix in each paper contains Matlab code.	[53,56]
Simple Map Data
The Global, Self-Consistent, Hierarchical, High-Resolution Geography Database (GSHHG) (coasts and political boundaries)	[359]	http://www.soest.hawaii.edu/pwessel/gshhg/
Shape File Format Definition	[79]	http://www.esri.com/library/whitepapers/pdfs/shapefile.pdf
M_Map: Free Software to Read Shapefiles in Matlab ( Cannot be Used in a Commercial Product!)		http://www2.ocgy.ubc.ca/ rich/map.html

All software and data except for M_Map appear to be usable in commercial products with nothing more than an attribution to the source. However, developers should verify the licensing terms with an attorney before use in any product.

6

David Frederic Crouse

2. CELESTIAL AND TERRESTRIAL COORDINATE SYSTEMS

Due to the complexity of defining global celestial and terrestrial coordinate systems, international stan¬ dards form the backbone of national standards. For example, since the advent of GPS, one of the most frequently encountered terrestrial coordinate systems is the United States Department of Defense’s (U.S. DoD’s) World Geodetic System 1984 (WGS 84) [69], While the WGS 84 standard is developed and main¬ tained by the United States’ National Geospatial-Intelligence Agency (NGA),1 1 it is in close agreement with international standards. The latest version of the WGS 84 standard used by U.S. DoD agrees with the align¬ ment and location of the axes in the latest (2008) version of the International Terrestrial Reference Frame (ITRF) [66,68], which is a realization of the International Terrestrial Reference System (ITRS) defined by the International Earth Rotation and Reference Systems Service (IERS). Differences between the WGS 84 and ITRF2008 arise due to differences in points of reference used to implement the standards, as discussed in [136, Ch. 2.11], and both standards are identical within the precision bounds of the ITRF definition.

Other countries maintain similar levels of compliance with international standard. The Russian Feder¬ ation’s GNSS system, GLONASS,12 whose interface is described in English in [263], uses the Parametry Zemli 1990 (PZ-90) coordinate system [348], The latest version of the PZ-90 coordinate system (PZ-90.11) is consistent to a centimeter-level precision with the ITRF2008 at epoch 2011 [349], The GNSS system run by China is the BeiDou Satellite Navigation System (BDS), whose interface is described in [48], The BDS uses the Chinese Geodetic Coordinate System 2000 (CGCS2000), which is consistent with an older (ITRF97) version of the ITRF [365], The axes of all GNSS systems coincide with those of the ITRF2008. The Galileo Terrestrial Reference Frame (GTRF), the coordinate system of Europe’s GNSS system (the interface is described in [83]) which is still under development, agrees with the ITRF2008 to the millime¬ ter level [6], Many countries arc harmonizing their GNSS systems through participation in groups like the United Nations Office for Outer Space Affairs, International Committee on Global Navigation Satellite Sys¬ tems.13 Other countries, such as Japan, with its Quasi-Zenith Satellite System (QZSS) [148], and India, with the GPS-Aided GEO Augmented Navigation (GAGAN) system [180,267], use GPS augmented with their own satellites.

The ITRF and WGS 84 standards arc realizations of ECEF coordinate systems: Their origins coincide with the center of mass of the Earth and their axes rotate with the Earth’s crust. The axes arc what arc termed “geographical axes” in [237] in that they arc defined by the relative position and motion of numerous reference stations on the Earth’s surface. The ITRF is a realization of the ITRS in that it defines specific points of reference, whereas the ITRS defines the coordinate system in a more abstract manner. According to Chapter 4.1 of the (older) 2003 version of the IERS conventions [225], 14

A Terrestrial Reference System (TRS) is a spatial reference system co-rotating with the Earth in its diurnal motion in space. In such a system, positions of points attached to the solid surface of the Earth have coordinates which undergo only small variations with time, due to

uThe NGA was formerly known as the National Imagery and Mapping Agency (NIMA), which itself formed in part from the Defense Mapping Agency (DMA). A brief history of the agency can be found at [246],

12In Russian: rjIOHACC=rjI06ajibHaa HABiiraiproHHan CnyTHHKOBaa CucTeMa or Global Navigation Satellite System.

13The harmonization of satellite navigation systems simplifies the fusion of signals from multiple systems for more robust localization. For example, Apple’s iPhone starting with the 4S model uses the United States’ GPS jointly with Russia’s GLONASS system [273] for improved localization reliability.

14The current (2010) version of the IERS conventions [257] does not redefine a TRS and a TRF; it simply does not provide as concise an explanation of the two systems.

Coordinate Systems for Target Tracking

7

geophysical effects (tectonic or tidal deformations). A Terrestrial Reference Frame (TRF) is a set of physical points with precisely determined coordinates in a specific coordinate system (Cartesian, geographic, mapping...) attached to a Terrestrial Reference System. Such a TRF is said to be a realization of the TRS.

On the other hand, the term “WGS 84 standard” refers both to an abstract standard described in [69] and to its implementation. As mentioned in [66], implementations of the WGS 84 in terms of reference points often have names like G1150, where ‘G’ refers to the use of GPS measurements to obtain coordinates and ‘1150’ is the GPS week number when the coordinates were implemented into the NGA’s GPS ephemeris estimation process.

When describing celestial bodies, the U.S. Navy’s astronomical almanac uses a realization of the Inter¬ national Celestial Reference System (ICRS) as its primary coordinate system [70, Ch. B, L], Both the ICRS and ITRS arc defined by the IERS in their conventions [257], The ICRS defines the alignment of the axes of a celestial coordinate system (that does not rotate with the Earth). However, the Geocentric Celestial Ref¬ erence System (GCRS) and the Barycentric Celestial Reference System (BCRS) place the ICRS axes either at the center of mass of the Earth, for the GCRS, or at the center of mass of the solar system for the BCRS. Both the BCRS and GCRS arc defined within a general relativistic framework, which necessitates the def¬ inition of specific relativistically consistent timescales. Section 3 maps out the international organizations that play fundamental roles in defining celestial and terrestrial timescales.

When describing the location of objects on or near the Earth, non-Cartesian coordinate systems, such as the familial- geodetic latitude and longitude, are often used. Section 4 describes a number of common coordinate systems for specifying the location of points on or near the Earth and how they relate to the approximate gravitational shape of the Earth. Note, however, that a very large number of alternate, local implementations of terrestrial coordinate systems exist in different countries for legal as well as practical purposes. Countries generally define legal boundaries in terms of geographic landmarks or in terms of geodetic latitude and longitude using their own local datums (manner of defining horizontal and vertical locations) to maintain a consistent set of boundaries in spite of plate tectonics. For example, certain reference stations in Australia moved approximately 2cm between 1996 and 1998 [329], and if national boundaries were defined solely in terms of coordinates in the ITRF2008, then the western pari of France and Spain would be slowly becoming ocean at a rate of centimeters per year- while Russia consumes the eastern pari of Europe as the Eurasian plate moves East, as can be seen from Fig. 1, which shows the locations and direction vectors of the stations defining the ITRF2008.

Thus, the European Union has the European Terrestrial Reference System (ETRS), a datum defined such that coordinates throughout the European Union do not change much over time. The ETRS is realized through the European Terrestrial Reference Frame 1989 (ETRF89) (see http : //www. epncb.oma.be ), which can be transformed to the ITRF [117]. However, individual countries within the European Union define local coordinates in terms of their own datums. For example, SAPOS15 provides coordinates within the ETRS for places in Germany [1 17].

Even in the United States, where the WGS 84 coordinate system is a DoD standard, the National Geode¬ tic Survey (NGS) (see the web site at http: / /www. ngs.noaa.gov ) defines the horizontal datum for

15SAPOS stands for “Satellitenpositionierungsdienst der deutschen Landesvermessung,” which means satellite positioning ser¬ vice of the German national survey.

8

David Frederic Crouse

ITRF2008 Stations at Epoch 2005.0

J_ I_ I_ I_ I_ I_ L

-150 -100 -50 0 50 100 150

Longitude

Fig. 1: Location points and direction vectors indicating the relative degree of motion in the local tangent plane (with the vertical defined by the WGS 84 reference ellipsoid) of the stations defining the ITRF2008 at the epoch of the data 2005.0 (the timescale was not specified; presumably TT was used). Red points and arrows refer to GNSS stations, green to very long baseline interferometry (VLBI) stations, dark blue to satellite laser-ranging (SLR) stations and dark blue to Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS) stations. Station motion is on the order of centimeters per year, presumably a Julian year in terrestrial time (3 1 , 557, 600 seconds per Julian year). The data for the outlines of the continents are from the Global Self-Consistent, Hierarchical, High-Resolution Geography Database (GSHHG) listed in Table 2 along with links to software for reading one type of file format of the data. (The software was tested, but not used to make this plot. Rather, the ASCII version of the data was read directly into Matlab and plotted, with care taken at the +180° boundary.)

the United States to be the North American Datum of 1983 (NAD 83). 16 Consequently, a point in latitude and longitude on the WGS 84 reference ellipsoid (as one might obtain using a GPS receiver) will not align with a point having the same latitude and longitude in the NAD 83 standard [306]. A discussion of the time dependence of datums around the world is in [301]. To convert between a number of different local and global datums as well as many different map projections, one can use the free Mensuration Services Pro¬ gram (MSP) Geographic Translator (GEOTRANS) that is offered with source code in C++ and Java from the NGA at http : //earth-info. nga.mil/GandG/ geotrans/ and is documented in [17, 18]. Al¬ ternatively, one can perform a number of datum conversions using the function in the NGS Geodetic Tool Kit, which can be downloaded from http : // www.ngs.noaa. gov/ TOOLS / and is described in a number of articles in Professional Surveyor starting with [368], which are all available on the NGS web site. Datum conversions are constantly updated and are not described in this report. The European Petroleum Survey Group (EPSG) of the International Association of Oil and Gas Producers (IOGP) maintains an online reg¬ istry of datums at http: / /www. epsg.org that assigns each datum a spatial reference system identifier

16The current version is NAD 83(2011).

Coordinate Systems for Target Tracking

9

(SRID). For example, the WGS 84 standard has identifiers 7030, 6326, and 4326, each addressing different aspects of the standard.

The IERS conventions [257] form the basis of the most important Cartesian celestial and terrestrial coordinate systems. A large number of organizations play a role in contributing to the celestial and terrestrial standards set by the IERS, as illustrated in Fig. 2. The centers of the International Gravity Field Service (IGFS) and the IERS Product Centers are all sources of geophysical data. The Commission on International Coordination of Space Techniques for Geodesy and Geodynamics (CSTG) is an example of one of many entities that work to standardize techniques between space geodetic organizations.

Fig. 2: A subset of the tangled web of international organizations behind many of the standards that form the basis of terrestrial and celestial coordinate systems. The IERS is the primary organization in internationally standardizing celestial and terrestrial co¬ ordinate systems, which are described in its conventions [257], The ICSU WDS is a multidisciplinary organization that supplanted the Federation of Astronomical and Geophysical Data Analysis Services (FAGS). The ISG used to be known as the International Geoid Service (IGeS). A number of professional organizations, such as the Institute of Navigation (ION), as well as many govern¬ mental organizations, including the NGA and the Naval Observatory in the United States, also play a role in establishing standard coordinate systems.

10

David Frederic Crouse

A number of coordinate systems defined in the IERS conventions [257] are related to each other in Fig. 3. Though other celestial and terrestrial coordinate systems exist, such as the International Astronomical Union’s (IAU’s) galactic coordinate system [31], most such systems are not as useful to engineers designing tracking algorithms. The primary coordinate system for tracking distant objects in space, such as plan¬ ets, meteors, or stars, is the BCRS. It is the ICRS with barycentric coordinate time (TCB) defined as the timescale. Whereas in relativity theory there is no “absolute rest” so all motion must be defined in a rela¬ tive manner, there does exist an absolute state of non-rotation as exemplified by Mach’s paradox [187, Ch. 3.1.5], The ICRS defines a non-rotating, approximately inertial coordinate system whose origin is at the center of mass of the solar system, but does not explicitly define a timescale. Inertial or pseudo-inertial co¬ ordinate systems are useful for expressing orbital or ballistic dynamics as they avoid the need for “phantom” acceleration, such as the Coriolis effect, as explained in Appendix A.

Fig. 3: The base coordinate system for tracking in the solar system is the BCRS, which is the ICRS, which has no explicit time component, with TCB added as time. The ICRS is realized using star catalogs as the ICRF in radio wavelengths and the HCRF in optical wavelengths. Alternatively, the UCAC can be used in place of the HCRF for more precise results. The pseudo-inertial reference frame for tracking near-Earth objects is the OCRS, which is a type of ECI coordinate system. The GCRS is related to the BCRS by an affine transform, placing its origin at the center of the Earth and using an appropriately relativistically scaled timescale, the TGB. The GTRS is a fixed-Earth model of which the common realization is the ITRS, whose definition is consistent with many models of the past. The ITRS is realized through the ITRF, which is essentially the same as the current version of the WGS 84 standard, except for the definition of the reference ellipsoid. Conversions between frames can be complicated and are generally performed with the aid of two intermediate reference frames, namely, the Celestial Intermediate Reference Frame (CIRS) and the Terrestrial Intermediate Reference System (TIRS) [187, Ch. 9] [257]. ECEF coordinate systems are generally realized using a timescale other than TCG.

On the other hand, the GCRS has its origin located at the center of mass of the Earth and uses the geocentric coordinate time (TCG) scale but the orientation of its axes matches the orientation of the ICRS axes. In much of the target tracking literature, when discussing tracking satellites around the Earth or ballistic missiles, a nebulously defined Earth-centered inertial (ECI) coordinate system is used. The GCRS as realized using observation-based references can be a solid definition of a practical ECI coordinate system for tracking.

Coordinate Systems for Target Tracking

11

The timescales associated with each of the coordinate systems arc inherently linked to the location of the reference clock when determining the simultaneity of events to an extremely high precision. This is because, as discussed in [145], special relativity theory tells us that events that appear simultaneous to one observer might appear non-simultaneous to another moving observer, even after accounting for the propagation time of light. The offset between the clocks of the observers depends on their relative velocities. Thus, by defining the reference “clock” for TCB to be at the center of mass (or inertia) of the solar system and the reference clock for TCG to be located at the center of mass of the Earth, one defines time scales that are linked to spatial coordinate systems in which events can be deemed simultaneous or not. The relativity of simultaneity is usually ignored, but can be significant depending on the application. For example, as described in [4], one must account for the relativity of simultaneity to get GPS localization accuracy down to below one meter. The necessity of defining high-precision coordinate systems simultaneous with timescales arises from the link between space and time in relativity and theory and is discussed in more detail in Section 3.

Figure 4 shows the alignment of the axes in the BCRS and GCRS coordinate systems, which is that of the ICRS. The “celestial sphere” just represents sets of directions in space. The idea comes from the fact that the stars in the sky somewhat look as if they are placed on a faraway sphere (they are all around). The alignment of the axes is defined with respect to a specific date (called an epoch). The reference epoch for the ICRS is called J2000.0, which refers to noon (1200 hours) in terrestrial time (TT) on 1 January 2000 in the Gregorian calendar (the standard calendar used throughout most of the world today). Terrestrial time is shortly before noon time coordinated universal time (UTC) on that same date.

Fig. 4: The BCRS and GCRS coordinate systems have the same alignment of axes, that of the ICRS. In the BCRS, the origin is at the barycenter of the solar system. In the GCRS, it follows the center of mass of the Earth. The x-axis (sometimes labeled T) points approximately in the direction of the mean vernal equinox at epoch, and the celestial equator is approximately the mean equator at epoch. The mean ecliptic at epoch intersects the celestial equator along the x-axis. Over time, the ecliptic can shift. The z-axis is chosen to be in the direction of the northern hemisphere, perpendicular to the celestial equator. The y-axis completes a right-handed coordinate system. Right ascension a (celestial longitude) and declination 8 (celestial latitude) form a spherical coordinate system. The angle of the ecliptic to the celestial plane e is known as the obliquity of the ecliptic.

12

David Frederic Crouse

The orientation of the axes in the ICRS is nominally based on the orientation of the mean equator and dynamical equinox at the reference epoch. The x- and y-axcs arc in the plane of the celestial equator with the A-axis (also known as the origin of right ascension) pointing toward the dynamical vernal equinox at epoch.17 The z-axis is orthogonal to the plane of the celestial equator and points in a general northerly direction. The celestial equator is the projection of the Earth's equator onto the celestial sphere. The ecliptic is a circle projected on the celestial sphere that coincides with the orbit of the Earth around the Sun and is tilted roughly 23.4° from the celestial equator along the A-axis [285, Ch. 3.53]. The ecliptic is the same as the apparent direction of the Sun in the sky as the Earth orbits the Sun. An equinox is when the ecliptic intersects the celestial equator. That is, the equinox occurs when the Sun appeal's to pass through these points in the sky. There are two equinoxes each year. The vernal equinox occurs when the Sun appeal's to be going from South to North in the sky. In reality, the ecliptic is not a perfect circle and the equator of the Earth wobbles with respect to the celestial sphere. Thus, the ICRS is nominally defined in terms of mean values of these quantities.

As noted in [354], practical systems align their axes based on the locations of stars, so the precise definitions of the mean equator and dynamical equinox at epoch (which are complicated) do not matter in practice. The International Celestial Reference Frame (ICRF) is a realization of the ICRS based on observations of the radio-emitting stars [257, Ch. 2.2], with the current version being the ICRF2, which is documented in [147] and can be downloaded from http: //hpiers. obspm.fr/icrs-pc/ . When using visible-light stars, the Hipparcos Celestial Reference Frame (HCRF) is a realization of the ICRF [257, Ch. 2.2], documentation and data for which can be found at the links in Table 2. However, only the most distant stars make good points of reference, as their angular displacements in the sky are minimal over time. For example. Fig. 5 illustrates the predicted motion of the stars near Barnard’s star over a period of 80 years. Barnard’s star has the highest proper motion in the ICRS, though it is not the angularly closest star next to the Sun.18

Since the Hipparcos catalog does not have a sufficiently high density of stars for very precise alignment in a small region, the U.S. Naval Observatory Charge-Coupled Device (CCD) Astrograph Catalog (UCAC) might often be preferred for high-precision applications [370]. The UCAC (the current version is UCAC4) is very large and cannot be downloaded, but can be obtained for free on digital versatile disc (DVD) as described in http://www.usno.navy.mil/USNO/astrometry/optical-IR-prod/ucac . If one simply wishes to look up the identity of a single star in the sky, the Set of Identifications, Measurements and Bibliography for Astronomical Data (SIMBAD) database can be used online, accessible through the Strasbourg Astronomical Data Center (CDS) at http : // cdsweb.u-strasbg.fr , and also described in [358],

Whereas the axes of the BCRS and GCRS can, in practice, be determined from distant stars, the origins of the respective coordinate systems are more difficult to find. The origin of the BCRS is the barycenter of the solar system. However, the center of mass of the solar system is not something that can be directly measured from Earth, rather it must be inferred based on observations of planetary motion. The mathematics behind such determinations are very complicated due to the necessary inclusion of relativity theory [187, Ch. 4, 5,

17The symbol T, which also marks the A-axis in Fig. 4, represents the zodiacal constellation Aries. In ancient Greek times, the vernal equinox pointed in the direction of Aries. Nowadays, it points toward the constellation Aquarius due to precession (precession is discussed in detail in Appendix E). A brief discussion on the history of the constellations and the zodiacal precession is given in [280].

l sBarnard's star is too dim to see with the naked eye.

Coordinate Systems for Target Tracking

13

*

5.5

O 5

» O

4.5b °

4

° O O O

3.5 o

°o, °

I

O O

O O

o

8

-92 -91.5 -91 -90.5 -90 -89.5 -89

Right Ascension (degrees)

Fig. 5: The right ascension and declination in the ICRS of the stars near Barnard’s star (in red) over a period of 80 years based on data in the Hipparcos catalog (Barnard’s star is catalog entry 87937). One can see that Barnard’s star has a large proper motion (the overlapping red circles), whereas the other stars barely move. The declination of the star is increasing as time progresses. The plot demonstrates how not all stars are poor for use as reference points in defining a non-rotating celestial coordinate system due to their large motion over the years. The function iauPmsafe in the IAU’s Standards of Fundamental Astronomy (SOFA) library was used to change the epoch of the data. The points in red are at the epoch of the Hipparcos catalog (1994.25 TT), and 20 Julian years prior and 20, 40, and 60 Julian years posterior to the epoch of the catalog (a Julian year being precisely 365.25 days).

6]. In practice, precomputed planetary ephemeris tables,19 or programs for efficient ephemeris calculation based on such tables, are generally used. Given a set of ephemerides and a time, one can relate the position of the Earth (and a terrestrial observer) to the barycenter.20 The IERS conventions use the development ephemerides 421 (DE421) of the National Aeronautics and Space Administration (NASA) Jet Propulsion Laboratory (JPL) as the ephemerides for the ICRS standard [257, Ch. 3]. However, the DE421 ephemerides have been supplanted by the newer DE430 ephemerides, which can be obtained along with the code to use them through the site for the Spacecraft Planet Instrument C-matrix Events (SPICE) toolkit listed in Table 2. Note that the BCRS coordinate system uses TCB as its time parameter, whereas the DE421/ DE430 ephemerides use barycentric dynamical time (TDB) as described in [187, Ch. 9.5].

As most individuals designing target-tracking algorithms are concerned with near-Earth objects, the GCRS and Geocentric Terrestrial Reference System (GTRS)21 are the most relevant coordinate systems. Both coordinate systems have their origin located at the center of the Earth and use TCG as their timescale, the rate that a theoretical clock located at the center of the Earth would tick if the mass of the Earth were not present.

19Ephemerides are tables of the position of celestial bodies as a function of time.

20Rather than relying on Earth-based detections, in the future, spacecraft might make use of pulsars for navigation. In [24], it was theoretically shown that mean position errors as low as 5 km might be possible when using three X-ray pulsars for navigation in space.

21The GTRS used to be known as the Conventional Terrestrial Reference System (CTRS) [257, Ch. 4.1.4],

14

David Frederic Crouse

The GCRS is the ECI coordinate system most relevant for target tracking. The GCRS is not a truly inertial frame of reference, because the origin of the system is in constant acceleration about the Sun. How¬ ever, the equivalence principle of general relativity theory [187, Ch. 3.1.3, 3.1.4] states that locally,22 the effects of gravity are indistinguishable from the effects of an accelerated coordinate system. Consequently, an observer in free fall would be unable to perform local experiments to tell that he is not in a co-moving inertial coordinate system. Consequently, the GCRS is often used as an approximately inertial frame for satellite tracking or for describing the motion of ballistic objects near the Earth. Since the axes of the GCRS are oriented in the same manner as the ICRS, stars in the form of the ICRF2 or other star catalog can be used to determine the alignment of the axes.

The GTRS is essentially the same as the GCRS, except it rotates with the Earth. The GTRS is an ECEF coordinate system for tracking. The orientation of the axes of the GTRS are determined by the ITRS, which is the GTRS without the time component. The ITRS is realized through the ITRF [257, Ch. 4], whereby the center of mass and the rotational angle of the Earth are determined using relativistic modeling and measurements from a network of observation stations through the technique centers of the IERS: the International Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS) Service (IDS), which uses Doppler techniques for celestial and terrestrial positioning; the International Laser Ranging Service (ILRS), which analyzes the orbits of the Laser Geodynamics Satellites (LAGEOS); the International GNSS Service (IGS), which provides a number of services based on global navigation satellites; and the International Very Long Baseline Interferometry (VLBI) Service for Geodesy and Astrometry (IVS), which precisely measures the rotation angle of the Earth based on distant quasars. When high-precision conversions between celestial and terrestrial coordinate systems are necessary, external Earth-orientation parameters as described in Table 3 are needed. Though research such as reported in [45,229] has provided insight into how changes in the Earth's rotation axis are correlated with other observable effects, Earth-orientation parameters cannot be predicted over long periods of time and must be obtained through observation.

Table 3: Earth-Orientation Parameters Necessary to Precisely Link ECEF and ECI Coordinate Systems

Orientation Parameter	Description
AT	AT = r-rr — RjTl ■ This quantity varies with changes to the Earth’s rotation rate.
LOD	Length of day. The derivative of the AT with respect to TAI; proportional to the Earth’s angular velocity.
dX.dY	Celestial pole offsets. Offsets of the modeled rotation axis to the true rotation axis in the GCRS.
dx,dy	Coordinates of the pole. Offsets of the modeled rotation axis of the Earth to the true rotation axis in the ITRS (i.e., with respect to the Earth’s crust).

Celestial pole offsets are usually less than 1 mas in each coordinate and the coordinates of the pole are normally less than 600 mas in each coordinate as can be ascertained by plotting historical data from the sources in Table 2. On the other hand, AT is steadily increasing and as of 1 April 2014 was approximately 67.4s according to the U.S. Naval Observatory’s long-term parameters at the source in Table 2. AT is sometimes defined as a difference with respect to TAI or UTC and/or with a flipped sign. NOTE: As per the IERS conventions, tide and libration effects should be added to the celestial pole offsets via the function

referenced in Table 2 before use.

"That is, in a sufficiently small region such that gravity appears to be uniform.

Coordinate Systems for Target Tracking

15

Table 2 lists the sources of Earth-orientation parameters used in this report. If one does not need better than arcsecond precision, then the celestial pole offsets, the coordinates of the pole, and the length of day (LOD) can be assumed to be zero. On the other hand, a rough estimate of the AT parameter is generally necessary. A low-precision estimate is to just use the cumulative number of leapseconds added to UTC. Section 3.2 describes LOD and AT in more detail as the parameters relate to various timescales. Note that some of the Earth-orientation parameters from the sources in Table 2 have been modified to have short-term tides removed so as to simplify interpolation. Such tidal effects must be added back using the software listed in the table.

Figure 6 shows how to convert between celestial and terrestrial coordinate systems using functions from the International Astronomical Union’s (IAU’s) standards of fundamental astronomy (SOFA) library. The equinox-based transformations, namely, the J2000.0 dynamical reference system as well as the mean-of- date (MOD) and true-of-date (TOD) reference systems, arc older celestial coordinate systems that must be taken into consideration when using older celestial data. For example, the U.S. Navy’s astronomical almanac [70] lists approximations for the locations of the Sun and the Moon in “apparent” coordinates, which arc synonymous with the “true-of-date” reference system. When converting velocity and acceleration measurements between the ITRS and GCRS coordinate systems, one must account for the rotation rate of the Earth, which offsets the velocity and introduces additional accelerations due to the Coriolis effect, as discussed in Appendix A.

Figure 7 shows how to convert from a number of outdated coordinate systems to the ITRS. The True Equator Mean Equinox (TEME) system is used in the Simplified General Perturbations 4 (SGP4) orbital propagation algorithm, which is presented with code in [142,339]. The TEME system is only of modern interest because the U.S. Strategic Command (USSTRATCOM)/the North American Aerospace Defense Command (NORAD) publishes/has published23 satellite ephemerides in the form of “two-line element” (TLE) sets24 that are given in the TEME coordinate system. Figure 7 also shows an older TOD coordinate system of which one should be aware to avoid confusion with the modern TOD reference system when read¬ ing older literature. One should avoid using the coordinate systems in Figure 7 in new tracking algorithms as they arc based on outdated theories.

If one wishes to simply convert between GCRS and ITRS coordinates, then one might also want to consider the free25 Naval Observatory Vector Astrometry Software (NOVAS) [337], which has the functions ter2cel and cel2ter to perform the conversions in one step. The NOVAS and SOFA routines for ECEF/ECI conversion arc compared in [233] and are shown to convert directions between the coordinate systems to within microarcsecond precision of each other. The microarcsecond difference is because the SOFA code includes some minor corrections to nutation that arc absent from the NOVAS code.

The origin of the ITRF2008, which is the latest realization of the ITRS, is not defined as the true center of gravity of the Earth. Rather, it uses a position averaged over time [257, Ch. 4]. Consequently, deformations of the Earth due to tides and other effects can shift the true center of gravity with respect to the reference frame, an effect termed “geocenter motion.” Ground-based stations that form points of reference for satellite navigation systems to establish their coordinate systems with respect to the ground move relative to the true

23The TLEs are available from USSTRATCOM at http : //www. space-track.org (registration required) and for certain satellites from http : //www. celestrak.com/NORAD/elements/ .

24The format of a TLE is described in [71].

25The NOVAS and SOFA libraries are free and, though one should check with a lawyer, it appears that the code can be incorpo¬ rated into commercial products with just an acknowledgement to the source.

16

David Frederic Crouse

center of mass of the Earth due to the motion of the atmosphere, tides, and tectonic motion [363]. Models for such motion are given in [257, Ch. 7], [364, Ch. 5.4], and [28, Ch. 2.3].

If all receivers have the same Cartesian bias in their locations with respect to the average center of mass of the Earth, then one can still track orbital debris to a high precision relative to the stations. Consequently, it makes sense to track satellites not directly in the ITRF coordinate system but in a system defined by satellites, which orbit the center of mass of the Earth, which, as noted in [363], only varies a few millimeters from its mean value. As noted in [364, Ch. 5.4], if one localizes a station using GPS satellites, then the coordinates are consistent with the changing center of mass of the Earth. The use of a dynamic center-of- mass coordinate system is discussed in [363].

Finally, it should be noted that a number of physical effects must be taken into account when computing the observed location of an object in the sky. These include delays due to the finite speed of light and parallax when changing coordinate systems, atmospheric and general relativistic refraction, and special relativistic aberration. The effects arc described in somewhat more detail in Appendix I, where it is detailed how one can transform planet locations from the coordinate system in which they come (the ICRS of the J2000.0 coordinate system) into observations for an observer near the Earth. Additionally, Appendix J quantifies the magnitude of the errors introduced by special relativistic aberration when a satellite observes other satellites, as the effect is often overlooked and can produce significant time-varying measurement biases.

Coordinate Systems for Target Tracking

17

Inputs

Vector in GCRS

*GCRS

Terrestrial Time and UT1 TT1,TT2, UT11, UT12 Polar Motion Coordinates, xp,yp

Celestial Pole Offsets dX, dY

GCRS

Geocentric Celestial Reference System

Multiply by Bias, Precession- Nutation Matrix C

i*cirs = CrccRS

4a) Find the Earth rotation angle for the given UT1 time 0 = iauEraOO (UT11, UT12) ;

5a) Form the rotation matrix about the 2 -axis for the Earth rotation angle

cos(0) sin (6) 0 R,3(0) = — sin(0) cos(6>) 0

0 0 1

CIRS

Celestial Intermediate Reference System

i*tirs = R.3(0)rc.Rs

la) Get the CIO locator s and the (x,y) coordinates of the CIP iauXys06a (TT1 , TT2 , &x, &y, &s) ; 2a) Add in the CIP Offsets x += dX; y += dY;

3a) Get the GCRS-to-CIRS matrix, C iauC2ixys (x, y, s, C) ;

Multiply by Earth Rotation Angle Matrix R,3(0)

TIRS

Terrestrial Intermediate Reference System

Equinox-Based Transformations

6b) Form the rotation matrix about the z -axis for GAST

cos(gast) sin(GASi) 0 R-3(gast)= — sin(GAST) cos(gast) 0 0 0 1

Multiply by Polar Motion Matrix W

fiTRS = WrTiRs

ITRS

International Terrestrial Reference System

lc) Get the Terrestrial Intennediate Origin (TIO) locator sp=iauSp00 (TTl , TT2 ) ;

2c) Get the polar motion matrix W iauPomOO (xp, yp, sp, W ) ;

Multiply by

GAST Matrix R,3(gast)

Multiply by Bias Matrix B

TJ2000.0 = BrGCRS

J2000.0

J2000.0 Dynamical Reference System

Multiply by Precession Matrix P

l*MOD = P*~J2000.0

MOD

Mean Of Date Reference System

Multiply by Nutation Matrix N

rTOD = NrMOD

TOD

True Of Date Reference System

lb) Get the bias, precession and uncorrected nutation matrices

iauPn06a (TTl, TT2, &A ip, & Ae, &e, B, P, rbp, N, rbpn) ; 2b) Get the pole coordinates and find dZ

dZ = - 4- dX- 4- dY

X		O'
Y	= B'P'N'	0
Z		1

3b) Get the equator of date coordinates and the celestial pole offsets

dip —

4b) Get the corrected nutation matrix

iauPn06 (TTl, TT2, Aip + dip, Ae + de, e, rb, rp, rbp, N, rbpn) ; 5b) Find Greenwich apparent sidereal time (GAST)

GAST=iauGst06 (UT11, UT12, TTl, TT2, NPB) ;

~dX		'dX‘
de	= PB	dY
_dZ		dZ_

Fig. 6: Transformations between the GCRS, which is an ECI coordinate system, and the ITRS, which is an ECEF coordinate system. The modern transformation uses steps la-5a to get to the TIRS through the CIRS. An older, equinox-based method goes through three intermediate coordinate systems to get to the TIRS, using steps lb through 6b. Steps lc and 2c are used to compute the matrix to go from the TIRS to the ITRS. The functions referenced for some of the transformations refer to C-function calls in the IAU’s SOFA library when using the IAU 2006/2000A precession-nutation model. Appendix E describes the origin and difference between precession and nutation. The full transformation using the modern celestial intermediate origin (ClO)-based model is rlTRS = WR3(0)CrccRS- When using the equinox-based method, the transformation is rjxRS = WR3(GAST)NPBrccRS- Terrestrial time and universal time 1 (UT1), which is a measure of the Earth's rotation angle, are assumed to be given as two-part Julian dates. To convert terrestrial time to universal time, the function iauTtutl (TTl, TT2 , deltaT, &UT11, &UT12) ; from the IAU’s SOFA library can be used with the Earth-orientation parameter deltaT, which is the difference between TT and UT1. The polar motion coordinates and the celestial pole offsets are assumed to be given in radians. The inclusion of the celestial pole offsets in the equinox-based transformation is based on the technique of [175]. Inverse transformations are straightforward as all matrices are rotation matrices and the transpose of a rotation matrix is its inverse. Combinations of functions for more direct transformations between coordinate systems are tabulated in [137]. Note: the conversions for UT1 above are IERS-compliant. However, the WGS 84 coordinate system omits certain terms when computing UT1 for GPS satellite ephemerides [69, Ch. 2.1].

18

David Frederic Crouse

la) Get Greenwich mean sidereal time using the

IAU 1982 model

GMSTl 982=iauGmst 82 (UTll, UT12);

2a) Form the rotation matrix about the 2-axis

C0S(GMST1982) sin(GMST1982) 0

R-3( GMSTl 982) =

— sin(GMST1982) COS(gMST1982) 0

0 0 1

Get the polar motion matrix W using the IAU 1980 model
	cos(xp)	sin(xp) sin(yp)	sin(xp) cos(yp)
W =	0	cos(yp)	- sin(yp)
	— sin(xp)	cos(xp) sin(yp)	cos(xp) cos(yp)

Fig. 7 : The transformation between the TEME, which is a non-standard coordinate system used with the two-line element sets pub¬ lished by NORAD for the old SGP4 orbital propagator, and the ITRS, as ascertained from [340], Also shown is the transformation between an outdated TOD system, which one should be aware of to avoid confusion with the more modern TOD in Fig. 6. Note that an outdated MOD system (related to the TOD system by a nutation correction) also exists. Relationships between older and modem coordinate systems are mapped out in [341], The function references refer to C-function calls in the IAU’s SOFA library. The polar motion coordinates are two of the Earth-orientation parameters. Terrestrial time and UT1, assumed given as two-part Julian dates, are related as described in the caption to Fig. 6. The TEME coordinate system shown is an “of date” system. A “mean” TEME coordinate system is also mentioned in the literature, but appears to not be in use.

Coordinate Systems for Target Tracking

19

3. TEMPORAL COORDINATE SYSTEMS 3.1 Defining Basic Temporal Coordinate Systems

The precise timing reference mandated for use by the United States (U.S.) Department of Defense (DoD) is the U.S. Naval Observatory’s (USNO’s) master clock [338]. However, the U.S. Navy does not single-handedly determine the time in the United States, but rather works in collaboration with a number of international organizations to determine standard timescales. Figure 8 maps out part of the complicated web of international organizations and national laboratories that determine the time around the world. Two good monographs for understanding the evolution and construction of the standard timescales maintained by these organizations are [226] and [16].

Fig. 8: A subset of the web of international organizations and national laboratories that contribute to the definition of interna¬ tional time standards. The BIPM is the primary organization responsible for determining UTC, though the ITU-R is the primary international legal authority for changes in the standards. The abbreviations for the labs correspond to those associated with the published deviations of national timescales from UTC and TAI given in ftp : / / f tp2.bipm.org/pub/tai/publication/ cirt /cirt.30 9 as of 9 October 2013. For example, NRL and USNO refer to the Naval Research Laboratory and the U.S. Naval Observatory in Washington, DC; DLR refers to the Deutsche Zentrum fur Luft- und Raumfahrt in Oberpfaffenhofen, Germany, and OP stands for the Observatoire de Paris in Paris, France. A list of the acronym definitions is given in Appendix M.

The International Telecommunication Union (ITU) is an international agency created by the United Na¬ tions (UN) that, among other things, coordinates the shared global use of the radio spectrum, manages geo¬ stationary satellite orbital positions, and plays a significant role in defining legally binding international time

20

David Frederic Crouse

standards [23, 154]. 26 Specifically, Working Party 7A (Time Signals and Frequency Standard Emissions) of Study Group 7 (Science Services) of the Radiocommunication Sector (ITU-R)27 of the ITU is responsible for work related to time standards. However, the ITU does not solely determine the time internationally. The Convention of the Metre, an international treaty signed in 1875 and amended in 1921, 28 established the International Bureau of Weights and Measures (BIPM),29 which is currently located in Sevres, France, under the authority of the International Committee for Weights and Measures (CIPM), which itself is under the authority of the General Conference on Weights and Measures (CGPM), which is formed of delegates of all of the signing countries.30 The CGPM is presided by the president of the French Academy of Sciences.

The BIPM manages a number of international timescales, the most notable being international atomic time (TAI), terrestrial time, and UTC [10]. 31 UTC is TAI offset by leap seconds; terrestrial time is something of a frequency-stabilized version of TAI, but is offset from TAI by a constant number of seconds that has no relation to leap seconds and does not change. Some disagreement on the legal authority of BIPM over UTC, compared to the ITU, exists [265]. The U.S. Naval Observatory’s clock maintains a representation of UTC, which is fed into the BIPM. BIPM combines the results with clocks from many different laboratories around the world [338]. Abbreviations for participating laboratories arc shown in Fig. 8. BIPM then combines the time measures and feeds the results back to the national laboratories. The representation of UTC kept at a particular laboratory is given by UTC(k) where k is the abbreviation of a laboratory name [226, Ch. 13.4]. For example, UTC(MIKES) refers to the time maintained by the Center for Metrology and Accreditation32 in Finland and UTC(USNO) refers to the time maintained by the U.S. Naval Observatory in Washington, DC.

The time feedback process between national laboratories and the BIPM is outlined in [16, Ch. 7.2], [226, Ch. 13.5], [10]. The details of the methods for transferring time between stations can vary, but generally involves the use of satellites to compare times between different laboratories. The realizations of TAI from different laboratories are combined at the BIMP to get the echelle atomique fibre (EAL).33 The EAL does not form the official realization of TAI at the BIPM because it can jump discontinuously as new time differences from the national laboratories comes in. Thus, a steering algorithm is used with various frequency standards at the BIPM to obtain a realization of TAI with no jumps. TT is a smoothed form of TAI plus a constant offset of 32. 184s [256].

While not directly involved in timekeeping, the International Organization for Standardization (ISO) plays a role in defining a number of aspects of how times and dates arc formatted. For example, the ISO 8601:2004 standard34 specifies how dates, times, and time intervals should be specified.

26Many of the ITU’s time and frequency standards can be downloaded free at http : / /www.itu.int/ rec/R-REC-TF/ .

27The ITU-R was formerly known as the International Radio Consultative Committee (CCIR).

28The text of the amended Convention of the Metre is available (in French) at http: / /www. bipm.org/utils/en/pdf/ metre_convention.pdf .

29Many of the acronyms do not correspond to the words that they are abbreviating, because the acronyms come from French. Such acronyms are tabulated in Appendix M.

30 As listed at http : / / www.bipm.org/en/convention/member_states / , the United States is a member state, having signed in 1878.

3 'Greenwich mean time (GMT) is an obsolete time standard that came before UTC and that no longer exists. Presently, the term GMT is generally used to mean UTC, especially in the United Kingdom.

32The acronym MIKES comes from Mittatekniikan keskus in Finnish.

33No English translation of the scale’s name is in wide use. The term “echelle atomique libre” means “free atomic scale.”

34The ISO 8601:2004 standard can be purchased for 134 Swiss francs (about 147 U.S. dollars) at http: //www. iso.org/ iso/ cat alogue_detail?csnumber =40874 . It is the opinion of the author that charging such high prices for fairly basic standards impedes their widespread adoption, thus hindering the mission of the ISO.

Coordinate Systems for Target Tracking

21

Figure 9 outlines the relationship between timescales that tend to be the most relevant for target-tracking algorithms. The timescales used internally in existing and developing GNSS, GPS (United States), Galileo (European Union), GLONASS (Russia), and BeiDou (China) are particularly relevant since such satellites are often used to synchronize clocks. BeiDou time (BDT), GPS time (GPST), and Galileo system time (GST) are all offset from TAI; the number of leap seconds needed to obtain UTC is also broadcast by the satellites. The times of the satellite systems are only notionally related to TAI by the given offsets as shown; the real times in the satellites can vary but are kept within certain error tolerances. For example, the BeiDou interface control document [48] states that BDT will remain within 100ns of UTC as given by the National Time Service Center of the Chinese Academy of Sciences (NTSC), written as UTC(NTSC). Similarly, GPS timescale is held within 90ns of UTC(USNO) [1 1] though offset by an integer number of seconds.

Fig. 9: The relationship between major modem timescales. Functions in the IAU’s SOFA library to perform the transformations are indicated, when available. The Earth-orientation parameter AT must be obtained from an outside source to compute UT1 (see Table 2). The conversion from TT to TDB requires UT1 as an additional input to iauDtdb. TDB is a timescale used by NASA JPL that does not use the SI second as its base unit, though its relativistic correctness is the same as TCB [313], TDB is also sometimes referred to as ephemeris time (ET). TT is not precisely an offset form of TAI; it is computed by smoothing TAI over a period of time. The U.S. Navy’s astronomical almanac uses TT for expressing ephemerides [70, Ch. B], GPST, GST, and BDT are broadcast with leap seconds information and can thus be easily transformed into UTC. GLONASS time matches the UTC offset for the time zone in which Moscow, Russia is located.

TT, TCG, and TCB have very specific definitions in terms of relativity theory. All are “coordinate” times, which has a certain meaning in terms of correctness to relativity theory and that should not be confused with a “coordinated” time, which is a standard time formed from multiple clocks. Special relativity theory plays a role in establishing a temporal coordinate system because moving clocks tick more slowly. General relativity theory also plays a role in standardizing timescales because clocks can tick at different rates depending upon the distribution of massive objects around them. Near the Earth, clocks tick faster at higher altitudes. Figure 10 illustrates how a clock on a satellite at varying altitudes is biased compared to TCG. The relativistic bias causes a frequency bias, when viewed by a TCG observer, that translates into a range-rate bias when performing detection. Overviews of the role of relativity theory in time synchronization are given in [16, Ch. 3] [226, Ch. 7,8] and a good overview of basic relativistic physics and its relation to standard coordinate systems and dynamic models is given in [187].

22

David Frederic Crouse

(a) Bias After One TCG Day

10.

15

20 £ QJ BC

c

a

25 OS

30

(b) Bias of a 2 GHz Clock

Fig. 10: Time dilation resulting from general and special relativity cause clocks on satellites to become biased compared to TCG, as plotted. The plots are relative to a TCG observer because the rate difference of a clock on a satellite and an (inertial) geocentric observer is time-invariant. With respect to an observer on the ground or another satellite, however, special relativity dictates that the bias varies with time since the relative velocities of the objects vary with time. Figure (a) is the time on a satellite clock minus the TCG time after one TCG day. A clock offset causes an offset in frequency. Figure (b) assumes that the satellite broadcasts a 2 GHz signal. The figures show the clock bias and the bias in the range-rate measurement due to ignoring relativity. The altitudes range from 100km to 36,000km, which is just above the altitude of a geosynchronous equatorial orbit. Appendix C provides more details in how the plots were computed.

TT is based on the rate that a clock located on the geoid, a hypothetical surface of equal gravitational potential of the rotating Earth defining mean sea level (MSL),35 would tick (Section 4 discusses the geoid in more detail). Precise clocks at different altitudes on the Earth must correct for altitude to remain syn¬ chronized to TT (or TAI/UTC) when reporting to the BIPM. Measurement models including the effects of relativity on radar and optical measurements have been previously studied [158, 184,323] [187, Ch. 7]36 and must often be considered when designing high-precision tracking algorithms for objects in space. Rel¬ ativity theory affects more than just the rate at which different clocks in a detection system tick. Indeed, as explained in [145], special relativity can paradoxically cause non-collocated events that are simultane¬ ous in one frame of reference to not be simultaneous in another frame of reference. Relativist corrections are essential for the accurate operation of GPS receivers [13] [364, Ch. 5.3]. For example, GPS signals are broadcast on 20.46MHz wide bands at 1575.42MHz and 1227.6MHz carrier frequencies, which are derived from a common 10.46MHz source [11, Sec. 3.3. 1.1]. The source and carrier frequencies are only nominally defined for an observer on the ground. The true rate of the driving clock inside the satellite is 10.22999999543 MHz, whereby the 10.46 MHz comes from using special relativity to deal with relativistic

35Note that the both the gravitational potential of the geoid as well as the handling of the permanent tide must be provided to have a unique definition of MSL. The necessity of specifying the permanent tide model is often overlooked and is discussed in more detail in Section 6.

36Those interested in visualizing the effects of special relativity alone to gain a more intuitive understanding might have fun playing with some of the games created as part of the Massachusetts Institute of Technology’s Open Relativity project [230],

Coordinate Systems for Target Tracking

23

time dilation relative to an average ground-based observer. Thus, satellite -based observers trying to use GPS need to make different adjustments.

TCG is a timescale defined as the rate a clock at the center of mass of the Earth would tick if the Earth were not present (and there were no Earth-rotation effects). Similarly TCB is a timescale defined as the rate a clock at the center of mass of the solar system would tick if all of the massive objects in the solar system were removed. TDB is a timescale similar to TCB that has been linearly transformed to remain closer, on average, to TT than TCB does. No clock directly measures TCG or TCB. Such coordinate timescales must be computed rather than measured.37 When tracking targets in orbit around the Earth, TCG is a good approximation to a fully relativistically correct timescale to use in the tracker, whereas when tracking objects far in the solar system (or for very high-precision tracking of targets near the Earth ), TCB is often more precise.

Coordinate timescales stand in contrast to proper timescales in relativity theory. Proper time is the time that would be read by a clock that is collocated with an observer [187, Ch. 2.2.4. 6] [226, Ch. 7.4]. Observers moving at different relative velocities will not agree on the simultaneity of events. This comes up when using high-precision ephemerides. For example, the NASA SPICE toolkit, which can be downloaded from the link in Table 2, can be used with various sets of ephemerides to determine the location of the planets at a particular time. The time unit used for the ephemerides is TDB, which is a barycentric coordinate time linearly related to TCB. One could use the function iauDtdb in the IAU’s SOFA library to convert from TT at the location of an Earthbound observer to TDB. The function iauDtdb uses the algorithm of [84] and requires the geocentric location of the observer because the Earth rotates and observers at different locations experience different proper times, meaning that the relative offsets between their clocks and TCB or TDB differ. Coordinate time is any unambiguous reference frame with respect to which various proper timescales can be compared. The fact that rotation is included in the definition of TT means that the definition of the relevant geoid includes the effects of general relativity (from gravity) and special relativity (due to the motion). Such a definition of the geoid including the rotation of the Earth is present in the WGS 84 standard [69], which is the coordinate system underlying GPS.

The link between time, motion, and position means that in high-precision applications, temporal and positional coordinate systems cannot be treated separately. This connection between time and space is one of the many reasons why Fig. 8 includes the IERS as well as the IAU, the IGS, and the International Union of Geodesy and Geophysics (IUGG), who all play important roles in the development of interna¬ tional positional coordinate systems, which arc discussed in Section 2. Many such organizations are united through the ICSU. Despite its widespread standardization, UTC itself is not directly appropriate for use in target-tracking algorithms because it is a discontinuous timescale. UTC is occasionally adjusted with “leap seconds” to keep it relatively in synch with a timescale called Universal Time 1 (UT1), which is a timescale based on the rotation angle of the Earth with respect to distant quasars in space [226, Ch. 2.7, 14, 17] and to the location of a theoretical mean Sun. The distant quasars define the ICRF set by the IERS [257, Ch. 2.2], which is discussed in more detail in Section 2. Note that the current version (G1762) of the WGS 84 standard [69, Ch. 2.1] is not fully consistent with the IERS conventions, because it omits librations and polar motion when computing UT1. This matters when using government -provided precise ephemeris data

37It would be nice to have a universal timescale that can be directly measured anywhere in the universe regardless of the observer’s motion and position in gravitational fields. However, there do not appear to be any physical processes that can be used to define such a scale. In other words “stardates,” as used in the Star Trek television series, are going to have to remain fictional for the time being.

24

David Frederic Crouse

for GPS satellites. All simulations in this present report were done using an IERS-compliant routine for computing UT 1 .

Figure 1 1 illustrates the relationship between common physical, solar, and sidereal timescales. The sidereal timescales are defined only in terms of the rotational position of the Earth with respect to the vernal equinox and thus arc not directly coupled to the position of the Sun on a daily basis. Note that the position of the vernal equinox with respect to distant stars slowly moves, so a sidereal day does not put an observer back in the same position with respect to the ICRS.38 Figure 6 in Section 2 makes use of the sidereal timescales when considering equinox-based transformations, which, in light of the modern IERS conventions, are largely obsolete. Local apparent solar time (LAT) can be viewed as the time that one would measure based solely on the location of the Sun in the sky at a particular longitude. Note that as documented in [176], the East longitudes used for obtaining local mean sidereal time and local apparent sidereal time must be in the TIRS. In [224, Sec. 4. 2-4. 5], UT1 is presented as a definition of “Greenwich mean solar time” and is local mean solar time at 0° longitude.39 As UT1 is defined in terms of a rotation in the TIRS, one can assume that the observer’s East longitude for obtaining local mean solar time (LMT) must also be in the TIRS.

Fig. 1 1 : The relationship between the modern physical time, TT, the modem solar time, UT 1 , and older solar and sidereal timescales. Functions in the IAU’s SOFA library to perform the transformations are indicated, when available, though the inputs do not nec¬ essarily match the flow chart. For actual times, the East longitude of the observer is added under the convention that 2n radians (360°) is equivalent to 24 hours (86,400 seconds). Note that the East longitude is given in TIRS coordinates, not in ITRS coor¬ dinates, sidereal times are often given in radians, and LMT and LAT are often given in hours (0 — 24) and are thus periodic. The Earth-orientation parameter AT must be obtained from an outside source to compute UT1 (see Table 2). The 12 hour offset for LMT is because LMT and LAT days, when expressed as angles/hours start at midnight, not noon as in UT1. The equation of time must be obtained from ephemerides. It is tabulated in the Navy’s astronomical almanac [70, Ch. C],

38The term for one rotation with respect to the ICRS is a “stellar day.” The IERS lists a stellar day as 86164.098903691 s, in com¬ parison to sidereal day duration of 86164.09053083288 s and a conventional solar day (without leap seconds) of precisely 86,400 s. A list of such constants is available from the IERS at http : / / hpiers.obspm. f r /eop-pc/models /cons tants.html .

39 Days in UT1 begin at noon, not at midnight. Figure 1 1 has LAT and LMT offset from UT1 by 12 hours so that a zero hour of the day corresponds to midnight, not noon. This is the standard convention. However, if LAT and LMT are represented as a full Julian date, it makes more sense to not add the 12 hours, as the other common date formats (UTC, UT1, TAI, and TT) all reference Julian dates from noon, keeping with older conventions.

Coordinate Systems for Target Tracking

25

The definition of LAT is based on the "apparent” location of the Sun. However, as defined in [332, Appen. C ], apparent places in astronomy take into account light-time and stellar aberration, but do not consider atmospheric refraction. LAT is not exactly the time one would read from a sundial, but it is close. One can find LAT (as an angle) by finding the local hour angle of the Sun in the TIRS and adding n radians. The local hour angle of the Sun is found by converting the direction of the Sun in the TIRS, accounting for aberration and light-time but not for atmospheric refraction, into spherical coordinates. The negative azimuth angle in spherical coordinates is the hour angle. The local hour angle is the hour angle of the Sun plus the East longitude of the receiver in the TIRS (not the ITRS). Conversion to a time of day occurs with the rule that 24 hours is 2k radians.

When a leap second is added, it occurs in all time zones around the world simultaneously (time zones arc just offset from UTC by various increments). Leap seconds arc added in the final minute of a particular day. That is, a minute might have 59s or 61 s rather than 60s. UTC cannot be set to always coincide with UT1 because UTC is meant to tick at a constant rate (minus intercalary [leap] seconds added to keep it in synch with UT1), whereas the duration of a day in UT1 varies because the instantaneous angular velocity of the Earth itself varies as a function of the positions of the Moon and other bodies in the solar system, as well as due to processes internal to the Earth. For example, the length of the day is estimated to have increased between 2.7 jits and 6.8 / is as a result of the 2004 Sumatra earthquake [1 14].

The accuracy and stability of the timescale needed in different parts of a multistatic tracking system vary. For updating target tracks, a relatively high degree of synchrony between sensors is required. In most cases, synchronizing the system to the clocks on the United States’ GPS satellites will provide sufficient accuracy for tracker updates. For example, consider the case of accurately tracking a meteor. In [218], it was noted that very small meteors (with masses on the order of micrograms) detected by the Arecibo radar were observed traveling at speeds up to 63km/s (%Mach 185), which is an order of magnitude faster than any aircraft or manmade satellite. In [226, Ch. 16.3], it is noted that timing accuracies as low as 20ns can be obtained using GPS satellites and timing accuracies accuracies below 1 ns can be achieved using two-way satellite time and frequency transfer techniques. For an object traveling 63km/s, a synchronization error of 20 ns corresponds to an offset of 1.3 mm in the position of the object (the amount the object moves over 20 ns); a timing error of Ins is a position offset of just 63 jUm. In both instances the bias is probably significantly smaller than the precision of the detecting radar.

When performing detection, the clock in the receiver needs to be very stable to measure Doppler during the measurement period. Frequency stability during such a short interval stands in contrast to the long¬ term time accuracy that was just discussed. In a bistatic system, frequency stability means that the clocks across sensors should be highly syntonous (have high short-term stability). For example, when transmitting a signal at a 2 GHz carrier, an offset in the frequency of the transmitting and receiving clocks of just 1 .05 Hz is the same frequency offset due to the Doppler shift of a target moving 16cm/s. Small Doppler biases also become important in Doppler-only localization applications such as in the IDS, which is a network of ground-based transmitters used by satellites for orbit determination.

When a time-synchronization accuracy only on the order of a few milliseconds is needed, as might be the case in an aviation radar network, then network time protocol (NTP) servers can be used for clock synchro¬ nization. Modern operating systems typically come with the software necessary to synchronize the clocks in computers to NTP servers over the internet. Many open-access servers, such as time-b.timefreq.bldrdoc.gov at the National Institute of Standards and Technology (NIST) in Boulder, Colorado, USA or ntp.neu.edu.cn at Northeastern University, Shenyang, Liaoning, China are available. Modern operating systems including

26

David Frederic Crouse

Mac OS X, Windows, and many versions of Unix and Linux come with software for querying network time servers.

In many applications the relative offset between two highly precise clocks is more important than the synchronization of either clock to an external timescale. For example, utilizing general relativity theory, a pair of modern atomic clocks have achieved sufficient accuracy that local variations in the geoid height can be measured to an accuracy of 1 cm (using general relativity theory) by comparing the rates of clocks in two different locations using a fiber-optic cable connection [32].

Table 4 lists a number of resources related to time transfer. Precise synchronization of a network of clocks over a wide area using global satellite systems is a complicated process. Basic training materials on the problem can be downloaded from the BIPM at the web site listed in Table 4. In practice, times and positioning obtained using GPS satellites can be improved with an augmented system. The ephemeris data (information on the orbital trajectories of the satellites) and clock corrections broadcast by each satellite are limited in precision. Users with internet connections can obtain more accurate data using the services of the IGS or the Global Differential GPS (GDGPS) System of the NASA JPL, where sub-nanosecond timing accuracies are possible. However, the GDGPS service is not guaranteed to be free. The standard for NTP communications is specified by the Internet Engineering Task Force through a number of standards listed in Table 4.

Table 4: Time Transfer Resources

Description	Internet Address
BIPM training materials for time standards and time transfer	http ://w w w . bipm . org/ws/ AllowedDocuments .j sp ? ws = TAI_TR AININ G
BIPM publications, including current biases of national laboratories’ clocks from TAI and UTC	http :// www . bipm . org/en/publications/scientif/tfg . html
ITU standards related to time transfer	http://w w w . itu . int/rec/R- REC- TF
GPSTk, Code for GPS localization and synchronization	http://www.gpstk.org/
Source code for GPS algorithms published in GPS Solutions Quarterly	http://www . ngs . noaa.gov/gps- toolbox/
Public network time protocol server lists	http :// support . ntp . org/bin/view/Servers/
Internet Engineering Task Force standards related to the network time protocol	https://tools.ietf.org/html/rfc5905 http://tools.ietf.org/html/rfc5906 http://tools.ietf.org/html/rfc5907 http://tools.ietf.org/html/rfc5908
Real-Time Products of the IGS	http://rts.igs.org/access/
Information on the GDGPS System	http://w w w . gdgps . net /
Free online service for refining estimates from raw GPS data in Canada	http://www.nrcan.gc.ca/earth-sciences/products-services/land-geodetic-survey/ gps-processing/online-global-gps-processing/54 1 5

Parts of the GDGPS service are considered the intellectual property of institutions such as the California Institute of Technology and are not guaranteed to be free. However, everything else listed is free and the real-time products of the IGS are a free alternative to the GDGPS service. Current versions of Mac OS X, Windows, and many versions of Linux come with software for

synchronization to an NTP server.

3.2 Time and Earth-Orientation Parameters

The differences ?uti ~ Lt- Tjti — Lai - Tti —Itai> or those quantities with the signs flipped are various definitions of the AT Earth-orientation parameter mentioned in Section 2. Being able to relate one’s local

Coordinate Systems for Target Tracking

27

time to UT1 is very important for localizing oneself using the stars and for converting measurements from Earthbound observers into the GCRS or ICRS, the coordinate systems in which one is most likely to track satellites and celestial objects. For example, the maximum allowable offset between UTC and UT1 is 0.9s [23]. Using the rotational velocity of the Earth of the WGS 84 coordinate system (given in Table 5 in Section 4.2), the maximum allowable time offset corresponds to an angular bias of 13.527 arcseconds. Considering an equatorial observer, using the WGS 84 Earth radius from the same table, the same angular bias when localizing oneself using stars corresponds to an offset on the surface of the Earth of 418.591 m.4()

Figure 12 shows a plot of AT over time and a plot showing how UTC has been updated to keep in synch with UT1. In [208], a description of how AT is estimated based on observational data is provided. Though a number of papers try to fit curves to historical AT observations for purposes of interpolation and prediction [1 18, 156, 157, 238], the models are not accurate for long-term predictions and observational data must be used. Based upon past prediction methods, the authors of [208] concluded that predictions over a four-year interval have a root mean squared (RMS) accuracy of 0.4 s, with a peak of 0.8 s. However, in addition to determining AT by observing distant stars, AT can be approximated based on observations of the motion of GNSS satellites [169].

0 1000 2000 3000 4000 5000 6000 7000 8000 0 1000 2000 3000 4000 5000 6000 7000 8000

Days since January 1 1992 (UTC) Days since January 1 1992 (UTC)

(a) The AT = TT — UT1 Term (b) UTC — UT1 with Leapsecond Jumps

Fig. 12: Plots of (a) the AT Earth-orientation parameter and (b) the difference between UTC and UT1, which is discontinuous due to the inclusion of leapseconds in UTC. The data was taken from the U.S. Naval Observatory’s finals2000A. data file available at the link in Table 2. Interpolation between points, including the addition of tidal terms, was performed using Matlab code based off of interp . f that the IERS provides at the link in Table 2. Though the UTC day duration is not quite uniform due to the presence of leapseconds, the scale of the chart is such that it is not noticeable.

40Longitudinal differences can be approximated by setting a clock for noon when the Sun is highest in the sky at a reference location and then observing when noon occurs after one has moved to a different location. This requires precise timekeeping, which was only possible after John Harrison invented the marine chronometer —a major improvement over some of the bizarre superstitious methods that some sailors tried (unsuccessfully) to use to keep time over distances. A monograph on the topic including discussions on why sailors thought they could tell time using injured dogs is [303],

28

David Frederic Crouse

The AT Earth-orientation parameter is directly linked to the LOD parameter. The LOD parameter is defined to be

Lop* A _ ^(Tjti ~ Hai) dtjAi

(1)

Consequently, by integrating over the LOD parameter in TAI, one obtains a version of AT in TAI time. Since a clock in TAI tics at the same rate as one in TT, one can replace TAI with TT in Eq. (1) without changing the value of LOD for an equivalent time in TT. The LOD parameter is related to the rotation rate of the Earth as

co = co

N

(2)

where co is the instantaneous angular velocity of the Earth's rotation in radians per second, Tjayi TAI = 86400 s is the nominal length of one day in TAI, and coN = 72921151.467064 x 10-12^ is the nominal angular ve¬ locity of the Earth at epoch 1820 based on the duration of the mean solar day being 7ja>, tAi- Note that 0JvTc|ay xai A 2 Ji, because the Earth’s rotation around the Sun means that the Earth must rotate more than 2k radians between noon on one day and noon (as defined by the peak height of the Sun in the sky) the next day. The specific numerical value used is from http : //hpiers. obspm.fr/eop-pc/earthor/ ut 1 lod/UT 1 .html . Though a very complicated explicit formula for LOD as a function of various celes¬ tial parameters has been derived [352], no accurate long-term predictions exist.

3.3 Units of Time

In general, to maximize the portability of data across sensor networks, it is good practice to always use SI units, which are defined in [39], because all countries except for the United States, Burma (Myanmar), and Liberia have legally adopted the SI standard [43]. The SI unit of time is the second and NASA SPICE software defines ET as seconds in TDB past the J2000.0 epoch.

However, ephemerides arc often listed in terms of Julian dates, which arc specified as a count of days, where a Julian day contains precisely 86400s and a Julian year is precisely 365.25 Julian days. For example, noon in TT on the first of January in the year 2000 according to the Gregorian calendar (the calendar used by most countries) is JD 2451545.0 TT, where JD indicates Julian date and TT notes that the timescale is terrestrial time. An explanation of Julian dates and the related Julian calendar is given in [285, Ch. 1.25,12.7, 12.8]. Julian days begin at noon rather than midnight. A modified Julian date (MJD) is a Julian date minus 2400000.5. By subtracting off such a large constant term, Julian dates can be represented to a slightly higher degree of precision using the same number of bits. For example, a single-precision floating-point number has 6 digits of precision [149], meaning that the time resolution of a normal Julian date is worse than one day and that of a modified Julian date is 1/10 of a day. On the other hand, a double-precision floating-point number has 15 digits of precision [149], allowing for 864 pis precision for a regular Julian date and 8.64 pi s precision with a modified Julian date. Julian dates are used in the Navy’s astronomical almanac [70, Ch. B, L]. Gregorian calendar years coupled with offsets in terms of Julian days are used in the TLE set data format used by NORAD for specifying orbital parameters [339]. The IAU SOFA software as well as the U.S. Navy’s NOVAS software both use two-part Julian dates. The date is split into two double-precision floating-point numbers. The date can be split in any manner. However, it is common for the first number to be an integer number of days and the second number to be the fraction of the day. This means that 24 hours can be expressed with approximately 15 digits of precision (better than nanosecond precision).

Coordinate Systems for Target Tracking

29

It is recommended that, when possible, data in a target-tracking system be recorded with a full date, preferably a Julian date, since it is standard, rather than providing the time of the data as a simple offset from when the system was turned on. The use of true dates allows users to later associate the measurements with other relevant information collected at the same time. For example, radar measurements with appropriate global timestamps can be associated with ADS-B data that can serve as “truth” for assessing performance.

3.4 Time for Celestial Observations

When determining the location of a celestial object in the sky at a particular time, the finite speed of light must be taken into account. One does not see the current location of the object. Rather, one sees where the object was when the light being observed left the object. Appendix I describes how the light-time correction is performed (with an approximate general relativistic correction) for an observer on the Earth using NASA’s SPICE toolkit, as the toolkit only provides estimates with respect to the center of the Earth. Uncertainty in the time light hits the target leads to uncertainty in the time to which a tracking algorithm should be predicted for a measurement update. Only a very small number of recursive tracking algorithms have been published that account for the light-time delay [97,231]. When given full (bistatic range and angle) measurements, one can use the apparent time the light hit the target. However, certain high bistatic angle arrangements between the transmitter and the receiver can be particularly sensitive to error when using the apparent time [52].

Note, however, that analogous speed of light delays for the effects of gravitational forces should not be used when integrating using common Newtonian or post-Newtonian dynamic models to determine the orbits of objects. The IERS conventions suggest a parameterized post-Newtonian model for orbital dynamics [257, Ch. 10]. The papers [216,344,345] debate whether the “speed of gravity” is infinite. When using Newton’s law of universal gravitation to compute the orbits of the planets, very large errors are introduced if one tries to limit the speed of gravity to the speed of light. However, the controversy is resolved in [40] where the general relativistic aberration of gravity explains the phenomenon: Orbits based on Newton’s law of universal gravitation assuming an infinite speed of gravity arc a good approximation to orbital dynamics using general relativity under which theory the speed of gravity is equal to the speed of light. More recent evidence considering the Shapiro time delay based on the bending of light from quasars around Jupiter confirms that under general relativity, the speed of gravity is finite [188]. Thus, when using a Newtonian or post-Newtonian gravitational model, gravitational forces should be taken as instantaneous. However, when using a full general relativistic model, gravitational effects travel at the speed of light.

30

David Frederic Crouse

4. COORDINATE SYSTEMS FOR DESCRIBING LOCATIONS ON AND NEAR THE EARTH

This section describes coordinate systems that are commonly used for describing points on or near the Earth. However, the results can be applied to other solar system bodies when given the appropriate data for the body. For example, whereas Earth gravitational models are discussed in this section regarding or¬ thometric heights, similar models are available for the Moon and other planets. For example, a recent spherical-harmonic gravitational model for the Moon called GL0660B is described in [186] and the associ¬ ated coefficients can be downloaded from the International Centre for Global Earth Models (ICGEM).

4.1 A Number of Common Spherical and Direction-Cosine Systems

A variety of non-Cartesian coordinate systems are used to describe the location of points with respect to near-Earth observers. Unfortunately, there is often little consensus on how many such coordinate systems are implemented. For example, Fig. 13 illustrates the various ways that a set of Cartesian axes is mapped to various coordinate systems for applications that might arise when performing tracking and navigation. When considering points on the Earth using spherical coordinates, one often uses longitude and geocentric latitude, as illustrated in Fig. 13. Note, however, that except when synthesizing magnetic and gravitational potential models, as discussed in Appendix F, points near the Earth arc generally given in longitude and geodetic latitude, which is an ellipsoidal coordinate system, discussed in Section 4.2.

z

Fig. 13: Various quantities related to spherical coordinate systems commonly used for target tracking and localization. In all instances, r is one-way range. When discussing points on the Earth, one often uses longitude A and geocentric latitude (j)c , where the z-axis is the axis of rotation of the Earth. When discussing the local horizon coordinate system, one often uses azimuth A measured clockwise from the y-axis (East) and an angle called the zenith distance f measured from the z-axis, which is the gravitational vertical. Alternatively, instead of the zenith angle, one sometimes uses the altitude a (which is another name for </>c in this system). When discussing measurements by a radar, one might use the r—u — v coordinate system, where the ;-axis is the pointing direction of the radar (boresight) and it is assumed that the z component is positive (in front of the radar). Alternatively, spherical coordinates for a radar often use azimuth 6 and elevation <p, where the ( x — z) -plane is level with the radar face. In some instances, 6 is measured clockwise as opposed to counterclockwise as shown, or a different axis is chosen to be the boresight.

Coordinate Systems for Target Tracking

31

The horizon coordinate system, which is mentioned in [332, Ch. 5.2.5. 1], is often used as a local coordinate system for astronomical observation. The z-axis is the gravitational vertical (the zenith), which points in the opposite direction of acceleration due to gravity. Such a system often specifies points using an azimuth angle A, which can be viewed similarly to a co-longitude, as it is measured from the y-axis clockwise in the (x — y) -plane, and an angle called a zenith distance £ measured from the gravitational vertical, or an angle called an altitude a (an elevation angle), measured from the (x — y) -plane.

In general, however, it is better to express directions using the zenith distance instead of the altitude angle as it eliminates possible confusion regarding the definition of the “horizon.” In the horizon coordinate system, the x and y axes are nominally located in the plane of the horizon. One understands the horizon to be the boundary between land and sky. However, refractive effects cause the apparent horizon to dip below the local tangent plane. Approximations to determine the dip angle as a function of the height of the observer are given in [332, Ch. 12.3.3.1], However, atmospheric turbulence can make observations near the horizon (or observations of the horizon itself) unreliable. For example, in [367], it is concluded that refraction within approximately 5° of the horizon is so variable that angular measurements of astronomical objects are limited to approximately one arcminute when considering visible light. Confusion when using altitude measurements can arise depending upon whether they are with respect to the apparent, refraction- corrupted horizon or with respect to the true geometric horizon. Though tracking and navigation algorithms such as those in [266] have made use of measurements of the horizon for the purposes of navigation and orientation estimation, horizon measurements are unreliable due to atmospheric turbulence. Moreover, the apparent location of the horizon that a sensor might detect can vary due to the presence of mountains versus ocean. Note that it is usually possible to see objects over the geometric horizon. Figure 14 illustrates this fact when considering a sunrise viewed from Hilo, Hawaii.

In addition to illustrating angular relations when discussing a spherical coordinate system for points on the Earth or astronomical observations, Fig. 13 also shows the angular relationship often used when discussing measurements by a radar. In the context of radar signal processing, the z-axis is often the pointing direction of the radar, and azimuth 6 and elevation (p are defined in completely different planes than are used when discussing longitude and geocentric latitude. Such a spherical coordinate system provides the basis of the r — u — v direction cosine system, which is described in the first tutorial [55] in the series [54,55,57],

For the purposes of this report, the spherical coordinate systems most relevant to target tracking are the longitude, geocentric latitude system, where a point consists of (r,A,0c), and the range-azimuth-elevation system ( r, 6. (p). The former is necessary for using gravitational- and magnetic-field models and the latter is commonly used to describe measurements by a radar.

Given a point (r, A , <j>c) in the longitude, geocentric latitude system at range r, longitude A , and geocentric latitude (j)c, the location in Cartesian coordinates is given by

x = rcos(A)cos(0c)	(3)
y = rsin(A)cos(0c)	(4)
Z = rsin(0c).	(5)

32

David Frederic Crouse

0.4

0.2

-0.2

o

©

c

o

-0.4

-0.8

Q

-1.2

-1.4

65.5

66

66.5

67

Degrees East of North

Fig. 14: A view of the rising Sun from Hilo, Hawaii (19.7056° geodetic latitude, —155.0858° longitude) at sea level (assumed to be the location of the ground, which is assumed to have zero WGS 84 ellipsoidal height), facing Northeast on 1 June 2013, at 15:42 UTC assuming an ambient temperature of 20° C, 70% humidity, and standard pressure (101, 325 Pa). The horizontal line is the horizon. The lower, blue circle is the outline of the Sun, ignoring atmospheric refraction, based on the DE430 ephemerides accessed using the SPICE toolkit available at the links in Table 2. The method of Appendix I for loading the Sun’s ephemerides was used along with the method of Appendix L for determining the outlining points of the Sun. The upper red circle is the apparent location of the Sun when using the numerical integration ray-tracing method of [138] and [332, Ch. 7.2] with the Sinclair atmospheric model of [295], An observer at sea level can see the Sun over the horizon. The look direction in degrees East of North and the vertical for determining the horizon are based on the ellipsoidal East-North-Up model described in Section 4.2 and not on a true horizon coordinate system.

Similarly, the conversion from Cartesian coordinates to spherical geodetic coordinates is

(6)

(7)

(8)

A = arctan2 (y,.r) = arcsin

where the function arctan2 represents a four-quadrant inverse tangent of y/x. In such a system, a local Cartesian set of orthogonal basis vectors can be defined at each point by differentiation of Eqs. (3), (4), and

Coordinate Systems for Target Tracking

33

(5) by r, X , and <j)c to get

dr

dX

=c 1U1

— rsin(A)cos(0c) rcos(A)cos(0c) 0

dr

d(j>c

dr

dr

=c iu2

=U3

— rcos(A) sin(0c) = — rsin(A)sin(0L.)

rcos(0c) cos(A)cos(0c) sin(A) sin(0c) sin(0c)

(9)

(10)

(11)

where m, uo, and U3 form an orthonormal basis and the constants c 1 and ci in Eqs. (9), (10), and (11) represent the fact that the vectors do not have unit magnitudes. This local coordinate system is the spher¬ ical equivalent to the local East-North-Up (ENU) coordinate system discussed in Section 4.2 when using ellipsoidal coordinates.

On the other hand, if a point is given in the range-azimuth-elevation system (r,6,(p) at range r, azimuth 6 and elevation <p, then the location in Cartesian coordinates is given by

v = rsin(0)cos(<p) y = rsin(tp) z = rcos(0)cos(<p)

so the inverse transformations are

r = \Jx1 +y 2 +z2 6 = arctan2 (x, z)

cp = arc sin

(12)

(13)

(14)

(15)

(16)

(17)

Though geocentric spherical coordinates are often used for physical Earth models, ellipsoidal coordi¬ nates, which use geodetic latitude, arc more commonly used for discussing the location of objects near the Earth. Section 4.2 discusses ellipsoidal coordinates and the reference ellipsoid, which is a good approxima¬ tion to the shape of the Earth.

4.2 Coordinates on the Reference Ellipsoid

The Earth is approximately oblate spheroidal in shape41 (an ellipsoidal shape), as illustrated in Fig. 15, which shows an ellipsoid with a graticule (lines of geodetic latitude [horizontal] and longitude [vertical] drawn). The hue shape of the Earth can be defined either in terms of the terrain or gravitationally, in terms of

41 The study of the shape of the Earth is termed geodesy; the use of computers to aid in geography is geomatics.

34

David Frederic Crouse

a hypothetical surface of equal gravitational potential called the geoid, which defines MSL. Understanding the geoid and geoidal heights is important for using terrain data and various physical models. Sections 5 and 6 discuss the gravitational shape of the Earth and the effects of the geoid and tides on coordinate-system determination.

z

Fig. 15: The Earth is approximately shaped like an oblate spheroid; that is, an ellipse rotated about the z-axis (a type of ellipsoid). The non-spheroidal nature of the Earth's shape has been exaggerated in the diagram. The horizontal lines represent lines of latitude; vertical lines represent lines of longitude. Longitude is an angle measured with respect to the.v-axis in the (x — y)-plane. It is analogous to right ascension in Fig. 4. This figure was also used in [58].

Geodetic latitude and longitude form the most common coordinate system used on maps and require the definition of a reference ellipsoid. A history of the development of the reference ellipsoid is given in [136, Ch. 2.11]. The WGS 84 standard [69, Ch. 3] and the IERS conventions [257, Ch. 4.2.6] define reference ellipsoids for the Earth to serve as horizontal datums. Datums serve as reference systems for positioning. The reference ellipsoid suggested in the IERS conventions is the Geodetic Reference System 1980 (GRS 80), which is a standard set by the IUGG [236]. Table 5 lists the defining parameters of the two aforementioned reference ellipsoids, as well as the parameters used specifically in ephemerides broadcast to GPS receivers and the parameters of two reference ellipsoids associated with the Earth Gravitation Model 2008 (EGM2008), which is used in gravitational models in this report.42 The non-geometric parameters GMfj and 0)5 are explained below. The subscript 6 on the parameters in the table is a symbol representing the Earth.43 Such subscripts are omitted when discussing reference ellipsoids in general.

42Though the EIGEN-6C4 model described in [102], which is available from the ICGEM's web site at http : // icgem.gf z- potsdara.de/ICGEM/ , is potentially more accurate than the EGM2008 model, the EGM2008 model is used in this report due to its widespread use and the availability of geoid models based on it.

43 An alternative symbol for the Earth is ©. A number of astronomical symbols exist, such as © for the Sun, <X for the Moon, d for Mars, and 9 for Venus. The IAU’s style manual [361] discourages the use of such symbols and recommends that they be defined if they are used. However, such symbols are commonly used as subscripts in mathematical expressions in technical articles.

Coordinate Systems for Target Tracking

35

Table 5: Parameters for Various Ellipsoidal Models

Quantity	WGS 84	WGS 84	GRS-80	EGM2008	EGM2008 for
GPS Navigation	(IERS Suggested)	Mean-Earth Ellipsoid	Geopotential Coefficients
GM&	3.9860050 x 1014	3.986004418 x 1014 s£	3.986005 x 1014		3.986004415 x 1014
as		6378137.0 m	6378137.0 m	6378136.58 m	6378136.3 m
1 /f&		298.257223563	298.257222101	298.257686
ujs	7.2921158553 x i()-5^ + 4.3 x 10-15 UT1f625f51M5 rad	7292115 X lQ-11^	7292115 x 10~14 — S

GM& is the universal gravitational constant times the mass of the Earth; a& is the semi-major axis of the Earth; f$ is the flattening factor of the ellipsoidal Earth; and 0)5 is the mean rotational velocity of the Earth. The WGS 84 reference is [69, Ch. 3]. The GPS version of it is used only when dealing with ephemerides broadcast by GPS satellites. The Geodetic Reference System 1980 is recommended by the IERS in [257, Sec. 4.2.6]. The parameters for the EGM2008 mean-Earth ellipsoid and the ellipsoid parameters for the EGM2008 geopotential coefficients are taken from the links for the EGM2008 model in Table 2. The values for the EGM2008 geopotential coefficients do not fully define an ellipsoid.

The reference ellipsoids are obtained by rotating a two-dimensional ellipse defined in the {x — y) -plane about the z-axis to make it three dimensional. Latitude and longitude define locations on the surface of a reference ellipsoid. Altitude (the height; not the angle), on the other hand, can be expressed in terms of a height above the ellipsoid, or in terms of another vertical datum, as discussed in Section 5. The WGS 84 standard [69, Ch. 3] defines the standard ellipse approximating the Earth shape (sometimes referred to as the “figure of the Earth”) to be centered at the center of mass of the Earth with the radius being the equatorial radius of the Earth a* and the inverse of the flattening factor being 1 //* as given in Table 5. For a generic ellipse or ellipsoid, the flattening factor is related to the equatorial radius a (distance from the center to a point on the equator), which is the same as the semi-major axis, and the polar radius b (distance from the center to the North or South pole), which is the same as the semi-minor axis, of the ellipse/ellipsoid by

/

a — b a

(18)

The difference between the equatorial radius and the polar radius of the reference ellipsoid is approximately U6 — b(, ~ 21.4 km, meaning that the Earth is nearly spherical. Points on a reference ellipsoid in three dimensions satisfy the equation

x2 + y2 z 2

a 2 +«2(1-/)2

(19)

Longitude is an angle measured counterclockwise (when looking down from the North pole) from the x- axis in the (jc — y) -plane, analogous to right ascension in Fig. 4.44 However, latitude is more complicated. To be able to properly locate points on a map and to use common gravitational models in tracking algorithms, one must differentiate between geocentric latitude (spherical coordinates) and geodetic latitude (ellipsoidal

44On other planets, planetocentric and planetographic longitudes can differ in direction. Planetocentric longitudes are measured counterclockwise from the x-axis when looking down from the North pole. On the other hand, planetographic longitudes are measured positively in the direction opposite the rotation direction of the planet [285, Ch. 4.22], The two coincide for the Earth.

36

David Frederic Crouse

coordinates), because both arc needed. Figure 16 shows the difference between geocentric latitude (j)c and geodetic latitude (j). Geocentric latitude is essentially the “elevation” in a spherical coordinate system. It is an angle measured up from the (jt — y) -plane (the equatorial plane). Geodetic latitude, on the other hand, for a particular point is obtained by drawing a line from the point perpendicularly intersecting the reference ellipse. The angle obtained due to the intersection of that line and the equatorial plane is the geodetic latitude.

z

Fig. 16: An illustration of how latitude is measured. When discussing the latitude of a point (the red circle), unless specified otherwise, one generally means the geodetic latitude <j>, which is taken with respect to a normal to the reference ellipsoid and generally does not intersect the origin. On the other hand, the geocentric latitude <pc, which is seldom used, is the angle measured with respect to the origin, h is the ellipsoidal height of the point. The horizontal axis is anywhere in the (.r — y) -plane. This diagram is also used in [58],

Unless otherwise specified, the latitude shown on maps is almost always geodetic latitude. While the definition of geodetic latitude is more complicated than that of geocentric latitude, before the advent of satellite geodesy using GPS and other means, geodetic latitude was significantly easier to measure than geocentric. A normal vector to a point on the ellipsoid of the Earth is approximately the gravitational vertical, which is what one perceives as “up,” so geodetic latitude can be approximately ascertained by measuring the positions of stars relative to the local vertical. However, high precision determination of the local “up” vector in a global coordinate system cannot rely on an ellipsoidal approximation. As mentioned in [87], gravitational deflections of the vertical from the ideal ellipsoidal model can be as high as 70 arcseconds (about 1.94 x 10~3 degrees) and are quantified in more detail in Section 5.4.

Figure 16 shows a point above the reference ellipsoid with its height h (the ellipsoidal height) measured with respect to the normal to the reference ellipsoid. The conversion of a point given in geodetic latitude,

Coordinate Systems for Target Tracking

37

longitude, and height (0, A,/z) coordinates to global Cartesian coordinates can be done as

x = (Ne + /j)cos(0)cos(A)	(20)
y = (Ne + /?.)cos(0) sin(A)	(21)
z = (Ne(l - e2) + h) sin (0)	(22)
a, a Ne — j - Y 1 — e1 (sin(0))2	(23)
e=y/2 f-p,	(24)

where 0 and A are the latitude and longitude, respectively, e is the first numerical eccentricity of the ellipsoid, and Ne is the radius of curvature in the meridian [285, Ch. 4.22]. On the other hand, the conversion from Cartesian coordinates is significantly more complicated and is given in Appendix B.

As described in [58], a set of local orthogonal basis vectors for defining a local ENU coordinate system can be obtained by differentiating Eqs. (20), (21), and (22) with respect to 0, A, and h to get

dr

dX =ClUl =

dr

d(j>

=C 2U2 =

— {Ne + h) cos(0 ) sin(A ) (Ne 4" It) cos(0 ) cos (A )

0

c°s(0 ) ~ (Ne + h) sin(0)^j cos(A)

c°s(0 ) ~ (Ne + h) sin(0)^ sin(A)

(Ne{l-e2)+h) cos(0) + (l -e2)^-sin(0)

dr

dh

=u3 =

cos(0)cos(A) cos(0)sin(A) , sin(0)

(25)

(26)

(27)

where

dNe nc2cos(0) sin(0)

= - r- Uo)

(l — e2 (sin(0))2^ 2

The constants c\ and ci in Eqs. (26) and (25) represent the fact that the vectors do not have unit magnitudes. To form basis vectors for the ENU coordinate system, note that ^ points up, 0 points in the direction of geographic North along the surface of the ellipse, and points East along the surface of the ellipse. The ENU coordinate system refers to the normalized forms of the bases in Eqs. (25) through (27).

The East vector m cannot be uniquely determined at the poles because has zero magnitude. To guarantee the existence of a local ENU coordinate system at all locations in ellipsoidal coordinates (0 , A , h), it is better to define m by the cross product

Ul = u2 X u3

(29)

38

David Frederic Crouse

to guarantee that it always exists. The longitude A determines the orientation of the North and East axes at the poles.

Table 5 provides more parameters than arc necessary to define the reference ellipsoid of the Earth. Specifically, some models define the universal gravitational constant times the mass of the Earth, GM&, and the angular velocity of the Earth's rotation, co$, in addition to the necessary semi-major axis as and inverse flattening factors I//5. Those two additional parameters suffice to define a full ellipsoidal gravitational approximation for the Earth. Taken together, the full set of parameters defining the WGS 84 and similar ellipsoids imply a surface of constant gravity potential (a geoid), where the potential is [136, Ch. 2.7]

U0

^ - arctan(e) + ^ co2a 2,

L-j J

(30)

where GM is the universal gravitational constant times the mass of the body in question and

* 1-/	(31)
7-1 A t —ae	(32)

arc the second e numerical eccentricity and the linear eccentricity E of the reference ellipsoid, respectively. The first numerical eccentricity e is given by Eq. (24). The gravity potential has units of watts per meter. As noted in [136, Ch. 2.7], when discussing gravity potentials other than those due to ellipsoidal approxi¬ mations, the symbol W is generally used instead of U . A subscript of 0 on U in Eq. (30) indicates that the value is taken on the surface of the ellipsoid.

Knowledge of the relationship between an ellipsoidal gravitational model of the Earth and the reference ellipsoid is important, as some sources do not explicitly provide the flattening factor of the reference ellip¬ soid. Rather, they provide a single spherical-harmonic gravitational coefficient from which the flattening factor must be deduced. This is the case in the appendix of [252], where instead of a flattening factor, the second-degree, fully normalized second zonal coefficient for a spherical-harmonic model of the Earth’s gravitational potential C2,o is given. Such models arc discussed in more detail in Appendix F. Here, the transformation of such a coefficient into the flattening factor for the reference ellipsoid is of interest.

The fully normalized second zonal coefficient for a spherical-harmonic model of the Earth's gravitational potential C2,o is related to the unnormalized second zonal coefficient C2.o, which is sometimes just called C2, as

C2 = C2, oV5. (33)

Similarly, C2 is related to the coefficient /2, which as in [202] is often referred to as a Jeffery constant, by

h = ~C2. (34)

Given /2, 00, and GM for a particular ellipsoidal-Earth model, the corresponding flattening factor of the Earth can be found through estimation of the first numerical eccentricity e of the reference ellipsoid using

Coordinate Systems for Target Tracking

39

the iterative relations from [236], which arc

(35)

(36)

(37)

where en is the first numerical eccentricity at iteration n. A reasonable starting value is eo = 0-8. One iterates Eqs. (35) through (37) until en no longer changes significantly. The first numerical eccentricity can be converted into the flattening factor for the ellipsoid using

f = l-VT^e>. (38)

V 1 - en

en+ 1 — A /3/2 + ~Tc

1 + arctan(<?„) —

o'-

( 4 \ ( co2a3el

V 15 y \2GMqn

When given with co and GM, reference ellipsoids directly define the gravitational potential and accelera¬ tion due to gravity anywhere on Earth. The direction of acceleration due to gravity under an ellipsoidal-Earth model is simply given by normalizing (27) to get 113. However, the magnitude of the gravitational accelera¬ tion under the ellipsoidal-Earth model can be useful for estimating the reading of a linear accelerometer that is on the ground, not accelerating, relative to the surface of the Earth. As derived in [136, Ch. 2-8], accel¬ eration due to gravity under an ellipsoidal model can be explicitly written in terms of ellipsoidal harmonic coordinates (j3,A,n), where /3 is the reduced latitude, A is the longitude, and u is the ellipsoidal harmonic semi-major axis. Appendix D provides the equations for converting to and from ellipsoidal harmonic co¬ ordinates. Given a point (j3,A,m), in ellipsoidal harmonic coordinates, and the linear eccentricity E of the reference ellipsoid, as well as associated values of GM and co, the acceleration vector g due to gravity seen by a stationary observer in an Earth-centered Earth-fixed frame of reference is

g = ypup + 7a ua + 7<ui

(39)

40

David Frederic Crouse

where up, u^, and u„ arc unit vectors in the directions of /3, A, and u as given in Appendix D and

» sm(/j)cos(p)

7a =0

Yu = -

1 / GM

(02a2Eqp \ / 1 .

w =t

w\n2+£'2 \qo(u2 + E2) I u2 + (E sin(/3))2

1

6

«2 + E2 u2 \

cb =3 1 + A3 1 “ A arctan -

40 =

4 =

-3

E3J

b 2

14-3 ^ I arctan El

u 2 \ i “ i \

1 + 3^ J arctan ( — ) — 3 — ) .

1

— sin^ ()3 ) — — H — (Q~u cos-(j8)

vv

(40)

(41)

(42)

(43)

(44)

(45)

(46)

Note that if the flattening factor is zero (the ellipsoid is actually a sphere) then E is zero and one will have 0/0 errors in the computation of various terms. In such an instance, one can use the asymptotic values. That is, in (40), substitute the value

4 b lim — = — 7 e-> o q0 iC

(47)

and in (42) substitute

lim

e-> o

Eqp

40

(48)

While the issue of a zero flattening factor would not arise when dealing with coordinate systems on the Earth, it can arise when dealing with coordinates on other celestial bodies. For example, in [274], a zero flattening factor is suggested for an ellipsoidal model of the Moon.

On the other hand, if the gravity acceleration is only needed on the surface of the reference ellipsoid go, then Somigliana’s formula can be used to easily find the magnitude of the gravitational acceleration vector. The direction of the acceleration vector is then the normalized version of the “up” vector in (27). Somigliana’s formula is [237, Ch. 2]

I go II =

qggcos2Q) +frgpsin20) \Ja2cos2((j)) +b 2 sin2(0)

(49)

where ge is the magnitude of acceleration due to gravity at the equator, gp is the magnitude of acceleration due to gravity at the poles, (p is geodetic latitude (the formula does not depend on longitude), a is the semi¬ major axis of the reference ellipsoid, and b is the semi-minor axis of the reference ellipsoid, which can be found by inverting (18) to get

b = «(!-/)•

(50)

Coordinate Systems for Target Tracking

41

If the gravitational potential due to the ellipsoidal-Earth model is needed everywhere in space, not just on the surface of the reference ellipsoid, one can use the formula from [136, Ch. 2-7] in terms of ellipsoidal harmonic coordinates

U

GM

= - arctan

C Cra2q

2qo

(, u 2 +£'2)cos2(j3).

(51)

Section 5 demonstrates the importance of gravitational potential and acceleration due to gravity in using data in different coordinate systems.

When considering the use of ellipsoidal latitude and longitude coordinates with common commercial products, one must be aware that a number of companies base their coordinates on a warped, incorrect map projection. The Mercator projection is one of the most commonly used map projections because it is conformal, meaning that angles are locally preserved, though distances are warped. However, many common online and mobile map applications use a warped Mercator projection called the “web Mercator projection” that was invented by Google, does not preserve angles, and on which latitude and longitude coordinates usually do not correctly correspond to the WGS 84 ellipsoidal coordinate system. The NGA has a presentation detailing why the web Mercator projection produces maps that are incompatible with the WGS 84 coordinate system, causing positional biases to be as high as 40 km in some areas [243].

Additional details and derivations related to the geometry of ellipsoidal-Earth models are given in [270, 271], including how parameters for the reference ellipsoid arc empirically found from measurements.

42

David Frederic Crouse

5. VERTICAL HEIGHT AND GRAVITY 5.1 Definitions of Height

Whereas altitudes (heights) in the IERS conventions (when not using a Cartesian coordinate system) arc defined as h with respect to the reference ellipsoid [257, Ch. 4.2.6], the WGS 84 standard defines vertical datums with respect to the geoid. The geoid is a surface of constant gravitational potential defining MSL, or a local measure of MSL. An older version of the WGS 84 standard used to make an exception when describing the depth of the ocean [67, Pgs. xii-xiv], where local sounding datums were suggested, but the current version suggests that all heights be reported at ellipsoidal heights [69, Ch. 10.6].

Figure 17 illustrates the differences between ellipsoidal and MSL (orthometric) elevations for a point on the surface of the Earth. The geoid is a hypothetical surface of constant gravitational potential45 chosen to approximately coincide with the height of the sea surface in the absence of tides. (This is one possible definition of MSL. Multiple definitions of MSL exist [75].) The distance between the geoid and the reference ellipsoid measured on a vertical from the reference ellipsoid is the geoid undulation N. The orthometric height H of a point is the distance measured from that point to the geoid by following the direction of gravity downward. This is the usual definition of “height” when one describes the height of a mountain, for example. The variable direction of gravity is referred to as the “curvature of the plumb line.” Though orthometric heights arc defined in a manner related to how acceleration due to gravity changes as one moves, surfaces of constant orthometric height are not surfaces of constant gravity or gravitational potential. The ellipsoidal height h is approximately the sum of the orthometric and geoidal heights.

Figure 1 8 plots (a) the EGM2008 geoid with magnified undulations and (b) the actual terrain according to the Earth20 1 2 model with altitudes magnified with respect to MSL. Based solely upon the distance from the center of the Earth, one might expect there to be no ocean near Europe based on the elevation in (b). However, since the geoid is high (red) near Europe, the level of the ocean by Europe is more distant from the center of the Earth than the level of the ocean in the Caribbean. Orthometric height determines how high something seems in terms of properties such as the thinness of the air. The air on the top of Mount Everest (the highest mountain in the world in terms of MSL) is thin because the top of Mount Everest is far above MSL. On the other hand, using the Earth2012 terrain data, one finds that the most distant point from the center of the Earth is not Mount Everest, but rather is located in the northern Andes mountain range in South America.46 This is because the Earth is roughly ellipsoidal in shape, so points near the equator tend to be farther from the center of the Earth.

Table 6 lists a number of different definitions of height/altitude that one might have to deal with when designing a target-tracking system. Absolute altitude is often the simplest height measure to use when tracking satellites and ballistic objects. Orthometric heights (MSL heights) arc seldom used as components of a tai'gct state in a tracking algorithm. However, they are often used in atmospheric models and terrain height models. For example, the Naval Research Laboratory Mass Spectrometer and Incoherent Scatter Radar Extended Model 2000 (NRLMSISE-00) empirical atmospheric model47 described in [74,261] is pa¬ rameterized in terms of orthometric heights. Similarly, the Global Multi-Resolution Terrain Elevation Data

45Contrary to what one might expect, the magnitude of acceleration due to gravity is generally not constant on a surface of constant gravitational potential.

46The farthest point front the Earth's center was the topic of the final question in the 2013 National Geographic Bee, with Chimborazo in Ecuador being the most distant peak from the Earth's center [14].

47The official Fortran release can be obtained from http://ccmc.gsfc.nasa.gov/modelweb/atmos/ nrlmsiseOO.html and an unofficial C-code release is at http : //www. brodo.de/ space/nrlmsise/ .

Coordinate Systems for Target Tracking

43

%

Fig. 17: A grossly exaggerated representation of the deviations between the reference ellipsoid, the geoid, and the terrain. The solid blue arrows represent the local gravitational vertical on the reference geoid (MSL) as determined by gravity. The elevation of a marked point on the terrain above the reference ellipsoid is h (ellipsoidal height), which is calculated from a normal to the ellipsoid. Similarly, the geoid undulation N is the height of a point on the geoid above the reference ellipsoid traced out by a normal vector. On the other hand, the orthometric height H is the height that would be measured from the point to the geoid by following the (curved) direction gradient from the point up or down until it intersected the geoid (dashed blue arrow). If the land were not there, it is the path that an ideal plumb line would follow. The distance N is often very close to N (both defined along line normal to the surface of the reference ellipsoid), so sometimes one approximates H % h — N.

2010 (GMTED2010) model provided free online48 by the U.S. Geological Survey (USGS) and the NGA provides terrain heights in terms of orthometric heights based on a geoid defined by the EGM96 gravita¬ tional model [61]. Terrain heights are important in target-tracking algorithms for predicting where ballistic objects will land (Will it go over the mountain or run into the side of the mountain?), and for determining whether terrain is obstructing the view of a target in a multistatic network. MSL heights are also used in the code released with the Enhanced Magnetic Model 2010 (EMM2010) (see Table 2 for documentation and a link to the code). The EMM2010 model is used in Sections 7 and 8.2, though implemented in a manner not using MSL heights. Pressure altitudes are commonly reported by transponders on civilian aircraft, including the common mode S transponder [150]. Geopotential heights are often used in atmospheric models. For example, the U.S. Standard Atmosphere 1976 [241] uses an approximate geopotential height.

Today, it would be significantly more convenient for terrain elevation data to be given in ellipsoidal heights when tracking in GTRS coordinates. However, prior to the use of satellites for measuring the ter¬ rain, the elevation of the terrain was determined through the direct or indirect measurement of changes in the gravitational potential on the ground [136, Ch. 4] [160], making elevations with respect to the geoid the ones directly measured. Researchers still make use of traditional leveling techniques for height deter¬ mination to validate and augment newer GPS-aided techniques, such as when Chinese researchers sought to more accurately determine the height of Mount Everest [44]. The use of gravitational potentials for al¬ titude measurements is unlikely to change for high-precision measurements, as modern atomic clocks have achieved sufficient accuracy that local variations in the geoid height can be measured to an accuracy of 1 cm by comparing the rates of clocks in two different locations [32].

48One can download the terrain elevation data from https : / / lta.cr.usgs.gov/earth_explorer .

44

David Frederic Crouse

(a) The Geoid, with Magnified Undulations (b) The Terrain, with Magnified Features

Fig. 18: In (a), the tide-free EGM2008 geoid with undulations magnified by a factor of 104 is displayed viewing North and South America. One can make out the Andes mountain range in South America. In (b), the Earth2012 terrain model with mean ocean, which is the distance of the terrain from the center of the Earth, is shown with the height above MSL (as given by the WGS 84 reference ellipsoid) exaggerated by a factor of 100. The spherical-harmonic coefficients used to make each plot were obtained from the sources in Table 2. The undulations and terrain heights were evaluated on a 2500 x 2500 grid in ellipsoidal latitude and longitude for the geoid and in spherical coordinates for the terrain.

5.2 Computing Orthometric Heights

A number of steps are necessary to convert between the different definitions of height in Table 6 and to utilize publicly available data models to reproduce the geoid and terrain plots in Fig. 18. A large amount of terrain data is only easily obtainable in terms of orthometric heights. For example, Fig. 19 plots the Digital Terrain Model 2006.0 (DTM2006.0), which is provided in terms of spherical-harmonic coefficients for orthometric heights with respect to the EGM2008 model. The coefficients can be downloaded from the link in Table 2. This section discusses the definition of the geoid and how geoid and orthometric heights are related.

A number of methods of estimating the geopotential at MSL are presented in [297], Some utilize satellite altimetry along with terrain density models and Earth geopotential models, and others simply try to average measures of the geopotential over the oceans. The definition of “mean seal level” that is suggested for use with the EGM2008 geopotential model is based on the geopotential implied by the EGM2008 mean-Earth ellipsoid, whose parameters are given in Table 5 and which is defined with the height-anomaly-to-geoid undulation coefficients at the link in Table 2. The documentation for the EGM2008 mean-Earth ellipsoid does not define the necessary GM& term or the rotation rate of the Earth (D& that are needed to use (30) to define the potential of the geoid. Thus, it is assumed that the values for the WGS 84 model from Table 5 are

Coordinate Systems for Target Tracking

45

Table 6: Common Height Measures Arising in Target Tracking Problems

Height Type	Description
Absolute Altitude	Distance from the center of the Earth
Orthometric Height	Height measured along the curved plumb line with respect to MSL
Pressure Altitude	Altitude based on a local pressure reading and a standard atmospheric model
Geoid Height Geoid Undulation	The distance between the reference ellipsoid and MSL (the geoid)
Geodetic Height Ellipsoidal Height	Height above the reference ellipsoid
Geopotential Height	Geopotential difference from surface, divided by a standard acceleration magnitude go

In many applications, only approximate values of orthometric height and geopotential height are used.

used. Using the EGM2008 mean-Earth ellipsoid with the additional parameters for the WGS 84 ellipsoid, (30) produces a gravity potential on the geoid of

m2

wEGM2008 = 6.263685 17 1456948 x 107 . (52)

s2

Note that the symbol Wo is used rather than Uq. The difference in notation shall be explained shortly. The gravity potential of the EGM2008 model is used to define the geoid in this report. However, different applications use different definitions. For example, the U.S. Department of Transportation’s Standard for Aeronautical Surveys and Related Products in the United States stipulates that orthometric heights shall be referenced to the North American Vertical Datum of 1988 (NAVD 88) [333],

Given a definition of the gravity potential of the geoid, to determine the height of the geoid above the reference ellipsoid using terrain and geopotential models, which is a necessary step in computing orthomet¬ ric heights, one must first differentiate between the different potentials that are listed in Table 7. The gravity potential W of the Earth at a particular point r = [rx,ry, r- \' in space under Newtonian mechanics is defined to be

Gravitational Centrifugal

Potential Potential

W= 'v^'T + '0(0' (53)

where V(r) is the gravitational potential at r, which is given by the mass distribution of the Earth, and Q(r) is the potential due to the rotation of the Earth. All of the potentials have units of m2/s2, which is the same as joules (energy) per unit mass in kilograms. Since the WGS 84 and ITRF standards set the z-axis of their ECEF reference frames to be the axis of rotation of the Earth, the centrifugal potential is

^(r) = \«>l(r2 + rf). (54)

The discussion in this section focuses on the potential of the Earth as expressed in Newtonian mechan¬ ics. Though more theoretically sound relativistic geodesy techniques exist as discussed in [187, Ch. 8],

46

David Frederic Crouse

Fig. 19: The DTM2006.0 terrain/bathymetric model plotted on a 2500 x 2500 uniform grid in latitude and longitude. The model provides terrain heights and bathymetric depths as orthometric heights. The tide-free EGM2008 geoid height was added to the orthometric heights (Helmert’s projection) to obtain ellipsoidal heights with respect to the WGS 84 reference ellipsoid to plot. The heights above the reference ellipsoid were exaggerated by a factor of 100 when plotted. At the current scale, an omission of the geoid height correction would not have noticeably changed the plot. However, to use the data to a high precision, one must consider the geoid height.

for example, they are complicated and Newtonian models appear to be the primary ones in use. The grav¬ ity potential W is the negative of the potential energy per unit mass of an object rotating with the Earth. Consequently, the gradient of the gravity potential is the gravity acceleration

g =VW(r)

=VV(r) + co|

rx

ry

0

where the gradient operator is defined in Cartesian coordinates as

V =

' _d_ _d_ d_y

drx ’ dry ’ drz

(55a)

(55b)

(56)

The ellipsoidal potentials U and Uq play an important role when one wishes to express geoid heights in WGS 84 coordinates even though the reference ellipsoid defining MSL for the EGM2008 model uses different values of the semi-major axis and flattening factors of the Earth. The ellipsoidal potential U (r)

Coordinate Systems for Target Tracking

47

Table 7: Commonly Used Symbols for the Different Earth Gravitational Potentials

Symbol	Name	Description
IV	Gravity Potential	Potential of an observer rotating with the Earth
V	Gravitational Potential	Potential of an observer not rotating with the Earth
a	Centrifugal Potential	Potential due to the Earth’s rotation
u	Ellipsoid Potential	Potential under an ellipsoidal-gravity model
Wq	Geoid Potential	Potential defining MSL
U0	Ellipsoid Surface Potential	Potential on the surface of a reference ellipsoid
w	Gravitational Acceleration	Acceleration due to gravity, not counting the Earth’s rotation.
g	Gravity Acceleration	Acceleration due to gravity, including the Earth’s rotation

Many papers do not differentiate between the gravity acceleration and the gravitational acceleration and use g for both. To avoid confusion, this report uses different symbols for each.

is the gravity potential at a point r under a particular ellipsoidal-Earth model, as given by (51), and Uq is the gravity potential on the surface of a specific ellipsoid as per (30). Usually, depending on the reference ellipsoid chosen, JJq A Wq.

At a point r chosen to be outside of all of the mass of the Earth, the gravitational potential V due to the Earth can be expressed in terms of spherical-harmonic coefficients [136, Ch. 1, 2] [29, Ch. 3.4]. 49 Though the spherical-harmonic representation is not technically valid at points within the mass of the Earth, it is usually assumed valid within the atmosphere. The gravitational potential of the Earth as represented by the NGA’s EGM2008 standard is given as a set of normalized spherical-harmonic coefficients. The same is true of the many other models listed by the ICGEM. The EGM2008 model is the recommended geopotential model in the IERS conventions [257, Ch. 6] and in the current version of the WGS 84 standard [69, Ch. 5] and it was shown to be the best available model based on laser satellite-tracking observations [307], though a newer model might be more accurate [102]. The mathematics of using spherical-harmonic coefficients can be difficult and are reviewed in Appendix F. It is worth noting that naive implementations of the harmonic equations can lead to a loss of precision [140] or can result in numerical overflow or underflow problems.

The geoid is always within the mass of the Earth as it is always within the atmosphere, and is often underground. In practice, the effects of the atmosphere are often neglected but the effects of terrain must be taken into account, necessitating information on both the height and density of the terrain. A commonly used technique for using gravitation and terrain models for determining the height of the geoid above a reference ellipsoid is that of Rapp [269]. Rapp's method assumes that the gravitational potential implied by the reference ellipsoid by (30) is the same as the geoid potential Wq. However, when using the geoid potential given by (52) for the EGM2008 model with the WGS 84 ellipsoid, Rapp’s assumption is false. Thus, a modified expression of Rapp’s method as given by Eq. 2-347 in [136, Ch. 2.17], as well as in [300], corrects for using a reference ellipsoid with an implied surface potential different from that of the geoid. The key equations for geoid height at point p, which might be expressed as (A. </y) in spherical coordinates

49Note that even if the gravitational potential everywhere outside of a massive body is known, one can not directly infer the mass distribution within the body. Additional information, such as seismic data, is necessary [136, Ch. 1.13].

48

David Frederic Crouse

or as (0, A) in ellipsoidal coordinates,50 when using the data provided with the EGM2008 model, are51

Np=£p+Kp r Tp Wo -Up

bP —

1 r\r r\/

(57)

(58)

where Cp is the height anomaly on the surface of the reference ellipsoid. Wo is the defining gravity potential of the geoid (for example, (52)), U p is the gravity potential on the surface of the reference ellipsoid at point p as computed using (30), yp is the magnitude of the acceleration due to gravity at point p on the surface of the reference ellipsoid using an ellipsoidal gravitational model via (39) or (49); that is,

Yp

A

go

(59)

The computation of go and Uq uses the value of GMs that is associated with the reference ellipsoid, not the value of GM that is used in the geopotential model, which might be different. The disturbing potential for a point on the surface of the reference ellipsoid is given by Tp, which is defined as

Tp = Wp - Up

(60)

where Wp is the gravity potential at point p as defined in (53). The disturbing potential Tp is often evaluated by modifying the fully normalized even spherical-harmonic coefficients Cn,n that are used to evaluate the gravitational potential V. Such coefficients are described in Appendix F. As described in [300], the modified spherical-harmonic coefficients to directly compute Tp are the same odd coefficients Sn_m as are used to compute V along with

rfp — r

n.m —

(jA/gHipsoid / flgHipsoid

GM* \ as

/^■ellipsoid

’-'n.m

(61)

where GMs and as are the the values of GM and the semi-major axis of the Earth as used by the gravitational potential model (e.g., the values for the EGM2008 model) and GMeuips0id and ^ellipsoid are the same as those used by the reference ellipsoid with which one wants the geoid undulation values (e.g., the values for the WGS 84 model).52 The value C*£psoid is a fully normalized spherical-harmonic coefficient for an ellipsoidal gravitational model, which, as can be ascertained from [136, Ch. 2.5, 2.7-2.9], is given by multiplying Cik.o in Eq. (F85) of Appendix F.4 by l/\/4k+\ for k ^ 0. Note that the values of C given in Appendix F.4 are only for even values of C„ m, hence the Cikfi notation. Value of C for other subscripts can be assumed to be zero.

Unlike as described in [300], the geoid height in (57) is not iteratively computed since the NGA has provided a set of spherical-harmonic coefficients to convert from height anomalies on the reference ellip¬ soid to geoid undulations. Kp is the 'C, — to — N coefficient that the NGA provides at a site listed in Table

50The lack of a radius or ellipsoidal height is because the geoid height does not depend on elevation. For all evaluations, the height is taken as 0, the surface of the reference ellipsoid.

5 ’When using the older EGM96 gravitational model, the value of Kp obtained using the coefficients provided by the NGA is scaled too large by a factor of 100. Thus. Eq. (57) would have Kp divided by 100 when using the older EGM96 model.

52In the NGA’s algorithm for evaluating geoid heights with the older EGM96 model, the F477 function at http : / /earth- inf o.nga. mil/ GandG/ wgs84 /gravitymod/ egm9 6/ egm96.html , it is approximated that the coefficients of are all 1 . This approximation is dropped in the reference implementation of the EGM2008 geoid.

Coordinate Systems for Target Tracking

49

2.53 The computation of Wp and Kp requires the evaluation of spherical-harmonic coefficients as described in Appendix F. To understand the origin of the model and the assumptions that must have gone into the Kp term, one is encouraged to consult [269] and [136, Ch. 2.17], A large number of undocumented assump¬ tions regarding the height and density of the land must have been used to obtain the spherical-harmonic coefficients for Kp.

Based on the code and documentation on the NGA’s web site, it appeal's that the NGA uses the approxi¬ mation that

^-Co.o— -0.41 (62)

7 ' Ve

where GM& is the value of GM as used by the gravitational potential model, Co o is a fully normalized spherical-harmonic coefficient as discussed in Appendix F, and r/? is the distance from the center of the Earth to point p on the surface of the reference ellipsoid, which can be found by taking the magnitude of the Cartesian point obtained by using (20), (21), and (22) to convert the ellipsoidal point (0,A, 0) (zero ellipsoidal height) from ellipsoidal to Cartesian coordinates. The 0.41 constant assumes that everything is expressed in SI units and represents an approximation to the second term in (58) plus the second term in (62).54 The NGA approximation was used when plotting Fig. 19. The difference between the NGA approximation and (58) when using the EGM2008 model is on the order of a few tenths of a centimeter, but is larger when using the older EGM96 gravitational model.

Geoid heights above the reference ellipsoid vary from from a high of approximately +84 m to a low of —107 m. The values can be determined either by evaluating Np directly as described or by using the precomputed geoid height values available from the NGA [240]. Consequently, the magnitudes of treating orthometric heights as ellipsoidal heights, such as when using terrain elevation data, is close to the same range.

The rigorous determination of orthometric heights is complicated by the presence of terrain and is de¬ scribed in [327]. When orthometric height H in Fig. 17 of a point is found by correctly following the curved plumb line from the geoid through any terrain up to the point, then one is discussing Pizzetti’s orthometric height. On the other hand, when one simply uses the approximation that

H = h — N (63)

one is said to use Helmert’s orthometric height. The value N is the geoid undulation at the latitude and longitude of the point, whereas N is the geoid undulation at wherever following the curved plumb line down from the point intersects the geoid. As discussed in [136, Ch. 5-2], the difference between Flelmert’s and Pizzetti’s orthometric heights is often negligibly small. As most applications using orthometric heights do not require high precision, and rigorously determining orthometric heights can be difficult, techniques for determining Pizzetti’s orthometric heights are not given in this report. In general, it is most convenient if those designing target-tracking algorithms entirely avoid the use of orthometric heights and simply express everything in terms of WGS 84 ellipsoidal or Cartesian coordinates.

53The NGA also provides precomputed tables of geoid undulations with interpolation software in Fortran.

54Under the older EGM96 model, the NGA used —0.53 instead of —0.41 . As documented in the latest WGS 84 standard [69, Sec. 6.2], the discrepancy is due to an error in an oscillator drift correction in satellite-measured altimeter data underlying the EGM96 model.

50

David Frederic Crouse

The discussion of orthometric heights thus far has assumed that “the” EGM2008 geopotential model is used for geoid determination. However, the unique definition of the geoid requires both a defining potential Wo, as in (53), and a specification of the tide model. As noted in [77], different countries use different tide models when defining orthometric heights. Orthometric heights arc always defined without consideration to the time-variant paid of a tidal model. However, time-invariant effects of the Sun and the Moon, permanent tides, also work to deform the Earth and change the gravity potential. The three permanent tide systems arc tide-free, zero-tide, and mean-tide. The various permanent tide systems are defined in Section 6 where the effects of time- variant as well as time-invariant tides arc analyzed.

The basic spherical-harmonic coefficients for the EGM2008 geopotential model (which lets one compute V in (53)) and for the Kp terms to go from height anomaly to geoid undulations are given in a tide-free geoidal model. The conversion of a geoid from a tide-free system to other tide systems requires the addition of extra terms as discussed in Section 6.

5.3 Considering Pressure Altitudes

A pressure altitude is an altitude that is computed based on an air-pressure measurement and an at¬ mospheric model. In the U.S., aircraft above 18,000ft (5484.6m) generally navigate based on pressure altitudes that have been calibrated to a standard setting for air pressure and temperature at MSL. No air- pressure model can perfectly match reality. However, when used in aviation, if all aircraft arc using the same incorrect model, they will experience the same model-based biases in altitude measurements and can thus maintain a safe vertical separation. The use of air-pressure measurements as reported by Mode C transponders in aircraft in the context of target tracking and sensor registration has been considered in [235]. In the U.S., the Federal Aviation Administration (FAA) stipulates [334] that altimeters in aircraft should be tested using the 1976 U.S. Standard Atmosphere of [241]. The U.S. Standard Atmosphere is also valid for civilian aviation outside of the U.S., because up to an altitude of 32km (which is above the ceiling of current commercial jets), the U.S. Standard Atmosphere is identical to the standard atmosphere of the In¬ ternational Civil Aviation Organization (ICAO), as noted in [8], where a number of atmospheric models are contrasted. Many commercial atmospheric pressure sensors, such as that described in [8], use the U.S. Standard Atmosphere to provide the user with an orthometric altitude reading.

For the purposes of target tracking using transponders, if one wants to use the altitude readings, one must convert the standard pressure altitude used by the aircraft into an altitude that the tracker can use. This generally means finding the actual air pressure measured by the aircraft from the reported altitude and then putting the air-pressure measurement into a more accurate model. This section reviews the pressure model of the U.S. Standard Atmosphere, so that one can extract air-pressure readings from transponder data. Due to the preponderance of high-fidelity atmospheric models,55 no specific model is considered for converting the aircraft’s air-pressure reading into an accurate height usable for tracking (such as an ellipsoidal height). However, using the U.S. Standard Atmosphere with an accurate local altimeter setting (rather than the standard one used above 18,000ft) can provide less-biased orthometric heights that can be converted into a format usable in a tracking algorithm.

To use the U.S. Standard Atmosphere 1976 to compute an orthometric height from a pressure altitude or vice versa, an intermediate conversion to a geopotential height is necessary, as the model below 86 km alti¬ tude is parameterized on geopotential heights. A geopotential height is essentially a gravitational potential

55Consider, for example, some of the models listed on NASA’s model web site at http://ccmc.gsfc.nasa.gov/ modelweb/ .

Coordinate Systems for Target Tracking

51

difference divided by a standard acceleration due to gravity. It can be defined as

„ AWP-Wo

— _

80

1

go

gdZ

(64a)

(64b)

where Wp is the gravity potential of (53) at the target and Wo is the gravity potential on the reference surface, typically on a specific definition of the geoid. The integral form of the geopotential height integrates from an orthometric height of 0 (on the geoid) to the points, which is at an orthometric height zp. The magnitude of acceleration due to gravity g depends on position. The integral is a path integral (there arc multiple ways to get from an orthometric height of 0 to zp), but the actual path taken does not matter, only the endpoints. The dZ is a differential orthometric height. The U.S. Standard Atmosphere takes go = 9.80665m/s2 as the reference acceleration due to gravity.

The standard atmosphere uses a low-precision approximation of the conversion from orthometric heights to geopotential heights. Using Newton’s law of universal gravitation, which one can find in most introduc¬ tory physics textbooks, one can express the magnitude of acceleration due to gravity g at a height of h above a spherical Earth with radius as

GM6 {r& +h)2

(65)

Assuming that go is the acceleration due to gravity when li = 0, multiplying both sides of (65) by the inverse of the right side of Eq. (65), one obtains

which can be solved to get

g

h2 + 2hrs ' GM6

_ gorl

( h + rs )2’

(66a)

(66b)

(67)

which no longer requires one to know GM&. The standard atmosphere takes the radius of the spherical Earth to be r$ = 6,356,766m. Substituting the approximation of acceleration due to gravity of (67), into the integral definition of geopotential height of (64b), one obtains the approximation

H ^ r6Zp ~ r6+zp

(68)

The standard atmosphere uses (68) to convert between orthometric height z,p of a point and geopotential height H. However, some air-pressure sensors that provide height readings do not perform the conversion. Since 86, 000 m orthometric height corresponds to 84, 852 m geopotential height, this can introduce a notable bias at higher altitudes.

52

David Frederic Crouse

Given the geopotential height 77 of a point, the standard air pressure for altitudes 86 km and below, which are relevant for general aviation, in units of pascals (N/m2) is

_ TM,b _

Tm,b + Lmj, (H — H/,)

( H-h exp

where R* = 8.3 1432 Nm/(molK) is a 1962 standard value for the universal gas constant and Mo = 28.9644 x 10~3 kg/mol is the mean molecular weight of the atmosphere at sea level. The values 7.y//, and 74/./; are constants that depend on the geopotential altitude 77 whose values are listed in Table 8 for different values of b. The value of b is chosen such that H\ > ^ 77. In some areas, the land may lie a significant distance below MSL. For example, the land around the Dead Sea between Israel and Jordan is between 700m and 730m below MSL, though it is rising by approximately 4mm per year due to the crust rebounding after water loss [244], The standard says that the atmospheric pressure model is valid down to 5 km below MSL. The pressure Pb is the air pressure at geopotential height 77/,. Given P/,, one can find Pb+\ by using (69) at b with 77 = 74+1. The model is parameterized on Pq and To, the air pressure and temperature at MSL. Local pressure and temperature settings that are above sea level can be extrapolated back down to sea level using the model and knowing the orthometric height of the observer.

gQM0

h

If LM,b A 0 If Lmm = 0

(69)

Table 8: Constants for Air Pressure and Orthometric Altitude Conversion

b	Hb{ m)		TM,b( K)
0	0	-6.5 x 10~3	288.150
1	11000	0	7m, o + Lm, o(77i — Ho) = 216.650
2	20000	1 x 10-3	7m, i = 216.650
3	32000	2.8 x 10“3	7m, 2 + 7/m,2(773 — H2) = 228.650
4	47000	0	7m, 3 + 7m, 3(774 — H3) = 270.650
5	51000	-2.8 x 10-3	7m, 4 = 270.650
6	71000	-2.0 x 10~3	7m, 5 + Lm,5{H6 — Hb) = 214.650
7	84852	0	7m, 6 + 7/m, 6 (H7 — Ho) = 186.946

The constants are for the U.S. Standard Atmosphere 1976 model up to 86km altitude. The numerical values for T^,b have not been rounded nor truncated.

Figure 20 plots altitude as a function of air pressure in the U.S. Standard Atmosphere 1976. The pressure Tb and temperature To at MSL, which define the model, were taken to be the standard pressure and tempera¬ ture used by aircraft over 18,000ft altitude in the United States, a pressure of 101, 325 Pa and a temperature of 288.15 K (15°C) [335].56

56The FAA’s aeronautical information manual lists the standard air pressure as 29.92 inches of mercury, which depending on the chosen temperature and altitude of the sensor, is not precisely 101, 325 Pa. However, 101, 325 Pa is the generally used conversion. The National Institute of Standards and Technology (NIST) lists multiple standard conversions from inches of mercury to pascals on their site at http: //physics. nist.gov/Pubs/SP811/appenB9.html . Inches of mercury should generally not be used to specify pressures as it is ambiguous and archaic.

Coordinate Systems for Target Tracking

53

Fig. 20: The orthometric height as a function of air pressure in the U.S. Standard Atmosphere 1976 using standard temperature and pressure at MSL. The pressure scale is logarithmic.

5.4 Gravitational Models and the Local Vertical

A frequently used local coordinate system in target-tracking systems is the ENU system. The coordinate system can be useful as one can approximate ellipsoidal “up” by the gravitational vertical and North using a gyrocompass, though more accurate results are theoretically possible using stars. As mentioned in Section 8.1.1, given at least two vectors in multiple coordinate systems, one can easily find the rotation matrix between the systems. Section 4.2 discusses the ENU coordinate system in terms of coordinates on the reference ellipsoid. However, if one assumes that the local ellipsoidal vertical coincides with the vertical implied by the reference ellipsoid, one loses precision. This section quantifies how much the non-ellipsoidal nature of the Earth can matter. The ability to relate the gravitational vertical direction to a global coordinate system can play an important role in the design of inertial-navigation devices. For example, in [193], the navigation error of an inertial measurement unit (IMU) operating using various gravitational models is investigated. It is shown that the best models offer horizontal -position accuracies better than 5 m after a period of 53 minutes. A similar analysis performed with respect to moving vehicles on land is in [191], where it is noted that the navigation errors can be reduced if the vehicle periodically stops and recalibrates its “motion-free” gravitational acceleration.

Using the EGM2008 zero-tide spherical-harmonic gravitational model and the method of Appendix F for utilizing such models, Fig. 21 plots the bias in the magnitude and direction of the gravity-acceleration vector as obtained using the WGS 84 reference ellipsoid versus using the zero-tide EGM2008 model. An altitude of 9km is used as it is above all terrain albeit still in the atmosphere. The angular distance between two vectors Vi and V2, determining the direction of the vertical in both models, was determined using the cross-product relation

6 = arctan(||vi x V2 1| , Vi • V2)

(70)

54

David Frederic Crouse

rather than the simpler dot product relation

a ( Vl'V2 ^

6 = arccos - — rpr; — - (71)

VllvilllNiy

as it was observed (most notably in Section 6 when considering the effects of tides) that the cross-product formulation is more robust to numerical precision limitations.

Gravitational Acceleration Difference from the Ellipsoidal Model The Deflection of the Vertical from the Ellipsoidal Model

Longitude Longitude

(a) Magnitude Bias (b) Direction Bias

Fig. 21: The difference in the magnitude (a) and the angle (b) between the gravity-acceleration vector g under the EGM2008 zero-tide model and the WGS 84 ellipsoidal model using the parameters of Table 5 and Eq. (39) for the acceleration vector under an ellipsoidal gravitational model for points 9 km above the reference ellipsoid. One milliGal (mGal) is equal to 1 x 10 5 m/s2. The points were evaluated on a 2500 x 2500 grid of points in latitude and longitude. As noted in [28, Ch. 3.8], the gravitational- acceleration differences of (a) will tend to bias the prediction of a low-Earth orbiting satellite by approximately 3.5km after one day.

The directional bias of the vertical direction according to the EGM2008 model in Fig. 21 is observed to exceed 50 arcseconds. This is line with with [87], which notes that the deflection of the vertical can be as high as 70 arcseconds. While probably insignificant when performing radar tracking of aircraft, such an angular bias can be significant when frying to perform high-precision optical tracking of space debris. The acceleration biases compared to using the ellipsoidal model in Fig. 21 appear to be quite small. However, they are many orders of magnitude larger than the highest reported 3D accelerometer accuracy, which is on the order of 9.8 x 10_13m/s2 [279]. 57 Also, as quantified in [28, Ch. 3.8], ignoring higher order terms of the gravitation model can lead to errors in the direction of motion of a low-Earth orbital object on the order of 3.5 km after one day, which is likely to be significant when tracking orbital debris. Both [28, Ch.

57The citation [279] refers only to the numbers reported in a very long abstract. The actual paper as well as the complete conference proceedings appear to be unobtainable.

Coordinate Systems for Target Tracking

55

3.8] and [234, Ch. 3.1] quantify the magnitudes of the effects of various forces on satellites. In [234, Ch. 3.1], it is noted that a constant radial gravitational-acceleration bias of 1 0 8 m/s2 changes the semi-major axis of the orbital ellipse of a geostationary satellite by approximately 1 m. Note, however, that acceleration contribution due to the Sun and Moon are on the order of 1 0 N m/s2 to 1 0 9 m/s2, as noted in [234, Ch. 3. 1] and verified for the Moon in Section 6.3.

If one wishes to use measurements of the gravitational vertical for orientation estimation, as an input in a dynamic model in a target-tracking algorithm, or in an inertial-navigation unit, it is generally desirable to have covariance matrices expressing the accuracy of the model. The EGM2008 model comes with standard- deviation terms for all of the coefficients in the model. Using the method of Appendix F.3, one can find a covariance matrix for the gravitational-acceleration vector implied by the EGM2008 model at any point. Figure 22 plots the average bias in the magnitude of acceleration due to gravity and the bias in the direction of the vertical implied by the standard deviations of the spherical-harmonic coefficients in the EGM2008 model. In [134], the vertical deflection error of the EGM2008 model is analyzed in the mountains of Bavaria and Switzerland. The RMS empirical accuracy of the model is observed to be as high as 4.19 arcseconds, which is larger than the magnitude of the error predicted in Fig. 22, though one would expect errors in the mountains to be higher than elsewhere. In [159], the magnitudes of various approximations used in the deflection of the vertical calculations, such as truncating the spherical-harmonic series or using different permanent tide system, arc quantified.

Mean Error in the Gravitational Acceleration Magnitude

-150 -100 -50 0 50 100 150

Longitude

(a) Acceleration Bias

Mean Error in the Vertical Direction

14

12

10

8

6

4

2

Ld

£

-150 -100 -50 0 50 100 150

Longitude

(b) Direction Bias

Fig. 22: The average bias of the gravity-acceleration vector (a) and the direction of the vertical (b) as implied by the standard- deviation terms that come with the zero-tide EGM2008 gravitational model evaluated on a 256 x 256 grid in latitude and longitude at 9km altitude. The method of Appendix F.3 was used to obtain covariance matrices for the acceleration vector due to gravity at each point. Then, Monte Carlo simulation of the errors assuming additive Gaussian noise was performed to determine the expected vertical bias. At 60° latitude (the latitude of southern Greenland), the acceleration bias is on the order of 17mGal and the vertical orientation bias is on the order of 2 arcseconds. The errors are larger than the commission errors plotted in [252], The results become worse near the poles. The plots only go up to +85° latitude, due to numerical issues discussed in Appendix F.3. The issues can be avoided by using quadruple precision instead of double-precision floating-point numbers.

56

David Frederic Crouse

The average biases in Fig. 22 were determined through Monte Carlo simulation. Given the covari¬ ance matrix of gravitational acceleration £yy, as well as the nominal gravity-acceleration vector based on the EGM2008 model g at a particular point, noisy samples of the gravity-acceleration vector were drawn assuming that they were normally distributed. That is,

i

Ssamp = § T Eyy W (72)

where w is a 3 x 1 standard ./F{0.1} random vector and the matrix square root is meant to be a lower- triangular Cholesky decomposition, which is chol (Sigma, ' lower' ) in Matlab and which can be implemented using the algorithm of [111, Ch. 4.2]. 58 The noise could be added directly to the gravity- acceleration g rather than the gravitational acceleration, since they differ by only an additive constant. The deflection of the vertical was computed using (71) with gsamp and g as vi and V2 for each random sample and the results were then averaged. The bias in gravity-acceleration was computed by averaging the magnitude

| |S — gsamp || •

The gravitational model considered in this section is time-invariant. Section 6 discusses time-invariant and time-variant tide models and quantifies the magnitude of the error in the vertical direction estimate due to time-varying tides. It is demonstrated that the vertical direction error due to time- variant tides is far less than the precision of the vertical in the EGM2008 model, as plotted in Fig. 22.

58 Strictly speaking, as noted in [111, Ch. 4.2.4], the Cholesky decomposition is not a matrix square root, though it is often called

one. Note that if £yy is the Cholesky decomposition of £yy , lhen Tyy = Tyy I Eyy I and usually £yy # EyyEyy .

Coordinate Systems for Target Tracking

57

6. THE EFFECTS OF TIDES ON COORDINATE SYSTEMS 6.1 Why Tides Can Matter

The oceans and crust of the Earth arc constantly in motion. Many arc familiar with the diurnal to semidiurnal ebb and flow of ocean tides, which can be quite large. For example, the largest observed mean range of tides is at Burntcoat Head, Minas Basin, Bay of Fundy, Nova Scotia, where the difference between mean high-tide and mean low-tide water levels is 38.4ft (% 11.7m) [42]. However, the crust of the Earth moves due to the effects of tides as well. While tides might initially appear to be rather irrelevant to those designing target-tracking algorithms, the motion of the ground due to tides can be larger than the precision of sensors for laser tracking of satellites. For example, laser detection and tracking of orbital debris has been demonstrated [3, 182], producing monostatic range detections of 0.7 m RMS accuracies for targets less than 0.3m2 [183], The IERS conventions have set an accuracy goal for satellite laser ranging to a millimeter- level accuracy [257, Ch. 9.1.3]. In contrast, the total variation in the position of the ground over time due to solid-Earth tides can be on the order of 30cm. This is shown in Fig. 23, where the offset of the ground due to solid-Earth tides using the models described in the caption is plotted. Consequently, one cannot simply localize a receiver in an ECEF coordinate system, such as the ITRF, once and never update the position. Rather, the location of the receiver must be updated constantly during the day.

Solid Tidal Offset without Ocean/Atmospheric Loading

-150 -100 -50 0 50 100 150

Longitude

Fig. 23: The magnitude of the Cartesian offset of the ground due to solid-Earth tides on 12 June 2012 at 14:00 UTC without the permanent tide, plotted on a 1000 x 1000 grid in latitude and longitude, using the model described in the IERS conventions [257], not including ocean and atmospheric loading effects. The IERS provides a reference implementation in Fortran called dehanttideinel, which can be obtained from a link in Table 2. As the direction changes over time, the total displacement over time is actually on the order of 30cm. The algorithms in the included subroutines step2diu and step21on, which provide tidal corrections due to mantle anelasticity in the diurnal and long-period bands respectively, are not completely documented in the IERS conventions. The necessary Moon and Sun direction vectors were obtained using the DE430 ephemerides, which were read using the cspice_spkezr function in NASA’s SPICE toolkit for Matlab (see Appendix I). Time-conversion routines from the IAU’s SOFA library were also used.

58

David Frederic Crouse

Due to the periodic tidal motion of the ground, when using a GNSS receiver to localize a sensor for tracking, one must know the specific algorithms for tidal compensation, if any. For example, the data pro¬ cessing routines in the European Space Agency’s book on GNSS data processing contain models for the motion of ground due to tides [320, Ch. 5.7], Solid-Earth tide estimation routines arc included in the soft¬ ware accompanying the book,59 allowing one to provide a user their location on the Earth in coordinates that (ideally) do not drift over time due to tidal effects. Similarly, when one examines the source code of the open-source GPS Toolkit, which is described in [120], 60 one sees that code for solid-Earth tide corrections is also present. However, when performing high-precision tracking of orbital objects, the actual time-varying ECEF location of the receiver is generally desired, in which case, one will not want to use a GPS receiver that provides coordinates that do not change over time due to tides. Minute variations in the displacement of the ground over moderate distances have proven significant in numerous other high-precision applications. For example, the motion of the ground due to tides matters when considering the calibration of large par¬ ticle accelerators. In such an application, relative displacements less than the size of a human hair can be significant [287].

When considering target tracking, tides affect more than just the location of a sensor in a global co¬ ordinate system. One must understand the basics of tidal models to properly use high-precision gravita¬ tional models, which arc provided in terms of varying permanent tide systems. Knowledge-based target identification can utilize a wide variety of parameters about the environment in which a target is being de¬ tected [109, Ch. 9]. For example, a ship detected in a region where the current water level is low must have a shallow draft or have run aground. If the ship is moving in the shallow region, it cannot be any type of ship with a deep draft.

The development of numerical tide models, as described in [171], is extremely complicated. Figure 24 provides an overview of many of the international organizations that play various roles in the development of tidal models. One designing target-tracking algorithms does not necessarily need to know how such models are developed, but must understand the terminology surrounding tide models and should know how significant tidal effects are. The following sections provide an overview of the terminology and models surrounding the effects of tides when considering the motion of the ground and sea as well as when using gravitational models.

6.2 Tides in Terms of Ground and Ocean Motion

Table 9 lists a number of terms that arise when discussing the motion of the ground, sea, and atmosphere due to tidal forces. The periodic motion of the ground due to tides is not the same thing as an earthquake. Whereas the magnitude 9.0 earthquake off the Pacific coast of Tohoku, Japan in 2011, which lasted ap¬ proximately 180s, had peak ground acceleration up to 3 times acceleration due to gravity [101, Fig. 3], solid-Earth tides arc the slow daily motions of the ground due to external tidal forces. The total motion of the sea surface and/or ground due to time-variant tides is often decomposed into many parts. Whereas Fig. 23 in the previous section showed the motion of the ground due to solid-Earth tides, the total motion of the ground due to tides involves a number of other tidal components. Table 0.1 in the IERS conventions quan¬ tifies the magnitudes of the effects of many tidal components beyond solid-Earth tides, with ocean loading being on the order of centimeters, atmospheric loading millimeters, the ocean pole tide on the order of 2 cm

59The book as well as the accompanying software CD can be downloaded from the European Space Agency’s Navipedia web site at http : / / www.navipedia.net/ index.php/ GNSS : Tools .

60The source code for the GPS Toolkit can be obtained on the web site http : / 7www.gpstk.org/ .

Coordinate Systems for Target Tracking

59

Fig. 24: A subset of the tangled web of international organizations that play a role in monitoring and modeling tides. This includes both gravitational models and models for the physical motion of the Earth. No one organization can be viewed as a definitive source of tidal information. Rather one must often draw data and models from multiple sources. Numerous other organizations, such as the International Hydrographic Organization (IHO), play various additional roles in mapping the sea surface and floor and in monitoring tides.

radially and in millimeters tangentially, and the ocean-pole -tide loading on the order of 2mm radially and less than 1mm tangentially [257, Ch. 0.2]. Additionally, sea levels can be affected by wind as well as the inverse barometric effect, which as discussed in [116,217]. The inverse barometric effect explains the difficult-to-model phenomenon that local sea levels rise when the local barometric pressure drops.

Global models for the displacement of points for ocean loading are very complex and include Some Pro¬ grams for Ocean-Tide Loading (SPOTL) available at http : //igppweb. ucsd.edu/ -agnew/Spotl/ spotlmain.html and the programs olfg and olgg available at the links on the web site http: // holt.oso.chalmers.se/loading/faq.html . Various tidal-loading models are compared in [254]. The IERS conventions describe models for pole tides and atmospheric loading [257, Ch. 7], Many tidal models, both for the motion of the ground as well as for changes in Earth geopotential models, decompose tidal components into sets of periodic values with periods that are proportional to the duration of lunar and solar effects [281]. The coefficients for the periodic terms are generally called Doodson parameters, which stem from Doodson’s harmonic development of the tidal generating potential [72].

60

David Frederic Crouse

Table 9: Commonly Used Terminology Relating to Time- Variant Tides that Move the Sea Surface and Ground

Term	Definition
Elastic Ocean Tide/ Geocentric Ocean Tide	The sum of the ocean tide and the radial load tides.
Ocean Tide	The vertical displacement of the ocean above the sea floor (which moves) due to time- variant tides.
Ocean Loading	The vertical displacement of the ground due to the weight of the ocean moving due to time- variant tides.
Atmospheric Loading	The vertical displacement of the ground due to the weight of the atmosphere moving due to time-variant tides.
Solid-Earth Tide	The motion of the ground due to tidal forces, not including the effects of ocean and atmospheric loading.
Radial-Load Tides	The combined effect of ocean and atmospheric loading as well as the solid-Earth tides.
Radiational Tide	The motion of the crust, oceans, and atmosphere due to heating from the Sun.
Pole Tide/ Polar Motion	A change in the direction of the Earth’s figure axis (the figure axis is defined in Appendix E).
Solid Earth Pole Tide	Motion of the ground due to the centrifugal effects of changes in the Earth’s figure axis.
Ocean-Pole-Tide Loading	An offset in the location of points due to the motion of the oceans resulting from the pole tide.
Dynamic Topography/ Dynamic Sea-Surface Topography/ Absolute Dynamic Topography/ Sea-Surface Topography/ Ocean-Surface Topography/	The true height of the sea surface with respect to the geoid; that is, the orthometric sea-surface height. It is the sum of the elastic ocean tide and the mean dynamic topography (and possibly the radiational tide). Sometimes, the terms are used to (ambiguously) mean the mean dynamic topography.
Mean Dynamic Topography/ Quasi-Stationary Sea Surface Topography (QSST)	The (nearly) time-invariant part of the dynamic topography. Does not include seasonal variations, non-permanent winds, eddies. Does include permanent currents, circulation linked to the Earth’s rotation, permanent winds.
Sea-Surface Height (SSH)	The height of the sea surface with respect to the reference ellipsoid. SSH = MSS + SLA
Mean Sea Surface (MSS)	The ellipsoidal height of the mean dynamic topography. Approximately, the geoid height plus the mean dynamic topography
Sea-Level Anomaly (SLA)/ Sea-Surface Height Anomaly (SSHA)	The height of the sea surface above the MSS.
Mean Sea-Level Anomaly	A mean sea level anomaly (sometimes just referred to as an SLA, lending confusion to the previous definition of the SLA), as defined by the NOAA [44], occurs when the 5-month running average of the interannual variation is at least 0.1 m greater than or less than the long-term trend.
Interannual Variation	The monthly MSL after the trend (increasing/decreasing over time) and the average seasonal cycle are removed.
Sea-State Bias/ Electromagnetic Bias	A bias in altimeter ranging to the sea surface due to wave troughs reflecting more energy than wave crests.
Tidal Aliasing	An aliasing problem in tidal estimation when the revisit period of an altimeter satellite is too low to observe high-frequency periodic tidal components.
Inverse Barometric Effect	The observation that tides are higher on the order of 1 cm when atmospheric pressure goes down 1 mbar.

The definitions have been inferred from common usage in many papers.

The tidal models discussed in the IERS conventions are specified in terms of a tide-free permanent tide model. Permanent tides are time-invariant theoretical tidal deformities of the Earth that arise due to theoretical modeling of how the Earth deforms in the presence of other bodies in the solar system (permanent tides are discussed in the following section). As noted in [320, Ch. 5.7.1], the tide model under which the ground displacements are taken can be changed by adding constant latitude-dependent offsets to the time- variant tidal offsets. Such offsets change the locations of points on the Earth if one attempts to take raw GNSS positions and reduce them to time-invariant locations, which is desirable when making or using maps, but which is not desirable when localizing a sensor for high-precision tracking. As sensors can often use filtered real-time GNSS signals to monitor their positions in a global ECEF coordinate system without directly using tidal models, the modeling of the motion of the solid ground due to tidal effects is not further considered in this report.

The geoid is often synonymous with the concept of MSL: a surface of equal gravitational potential to which the Earth’s oceans would rest in the absence of time-variant tides, under a certain permanent tide definition. However, the true definition of MSL and the ocean sea surface is somewhat more complicated

Coordinate Systems for Target Tracking

61

due to semi-permanent currents in the oceans. Just as swirling water around in a cup makes the water at the sides of the cup higher than the water in the center of the cup, the motion of water in the oceans causes some parts of the ocean to be constantly higher than others after compensating for tidal effects. The quasi¬ stationary sea-surface topography (QSST) is the mean distance between the actual surface of the sea and the geoid after time-variant tides have been removed. As plotted in [196], the maximum difference between the high and low values of the QSST worldwide is approximately 2m. In [252], it is plotted as approximately 3 nr. To put this in context, the World Meteorological Organization defines sea-state codes on a scale from 1 to 9, where waves having a height of 2 nr fall into sea states 4 and 5, which is moderate to rough seas [362, Pg. 326]. On the Beaufort scale of wind, such waves are typical when the wind 10m above the ground is about 8m/s to 10.7m/s. Consequently, without detailed, high-resolution, real-time oceanic wave measurements, it is unlikely that QSST will matter when predicting ahead trajectories of ships at sea. Indeed, in extreme conditions, the difference between the geoid and the mean sea surface can be very large. For example, the highest reported wave height appeal's to be 112ft 34.14m) reported by the USS Rampao on 7 February

1933 during winds of 60 kt (% 31 m/s) [360].61

6.3 Tides in Terms of the Gravitational Field

Table 10 defines a number of terms that are important for describing the effects of tides when discussing the effects of tides with respect to gravitational models. The gravitational potential models discussed in Sections 5 and 5.4 are defined in terms of various permanent tide systems. Permanent tides are time-invariant distortions of the Earth's shape and gravitational potential resulting from the presence of other bodies in the solar system. Gravitational potential models, and thus geoids, are generally defined in terms of one of three permanent tide models, tide-free, zero-tide, and mean-tide systems, onto which time-variant components can be added. This section explains the difference between the tidal models, and demonstrates that the direct tidal contribution as well as all time-variant tidal contributions can usually be ignored in algorithms for determining the vertical direction of an observer on the ground, as the effects arc an order of magnitude smaller than the precision of the zero-tide EGM2008 model.

In a tide-free system, the gravitational potential V of a point in an ECI coordinate system is computed based on the shape that the Earth would be if the Earth were the only object in space (the external potential [direct tidal potential] has been removed) and all time-variant and time-invariant tidal effects were removed. The tide-free model can only be computed based upon models of the elasticity of the Earth; it cannot be directly observed. The zero-tide system computes the gravitational potential of an observer in an ECI coordinate system based on the mass distribution of the true shape of the Earth after time-variant components and the external potential have been removed. Consequently, the zero-tide model is the best choice to use when one wishes to consider the gravitational effects of the Earth independent of other bodies in the solar system without regal'd to time-variant tides.

The mean tidal model is the zero-tide model with the direct tidal potential added. The direct tidal potential is the apparent change in potential experienced by an observer in an ECI or ECEF coordinate system due to the physical presence of bodies other than the Earth. To understand this, consider a coordinate system whose origin is at the Newtonian center of mass of a system containing N bodies. In an arbitrary

61 That is the highest reported wave not caused by a seismic disturbance. The highest tsunami wave appears to have been a 1720ft (ss 524m) wave in Lituya Bay on the northeast shore of the Gulf on Alaska on 9 July 1958 caused by a rockslide into the bay resulting from an earthquake [228].

62

David Frederic Crouse

Table 10: Definitions of the Different Tidal Potentials, Tide Systems, and Tides

Term	Definition
Non-Tidal Potential	The potential of the Earth if no other bodies existed in space. Changes over time are due to Earth mechanics, such as volcanism.
Indirect Tidal Potential	The potential due to the deformation of the Earth’s crust by other celestial bodies (Sun, Moon, etc). There are time-invariant as well as time-variant parts.
Direct/External Tidal Potential	The time-variant and invariant gravitational potential due entirely to other celestial bodies, almost always computed with respect to an observer accelerating with the Earth through space.
Centrifugal Potential	The apparent change in gravitational potential experienced by a stationary observer on a rotating Earth compared to an observer on a theoretical non-rotating Earth.
Tide-Free/Non-Tidal System	A system influenced only by the non-tidal potential.
Zero-Tide System	A system influenced only by the non-tidal potential and the permanent indirect-tidal potential.
Mean-Tide System	A system influenced only by the non-tidal potential and the permanent direct-and-indirect-tidal potentials.
Permanent Tide	Contributions to the tidal potential (or the shape of the Earth) that are time-invariant.
Solid-Earth Tide	Indirect tidal-potential contributions due to the deformation of the non-oceanic (or any water) and atmospheric parts of the Earth. Sometimes, direct tidal potentials are included.
Ocean Tide	The indirect tidal-potential contribution due to the movement of water in the oceans.
Pole Tide/ Polar Motion	Changes in the orientation of mass within the Earth leading to a change in the Earth’s figure axis, which changes acceleration due to gravity with respect to fixed points on the crust.
Solid Earth Pole Tide	A time-variant solid-Earth tidal potential arising from the motion of the Earth’s figure axis with respect to the ITRS (polar motion), presumably due to motions in the Earth’s core and mantle with respect to the crust.
Ocean Pole Tide	A time-variant tidal potential due to the motion of the ocean as a result of the centrifugal effect of polar motion.

The definitions are based on information in [257,298,299].

coordinate system, the location of the center of mass rcenter is

Tcenter

(73)

where Mj is the mass of the ith body whose center of mass is located at r j. A coordinate system defined such that the center of mass is the origin has rcenter = 0 over all time. Consequently, the location of the first object in the center-of-mass coordinate system can be expressed as

<74>

Now, consider moving the origin of the coordinate system to the center of mass of the first object. A point in the center-of mass coordinate system p can be expressed in the coordinate system at the center of mass of the first object p by the affine transformation

P = ri+p.

(75)

Using (74), taking the second derivative of (75) with respect to time, and solving for p, one finds that an acceleration can be transformed from the center-of-mass system to the center of mass of Object l’s system using

(76)

Coordinate Systems for Target Tracking

63

The acceleration due to the direct tidal potential (both permanent and time-variant) on Object 1 is p minus the contribution to gravitational acceleration due to Object 1 itself. The direct tidal potential is thus the potential whose gradient (or negative gradient, depending on the definition) gives the acceleration due to the direct tidal potential.

The significance of the transformation from an inertial center-of-mass coordinate system and its relation to tidal potentials can be easily understood when considering the ideal Earth-Moon system on its own, as illustrated in Fig. 25, where the Earth is approximated as a sphere having radius r* and the center of mass of the Moon is a distance of Ra from the center of mass of the Earth. The example here is inspired by a similar example in [171, Ch. 6.1]. Points A, B, and C are all on the surface of the spherical Earth. Using the law of cosines and the law of sines, one can find that

rbm =Rl + rl -Rdf 6 cos(e)

a =arcsin f — — sin(0) ) .

\rbm )

The force per unit mass (acceleration) experienced by a small object at point p on the surface of the spherical

(77)

(78)

y

Fig. 25: The geometry of the Earth and the Moon used in the discussion on the direct tidal potential. When the gravitational force of the Moon on the Earth is transformed into an ECI coordinate system, the transformed forces point roughly towards the center of the Earth at the poles (along the v-axis) and away from the center of the Earth at the equator (along the x-axis) resulting in an equatorial bulge and explaining why there are two high tides per day on most of the Earth. The gravitational potential associated with the force of the Moon and all other solar system bodies in ECI coordinates is the direct potential. The time-invariant part is the permanent direct potential. The relative radii of the Earth and Moon are to scale, but the Earth-Moon distance is not to scale.

Earth, but not rotating with the Earth, is

ap=W6 + Wc (79a)

(_ uxcos(0) — uvsin(0)) + (uxcos(a) — uvsin(a)) (79b)

r6 Rbm

where simple uniform-mass-distribution approximations were used to model the gravitational effects of the Earth and Moon. GM& and GMc are the universal gravitational constant times the mass of the Earth and the

64

David Frederic Crouse

mass of the Moon, respectively, r$ is the radius of the Earth, and uA and uv arc unit vectors in the directions of the x and y axes, respectively. Noting that the acceleration of the Moon due to the Earth is

ac = -

GM6

Rl

Uv

(80)

where R& is the distance between the centers of mass of the Earth and Moon, the acceleration of the point transformed into a coordinate system centered at the center of mass of the Earth using (76) is

Up Up o (81)

The direct tidal contribution to the acceleration in (81) is all accelerations resulting from forces that arc not due directly to the physical mass of the Earth. That means that one must subtract off VV$ to get the direct tidal acceleration contribution

~ direct

ap

GM& ■ t w GMd

(uvcos(a) — uvsm(a)) - — ^-ux.

kbm K<z

(82)

In the examples below, the following values arc used:

GM6 =3.986004418 x 1014 -= (83)

s2

H =0.0123000371 (84)

GMd =hGM6 (85)

r6 =6378136. 6m (86)

Rd =3.994833385218484 x 108m. (87)

All of the values except for that of R arc taken from the IERS conventions [257, Ch. 1], The value of R<t was determined using the DE430 ephemerides read using the cspice.spkezr function in NASA’s SPICE toolkit for Matlab for the date of 12 June 2012 at 14:00 UTC (see Appendix I). At point A in Fig. 25, when 0 = 0, Eq. (82) simplifies to

<reC1e=o GMc((^_r6)2

» 1.005 x 10_6^ux.

s2

Similarly, at point C in Fig. 25, when 6 = n, (82) simplifies to

~ direct

ap

0 = TI

-GMt ((ATPp- 4)

9.580 x 10~7^ux.

s2

Ux

(88a)

(88b)

(89a)

(89b)

Coordinate Systems for Target Tracking

65

On the other hand, when 6 = |, then

£ direct I *P lfl=§

GMd

Rl+n

r 6

4

GMtr

R2

Ur

- 1.174 x 10“ V- 4.903 x 10“ 'uv

(90a)

(90b)

From (88b), (89b), and (90b), it can be seen that the total direct tidal potential of the Moon produces a force per unit mass on the order of 10 6 m/s2 around the globe. As the Moon orbits the Earth, the directions of the accelerations will change. However, since the Moon stays in roughly the same plane, the acceleration of a particle near the poles contains a time-invariant portion in the y direction, thus compressing the Earth. This time-invariant contribution is the direct permanent tidal potential.

The direct permanent tidal potential due to the Moon is on the order of 10 6 m/s2. In comparison, the bias of acceleration due to gravity in the tide-free EGM2008 model in Fig. 22 is on the order of 10-5m/s2. As shown in a plot in [234, Ch. 3.1], the acceleration experienced in an ECI coordinate system due to the Sun is less than that due to the Moon. Consequently, when using a gravitational-acceleration measurement to determine the orientation of an observer on or near the surface of the Earth, one can probably ignore the direct tidal potential. Consequently, the zero-tide permanent-tide model is often sufficient.