Solving Molodensky’s series by fast Fourier transform techniques

The use of the fast Fourier transform algorithm in the evaluation of the Molodensky series terms is demonstrated in this paper. The solution by analytical continuation to point level has been reformulated to obtain convolution integrals in planar approximation which can be efficiently evaluated in the frequency domain. Preliminary results show that the solution by Faye anomalies is not sufficient for highly accurate deflections of the vertical and height anomalies. The Molodensky solution up to at least the second-order term must be carried out. Part of the unrecovered deflection and height anomaly signal appears to be due to density variations, verifying the essential role of density modelling. A remove-restore technique for the terrain effects can improve the convergence of the series and minimize the interpolation errors. Paper presented at theI Hotine-Marussi Symposium on Mathematical Geodesy, Rome, June 3–6, 1985.     

The Abel-Poisson kernel and the Abel-Poisson integral in a moving tangent space

The upward-downward continuation of a harmonic function like the gravitational potential is conventionally based on the direct-inverse Abel-Poisson integral with respect to a sphere of reference. Here we aim at an error estimation of the “planar approximation” of the Abel-Poisson kernel, which is often used due to its convolution form. Such a convolution form is a prerequisite to applying fast Fourier transformation techniques. By means of an oblique azimuthal map projection / projection onto the local tangent plane at an evaluation point of the reference sphere of type “equiareal” we arrive at a rigorous transformation of the Abel-Poisson kernel/Abel-Poisson integral in a convolution form. As soon as we expand the “equiareal” Abel-Poisson kernel/Abel-Poisson integral we gain the “planar approximation”. The differences between the exact Abel-Poisson kernel of type “equiareal” and the “planar approximation” are plotted and tabulated. Six configurations are studied in detail in order to document the error budget, which varies from 0.1% for points at a spherical height H=10km above the terrestrial reference sphere up to 98% for points at a spherical height H = 6.3×10⁶km. Received: 18 March 1997 / Accepted: 19 January 1998     

A solution to the downward continuation effect on the geoid determined by Stokes' formula

 The analytical continuation of the surface gravity anomaly to sea level is a necessary correction in the application of Stokes' formula for geoid estimation. This process is frequently performed by the inversion of Poisson's integral formula for a sphere. Unfortunately, this integral equation corresponds to an improperly posed problem, and the solution is both numerically unstable, unless it is well smoothed, and tedious to compute. A solution that avoids the intermediate step of downward continuation of the gravity anomaly is presented. Instead the effect on the geoid as provided by Stokes' formula is studied directly. The practical solution is partly presented in terms of a truncated Taylor series and partly as a truncated series of spherical harmonics. Some simple numerical estimates show that the solution mostly meets the requests of a 1-cm geoid model, but the truncation error of the far zone must be studied more precisely for high altitudes of the computation point. In addition, it should be emphasized that the derived solution is more computer efficient than the detour by Poisson's integral. Received: 6 February 2002 / Accepted: 18 November 2002 Acknowledgements. Jonas ?gren carried out the numerical calculations and gave some critical and constructive remarks on a draft version of the paper. This support is cordially acknowledged. Also, the thorough work performed by one unknown reviewer is very much appreciated.     

Accounting for topography in the calculation of quasigeoidal heights and plumb-line deflections from gravity anomalies

A calculation of quasigeoidal heights and plumb-line deflections according to Molodensky formulae was carried out under elimination of the effect of topography from gravity anomalies. After the masses of topography had been removed a smoothed-out surface passing through astronomical and gravity stations was considered as representing the physical surface of the Earth. Thus it has been practically rendered possible to use the first-approximation formulae of Molodensky, and, in many cases, also the “zero-approximation” formulae analogous to the formulae of Stokes and Vening-Meinesz. The effect of the restored masses of topography was then added to the quantities found; the said effect was expressed as the effect of topography condensed on the normal equipotential surface passing through the point under investigation, plus a correction for condensation. Following some transformations, the resulting formulae (13) and (18) were obtained which formulae differ in their “zero-approximation” (15) and (20) from traditional formulas in that they contain terrait reductions added to free-air anomalies. Moreover, in the calculation of plumb-line deflections directly in mountain regions a correction for differing effects of topography before and after its condensation is to be introduced. A tentative expansion of terrain reduction in terms of spherical harmonics up to the third order is given; it can be seen therefrom that the Stokes series in its usual form is subject to a mean arror about 15–20%. It is also shown that the expansion of free-air anomalies in terms of spherical functions contains a first-order harmonic with a mean values about ±0.3 mgl. The said harmonic practically disappears in the expansion of the sum of free-air anomalies and terrain reductions.     

On some problems of the downward continuation of the 5′×5′ mean Helmert gravity disturbance

This research deals with some theoretical and numerical problems of the downward continuation of mean Helmert gravity disturbances. We prove that the downward continuation of the disturbing potential is much smoother, as well as two orders of magnitude smaller than that of the gravity anomaly, and we give the expression in spectral form for calculating the disturbing potential term. Numerical results show that for calculating truncation errors the first 180^∘ of a global potential model suffice. We also discuss the theoretical convergence problem of the iterative scheme. We prove that the 5^′×5^′ mean iterative scheme is convergent and the convergence speed depends on the topographic height; for Canada, to achieve an accuracy of 0.01 mGal, at most 80 iterations are needed. The comparison of the “mean” and “point” schemes shows that the mean scheme should give a more reasonable and reliable solution, while the point scheme brings a large error to the solution. Received: 19 August 1996 / Accepted: 4 February 1998     

Global height datum unification: a new approach in gravity potential space

The problem of “global height datum unification” is solved in the gravity potential space based on: (1) high-resolution local gravity field modeling, (2) geocentric coordinates of the reference benchmark, and (3) a known value of the geoid’s potential. The high-resolution local gravity field model is derived based on a solution of the fixed-free two-boundary-value problem of the Earth’s gravity field using (a) potential difference values (from precise leveling), (b) modulus of the gravity vector (from gravimetry), (c) astronomical longitude and latitude (from geodetic astronomy and/or combination of (GNSS) Global Navigation Satellite System observations with total station measurements), (d) and satellite altimetry. Knowing the height of the reference benchmark in the national height system and its geocentric GNSS coordinates, and using the derived high-resolution local gravity field model, the gravity potential value of the zero point of the height system is computed. The difference between the derived gravity potential value of the zero point of the height system and the geoid’s potential value is computed. This potential difference gives the offset of the zero point of the height system from geoid in the “potential space”, which is transferred into “geometry space” using the transformation formula derived in this paper. The method was applied to the computation of the offset of the zero point of the Iranian height datum from the geoid’s potential value W ₀=62636855.8 m²/s². According to the geometry space computations, the height datum of Iran is 0.09 m below the geoid.     

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.     

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     

航空重力测量数据向下延拓的正则化算法及其谱分解

基于Poisson积分方程,提出了以地面重力观测值作为控制并顾及外区影响的向下解析延拓数学模型,推导了向下解析延拓的谱分解式,在频域内分析了造成向下延拓结果不稳定的原因,进而给出了向下延拓的正则化算法,并讨论了向下延拓中的地形影响.通过对我国首次航空重力测量试验数据的处理表明,提出的方法可获得稳定、精确的向下延拓结果.     

A general method for an explicit determination of the shape of the earth from gravimetric data

Summary The principal formulae of the geophysical geodesy are based on the famous explicit expression of Stokes (1849). Up to now, there has been no method for a computation of the corresponding explicit expression of a (non-spherical) surface with masses outside the geoid. In this paper there is a solution of this problem. Another paper on this subject was presented to “Nordiska geodetm?tet” in Copenhagen May 1959 (in Swedish).     

Far-zone effects for different topographic-compensation models based on a spherical harmonic expansion of the topography

The determination of the gravimetric geoid is based on the magnitude of gravity observed at the surface of the Earth or at airborne altitude. To apply the Stokes’s or Hotine’s formulae at the geoid, the potential outside the geoid must be harmonic and the observed gravity must be reduced to the geoid. For this reason, the topographic (and atmospheric) masses outside the geoid must be “condensed” or “shifted” inside the geoid so that the disturbing gravity potential T fulfills Laplace’s equation everywhere outside the geoid. The gravitational effects of the topographic-compensation masses can also be used to subtract these high-frequent gravity signals from the airborne observations and to simplify the downward continuation procedures. The effects of the topographic-compensation masses can be calculated by numerical integration based on a digital terrain model or by representing the topographic masses by a spherical harmonic expansion. To reduce the computation time in the former case, the integration over the Earth can be divided into two parts: a spherical cap around the computation point, called the near zone, and the rest of the world, called the far zone. The latter one can be also represented by a global spherical harmonic expansion. This can be performed by a Molodenskii-type spectral approach. This article extends the original approach derived in Novák et al. (J Geod 75(9–10):491–504, 2001), which is restricted to determine the far-zone effects for Helmert’s second method of condensation for ground gravimetry. Here formulae for the far-zone effects of the global topography on gravity and geoidal heights for Helmert’s first method of condensation as well as for the Airy-Heiskanen model are presented and some improvements given. Furthermore, this approach is generalized for determining the far-zone effects at aeroplane altitudes. Numerical results for a part of the Canadian Rocky Mountains are presented to illustrate the size and distributions of these effects.     

Mixed Integer-Real Valued Adjustment (IRA) Problems: GPS Initial Cycle Ambiguity Resolution by Means of the LLL Algorithm

In order to achieve to GPS solutions of first-order accuracy and integrity, carrier phase observations as well as pseudorange observations have to be adjusted with respect to a linear/linearized model. Here the problem of mixed integer-real valued parameter adjustment (IRA) is met. Indeed, integer cycle ambiguity unknowns have to be estimated and tested. At first we review the three concepts to deal with IRA: (i) DDD or triple difference observations are produced by a properly chosen difference operator and choice of basis, namely being free of integer-valued unknowns (ii) The real-valued unknown parameters are eliminated by a Gauss elimination step while the remaining integer-valued unknown parameters (initial cycle ambiguities) are determined by Quadratic Programming and (iii) a RA substitute model is firstly implemented (real-valued estimates of initial cycle ambiguities) and secondly a minimum distance map is designed which operates on the real-valued approximation of integers with respect to the integer data in a lattice. This is the place where the integer Gram-Schmidt orthogonalization by means of the LLL algorithm (modified LLL algorithm) is applied being illustrated by four examples. In particular, we prove that in general it is impossible to transform an oblique base of a lattice to an orthogonal base by Gram-Schmidt orthogonalization where its matrix enties are integer. The volume preserving Gram-Schmidt orthogonalization operator constraint to integer entries produces “almost orthogonal” bases which, in turn, can be used to produce the integer-valued unknown parameters (initial cycle ambiguities) from the LLL algorithm (modified LLL algorithm). Systematic errors generated by “almost orthogonal” lattice bases are quantified by A. K. Lenstra et al. (1982) as well as M. Pohst (1987). The solution point of Integer Least Squares generated by the LLL algorithm is = (L')⁻¹[L'◯] ∈ ℤ ^m where L is the lower triangular Gram-Schmidt matrix rounded to nearest integers, [L], and = [L'◯] are the nearest integers of L'◯, ◯ being the real valued approximation of z ∈ ℤ ^m , the m-dimensional lattice space Λ. Indeed due to “almost orthogonality” of the integer Gram-Schmidt procedure, the solution point is only suboptimal, only close to “least squares.” ? 2000 John Wiley & Sons, Inc.     

应用格林积分直接以地面边值确定外部扰动重力场

以地面为边界的Molodensky问题通常得到的是级数形式的解,高阶项体现了地面边值到某一光滑面上的改正,在应用中不仅遇到计算的复杂性、稳定性问题,也存在对数据密集要求的困难。本文从推求外部扰动重力场的应用出发,将格林公式用于数字地形面,在忽略水平分量的积分影响的情况下,得到以地面重力异常和高程异常差为边值确定外部扰动位的表达式,其核函数分别为距离倒数和Poisson核。该方法不需要对地面数据进行延拓处理,且核函数形式简洁,适于外部扰动重力场的随机计算。     

The effect of topography on the determination of the geoid using analytical downward continuation

The method of analytical downward continuation has been used for solving Molodensky’s problem. This method can also be used to reduce the surface free air anomaly to the ellipsoid for the determination of the coefficients of the spherical harmonic expansion of the geopotential. In the reduction of airborne or satellite gradiometry data, if the sea level is chosen as reference surface, we will encounter the problem of the analytical downward continuation of the disturbing potential into the earth, too. The goal of this paper is to find out the topographic effect of solving Stoke’sboundary value problem (determination of the geoid) by using the method of analytical downward continuation. It is shown that the disturbing potential obtained by using the analytical downward continuation is different from the true disturbing potential on the sea level mostly by a −2πGρh 2/p. This correction is important and it is very easy to compute and add to the final results. A terrain effect (effect of the topography from the Bouguer plate) is found to be much smaller than the correction of the Bouguer plate and can be neglected in most cases. It is also shown that the geoid determined by using the Helmert’s second condensation (including the indirect effect) and using the analytical downward continuation procedure (including the topographic effect) are identical. They are different procedures and may be used in different environments, e.g., the analytical downward continuation procedure is also more convenient for processing the aerial gravity gradient data. A numerical test was completed in a rough mountain area, 35°<ϕ<38°, 240°<λ<243°. A digital height model in 30″×30″ point value was used. The test indicated that the terrain effect in the test area has theRMS value ±0.2−0.3 cm for geoid. The topographic effect on the deflections of the vertical is around1 arc second.     

A new series solution of Molodensky's problem

A complete series solution of Molodensky's boundary-value problem is derived using, instead of an integral equation, analytical continuation by means of power series. This solution is shown to be equivalent, term by term, to the Molodensky-Brovar series, but is simpler and practically more convenient.

---

Zusammenfassung Eine vollst?ndige Reihenl?sung des Problems von Molodensky wird hergeleitet, wobei anstatt einer Integralgleichung die analytische Fortsetzung mittels Potenzreihen zugrunde gelegt wird. Es zeigt sich, dass diese L?sung gliedweise ?quivalent zur Reihe von Molodensky-Brovar ist, aber sie ist einfacher und praktisch brauchbarer.

---

Résumé On déduit une série qui donne une solution complète du problème de Molodensky, en utilisant, au lieu d'une équation intégrale, la continuation analytique par une série de puissances. Il s'ensuit que cette solution est équivalente à la série de Molodensky-Brovar, mais elle est plus simple et plus pratique.

     

顾及地形效应的地面重力向上延拓模型分析与检验

重力向上延拓在外部重力场逼近和航空重力测量数据质量评估中具有重要应用。本文深入分析研究了6种向上延拓计算模型的技术特点和适用条件,提出了应用超高阶位模型加地形改正、点质量方法结合移去-恢复技术实现“先向下后向上延拓”计算的实施策略,探讨了计算过程特别是前端向下延拓过程的稳定性问题。通过实际数值计算,定量评估了地形质量对不同高度向上延拓结果的影响,对比分析了不同向上延拓模型顾及地形效应的实际效果,同时对向上延拓模型计算精度进行了估计。在地形变化比较激烈的山区,地形质量对向上延拓结果的影响最大可达几十个mGal（10-5m·s-2）,当计算高度为10 km时,该项影响超过3 mGal;向上延拓计算模型误差（不含数据误差影响）一般不超过1 mGal;基于超高阶位模型和地形改正信息实施向下延拓过渡的布阿桑（Poisson）积分向上延拓模型,具有计算过程简便、计算结果稳定可靠等优点。     

＂Digital Region＂ in the Context of the ＂Grid Computing＂

This paper is to construct a “digital local, regional, region“ information framework based on the technology of “SIG“ and its significance and application to the regional sustainable development evaluation system. First, the concept of the “grid computing“ and “SIG“ is interpreted and discussed, then the relationship between the “grid computing“ and “digital region“ is analyzed, and the framework of the “digital region“ is put forward. Finally, the significance and application of “grid computing“ to the “region sustainable development evaluation system“ are discussed.     

The topographic bias by analytical continuation in physical geodesy

This study emphasizes that the harmonic downward continuation of an external representation of the Earth’s gravity potential to sea level through the topographic masses implies a topographic bias. It is shown that the bias is only dependent on the topographic density along the geocentric radius at the computation point. The bias corresponds to the combined topographic geoid effect, i.e., the sum of the direct and indirect topographic effects. For a laterally variable topographic density function, the combined geoid effect is proportional to terms of powers two and three of the topographic height, while all higher order terms vanish. The result is useful in geoid determination by analytical continuation, e.g., from an Earth gravity model, Stokes’s formula or a combination thereof.     

重力异常向上延拓严密改化模型及向下延拓应用

重力异常向上延拓全球积分模型在航空重力测量数据质量评估和向下延拓迭代计算等领域具有广泛的应用。为了消除积分核函数奇异性影响,需要对该模型进行基于积分恒等式的移去-恢复转换及全球积分域的分区改化处理。在此过程中,传统改化处理方法往往忽略了全球积分过渡到局域积分引起的积分恒等式偏差影响,从而导致不必要的计算模型误差,最终影响向上延拓计算结果的可靠性,甚至影响向下延拓迭代解算结果的稳定性。针对此问题,本文开展了重力异常向上延拓积分模型改化及向下延拓应用分析研究,依据实测数据保障条件和积分恒等式适用条件要求,导出了重力异常向上延拓积分模型的分步改化公式,提出了补偿传统改化模型缺陷的修正公式,并将最终的严密改化模型应用于重力异常向下延拓迭代解算。使用超高阶地球位模型EGM2008作为标准位场开展数值计算检验,分别对重力异常向上延拓分步改化模型的计算精度及在向下延拓迭代解算中的应用效果进行了检核评估,验证了采用严密改化模型的必要性和有效性。     

Short Note: On the Relativistic Doppler Effect for Precise Velocity Determination using GPS

The Doppler effect is the apparent shift in frequency of an electromagnetic signal that is received by an observer moving relative to the source of the signal. The Doppler frequency shift relates directly to the relative speed between the receiver and the transmitter, and has thus been widely used in velocity determination. A GPS receiver-satellite pair is in the Earth’s gravity field and GPS signals travel at the speed of light, hence both Einstein’s special and general relativity theories apply. This paper establishes the relationship between a Doppler shift and a user’s ground velocity by taking both the special and general relativistic effects into consideration. A unified Doppler shift model is developed, which accommodates both the classical Doppler effect and the relativistic Doppler effect under special and general relativities. By identifying the relativistic correction terms in the model, a highly accurate GPS Doppler shift observation equation is presented. It is demonstrated that in the GPS “frequency” or “velocity” domain, the relativistic effect from satellite motion changes the receiver-satellite line-of-sight direction, and the measured Doppler shift has correction terms due to the relativistic effects of the receiver potential difference from the geoid, the orbit eccentricity, and the rotation of the Earth.