The difference between an ellipsoidal and a spherical spectrum of gravimetric data lies with the different surfaces on which the data reside. Therefore, all “spherical” harmonic coefficient sets determined from terrestrial data should more correctly be called ellipsoidal spectra. A spherical spectrum, on the other hand, is obtained from data provided, for example, by a satellite in a circular polar orbit. This paper derives the exact forward and reverse transformations from ellipsoidal to spherical spectra with practical computations in mind.

An estimator was recently developed by R. Rummel that captures most advantages of least-squares collocation and the integral model approach without many of their disadvantages. It was developed as an alternative to the inverse Stokes integral for processing large, densely spaced satellite altimetry data sets by accounting also for measurement errors. This idea is extended here to heterogeneous gravity gradient data sets at altitude. The appropriate kernel function is expressed as a Fourier transform of its spectrum which in the flat-earth approximation is fitted by pieces of powers of spatial frequency and thereby can be integrated analytically. The estimator is applied to the computation of gravity disturbance differences referenced to a least-squares collocation surface defined by a set of tie-point gravity disturbances.

The Brillouin sphere is defined as the smallest sphere, centered at the origin of the geocentric coordinate system, that incorporates all the condensed matter composing the planet. The Brillouin sphere touches the Earth at a single point, and the radial line that begins at the origin and passes through that point is called the singular radial line. For about 60 years there has been a persistent anxiety about whether or not a spherical harmonic (SH) expansion of the external gravitational potential, V, will converge beneath the Brillouin sphere. Recently, it was proven that the probability of such convergence is zero. One of these proofs provided an asymptotic relation, called Costin's formula, for the upper bound, E-N, on the absolute value of the prediction error, e(N), of a SH series model, V-N(theta,lambda,r), truncated at some maximum degree, N = n(max). When the SH series is restricted to (or projected onto) a particular radial line, it reduces to a Taylor series (TS) in 1/r. Costin's formula is E-N similar or equal to BN-b(R/r)(N), where R is the radius of the Brillouin sphere. This formula depends on two positive parameters: b, which controls the decay of error amplitude as a function of N when r is fixed, and a scale factor B. We show here that Costin's formula derives from a similar asymptotic relation for the upper bound, A(n) on the absolute value of the TS coefficients, a(n), for the same radial line. This formula, A(n) similar or equal to Kn(-k), depends on degree, n, and two positive parameters, k and K, that are analogous to b and B. We use synthetic planets, for which we can compute the potential, V, and also the radial component of gravitational acceleration, g(r) = partial derivative V/partial derivative r, to hundreds of significant digits, to validate both of these asymptotic formulas. Let superscript V refer to asymptotic parameters associated with the coefficients and prediction errors for gravitational potential, and superscript g to the coefficients and predictions errors associated with g(r). For polyhedral planets of uniform density we show that b(V) = k(V) = 7/2 and b(g) = k(g) = 5/2 almost everywhere. We show that the frequency of oscillation (around zero) of the TS coefficients and the series prediction errors, for a given radial line, is controlled by the geocentric angle, alpha, between that radial line and the singular radial line. We also derive useful identities connecting K-V, B-V, K-g, and B-g. These identities are expressed in terms of quotients of the various scale factors. The only other quantities involved in these identities are alpha and R. The phenomenology of 'series divergence' and prediction error (when r < R) can be described as a function of the truncation degree, N, or the depth, d, beneath the Brillouin sphere. For a fixed r <= R, as N increases from very low values, the upper error bound E-N shrinks until it reaches its minimum (best) value when N reaches some particular or optimum value, N-opt. When N > N-opt, prediction error grows as N continues to increase. Eventually, when N >> N-opt, prediction errors increase exponentially with rising N. If we fix the value of N and allow R/r to vary, then we find that prediction error in free space beneath the Brillouin sphere increases exponentially with depth, d, beneath the Brillouin sphere. Because b(g) = b(V) - 1 everywhere, divergence driven prediction error intensifies more rapidly for g(r) than for V, both in terms of its dependence on N and d. If we fix both N and d, and focus on the 'lateral' variations in prediction error, we observe that divergence and prediction error tend to increase (as does B) as we approach high-amplitude topography.

This note presents formulas to evaluate a spherical harmonic model of Earth’s gravitational potential for essential gravimetric quantities without spherical and linear approximation. Typically, 10–13 significant digits of numerical accuracy for such computations are obtained over the globe using EGM2008 with FORTRAN 77 code that is also provided.

We address a claim frequently raised by one author that determination of Earth's gravity field by satellite observations using traditional methods is mathematically flawed. Specifically, attention is drawn to the practice of setting to zero the initial sensitivity matrix in the variational equations for force model parameters. It is asserted that this would lead to mathematical contradictions. In this paper we establish necessary and sufficient conditions for the initial sensitivity matrix to be zero-conditions that are well founded and accepted worldwide in classical satellite-based determinations of the gravity field. To claim otherwise is shown to be without basis. We inspect a proposed counterexample, and find it, too, requires zero initialization of the sensitivity matrix. In addition, we review proofs and derivations from a classic textbook that also confirm zero initialization. In a numerical exercise, perfect, synthetic data for the central force problem are processed with standard procedures, and results confirm the validity of zero initialization of the sensitivity matrix.

This book covers all aspects of inertial navigation systems (INS), including the sensor technology and the estimation of instrument errors, as well as their integration with the Global Positioning System (GPS) for geodetic applications. Complete mathematical derivations are given. Both stabilized and strapdown mechanizations are treated in detail. Derived algorithms to process sensor data and a comprehensive explanation of the error dynamics provide not only an analytical understanding but also a practical implementation of the concepts. A self-contained description of GPS, with emphasis on kinematic applications, is one of the highlights in this book. The text is of interest to geodesists, including surveyors, mappers, and photogrammetrists; to engineers in aviation, navigation, guidance, transportation, and robotics; and to scientists involved in aerogeophysics and remote sensing.

When thinking of gravity in geodesy and geophysics, one usually thinks of its magnitude, often referred to a reference field, the normal gravity. It is, after all, the free-air gravity anomaly that plays the significant role in terrestrial data bases that lead to Earth Gravitational Models (such as EGM96 or EGM2008) for a multitude of geodetic and geophysical applications. It is the Bouguer anomaly that geologists and exploration geophysicists use to infer deep crustal density anomalies. Yet, it was also Pierre Bouguer (1698-1758) who, using the measured direction of gravity, was the first to endeavor a determination of Earth’s mean density (to “weigh the Earth”), that is, by observing the deflection of the vertical due to Mount Chimborazo in Ecuador. Bouguer’s results, moreover, sowed initial seeds for the theories of isostasy. With these auspicious beginnings, the deflection of the vertical has been an important, if not illustrious, player in geodetic history that continues to the present day. Neglecting the vertical deflection in fundamental surveying campaigns in the mid to late 18th century (e.g., Lacaille in South Africa and Méchain and Delambre in France) led to errors in the perceived shape of the Earth, as well as its scale that influenced the definition of the length of a meter. The vertical deflection, though generally excluded from modern EGM developments, nevertheless forms a valuable resource to validate such models. It is also the vertical deflection that is indispensable for precision autonomous navigation (i.e., without external aids such as GPS) using inertial measurement units. It is the deflection of the vertical that, measured solely along horizontal lines, would readily provide geoid undulation profiles, essential for the modernization of height systems (i.e., vertical geodetic control) without the laborious and traditional methods of spirit leveling. But, measuring the deflection of the vertical is itself an arduous undertaking and this has essentially contributed to its neglect and/or underusage. Even Vening-Meinesz’s formulas of convolution with gravity anomalies do not greatly facilitate its determination. This presentation offers a review of the many roles the vertical deflection has, or could have, played over the centuries, how it has been measured or computed, and how gravity gradiometry might eventually awaken its full potential.

SUMMARY The geological setting of southwestern Oklahoma and northeastern Texas is an ideal example of an aulacogen, the result of the tectonic evolution of a failed rift of the North American continent during the Palaeozoic era (540–360 Ma). The Wichita Province forms the uplifted basement portion of this Southern Oklahoma Aulacogen (SOA). The major fault zones to its north and south are clearly evident in gravity gradient maps produced by the recently constructed Earth Gravitational Model 2008 (EGM2008). Fault parameters, such as the dip angle, location and density contrasts have been estimated from profiles of seismic data and local gravimetry in the 1990s. On the other hand, gravitational gradients that are derived from EGM2008 and then combined to form the differential field curvature are particularly indicative of linear structures such as dip-slip faults. They are used here exclusively, that is, without additional geophysical constraints, in an optimal, least-squares estimation based on the Monte Carlo technique of simulated annealing to determine dip angle and location parameters of the major faults that border the Wichita Uplift region. Results show that these faults have small dip angles, in basic agreement with the low-angle faults inferred from seismic studies. The EGM2008 gradients also appear in some cases to provide an improved map of the major faults in the region, thus offering a strong constraint on their location.

In 2008 P. Xu (Celest Mech Dyn Astron, 100:231–249) proposed a strictly kinematic perturbation method for determining the Earth’s gravitational field from continuous satellite tracking. The main idea is to process orbital arcs of arbitrary length, thus minimizing superfluous parameter estimation associated with stitching together short-arc solutions, and at the same time formulating the problem in terms of standard linear parameter estimation. While the original formulation appears mathematically robust, its nested quadruple along-track integrations are computationally challenging. We reduce the formulation to double integrals and show that the method is numerically not feasible as originally envisaged. On the other hand, by abandoning the rigorous Gauss-Markov formalism, we show the numerical feasibility of processing multiple-day orbital arcs. The methodology lends itself to high-low and low-low satellite-to-satellite tracking, or combinations thereof, as for GRACE-like systems.

Inferring seafloor topography from gravity anomaly currently is the dominant method to obtain a global view of the oceans. Standard techniques rely on an approximate, linear relationship between topography and gravity, which is valid only if the local topography is smooth compared with the regional topography, so the estimation accuracy in the very rugged areas is low. Current methods can be improved by removing the linear approximation and estimating the topography through simulated annealing and by using gravity gradiometry that is more sensitive to topography at short wavelengths than the gravity anomaly. Simulated annealing is a global optimization technique that can process nonlinear inverse problems. It is developed to estimate the seafloor depths by minimizing the difference between the observed and forward-computed vertical gravity gradients. The method is tested on altimetry-derived gravity gradients in a 2(degrees)x2(degrees) area of rugged seafloor topography in the West Pacific Ocean and results in estimates with a root-mean-square error of +/- 236 m. Compared to estimates from an existing model obtained by standard techniques this represents an accuracy improvement of 22%.

Gravity gradiometry on a moving platform, whether ground or airborne, has the potential to offer an efficient and accurate determination of the deflection of the vertical by simple line integration. A significant error in this process is a trend error that results from the integration of systematic gradient errors. Using an airborne full-tensor gradiometry data set of regularly spaced and intersecting tracks over a 10 km square region and the USDOV2012 vertical deflection model to calibrate these long wavelength errors, it is shown that the gradient-derived deflections agree with the USDOV2012 model at the level of 0.6–0.9 arcsec. Moreover, it is shown by graphical inspection that these differences represent high-frequency signal rather than error. Another data processing technique is examined using only (simulated) single-gradiometer instrument data, i.e., the local differential curvature components, \((\varGamma _{{ yy}} - \varGamma _{{ xx}})/2\) and \(\varGamma _{{ xy}}\), of the gravity field. While in theory these data can yield deflection components using two parallel data tracks, the results in the tested case are unsatisfactory due to implicit additional cross-track integration errors that accumulate systematically. The analysis thus demonstrates the importance of using the individual horizontal gradient components, \(\varGamma _{{ xx}}\), \(\varGamma _{{ yy}}\), to derive the deflection of the vertical.

The energy balance approach to geopotential modeling from satellite-to-satellite tracking offers the advantage of emulating an in-situ measurement system, analogous to satellite-borne gravity gradiometry, but in terms of differences in potential. Although the theoretical background is well known, based on the conservation of energy in celestial mechanics, its application to geopotential determination from satellite tracking had to await the capability for accurate, independent, kinematic orbit determination. This is now possible with Global Navigation Satellite Systems (GNSS), such as GPS. The precision tracking between two co-orbiting satellites brings additional short-wavelength information to the in-situ measurement. These tutorial notes derive the exact relationships between the gravitational potential and the orbital state vectors from the energy balance perspective, both in the inertial frame and in the Earth-fixed frame. Also derived are the equations that specifically incorporate the range-rate between co-orbiting satellites. Particular attention is given to the rotation potential and the temporal dependence of the potential due to various sources, including tidal variations, Earth's orientation and deformation, and terrestrial mass fluxes. A detailed analysis of magnitudes then leads to possibly acceptable approximations. It also shows for a satellite configuration such as GRACE (Gravity Recovery and Climate Experiment) that the radial component of the relative velocity between two co-orbiting satellites is as important in magnitude as the along-track component. Thus, the measured inter-satellite range-rate, for example, cannot be used alone to determine the potential difference. However, it is also shown that the short-wavelength content of the potential resides more in the along-track component than in the radial component, which demonstrates the significance of the range-rate measurement. A simple error analysis of the system identifies the requirements for state-vector accuracies (both for position and velocity) in relation to the range-rate accuracy. The observational equations for satellite-to-satellite tracking are derived and related to both global and regional geopotential modeling. The discussion concludes with a sample of published case studies that have demonstrated the energy balance approach and achieved tangible results in large-scale hydrological mass flux monitoring.

The inverse problem of estimating parameters (e.g., size, depth) of subsurface structures can be considered as an optimization problem where the parameters of a constructed forward model are estimated from observations collected on or above the Earth’s surface by minimizing the difference between the predicted model and the observations. Traditional solutions based on gradient-based approaches applied to nonlinear and non-unique problems basically depend on the initial conditions and may not always converge to the global minimum of the cost function if the starting model is far away from the true model. Alternatives to these straightforward approaches are innovative methods such as random search techniques that operate directly on the nonlinear models. This study compares a Monte-Carlo optimization method called Simulated Annealing (SA) to the Least-Squares Solution (LESS) within the general Gauss-Helmert formulation to estimate the parameters of a dip-slip fault from gravity gradient measurements as might be collected on profiles surveyed by an airborne system. It is shown that the SA algorithm is a more robust technique with respect to initial conditions in that it proceeds more comprehensively in parameter space and converges to their true values and thus the global minimum of the cost function. The SA algorithm is able to estimate the parameters of the fault as well as or better than LESS, and in the presence of significant background geologic and observation noise.

Two approaches have been formulated to compute the gravitational potential difference using low–low satellite-to-satellite tracking data based on energy integral: one in the geocentric inertial reference system, and the other in the terrestrial reference system. The focus of this work is on the approach in the geocentric inertial reference system, where a potential rotation term appears in addition to the potential term. In former formulations, the contribution of the time-variable components of the gravitational potential to the potential term was included, but their contribution to the potential rotation term was neglected. In this work, an improvement to the former formulations is made by reformulating the potential rotation term to include the contribution of the time-variable components of the gravitational potential. A simulation shows that our more accurate formulation of the potential rotation term is necessary to achieve the accuracy for recovering the temporal variation of the Earth’s gravity field, such as for use to the Gravity Recovery And Climate Experiment GRACE observation data based on this approach.