Optimized formulas for the gravitational field of a tesseroid

Various tasks in geodesy, geophysics, and related geosciences require precise information on the impact of mass distributions on gravity field-related quantities, such as the gravitational potential and its partial derivatives. Using forward modeling based on Newton’s integral, mass distributions are generally decomposed into regular elementary bodies. In classical approaches, prisms or point mass approximations are mostly utilized. Considering the effect of the sphericity of the Earth, alternative mass modeling methods based on tesseroid bodies (spherical prisms) should be taken into account, particularly in regional and global applications. Expressions for the gravitational field of a point mass are relatively simple when formulated in Cartesian coordinates. In the case of integrating over a tesseroid volume bounded by geocentric spherical coordinates, it will be shown that it is also beneficial to represent the integral kernel in terms of Cartesian coordinates. This considerably simplifies the determination of the tesseroid’s potential derivatives in comparison with previously published methodologies that make use of integral kernels expressed in spherical coordinates. Based on this idea, optimized formulas for the gravitational potential of a homogeneous tesseroid and its derivatives up to second-order are elaborated in this paper. These new formulas do not suffer from the polar singularity of the spherical coordinate system and can, therefore, be evaluated for any position on the globe. Since integrals over tesseroid volumes cannot be solved analytically, the numerical evaluation is achieved by means of expanding the integral kernel in a Taylor series with fourth-order error in the spatial coordinates of the integration point. As the structure of the Cartesian integral kernel is substantially simplified, Taylor coefficients can be represented in a compact and computationally attractive form. Thus, the use of the optimized tesseroid formulas particularly benefits from a significant decrease in computation time by about 45 % compared to previously used algorithms. In order to show the computational efficiency and to validate the mathematical derivations, the new tesseroid formulas are applied to two realistic numerical experiments and are compared to previously published tesseroid methods and the conventional prism approach.     

Harmonic analysis of the Earth's gravitational field by means of semi-continuous ephemerides of a low Earth orbiting GPS-tracked satellite. Case study: CHAMP

An algorithm for the determination of the spherical harmonic coefficients of the terrestrial gravitational field representation from the analysis of a kinematic orbit solution of a low earth orbiting GPS-tracked satellite is presented and examined. A gain in accuracy is expected since the kinematic orbit of a LEO satellite can nowadays be determined with very high precision, in the range of a few centimeters. In particular, advantage is taken of Newton's Law of Motion, which balances the acceleration vector with respect to an inertial frame of reference (IRF) and the gradient of the gravitational potential. By means of triple differences, and in particular higher-order differences (seven-point scheme, nine-point scheme), based upon Newton's interpolation formula, the local acceleration vector is estimated from relative GPS position time series. The gradient of the gravitational potential is conventionally given in a body-fixed frame of reference (BRF) where it is nearly time independent or stationary. Accordingly, the gradient of the gravitational potential has to be transformed from spherical BRF to Cartesian IRF. Such a transformation is possible by differentiating the gravitational potential, given as a spherical harmonics series expansion, with respect to Cartesian coordinates by means of the chain rule, and expressing zero- and first-order Ferrer's associated Legendre functions in terms of Cartesian coordinates. Subsequently, the BRF Cartesian coordinates are transformed into IRF Cartesian coordinates by means of the polar motion matrix, the precession–nutation matrices and the Greenwich sidereal time angle (GAST). In such a way a spherical harmonic representation of the terrestrial gravitational field intensity with respect to an IRF is achieved. Numerical tests of a resulting Gauss–Markov model document not only the quality and the high resolution of such a space gravity spectroscopy, but also the problems resulting from noise amplification in the acceleration determination process.     

Effects of topographic–isostatic masses on gravitational functionals at the Earth’s surface and at airborne and satellite altitudes

Topographic–isostatic masses represent an important source of gravity field information, especially in the high-frequency band, even if the detailed mass-density distribution inside the topographic masses is unknown. If this information is used within a remove-restore procedure, then the instability problems in downward continuation of gravity observations from aircraft or satellite altitudes can be reduced. In this article, integral formulae are derived for determination of gravitational effects of topographic–isostatic masses on the first- and second-order derivatives of the gravitational potential for three topographic–isostatic models. The application of these formulas is useful for airborne gravimetry/gradiometry and satellite gravity gradiometry. The formulas are presented in spherical approximation by separating the 3D integration in an analytical integration in the radial direction and 2D integration over the mean sphere. Therefore, spherical volume elements can be considered as being approximated by mass-lines located at the centre of the discretization compartments (the mass of the tesseroid is condensed mathematically along its vertical axis). The errors of this approximation are investigated for the second-order derivatives of the topographic–isostatic gravitational potential in the vicinity of the Earth’s surface. The formulas are then applied to various scenarios of airborne gravimetry/gradiometry and satellite gradiometry. The components of the gravitational vector at aircraft altitudes of 4 and 10 km have been determined, as well as the gravitational tensor components at a satellite altitude of 250 km envisaged for the forthcoming GOCE (gravity field and steady-state ocean-circulation explorer) mission. The numerical computations are based on digital elevation models with a 5-arc-minute resolution for satellite gravity gradiometry and 1-arc-minute resolution for airborne gravity/gradiometry.     

利用卫星引力梯度确定地球重力场的张量不变方法研究

应用张量不变理论对利用卫星重力梯度数据确定地球重力场的方法进行了研究,对张量不变观测方程的线性化处理、非全张量观测值的数据处理策略以及采用白噪声特性下的梯度观测值恢复地球重力场的精度等进行了数值分析.结果表明,张量不变解实现了不同观测值的联合求解,基于先验重力场模型的线性化方法在实现张量不变观测模型线性化处理的同时,提...     

Recursive computation of finite difference of associated Legendre functions

The existing methods to compute the definite integral of associated Legendre function (ALF) with respect to the argument suffer from a loss of significant figures independently of the latitude. This is caused by the subtraction of similar quantities in the additional term of their recurrence formulas, especially the finite difference of their values between two endpoints of the integration interval. In order to resolve the problem, we develop a recursive algorithm to compute their finite difference. Also, we modify the algorithm to evaluate their definite integrals assuming that their values at one endpoint are known. We numerically confirm a significant increase in computing precision of the integral by the new method. When the interval is one arc minute, for example, the gain amounts to 2–4 digits for the degree of harmonics in the range 2 ≤ n ≤ 2,048. This improvement in precision is achieved at a negligible increase in CPU time, say less than 5%.     

Satellite gradiometry using a satellite pair

The GRACE mission has substantiated the low–low satellite-to-satellite tracking (LL-SST) concept. The LL-SST configuration can be combined with the previously realized high–low SST concept in the CHAMP mission to provide a much higher accuracy. The line of sight (LOS) acceleration difference between the GRACE satellite pair, the simplest form of the combined observable, is mostly used for mapping the global gravity field of the Earth in terms of spherical harmonic coefficients. As an alternative observable, a linear combination of the gravitational gradient tensor components is proposed. Being a one-point function and having a direct relation with the field geometry (curvature of the field at the point) are two noteworthy achievements of the alternative formulation. In addition, using an observation quantity that is related to the second-instead of the first-order derivatives of the gravitational potential amplifies the high-frequency part of the signal. Since the transition from the first- to the second-order derivatives includes the application of a finite-differences scheme, the high-frequency part of the noise is also amplified. Nevertheless, due to the different spectral behaviour of signal and noise, in the end the second-order approach leads to improved gravitational field resolution. Mathematical formulae for the gradiometry approach, for both linear and higher-degree approximations, are derived. The proposed approach is implemented for recovery of the global gravitational field and the results are compared with those of LOS acceleration differences. Moreover, LOS acceleration difference residuals are calculated, which are at the level of a few tenths of mGal. Error analysis shows that the residuals of the estimated degree variances are less than 10^–3. Furthermore, the gravity anomaly residuals are less than 2 mGal for most points on the Earth.     

A comparison of the tesseroid,prism and point-mass approaches for mass reductions in gravity field modelling

The calculation of topographic (and iso- static) reductions is one of the most time-consuming operations in gravity field modelling. For this calculation, the topographic surface of the Earth is often divided with respect to geographical or map-grid lines, and the topographic heights are averaged over the respective grid elements. The bodies bounded by surfaces of constant (ellipsoidal) heights and geographical grid lines are denoted as tesseroids. Usually these ellipsoidal (or spherical) tesseroids are replaced by “equivalent” vertical rectangular prisms of the same mass. This approximation is motivated by the fact that the volume integrals for the calculation of the potential and its derivatives can be exactly solved for rectangular prisms, but not for the tesseroids. In this paper, an approximate solution of the spherical tesseroid integrals is provided based on series expansions including third-order terms. By choosing the geometrical centre of the tesseroid as the Taylor expansion point, the number of non-vanishing series terms can be greatly reduced. The zero-order term is equivalent to the point-mass formula. Test computations show the high numerical efficiency of the tesseroid method versus the prism approach, both regarding computation time and accuracy. Since the approximation errors due to the truncation of the Taylor series decrease very quickly with increasing distance of the tesseroid from the computation point, only the elements in the direct vicinity of the computation point have to be separately evaluated, e.g. by the prism formulas. The results are also compared with the point-mass formula. Further potential refinements of the tesseroid approach, such as considering ellipsoidal tesseroids, are indicated.     

A comparison of different mass elements for use in gravity gradiometry

Topographic and isostatic mass anomalies affect the external gravity field of the Earth. Therefore, these effects also exist in the gravity gradients observed, e.g., by the satellite gravity gradiometry mission GOCE (Gravity and Steady-State Ocean Circulation Experiment). The downward continuation of the gravitational signals is rather difficult because of the high-frequency behaviour of the combined topographic and isostatic effects. Thus, it is preferable to smooth the gravity field by some topographic-isostatic reduction. In this paper the focus is on the modelling of masses in the space domain, which can be subdivided into different mass elements and evaluated with analytical, semi-analytical and numerical methods. Five alternative mass elements are reviewed and discussed: the tesseroid, the point mass, the prism, the mass layer and the mass line. The formulae for the potential, the attraction components and the Marussi tensor of second-order potential derivatives are provided. The formulae for different mass elements and computation methods are checked by assuming a synthetic topography of constant height over a spherical cap and the position of the computation point on the polar axis. For this special situation an exact analytical solution for the tesseroid exists and a comparison between the analytical solution of a spherical cap and the modelling of different mass elements is possible. A comparison of the computation times shows that modelling by tesseroids with different methods produces the most accurate results in an acceptable computation time. As a numerical example, the Marussi tensor of the topographic effect is computed globally using tesseroids calculated by Gauss–Legendre cubature (3D) on the basis of a digital height model. The order of magnitude in the radial-radial component is about  ± 8 E.U. Electronic supplementary material  The online version of this article (doi:) contains supplementary material, which is available to authorized users.     

Gravitational gradients by tensor analysis with application to spherical coordinates

This contribution deals with the derivation of explicit expressions of the gradients of first, second and third order of the gravitational potential. This is accomplished in the framework of tensor analysis which naturally allows to apply general formulae to the specific coordinate systems in use in geodesy. In particular it is recalled here that when the potential field is expressed in general coordinates on a 3D manifold, the gradient operation leads to the definition of the covariant derivative and that the covariant derivative of a tensor can be obtained by application of a simple rule. When applied to the gravitational potential or to any of its gradients, the rule straightforwardly provides the expressions of the higher-order gradients. It is also shown that the tensor approach offers a clear distinction between natural and physical components of the gradients. Two fundamental reference systems—a global, bodycentric system and a local, topocentric system, both body-fixed—are introduced and transformation rules are derived to convert quantities between the two systems. The results include explicit expressions for the gradients of the first three orders in both reference systems.     

A new methodology to compute the gravitational contribution of a spherical tesseroid based on the analytical solution of a sector of a spherical zonal band

A new methodology for computing the gravitational effect of a spherical tesseroid has been devised and implemented. The methodology is based on the rotation from the global Earth-Centred Rotational reference frame to the local Earth-Centred P-Rotational reference frame, referred to the computation point P, and it requires knowledge of the height and the angular extension of each topographic column. After rotation, the gravitational effect of the tesseroid is computed via the effect of a sector of the spherical zonal band. In this respect, two possible procedures for handling the rotated tesseroids have been proposed and tested. The results obtained with the devised methodology are in good agreement with those derived by applying other existing methodologies.     

Evaluation of gravitational curvatures of a tesseroid in spherical integral kernels

Proper understanding of how the Earth’s mass distributions and redistributions influence the Earth’s gravity field-related functionals is crucial for numerous applications in geodesy, geophysics and related geosciences. Calculations of the gravitational curvatures (GC) have been proposed in geodesy in recent years. In view of future satellite missions, the sixth-order developments of the gradients are becoming requisite. In this paper, a set of 3D integral GC formulas of a tesseroid mass body have been provided by spherical integral kernels in the spatial domain. Based on the Taylor series expansion approach, the numerical expressions of the 3D GC formulas are provided up to sixth order. Moreover, numerical experiments demonstrate the correctness of the 3D Taylor series approach for the GC formulas with order as high as sixth order. Analogous to other gravitational effects (e.g., gravitational potential, gravity vector, gravity gradient tensor), numerically it is found that there exist the very-near-area problem and polar singularity problem in the GC east–east–radial, north–north–radial and radial–radial–radial components in spatial domain, and compared to the other gravitational effects, the relative approximation errors of the GC components are larger due to not only the influence of the geocentric distance but also the influence of the latitude. This study shows that the magnitude of each term for the nonzero GC functionals by a grid resolution 15\(^{{\prime } }\,\times \) 15\(^{{\prime }}\) at GOCE satellite height can reach of about 10\(^{-16}\) m\(^{-1}\) s\(^{2}\) for zero order, 10\(^{-24 }\) or 10\(^{-23}\) m\(^{-1}\) s\(^{2}\) for second order, 10\(^{-29}\) m\(^{-1}\) s\(^{2}\) for fourth order and 10\(^{-35}\) or 10\(^{-34}\) m\(^{-1}\) s\(^{2}\) for sixth order, respectively.     

Integral transformations of gradiometric data onto a GRACE type of observable

Integral transformations of gravitational gradients onto a Gravity Recovery And Climate Experiment (GRACE) type of observable are derived in this article. The gravitational gradients represent components of the gravitational tensor in the local north-oriented frame. The GRACE type of observable corresponds to a difference between two gravitational vectors as projected onto the line of sight between the two GRACE satellites. In total, three integral transformations relating vertical–vertical, vertical–horizontal and horizontal–horizontal gravitational gradients with the GRACE type of observable are provided. Spectral and closed forms of corresponding isotropic kernels are derived for each transformation. Special cases show that the integral transformations are general and relate gravitational gradients to many other quantities of the gravitational field, such as the gravitational vector, and its radial and tangential components. Correctness of the mathematical derivations is validated in a closed-loop simulation using synthetic data.     

利用重力异常数据构建重力梯度场的方法比较

为实现大范围、高精度基准重力梯度数据库的构建,考虑到重力梯度场对地形质量的敏感效应,一般利用恒密度数字高程模型来求取重力梯度值,从而忽略了地形密度变化以及水准面以下密度异常对重力梯度的影响。根据重力位理论中求解边值问题的数值应用方法,直接利用重力异常数据求取重力梯度场,弥补了密度变化和密度异常在重力梯度上的反映。根据模型算例和实测重力异常数据求取了剖面重力梯度值,结果表明,限于重力数据空间分辨率的影响,利用重力异常数据可恢复中长波段重力梯度场。该方法与地形数据求取重力梯度和卫星重力梯度测量等方法技术相结合,对重力梯度数据库的建设具有实际应用价值。     

Methodology and use of tensor invariants for satellite gravity gradiometry

Although its use is widespread in several other scientific disciplines, the theory of tensor invariants is only marginally adopted in gravity field modeling. We aim to close this gap by developing and applying the invariants approach for geopotential recovery. Gravitational tensor invariants are deduced from products of second-order derivatives of the gravitational potential. The benefit of the method presented arises from its independence of the gradiometer instrument’s orientation in space. Thus, we refrain from the classical methods for satellite gravity gradiometry analysis, i.e., in terms of individual gravity gradients, in favor of the alternative invariants approach. The invariants approach requires a tailored processing strategy. Firstly, the non-linear functionals with regard to the potential series expansion in spherical harmonics necessitates the linearization and iterative solution of the resulting least-squares problem. From the computational point of view, efficient linearization by means of perturbation theory has been adopted. It only requires the computation of reference gravity gradients. Secondly, the deduced pseudo-observations are composed of all the gravitational tensor elements, all of which require a comparable level of accuracy. Additionally, implementation of the invariants method for large data sets is a challenging task. We show the fundamentals of tensor invariants theory adapted to satellite gradiometry. With regard to the GOCE (Gravity field and steady-state Ocean Circulation Explorer) satellite gradiometry mission, we demonstrate that the iterative parameter estimation process converges within only two iterations. Additionally, for the GOCE configuration, we show the invariants approach to be insensitive to the synthesis of unobserved gravity gradients.     

The spacetime gravitational field of a deformable body

The high-resolution analysis of orbit perturbations of terrestrial artificial satellites has documented that the eigengravitation of a massive body like the Earth changes in time, namely with periodic and aperiodic constituents. For the space-time variation of the gravitational field the action of internal and external volume as well as surface forces on a deformable massive body are responsible. Free of any assumption on the symmetry of the constitution of the deformable body we review the incremental spatial (“Eulerian”) and material (“Lagrangean”) gravitational field equations, in particular the source terms (two constituents: the divergence of the displacement field as well as the projection of the displacement field onto the gradient of the reference mass density function) and the `jump conditions' at the boundary surface of the body as well as at internal interfaces both in linear approximation. A spherical harmonic expansion in terms of multipoles of the incremental Eulerian gravitational potential is presented. Three types of spherical multipoles are identified, namely the dilatation multipoles, the transport displacement multipoles and those multipoles which are generated by mass condensation onto the boundary reference surface or internal interfaces. The degree-one term has been identified as non-zero, thus as a “dipole moment” being responsible for the varying position of the deformable body's mass centre. Finally, for those deformable bodies which enjoy a spherically symmetric constitution, emphasis is on the functional relation between Green functions, namely between Fourier-/ Laplace-transformed volume versus surface Love-Shida functions (h(r),l(r) versus h ^′(r),l ^′(r)) and Love functions k(r) versus k ^′(r). The functional relation is numerically tested for an active tidal force/potential and an active loading force/potential, proving an excellent agreement with experimental results. Received: December 1995 / Accepted: 1 February 1997     

卫星重力梯度数据解算位系数的最小二乘配置法

卫星重力梯度测量在恢复地球重力场的研究中已经得到了广泛应用。本文通过空间扰动位协方差函数特性,得出卫星重力梯度数据与引力位系数的相关协方差函数。利用最小二乘配置法,最终推导出由重力梯度数据直接解算引力位系数的函数表达式,并简要分析其实用性。     

The optimum expression for the gravitational potential of polyhedral bodies having a linearly varying density distribution

When topography is represented by a simple regular grid digital elevation model, the analytical rectangular prism approach is often used for a precise gravity field modelling at the vicinity of the computation point. However, when the topographical surface is represented more realistically, for instance by a triangular irregular network (TIN) model, the analytical integration using arbitrary polyhedral bodies (the analytical line integral approach) can be implemented directly without additional data pre-processing (gridding or interpolation). The analytical line integral approach can also facilitate 3-D density models created for complex geometrical bodies. For the forward modelling of the gravitational field generated by the geological structures with variable densities, the analytical integration can be carried out using polyhedral bodies with a varying density. The optimal expression for the gravitational attraction vector generated by an arbitrary polyhedral body having a linearly varying density is known. In this article, the corresponding optimal expression for the gravitational potential is derived by means of line integrals after applying the Gauss divergence theorem.     

基于χ2准则的磁梯度张量3D聚焦反演方法

针对铁磁性物质反演中正则化参数自适应选择的问题,提出了基于χ2准则的磁梯度张量3D聚焦反演方法。利用深度加权矩阵和最小支撑矩阵对经典Tikhonov正则化理论框架下的反演模型进行约束得到目标函数,避免了由于反演参数多于采集点数而导致反演解的多解性,并有效解决了核函数随深度增大而快速衰减的问题。通过对目标函数进行迭代奇异值分解获得最佳物性参数,并根据χ2准则自适应地确定目标函数在迭代过程中的正则化参数,提高了迭代速度和求解精度。仿真和实验结果表明:该方法能准确还原磁性异常体的轮廓形态,具有较好的模型分辨率。     

On the evaluation of the gravity effects of polyhedral bodies and a consistent treatment of related singularities

We show that the singularities which can affect the computation of the gravity effects (potential, gravity and tensor gradient fields) can be systematically addressed by invoking distribution theory and suitable formulas of differential calculus. Thus, differently from previous contributions on the subject, the use of a-posteriori corrections of the formulas derived in absence of singularities can be ruled out. The general approach presented in the paper is further specialized to the case of polyhedral bodies and detailed for a rectangular prism having a constant mass density. With reference to this last case, we derive novel expressions for the related gravitational field, as well as for its first and second derivative, at an observation point coincident with a prism vertex and show that they turn out to be more compact than the ones reported in the specialized literature.     

Simultaneous determination of gravitational acceleration and ground vibrations by free fall experiments

We investigated how well we can simultaneously determine the gravitational acceleration and ground vibrations with many free falls at short intervals by using synthetic data which contain quadratic functions, sinusoids and white noise in order to improve the accuracy of absolute gravimetry and to realize the measurements without a seismometer. As a result of simulations, lower white noise, longer dropping time and longer duration of a series of free falls improve the accuracy in determination of g. The sampling internvals and the dropping intervals, on the other hand, hardly affect the accuracy so long as they are shorter than Nyquist intervals. Both the gravitational acceleration and the waveform of the long period vibrations can be well determined with a series of free falls of about only 2 cm dropping distance without deterioration in the accuracy with this procedure.