Fourier transform summation of Legendre series and D-functions

The relation between D- and d-functions, spherical harmonic functions and Legendre functions is reviewed. Dmatrices and irreducible representations of the rotation group O(3) and SU(2) group are briefly reviewed. Two new recursive methods for calculations of D-matrices are presented. Legendre functions are evaluated as part of this scheme. Vector spherical harmonics in the form af generalized spherical harmonics are also included as well as derivatives of the spherical harmonics. The special dmatrices evaluated for argument equal to/2 offer a simple method of calculating the Fourier coefficients of Legendre functions, derivatives of Legendre functions and vector spherical harmonics. Summation of a Legendre series or a full synthesis on the unit sphere of a field can then be performed by transforming the spherical harmonic coefficients to Fourier coefficients and making the summation by an inverse FFT (Fast Fourier Transform). The procedure is general and can also be applied to evaluate derivatives of a field and components of vector and tensor fields.     

A method for computing the coefficients in the product-sum formula of associated Legendre functions

The product of two associated Legendre functions can be represented by a finite series in associated Legendre functions with unique coefficients. In this study a method is proposed to compute the coefficients in this product-sum formula. The method is of recursive nature and is based on the straightforward polynomial form of the associated Legendre function's factor. The method is verified through the computation of integrals of products of two associated Legendre functions over a given interval and the computation of integrals of products of two Legendre polynomials over [0,1]. These coefficients are basically constant and can be used in any future related applications. A table containing the coefficients up to degree 5 is given for ready reference.     

On the computation and approximation of ultra-high-degree spherical harmonic series

Spherical harmonic series, commonly used to represent the Earth’s gravitational field, are now routinely expanded to ultra-high degree (> 2,000), where the computations of the associated Legendre functions exhibit extremely large ranges (thousands of orders) of magnitudes with varying latitude. We show that in the degree-and-order domain, (ℓ,m), of these functions (with full ortho-normalization), their rather stable oscillatory behavior is distinctly separated from a region of very strong attenuation by a simple linear relationship: , where θ is the polar angle. Derivatives and integrals of associated Legendre functions have these same characteristics. This leads to an operational approach to the computation of spherical harmonic series, including derivatives and integrals of such series, that neglects the numerically insignificant functions on the basis of the above empirical relationship and obviates any concern about their broad range of magnitudes in the recursion formulas that are used to compute them. Tests with a simulated gravitational field show that the errors in so doing can be made less than the data noise at all latitudes and up to expansion degree of at least 10,800. Neglecting numerically insignificant terms in the spherical harmonic series also offers a computational savings of at least one third.     

Transformation between surface spherical harmonic expansion of arbitrary high degree and order and double Fourier series on sphere

In order to accelerate the spherical harmonic synthesis and/or analysis of arbitrary function on the unit sphere, we developed a pair of procedures to transform between a truncated spherical harmonic expansion and the corresponding two-dimensional Fourier series. First, we obtained an analytic expression of the sine/cosine series coefficient of the \(4 \pi \) fully normalized associated Legendre function in terms of the rectangle values of the Wigner d function. Then, we elaborated the existing method to transform the coefficients of the surface spherical harmonic expansion to those of the double Fourier series so as to be capable with arbitrary high degree and order. Next, we created a new method to transform inversely a given double Fourier series to the corresponding surface spherical harmonic expansion. The key of the new method is a couple of new recurrence formulas to compute the inverse transformation coefficients: a decreasing-order, fixed-degree, and fixed-wavenumber three-term formula for general terms, and an increasing-degree-and-order and fixed-wavenumber two-term formula for diagonal terms. Meanwhile, the two seed values are analytically prepared. Both of the forward and inverse transformation procedures are confirmed to be sufficiently accurate and applicable to an extremely high degree/order/wavenumber as \(2^{30}\,{\approx }\,10^9\). The developed procedures will be useful not only in the synthesis and analysis of the spherical harmonic expansion of arbitrary high degree and order, but also in the evaluation of the derivatives and integrals of the spherical harmonic expansion.     

Global spherical harmonic computation by two-dimensional Fourier methods

A method is presented for performing global spherical harmonic computation by two-dimensional Fourier transformations. The method goes back to old literature (Schuster 1902) and tackles the problem of non-orthogonality of Legendre-functions, when discretized on an equi-angular grid. Both analysis and synthesis relations are presented, which link the spherical harmonic spectrum to a two-dimensional Fourier spectrum. As an alternative, certain functions of co-latitude are introduced, which are orthogonal to discretized Legendre functions. Several independent Fourier approaches for spherical harmonic computation fit into our general scheme.     

Spherical harmonic modelling to ultra-high degree of Bouguer and isostatic anomalies

The availability of high-resolution global digital elevation data sets has raised a growing interest in the feasibility of obtaining their spherical harmonic representation at matching resolution, and from there in the modelling of induced gravity perturbations. We have therefore estimated spherical Bouguer and Airy isostatic anomalies whose spherical harmonic models are derived from the Earth’s topography harmonic expansion. These spherical anomalies differ from the classical planar ones and may be used in the context of new applications. We succeeded in meeting a number of challenges to build spherical harmonic models with no theoretical limitation on the resolution. A specific algorithm was developed to enable the computation of associated Legendre functions to any degree and order. It was successfully tested up to degree 32,400. All analyses and syntheses were performed, in 64 bits arithmetic and with semi-empirical control of the significant terms to prevent from calculus underflows and overflows, according to IEEE limitations, also in preserving the speed of a specific regular grid processing scheme. Finally, the continuation from the reference ellipsoid’s surface to the Earth’s surface was performed by high-order Taylor expansion with all grids of required partial derivatives being computed in parallel. The main application was the production of a 1′ × 1′ equiangular global Bouguer anomaly grid which was computed by spherical harmonic analysis of the Earth’s topography–bathymetry ETOPO1 data set up to degree and order 10,800, taking into account the precise boundaries and densities of major lakes and inner seas, with their own altitude, polar caps with bedrock information, and land areas below sea level. The harmonic coefficients for each entity were derived by analyzing the corresponding ETOPO1 part, and free surface data when required, at one arc minute resolution. The following approximations were made: the land, ocean and ice cap gravity spherical harmonic coefficients were computed up to the third degree of the altitude, and the harmonics of the other, smaller parts up to the second degree. Their sum constitutes what we call ETOPG1, the Earth’s TOPography derived Gravity model at 1′ resolution (half-wavelength). The EGM2008 gravity field model and ETOPG1 were then used to rigorously compute 1′ × 1′ point values of surface gravity anomalies and disturbances, respectively, worldwide, at the real Earth’s surface, i.e. at the lower limit of the atmosphere. The disturbance grid is the most interesting product of this study and can be used in various contexts. The surface gravity anomaly grid is an accurate product associated with EGM2008 and ETOPO1, but its gravity information contents are those of EGM2008. Our method was validated by comparison with a direct numerical integration approach applied to a test area in Morocco–South of Spain (Kuhn, private communication 2011) and the agreement was satisfactory. Finally isostatic corrections according to the Airy model, but in spherical geometry, with harmonic coefficients derived from the sets of the ETOPO1 different parts, were computed with a uniform depth of compensation of 30?km. The new world Bouguer and isostatic gravity maps and grids here produced will be made available through the Commission for the Geological Map of the World. Since gravity values are those of the EGM2008 model, geophysical interpretation from these products should not be done for spatial scales below 5 arc minutes (half-wavelength).     

Improved convergence rates for the truncation error in gravimetric geoid determination

 When Stokes's integral is used over a spherical cap to compute a gravimetric estimate of the geoid, a truncation error results due to the neglect of gravity data over the remainder of the Earth. Associated with the truncation error is an error kernel defined over these two complementary regions. An important observation is that the rate of decay of the coefficients of the series expansion for the truncation error in terms of Legendre polynomials is determined by the smoothness properties of the error kernel. Previously published deterministic modifications of Stokes's integration kernel involve either a discontinuity in the error kernel or its first derivative at the spherical cap radius. These kernels are generalised and extended by constructing error kernels whose derivatives at the spherical cap radius are continuous up to an arbitrary order. This construction is achieved by smoothly continuing the error kernel function into the spherical cap using a suitable degree polynomial. Accordingly, an improved rate of convergence of the spectral series representation of the truncation error is obtained. Received: 21 April 1998 / Accepted: 4 October 1999     

Recurrence relations for integrals of Associated Legendre functions

Recurrence relations for the evaluation of the integrals of associated Legendre functions over an arbitrary interval within (0°, 90°) have been derived which yield sufficiently accurate results throughout the entire range of their possible applications. These recurrence relations have been used to compute integrals up to degree 100 and similar computations can be carried out without any difficulty up to a degree as high as the memory in a computer permits. The computed values have been tested with independent check formulae, also derived in this work; the corresponding relative errors never exceed 10⁻²³ in magnitude. Contribution from the Earth Physics Branch No. 719     

On computing ellipsoidal harmonics using Jekeli’s renormalization

Gravity data observed on or reduced to the ellipsoid are preferably represented using ellipsoidal harmonics instead of spherical harmonics. Ellipsoidal harmonics, however, are difficult to use in practice because the computation of the associated Legendre functions of the second kind that occur in the ellipsoidal harmonic expansions is not straightforward. Jekeli’s renormalization simplifies the computation of the associated Legendre functions. We extended the direct computation of these functions—as well as that of their ratio—up to the second derivatives and minimized the number of required recurrences by a suitable hypergeometric transformation. Compared with the original Jekeli’s renormalization the associated Legendre differential equation is fulfilled up to much higher degrees and orders for our optimized recurrences. The derived functions were tested by comparing functionals of the gravitational potential computed with both ellipsoidal and spherical harmonic syntheses. As an input, the high resolution global gravity field model EGM2008 was used. The relative agreement we found between the results of ellipsoidal and spherical syntheses is 10^?14, 10^?12 and 10^?8 for the potential and its first and second derivatives, respectively. Using the original renormalization, this agreement is 10^?12, 10^?8 and 10^?5, respectively. In addition, our optimized recurrences require less computation time as the number of required terms for the hypergeometric functions is less.     

The geodetic approximations in the conversion of geoid height to gravity anomaly by Fourier transform

Six sources of error in the use of Fourier methods for the conversion of geoid heights to gravity anomalies are considered. The errors due to spherical approximation are unimportant. The errors due to approximations in Stokes' integral may be eliminated by use of the gravity coating rather than the gravity anomaly. The chord-to-arc error and the truncation error may be reduced by using a reference field. Tapering of the edges of the measurement window reduces the truncation error. The long-wavelength components of the high degree spherical harmonics cause small offsets in the resulting gravity anomalies. The errors due to the plane approximation can be reduced by appropriate choice of map projection and area of integration.     

On the spherical and spheroidal harmonic expansion of the gravitational potential of the topographic masses

This paper is devoted to the spherical and spheroidal harmonic expansion of the gravitational potential of the topographic masses in the most rigorous way. Such an expansion can be used to compute gravimetric topographic effects for geodetic and geophysical applications. It can also be used to augment a global gravity model to a much higher resolution of the gravitational potential of the topography. A formulation for a spherical harmonic expansion is developed without the spherical approximation. Then, formulas for the spheroidal harmonic expansion are derived. For the latter, Legendre’s functions of the first and second kinds with imaginary variable are expanded in Laurent series. They are then scaled into two real power series of the second eccentricity of the reference ellipsoid. Using these series, formulas for computing the spheroidal harmonic coefficients are reduced to surface harmonic analysis. Two numerical examples are presented. The first is a spherical harmonic expansion to degree and order 2700 by taking advantage of existing software. It demonstrates that rigorous spherical harmonic expansion is possible, but the computed potential on the geoid shows noticeable error pattern at Polar Regions due to the downward continuation from the bounding sphere to the geoid. The second numerical example is the spheroidal expansion to degree and order 180 for the exterior space. The power series of the second eccentricity of the reference ellipsoid is truncated at the eighth order leading to omission errors of 25 nm (RMS) for land areas, with extreme values around 0.5 mm to geoid height. The results show that the ellipsoidal correction is 1.65 m (RMS) over land areas, with maximum value of 13.19 m in the Andes. It shows also that the correction resembles the topography closely, implying that the ellipsoidal correction is rich in all frequencies of the gravity field and not only long wavelength as it is commonly assumed.     

球冠谐分析中非整阶Legendre函数的性质及其计算

局部重力场的谱方法是当前重力学的研究方向,该方法的核心问题是如何构造合适的谱函数以及如何对谱函数实施快速、有效的计算。当所研究的区域近似一个球冠时,迂冠谐函数是该区域对应的谱函数,它由非整阶勒让德（Ｌｅｇｅｎｄｒｅ）函数和三角函数组成,显然非整阶勒让德函数的构造和计算是研究球冠谐函数的关键。本文研究了非高小阶勒让德函数的性质和实用计算方法,包括如何对非整阶勒让德函数实施规格化处理。     

全张量重力梯度数据的谱表示方法

在文献「１」的基础上,进一步研究全张量重力梯度数据的全局和局部分量的广义球谐谱表示和轨道根数据表示,并给出了广义球谐函数与球谐函数这间的关系,从理论上得到了全张量重力梯度数据的描述方法和由全张量重力梯度网格数据恢复全球重力位谱系数的基本公式,本文对全张量重力梯度数据的谱表示和谱分析所做的工作,对由重力梯度张量的部分分量恢复全球重力位普系数有一定的参考价值。     

Computation of spherical harmonic coefficients and their error estimates using least-squares collocation

 Equations expressing the covariances between spherical harmonic coefficients and linear functionals applied on the anomalous gravity potential, T, are derived. The functionals are the evaluation functionals, and those associated with first- and second-order derivatives of T. These equations form the basis for the prediction of spherical harmonic coefficients using least-squares collocation (LSC). The equations were implemented in the GRAVSOFT program GEOCOL. Initially, tests using EGM96 were performed using global and regional sets of geoid heights, gravity anomalies and second-order vertical gravity gradients at ground level and at altitude. The global tests confirm that coefficients may be estimated consistently using LSC while the error estimates are much too large for the lower-order coefficients. The validity of an error estimate calculated using LSC with an isotropic covariance function is based on a hypothesis that the coefficients of a specific degree all belong to the same normal distribution. However, the coefficients of lower degree do not fulfil this, and this seems to be the reason for the too-pessimistic error estimates. In order to test this the coefficients of EGM96 were perturbed, so that the pertubations for a specific degree all belonged to a normal distribution with the variance equal to the mean error variance of the coefficients. The pertubations were used to generate residual geoid heights, gravity anomalies and second-order vertical gravity gradients. These data were then used to calculate estimates of the perturbed coefficients as well as error estimates of the quantities, which now have a very good agreement with the errors computed from the simulated observed minus calculated coefficients. Tests with regionally distributed data showed that long-wavelength information is lost, but also that it seems to be recovered for specific coefficients depending on where the data are located. Received: 3 February 2000 / Accepted: 23 October 2000     

椭球谐和球谐系数之间一个简单的转换关系

本文推导的椭球谐系数和球谐系数相互之间转换关系的核心思想是在ε~2量级下利用Legendre函数的正交性,从球谐系数求解的积分表示出发,将积分中的椭球坐标变量与球坐标变量相互转换,从而得出椭球谐系数与球谐系数之间的转换关系。本文导出的转换关系有以下优点:①对于第二类Legendre函数的计算采用Laurent级数表示,使计算第二类Legendre函数更为简单;②保留了ε~2量级下,导出的转换关系相比文献[2]的形式更简单,满足物理大地测量边值问题线性化的要求;③顾及了余纬和归化余纬的区别。     

Efficient propagation of error covariance matrices of gravitational models: application to GRACE and GOCE

We have applied efficient methods for computing variances and covariances of functions of a global gravity field model expanded in spherical harmonics, using the full variance–covariance matrix of the coefficients. Examples are given with recent models derived from GRACE (up to degree and order 150), and with simulated GOCE derived solutions (up to degree and order 200).     

Construction of spherical harmonic series for the potential derivatives of arbitrary orders in the geocentric Earth-fixed reference frame

The derivatives of the Earth gravitational potential are considered in the global Cartesian Earth-fixed reference frame. Spherical harmonic series are constructed for the potential derivatives of the first and second orders on the basis of a general expression of Cunningham (Celest Mech 2:207–216, 1970) for arbitrary order derivatives of a spherical harmonic. A common structure of the series for the potential and its first- and second-order derivatives allows to develop a general procedure for constructing similar series for the potential derivatives of arbitrary orders. The coefficients of the derivatives are defined by means of recurrence relations in which a coefficient of a certain order derivative is a linear function of two coefficients of a preceding order derivative. The coefficients of the second-order derivatives are also presented as explicit functions of three coefficients of the potential. On the basis of the geopotential model EGM2008, the spherical harmonic coefficients are calculated for the first-, second-, and some third-order derivatives of the disturbing potential T, representing the full potential V, after eliminating from it the zero- and first-degree harmonics. The coefficients of two lowest degrees in the series for the derivatives of T are presented. The corresponding degree variances are estimated. The obtained results can be applied for solving various problems of satellite geodesy and celestial mechanics.     

Comparison of undulation difference accuracies using gravity anomalies and gravity disturbances

Errors are considered in the outer zone contribution to oceanic undulation differences as obtained from a set of potential coefficients complete to degree 180. It is assumed that the gravity data of the inner zone (a spherical cap), consisting of either gravity anomalies or gravity disturbances, has negligible error. This implies that error estimates of the total undulation difference are analyzed. If the potential coefficients are derived from a global field of 1°×1° mean anomalies accurate to ε_Δg=10 mgal, then for a cap radius of 10°, the undulation difference error (for separations between 100 km and 2000 km) ranges from 13 cm to 55 cm in the gravity anomaly case and from 6 cm to 36 cm in the gravity disturbance case. If ε_Δg is reduced to 1 mgal, these errors in both cases are less than 10 cm. In the absence of a spherical cap, both cases yield identical error estimates: about 68 cm if ε_Δg=1 mgal (for most separations) and ranging from 93 cm to 160 cm if ε_Δg=10 mgal. Introducing a perfect 30-degree reference field, the latter errors are reduced to about 110 cm for most separations.     

Earth2014: 1 arc-min shape,topography, bedrock and ice-sheet models – Available as gridded data and degree-10,800 spherical harmonics

Since the release of the ETOPO1 global Earth topography model through the US NOAA in 2009, new or significantly improved topographic data sets have become available over Antarctica, Greenland and parts of the oceans. Here, we present a suite of new 1′ (arc-min) models of Earth’s topography, bedrock and ice-sheets constructed as a composite from up-to-date topography models: Earth2014. Our model suite relies on SRTM30_PLUS v9 bathymetry for the base layer, merged with SRTM v4.1 topography over the continents, Bedmap2 over Antarctica and the new Greenland bedrock topography (GBT v3). As such, Earth2014 provides substantially improved information of bedrock and topography over Earth’s major ice sheets, and more recent bathymetric depth data over the oceans, all merged into readily usable global grids. To satisfy multiple applications of global elevation data, Earth2014 provides different representations of Earth’s relief. These are grids of (1) the physical surface, (2) bedrock (Earth’s relief without water and ice masses), (3) bedrock and ice (Earth without water masses), (4) ice sheet thicknesses, (5) rock-equivalent topography (ice and water masses condensed to layers of rock) as mass representation. These models have been transformed into ultra-high degree spherical harmonics, yielding degree 10,800 series expansions of the Earth2014 grids as input for spectral modelling techniques. As further variants, planetary shape models were constructed, providing distances between relief points and the geocenter. The paper describes the input data sets, the development procedures applied, the resulting gridded and spectral representations of Earth2014, external validation results and possible applications. The Earth2014 model suite is freely available via http://ddfe.curtin.edu.au/models/Earth2014/.     

The relation between degree-2160 spectral models of Earth’s gravitational and topographic potential: a guide on global correlation measures and their dependency on approximation effects

Comparisons between high-degree models of the Earth’s topographic and gravitational potential may give insight into the quality and resolution of the source data sets, provide feedback on the modelling techniques and help to better understand the gravity field composition. Degree correlations (cross-correlation coefficients) or reduction rates (quantifying the amount of topographic signal contained in the gravitational potential) are indicators used in a number of contemporary studies. However, depending on the modelling techniques and underlying levels of approximation, the correlation at high degrees may vary significantly, as do the conclusions drawn. The present paper addresses this problem by attempting to provide a guide on global correlation measures with particular emphasis on approximation effects and variants of topographic potential modelling. We investigate and discuss the impact of different effects (e.g., truncation of series expansions of the topographic potential, mass compression, ellipsoidal versus spherical approximation, ellipsoidal harmonic coefficient versus spherical harmonic coefficient (SHC) representation) on correlation measures. Our study demonstrates that the correlation coefficients are realistic only when the model’s harmonic coefficients of a given degree are largely independent of the coefficients of other degrees, permitting degree-wise evaluations. This is the case, e.g., when both models are represented in terms of SHCs and spherical approximation (i.e. spherical arrangement of field-generating masses). Alternatively, a representation in ellipsoidal harmonics can be combined with ellipsoidal approximation. The usual ellipsoidal approximation level (i.e. ellipsoidal mass arrangement) is shown to bias correlation coefficients when SHCs are used. Importantly, gravity models from the International Centre for Global Earth Models (ICGEM) are inherently based on this approximation level. A transformation is presented that enables a transformation of ICGEM geopotential models from ellipsoidal to spherical approximation. The transformation is applied to generate a spherical transform of EGM2008 (sphEGM2008) that can meaningfully be correlated degree-wise with the topographic potential. We exploit this new technique and compare a number of models of topographic potential constituents (e.g., potential implied by land topography, ocean water masses) based on the Earth2014 global relief model and a mass-layer forward modelling technique with sphEGM2008. Different to previous findings, our results show very significant short-scale correlation between Earth’s gravitational potential and the potential generated by Earth’s land topography (correlation +0.92, and 60% of EGM2008 signals are delivered through the forward modelling). Our tests reveal that the potential generated by Earth’s oceans water masses is largely unrelated to the geopotential at short scales, suggesting that altimetry-derived gravity and/or bathymetric data sets are significantly underpowered at 5 arc-min scales. We further decompose the topographic potential into the Bouguer shell and terrain correction and show that they are responsible for about 20 and 25% of EGM2008 short-scale signals, respectively. As a general conclusion, the paper shows the importance of using compatible models in topographic/gravitational potential comparisons and recommends the use of SHCs together with spherical approximation or EHCs with ellipsoidal approximation in order to avoid biases in the correlation measures.