Formation of surface spherical harmonic normal matrices and application to high-degree geopotential modeling

 The structure of normal matrices occurring in the problem of weighted least-squares spherical harmonic analysis of measurements scattered on a sphere with random noises is investigated. Efficient algorithms for the formation of the normal matrices are derived using fundamental relations inherent to the products of two surface spherical harmonic functions. The whole elements of a normal matrix complete to spherical harmonic degree L are recursively obtained from its first row or first column extended to degree 2L with only O(L ⁴) computational operations. Applications of the algorithms to the formation of surface normal matrices from geoid undulations and surface gravity anomalies are discussed in connection with the high-degree geopotential modeling. Received: 22 March 1999 / Accepted: 23 December 1999     

The direct topographical correction in gravimetric geoid determination by the Stokes–Helmert method

 The direct topographical correction is composed of both local effects and long-wavelength contributions. This implies that the classical integral formula for determining the direct effect may have some numerical problems in representing these different signals. On the other hand, a representation by a set of harmonic coefficients of the topography to, say, degree and order 360 will omit significant short-wavelength signals. A new formula is derived by combining the classical formula and a set of spherical harmonics. Finally, the results of this solution are compared with the Moritz topographical correction in a test area. Received: 27 July 1998 / Accepted: 29 March 2000     

SHORT NOTE: Extending simplified high-degree synthesis methods to second latitudinal derivatives of geopotential

 Simplified techniques for high-degree spherical harmonic synthesis are extended to include gravitational potential second derivatives with respect to latitude. Received: 23 July 2001 / Accepted: 12 April 2002 Acknowledgement. The authors would like to thank Christian Tscherning for recommending Laplace's equation as an accuracy test. Our use of Legendre's differential equation, as the most direct means for extending our simplified synthesis methods to second-order derivatives, was a direct result of this suggestion. Correspondence to: S. A. Holmes     

The temporal variation of the spherical and Cartesian multipoles of the gravity field: the generalized MacCullagh representation

 The Cartesian moments of the mass density of a gravitating body and the spherical harmonic coefficients of its gravitational field are related in a peculiar way. In particular, the products of inertia can be expressed by the spherical harmonic coefficients of the gravitational potential as was derived by MacCullagh for a rigid body. Here the MacCullagh formulae are extended to a deformable body which is restricted to radial symmetry in order to apply the Love–Shida hypothesis. The mass conservation law allows a representation of the incremental mass density by the respective excitation function. A representation of an arbitrary Cartesian monome is always possible by sums of solid spherical harmonics multiplied by powers of the radius. Introducing these representations into the definition of the Cartesian moments, an extension of the MacCullagh formulae is obtained. In particular, for excitation functions with a vanishing harmonic coefficient of degree zero, the (diagonal) incremental moments of inertia also can be represented by the excitation coefficients. Four types of excitation functions are considered, namely: (1) tidal excitation; (2) loading potential; (3) centrifugal potential; and (4) transverse surface stress. One application of the results could be model computation of the length-of-day variations and polar motion, which depend on the moments of inertia. Received: 27 July 1999 / Accepted: 24 May 2000     

Two different methods of geoidal height determinations using a spherical harmonic representation of the geopotential, topographic corrections and the height anomaly–geoidal height difference

 It is suggested that a spherical harmonic representation of the geoidal heights using global Earth gravity models (EGM) might be accurate enough for many applications, although we know that some short-wavelength signals are missing in a potential coefficient model. A `direct' method of geoidal height determination from a global Earth gravity model coefficient alone and an `indirect' approach of geoidal height determination through height anomaly computed from a global gravity model are investigated. In both methods, suitable correction terms are applied. The results of computations in two test areas show that the direct and indirect approaches of geoid height determination yield good agreement with the classical gravimetric geoidal heights which are determined from Stokes' formula. Surprisingly, the results of the indirect method of geoidal height determination yield better agreement with the global positioning system (GPS)-levelling derived geoid heights, which are used to demonstrate such improvements, than the results of gravimetric geoid heights at to the same GPS stations. It has been demonstrated that the application of correction terms in both methods improves the agreement of geoidal heights at GPS-levelling stations. It is also found that the correction terms in the direct method of geoidal height determination are mostly similar to the correction terms used for the indirect determination of geoidal heights from height anomalies. Received: 26 July 2001 / Accepted: 21 February 2002     

Short-arc analysis of intersatellite tracking data in a gravity mapping mission

 A technique for the analysis of low–low intersatellite range-rate data in a gravity mapping mission is explored. The technique is based on standard tracking data analysis for orbit determination but uses a spherical coordinate representation of the 12 epoch state parameters describing the baseline between the two satellites. This representation of the state parameters is exploited to allow the intersatellite range-rate analysis to benefit from information provided by other tracking data types without large simultaneous multiple-data-type solutions. The technique appears especially valuable for estimating gravity from short arcs (e.g. less than 15 minutes) of data. Gravity recovery simulations which use short arcs are compared with those using arcs a day in length. For a high-inclination orbit, the short-arc analysis recovers low-order gravity coefficients remarkably well, although higher-order terms, especially sectorial terms, are less accurate. Simulations suggest that either long or short arcs of the Gravity Recovery and Climate Experiment (GRACE) data are likely to improve parts of the geopotential spectrum by orders of magnitude. Received: 26 June 2001 / Accepted: 21 January 2002     

How to handle topography in practical geoid determination: three examples

 Three different methods of handling topography in geoid determination were investigated. The first two methods employ the residual terrain model (RTM) remove–restore technique, yielding the quasigeoid, whereas the third method uses the classical Helmert condensation method, yielding the geoid. All three methods were used with the geopotential model Earth Gravity Model (1996) (EGM96) as a reference, and the results were compared to precise global positioning system (GPS) levelling networks in Scandinavia. An investigation of the Helmert method, focusing on the different types of indirect effects and their effects on the geoid, was also carried out. The three different methods used produce almost identical results at the 5-cm level, when compared to the GPS levelling networks. However, small systematic differences existed. Received: 18 March 1999 / Accepted: 21 March 2000     

New analytical and numerical approaches for geopotential modeling

 The standard analytical approach which is applied for constructing geopotential models OSU86 and earlier ones, is based on reducing the boundary value equation to a sphere enveloping the Earth and then solving it directly with respect to the potential coefficients _n,m . In an alternative procedure, developed by Jekeli and used for constructing the models OSU91 and EGM96, at first an ellipsoidal harmonic series is developed for the geopotential and then its coefficients _n,m ^e are transformed to the unknown _n,m . The second solution is more exact, but much more complicated. The standard procedure is modified and a new simple integral formula is derived for evaluating the potential coefficients. The efficiency of the standard and new procedures is studied numerically. In these solutions the same input data are used as for constructing high-degree parts of the EGM96 models. From two sets of _n,m (n≤360,|m|≤n), derived by the standard and new approaches, different spectral characteristics of the gravity anomaly and the geoid undulation are estimated and then compared with similar characteristics evaluated by Jekeli's approach (`etalon' solution). The new solution appears to be very close to Jekeli's, as opposed to the standard solution. The discrepancies between all the characteristics of the new and `etalon' solutions are smaller than the corresponding discrepancies between two versions of the final geopotential model EGM96, one of them (HDM190) constructed by the block-diagonal least squares (LS) adjustment and the other one (V068) by using Jekeli's approach. On the basis of the derived analytical solution a new simple mathematical model is developed to apply the LS technique for evaluating geopotential coefficients. Received: 12 December 2000 / Accepted: 21 June 2001     

A formula for gravimetric terrain corrections using powers of topographic height

Journal of Geodesy - ?The application of Stokes' formula to create geoid undulations requires no masses outside the geoid. However, due to the existence of the topography, terrain...     

Regularization of geopotential determination from satellite data by variance components

 Different types of present or future satellite data have to be combined by applying appropriate weighting for the determination of the gravity field of the Earth, for instance GPS observations for CHAMP with satellite to satellite tracking for the coming mission GRACE as well as gradiometer measurements for GOCE. In addition, the estimate of the geopotential has to be smoothed or regularized because of the inversion problem. It is proposed to solve these two tasks by Bayesian inference on variance components. The estimates of the variance components are computed by a stochastic estimator of the traces of matrices connected with the inverse of the matrix of normal equations, thus leading to a new method for determining variance components for large linear systems. The posterior density function for the variance components, weighting factors and regularization parameters are given in order to compute the confidence intervals for these quantities. Test computations with simulated gradiometer observations for GOCE and satellite to satellite tracking for GRACE show the validity of the approach. Received: 5 June 2001 / Accepted: 28 November 2001     

Spherical harmonic analysis and synthesis for global multiresolution applications

 On the Earth and in its neighborhood, spherical harmonic analysis and synthesis are standard mathematical procedures for scalar, vector and tensor fields. However, with the advent of multiresolution applications, additional considerations about convolution filtering with decimation and dilation are required. As global applications often imply discrete observations on regular grids, computational challenges arise and conflicting claims about spherical harmonic transforms have recently appeared in the literature. Following an overview of general multiresolution analysis and synthesis, spherical harmonic transforms are discussed for discrete global computations. For the necessary multi-rate filtering operations, spherical convolutions along with decimations and dilations are discussed, with practical examples of applications. Concluding remarks are then included for general applications, with some discussion of the computational complexity involved and the ongoing investigations in research centers. Received: 13 November 2000 / Accepted: 12 June 2001     

Assessment of three numerical solution strategies for gravity field recovery from GOCE satellite gravity gradiometry implemented on a parallel platform

 The recovery of a full set of gravity field parameters from satellite gravity gradiometry (SGG) is a huge numerical and computational task. In practice, parallel computing has to be applied to estimate the more than 90 000 harmonic coefficients parameterizing the Earth's gravity field up to a maximum spherical harmonic degree of 300. Three independent solution strategies (preconditioned conjugate gradient method, semi-analytic approach, and distributed non-approximative adjustment), which are based on different concepts, are assessed and compared both theoretically and on the basis of a realistic-as-possible numerical simulation regarding the accuracy of the results, as well as the computational effort. Special concern is given to the correct treatment of the coloured noise characteristics of the gradiometer. The numerical simulations show that the three methods deliver nearly identical results—even in the case of large data gaps in the observation time series. The newly proposed distributed non-approximative adjustment approach, which is the only one of the three methods that solves the inverse problem in a strict sense, also turns out to be a feasible method for practical applications. Received: 17 December 2001 / Accepted: 17 July 2002 Acknowledgments. We would like to thank Prof. W.-D. Schuh, Institute of Theoretical Geodesy, University of Bonn, for providing us with the serial version of the PCGMA algorithm, which forms the basis for the parallel PCGMA package developed at our institute. This study was partially performed in the course of the GOCE project `From E?tv?s to mGal+', funded by the European Space Agency (ESA) under contract No. 14287/00/NL/DC. Correspondence to: R. Pail     

The gravitational potential and its derivatives for the prism

 As a simple building block, the right rectangular parallelepiped (prism) has an important role mostly in local gravity field modelling studies when the so called flat-Earth approximation is sufficient. Its primary (methodological) advantage follows from the simplicity of the rigorous and consistent analytical forms describing the different gravitation-related quantities. The analytical forms provide numerical values for these quantities which satisfy the functional connections existing between these quantities at the level of numerical precision applied. Closed expressions for the gravitational potential of the prism and its derivatives (up to the third order) are listed for easy reference. Received: 18 August 1999 / Accepted: 15 June 2000     

Green's function solution to spherical gradiometric boundary-value problems

 Three independent gradiometric boundary-value problems (BVPs) with three types of gradiometric data, {Γ _rr }, {Γ _{r θ},Γ _{r λ}} and {Γ_θθ−Γ_λλ,Γ_θλ}, prescribed on a sphere are solved to determine the gravitational potential on and outside the sphere. The existence and uniqueness conditions on the solutions are formulated showing that the zero- and the first-degree spherical harmonics are to be removed from {Γ _{r θ},Γ _{r λ}} and {Γ_θθ−Γ_λλ,Γ_θλ}, respectively. The solutions to the gradiometric BVPs are presented in terms of Green's functions, which are expressed in both spectral and closed spatial forms. The logarithmic singularity of the Green's function at the point ψ=0 is investigated for the component Γ _rr . The other two Green's functions are finite at this point. Comparisons to the paper by van Gelderen and Rummel [Journal of Geodesy (2001) 75: 1–11] show that the presented solution refines the former solution. Received: 3 October 2001 / Accepted: 4 October 2002     

Geoid, topography, and the Bouguer plate or shell

 Topography plays an important role in solving many geodetic and geophysical problems. In the evaluation of a topographical effect, a planar model, a spherical model or an even more sophisticated model can be used. In most applications, the planar model is considered appropriate: recall the evaluation of gravity reductions of the free-air, Poincaré–Prey or Bouguer kind. For some applications, such as the evaluation of topographical effects in gravimetric geoid computations, it is preferable or even necessary to use at least the spherical model of topography. In modelling the topographical effect, the bulk of the effect comes from the Bouguer plate, in the case of the planar model, or from the Bouguer shell, in the case of the spherical model. The difference between the effects of the Bouguer plate and the Bouguer shell is studied, while the effect of the rest of topography, the terrain, is discussed elsewhere. It is argued that the classical Bouguer plate gravity reduction should be considered as a mathematical construction with unclear physical meaning. It is shown that if the reduction is understood to be reducing observed gravity onto the geoid through the Bouguer plate/shell then both models give practically identical answers, as associated with Poincaré's and Prey's work. It is shown why only the spherical model should be used in the evaluation of topographical effects in the Stokes–Helmert solution of Stokes' boundary-value problem. The reason for this is that the Bouguer plate model does not allow for a physically acceptable condensation scheme for the topography. Received: 24 December 1999 / Accepted: 11 December 2000     

Geoid determination with density hypotheses from isostatic models and geological information

 Geoid determination by Stokes's formula requires a complete knowledge of the topographical mass density distribution in order to perform gravity reductions to the geoid boundary. However, deeper masses are also of interest, in order to produce a smooth field of gravity anomalies which will improve results from interpolation procedures. Until now, in most cases a constant mass density has been considered, which is a very rough approximation of reality. The influence on the geoid height coming from different mass density hypotheses given by the isostatic models of Pratt/Hayford, Airy/Heiskanen and Vening Meinesz is studied. Apart from a constant mass density value, additional density information deduced from geological maps and thick sedimentary layers is considered. An overview of how mass density distributions act within Stokes's theory is given. The isostatic models are considered in spherical and planar approximation, as well as with constant and lateral variable mass density of the topographical and deeper masses. Numerical results in a test area in south-west Germany show that the differences in the geoid height due to different density hypotheses can reach a magnitude of more than 1 decimetre, which is not negligible in a precise geoid determination with centimetre accuracy. Received: 7 January 2002 / Accepted: 20 September 2002 M. Kuhn now at: Western Australian Centre for Geodesy, Curtin University of Technology, GPO Box U1987, Perth, WA 6845, Australia Acknowledgements. The author would gratefully thank Prof. Dr.-Ing. B. Heck, who was the supervisor of my PhD thesis, and the second examiner Prof. Dr.-Ing. K.H. Ilk, as well as all other colleagues for their support of this work. Particular thanks go to the Landesvermessungsamt Baden–Württemberg (Survey Department of Baden–Württemberg), Bureau Gravimetrique International (BGI, France) for providing the gravity data and the Geologisches Landesamt Baden–Württemberg (Geological Department of Baden–Württemberg) for providing data and maps of the sediment layers within the Rhine Valley. Grateful thanks goes to Prof. W.E. Featherstone and the reviewers Prof. S.D. Pagiatakis, Dr. U. Marti as well as an unknown reviewer for their helpful comments on this paper.     

Some modifications of Stokes' formula that account for truncation and potential coefficient errors

 Stokes' formula from 1849 is still the basis for the gravimetric determination of the geoid. The modification of the formula, originating with Molodensky, aims at reducing the truncation error outside a spherical cap of integration. This goal is still prevalent among various modifications. In contrast to these approaches, some least-squares types of modification that aim at reducing the truncation error, as well as the error stemming from the potential coefficients, are demonstrated. The least-squares estimators are provided in the two cases that (1) Stokes' kernel is a priori modified (e.g. according to Molodensky's approach) and (2) Stokes' kernel is optimally modified to minimize the global mean square error. Meissl-type modifications are also studied. In addition, the use of a higher than second-degree reference field versus the original (Pizzetti-type) reference field is discussed, and it is concluded that the former choice of reference field implies increased computer labour to achieve the same result as with the original reference field. Received: 14 December 1998 / Accepted: 4 October 1999     

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     

Spatially restricted data distributions on the sphere: the method of orthonormalized functions and applications

 In many geoscientific applications data are irregularly distributed and not globally available, e.g. caps around the poles which are uncovered due to non-polar satellite orbits, or signals being defined solely on bounded regions on the globe. Starting from a sequence of base functions with global support, which in the present case is composed of spherical harmonics being initially non-orthogonal on a bounded subdomain, a set of functions is generated that constitutes an orthonormal basis. Different approaches to realize this transformation are studied and compared with respect to numerical stability and computational effort, and the corresponding effects on the coefficient recovery are investigated. A number of synthetic tests demonstrate the applicability, the benefit, but also the limitations, of this method. Received: 24 March 2000 / Accepted: 9 October 2000     

The ERS-2 altimetric bias and gravity field enhancement using dual crossovers between ERS and TOPEX/Poseidon

 Dual satellite crossovers (DXO) between the two European Remote Sensing satellites ERS-1 and ERS-2 and TOPEX/Poseidon are used to (1) refine the Earth's gravity field and (2) extend the study of the ERS-2 altimetric range stability to cover the first four years of its operation. The enhanced gravity field model, AGM-98, is validated by several methodologies and will be shown to provide, in particular, low geographically correlated orbital error for ERS-2. For the ERS-2 altimetric range study, TOPEX/Poseidon is first calibrated through comparison against in situ tide gauge data. A time series of the ERS-2 altimeter bias has been recovered along with other geophysical correction terms using tables for bias jumps in the range measurements at the single point target response (SPTR) events. On utilising the original version of the SPTR tables the overall bias drift is seen to be 2.6±1.0 mm/yr with an RMS of fit of 12.2 mm but with discontinuities at the centimetre level at the SPTR events. On utilising the recently released revised tables, SPTR2000, the drift is better defined at 2.4±0.6 mm/yr with the RMS of fit reduced to 3.7 mm. Investigations identify the sea-state bias as a source of error with corrections affecting the overall drift by close to 1.2 mm/yr. Received: 25 May 2000 / Accepted: 24 January 2001