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.     

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.     

Ellipsoidal terrain correction based on multi-cylindrical equal-area map projection of the reference ellipsoid

An operational algorithm for computation of terrain correction (or local gravity field modeling) based on application of closed-form solution of the Newton integral in terms of Cartesian coordinates in multi-cylindrical equal-area map projection of the reference ellipsoid is presented. Multi-cylindrical equal-area map projection of the reference ellipsoid has been derived and is described in detail for the first time. Ellipsoidal mass elements with various sizes on the surface of the reference ellipsoid are selected and the gravitational potential and vector of gravitational intensity (i.e. gravitational acceleration) of the mass elements are computed via numerical solution of the Newton integral in terms of geodetic coordinates {,,h}. Four base- edge points of the ellipsoidal mass elements are transformed into a multi-cylindrical equal-area map projection surface to build Cartesian mass elements by associating the height of the corresponding ellipsoidal mass elements to the transformed area elements. Using the closed-form solution of the Newton integral in terms of Cartesian coordinates, the gravitational potential and vector of gravitational intensity of the transformed Cartesian mass elements are computed and compared with those of the numerical solution of the Newton integral for the ellipsoidal mass elements in terms of geodetic coordinates. Numerical tests indicate that the difference between the two computations, i.e. numerical solution of the Newton integral for ellipsoidal mass elements in terms of geodetic coordinates and closed-form solution of the Newton integral in terms of Cartesian coordinates, in a multi-cylindrical equal-area map projection, is less than 1.6×10^–8 m²/s² for a mass element with a cross section area of 10×10 m and a height of 10,000 m. For a mass element with a cross section area of 1×1 km and a height of 10,000 m the difference is less than 1.5×10^–4m²/s². Since 1.5× 10^–4 m²/s² is equivalent to 1.5×10^–5m in the vertical direction, it can be concluded that a method for terrain correction (or local gravity field modeling) based on closed-form solution of the Newton integral in terms of Cartesian coordinates of a multi-cylindrical equal-area map projection of the reference ellipsoid has been developed which has the accuracy of terrain correction (or local gravity field modeling) based on the Newton integral in terms of ellipsoidal coordinates.Acknowledgments. This research has been financially supported by the University of Tehran based on grant number 621/4/859. This support is gratefully acknowledged. The authors are also grateful for the comments and corrections made to the initial version of the paper by Dr. S. Petrovic from GFZ Potsdam and the other two anonymous reviewers. Their comments helped to improve the structure of the paper significantly.     

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.     

Upward continuation of Markov type anomalous gravity potential models

Linear gravity field state space models are still a useful tool to model the anomalous gravity field in vector gravimetry, airborne gravimetry, inertial geodesy and navigation. This paper deals with an idea ofJordan and Heller (1978) to solve analytically the upward continuation problem of Markov gravity models.In contrary to the standard Markov shaping filter approach the height dependency of the covariance function, i.e. variance factor and correlation length as function of height, is strictly introduced in state space and not neglected. Using some basic integral transforms, a general upward continuation integral is derived for the n-th order Markov process. The upward continuation integral is solved for the special and practically important case of 2nd order Markov process in very detail. This leads to the introduction of the special sine and cosine integral functions into the the mathematical covariance model. The features of the covariance model are analyzed analytically and the height dependency is discussed numerically.     

Downward continuation and geoid determination based on band-limited airborne gravity data

 The downward continuation of the harmonic disturbing gravity potential, derived at flight level from discrete observations of airborne gravity by the spherical Hotine integral, to the geoid is discussed. The initial-boundary-value approach, based on both the direct and inverse solution to Dirichlet's problem of potential theory, is used. Evaluation of the discretized Fredholm integral equation of the first kind and its inverse is numerically tested using synthetic airborne gravity data. Characteristics of the synthetic gravity data correspond to typical airborne data used for geoid determination today and in the foreseeable future: discrete gravity observations at a mean flight height of 2 to 6 km above mean sea level with minimum spatial resolution of 2.5 arcmin and a noise level of 1.5 mGal. Numerical results for both approaches are presented and discussed. The direct approach can successfully be used for the downward continuation of airborne potential without any numerical instabilities associated with the inverse approach. In addition to these two-step approaches, a one-step procedure is also discussed. This procedure is based on a direct relationship between gravity disturbances at flight level and the disturbing gravity potential at sea level. This procedure provided the best results in terms of accuracy, stability and numerical efficiency. As a general result, numerically stable downward continuation of airborne gravity data can be seen as another advantage of airborne gravimetry in the field of geoid determination. Received: 6 June 2001 / Accepted: 3 January 2002     

Wavelet Modeling of Regional and Temporal Variations of the Earth’s Gravitational Potential Observed by GRACE

This work is dedicated to the wavelet modeling of regional and temporal variations of the Earth’s gravitational potential observed by the GRACE (gravity recovery and climate experiment) satellite mission. In the first part, all required mathematical tools and methods involving spherical wavelets are provided. Then, we apply our method to monthly GRACE gravity fields. A strong seasonal signal can be identified which is restricted to areas where large-scale redistributions of continental water mass are expected. This assumption is analyzed and verified by comparing the time-series of regionally obtained wavelet coefficients of the gravitational signal originating from hydrology models and the gravitational potential observed by GRACE. The results are in good agreement with previous studies and illustrate that wavelets are an appropriate tool to investigate regional effects in the Earth’s gravitational field. Electronic Supplementary Material Supplementary material is available for this article at     

Accurate computation of gravitational field of a tesseroid

We developed an accurate method to compute the gravitational field of a tesseroid. The method numerically integrates a surface integral representation of the gravitational potential of the tesseroid by conditionally splitting its line integration intervals and by using the double exponential quadrature rule. Then, it evaluates the gravitational acceleration vector and the gravity gradient tensor by numerically differentiating the numerically integrated potential. The numerical differentiation is conducted by appropriately switching the central and the single-sided second-order difference formulas with a suitable choice of the test argument displacement. If necessary, the new method is extended to the case of a general tesseroid with the variable density profile, the variable surface height functions, and/or the variable intervals in longitude or in latitude. The new method is capable of computing the gravitational field of the tesseroid independently on the location of the evaluation point, namely whether outside, near the surface of, on the surface of, or inside the tesseroid. The achievable precision is 14–15 digits for the potential, 9–11 digits for the acceleration vector, and 6–8 digits for the gradient tensor in the double precision environment. The correct digits are roughly doubled if employing the quadruple precision computation. The new method provides a reliable procedure to compute the topographic gravitational field, especially that near, on, and below the surface. Also, it could potentially serve as a sure reference to complement and elaborate the existing approaches using the Gauss–Legendre quadrature or other standard methods of numerical integration.     

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.     

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

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

An integrated wavelet concept of physical geodesy

For the determination of the earth's gravity field many types of observations are nowadays available, including terrestrial gravimetry, airborne gravimetry, satellite-to-satellite tracking, satellite gradio-metry, etc. The mathematical connection between these observables on the one hand and gravity field and shape of the earth on the other is called the integrated concept of physical geodesy. In this paper harmonic wavelets are introduced by which the gravitational part of the gravity field can be approximated progressively better and better, reflecting an increasing flow of observations. An integrated concept of physical geodesy in terms of harmonic wavelets is presented. Essential tools for approximation are integration formulas relating an integral over an internal sphere to suitable linear combinations of observation functionals, i.e. linear functionals representing the geodetic observables. A scale discrete version of multiresolution is described for approximating the gravitational potential outside and on the earth's surface. Furthermore, an exact fully discrete wavelet approximation is developed for the case of band-limited wavelets. A method for combined global outer harmonic and local harmonic wavelet modelling is proposed corresponding to realistic earth's models. As examples, the role of wavelets is discussed for the classical Stokes problem, the oblique derivative problem, satellite-to-satellite tracking, satellite gravity gradiometry and combined satellite-to-satellite tracking and gradiometry. Received: 28 February 1997 / Accepted: 17 November 1997     

固体潮对地球重力场时变特征影响的潮波公式

精密和详细测定地球重力场及其随时间变化,是目前卫星重力测量的主要课题.基于目前高精度的固体潮展开,研究固体潮对地球重力场时变特征的影响.不同于IERS2000推荐的固体潮对重力影响的理论模型,独立于天体历书而基于精密引潮力位展开直接给出在毫微伽精度下的潮波公式,并考虑到四阶潮汐的效应.同时对重力数据归算中的永久性潮汐的处理进行总结和说明.本文的工作可为高精度的地球重力场的研究提供理论依据和参考.     

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.     

ITG-CHAMP01: a CHAMP gravity field model from short kinematic arcs over a one-year observation period

Global gravity field models have been determined based on kinematic orbits covering an observation period of one year beginning from March 2002. Three different models have been derived up to a maximum degree of n=90 of a spherical harmonic expansion of the gravitational potential. One version, ITG-CHAMP01E, has been regularized beginning from degree n=40 upwards, based on the potential coefficients of the gravity field model EGM96. A second model, ITG-CHAMP01K, has been determined based on Kaulas rule of thumb, also beginning from degree n=40. A third version, ITG-CHAMP01S, has been determined without any regularization. The physical model of the gravity field recovery technique is based on Newtons equation of motion, formulated as a boundary value problem in the form of a Fredholm-type integral equation. The observation equations are formulated in the space domain by dividing the one-year orbit into short sections of approximately 30-minute arcs. For every short arc, a variance factor has been determined by an iterative computation procedure. The three gravity field models have been validated based on various criteria, and demonstrate the quality of not only the gravity field recovery technique but also the kinematically determined orbits.     

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.     

The spherical horizontal and spherical vertical boundary value problem – vertical deflections and geoidal undulations – the completed Meissl diagram

 In a comparison of the solution of the spherical horizontal and vertical boundary value problems of physical geodesy it is aimed to construct downward continuation operators for vertical deflections (surface gradient of the incremental gravitational potential) and for gravity disturbances (vertical derivative of the incremental gravitational potential) from points on the Earth's topographic surface or of the three-dimensional (3-D) Euclidean space nearby down to the international reference sphere (IRS). First the horizontal and vertical components of the gravity vector, namely spherical vertical deflections and spherical gravity disturbances, are set up. Second, the horizontal and vertical boundary value problem in spherical gravity and geometry space is considered. The incremental gravity vector is represented in terms of vector spherical harmonics. The solution of horizontal spherical boundary problem in terms of the horizontal vector-valued Green function converts vertical deflections given on the IRS to the incremental gravitational potential external in the 3-D Euclidean space. The horizontal Green functions specialized to evaluation and source points on the IRS coincide with the Stokes kernel for vertical deflections. Third, the vertical spherical boundary value problem is solved in terms of the vertical scalar-valued Green function. Fourth, the operators for upward continuation of vertical deflections given on the IRS to vertical deflections in its external 3-D Euclidean space are constructed. Fifth, the operators for upward continuation of incremental gravity given on the IRS to incremental gravity to the external 3-D Euclidean space are generated. Finally, Meissl-type diagrams for upward continuation and regularized downward continuation of horizontal and vertical gravity data, namely vertical deflection and incremental gravity, are produced. Received: 10 May 2000 / Accepted: 26 February 2001     

引潮力引起的地球引力场时空变化的理论研究

地球重力场的时空结构与分布特性无论在基础理论研究,还是在地理空间信息建设中都具有重要的意义。地球表面上的观测仪器检测到的只是某一个质点的重力变化(拉格朗日变化),但是理论研究却是基于欧拉引力和引力位场方程进行的。本文基于连续介质力学的基本理论,在改正了文献[5]的一处原则性错误之后,推导出了引力和引力位的Lagrange和Euler增量表达式。本文工作对于高精度地球重力场的时空变化的理论研究具有参考价值。     

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

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

由星载GPS相位数据确定地球重力场模型若干问题研究

人造地球卫星在地球引力场中运动,可以探测地球重力场的长波信息。随着GPS技术的发展,星载GPS技术日趋成熟,因此由星载GPS相位数据确定地球重力场模型是当前国际地学研究的热点之一。本文给出了确定地球重力场模型中的星载GPS星地相位双差观测量,阐述了Cowell II数值轨道积分公式,导出了参数估计中星地双差观测量的偏导数,利用分块Bayes最小二乘参数估计地球引力场位系数等有关参数。