Band-limited functions on a bounded spherical domain: the Slepian problem on the sphere

The Slepian problem consists of determining a sequence of functions that constitute an orthonormal basis of a subset of ℝ (or ℝ²) concentrating the maximum information in the subspace of square integrable functions with a band-limited spectrum. The same problem can be stated and solved on the sphere. The relation between the new basis and the ordinary spherical harmonic basis can be explicitly written and numerically studied. The new base functions are orthogonal on both the subspace and the whole sphere. Numerical tests show the applicability of the Slepian approach with regard to solvability and stability in the case of polar data gaps, even in the presence of aliasing. This tool turns out to be a natural solution to the polar gap problem in satellite geodesy. It enables capture of the maximum amount of information from non-polar gravity field missions. Received: 10 June 1998 / Accepted: 20 May 1999     

Construction of Green's function for the Stokes boundary-value problem with ellipsoidal corrections in the boundary condition

Green's function for the boundary-value problem of Stokes's type with ellipsoidal corrections in the boundary condition for anomalous gravity is constructed in a closed form. The `spherical-ellipsoidal' Stokes function describing the effect of two ellipsoidal correcting terms occurring in the boundary condition for anomalous gravity is expressed in O(e <sup>2</sup> <sub>0</sub>)-approximation as a finite sum of elementary functions analytically representing the behaviour of the integration kernel at the singular point ψ=0. We show that the `spherical-ellipsoidal' Stokes function has only a logarithmic singularity in the vicinity of its singular point. The constructed Green function enables us to avoid applying an iterative approach to solve Stokes's boundary-value problem with ellipsoidal correction terms involved in the boundary condition for anomalous gravity. A new Green-function approach is more convenient from the numerical point of view since the solution of the boundary-value problem is determined in one step by computing a Stokes-type integral. The question of the convergence of an iterative scheme recommended so far to solve this boundary-value problem is thus irrelevant. Received: 5 June 1997 / Accepted: 20 February 1998     

Numerical calculation of weakly singular surface integrals

We consider semi-analytical and purely numerical integration methods for weakly singular integrals with point singularities on curved smooth surfaces. The methods can be applied to many practical computations in Geodesy, e.g. terrain corrections, Stokes' and Hotines' integral, surface potentials, and the solution of geodetic boundary value problems using integral equations. Current numerical integration techniques are reviewed. The most important semi-analytical and purely numerical techniques are described. Test calcualtions are done and the techniques are compared as regards accuracy and computational efficiency. Semi-analytical methods, which are based on some regularizing parameter transformations, are superior to purely numerical techniques. The best choice are modified polar coordinates defined in the parameter domain with the singularity as pole. Triangular coordinates show similar performance if carefully tuned. Extrapolation techniques and adaptive subdivision techniques behave poorly as regards accuracy and numerical efficiency. Standard integration techniques, which ignore the singularity, completely fail.     

On the approximation of external gravitational potential with closed systems of (trial) functions

Summary Let S be the (regular) boundary-surface of an exterior regionE _e in Euclidean space ℜ³ (for instance: sphere, ellipsoid, geoid, earth's surface). Denote by {φ_n} a countable, linearly independent system of trial functions (e.g., solid spherical harmonics or certain singularity functions) which are harmonic in some domain containingE _e ∪ S. It is the purpose of this paper to show that the restrictions {ϕ_n} of the functions {φ_n} onS form a closed system in the spaceC (S), i.e. any functionf, defined and continuous onS, can be approximated uniformly by a linear combination of the functions ϕ_n. Consequences of this result are versions of Runge and Keldysh-Lavrentiev theorems adapted to the chosen system {φ_n} and the mathematical justification of the use of trial functions in numerical (especially: collocational) procedures.     

Some cubature formulas for singular integrals in physical geodesy

In terms of a change in variable and Newton-Cotes integration formulas, some cubatures for singular integrals in physical geodesy are formulated in detail. Numerical examples are attached to illustrate their feasibility and their accuracy. Received: 10 July 1995 / Accepted: 4 December 1996     

Construction of Green's function to the external Dirichlet boundary-value problem for the Laplace equation on an ellipsoid of revolution

Green's function to the external Dirichlet boundary-value problem for the Laplace equation with data distributed on an ellipsoid of revolution has been constructed in a closed form. The ellipsoidal Poisson kernel describing the effect of the ellipticity of the boundary on the solution of the investigated boundary-value problem has been expressed as a finite sum of elementary functions which describe analytically the behaviour of the ellipsoidal Poisson kernel at the singular point ψ = 0. We have shown that the degree of singularity of the ellipsoidal Poisson kernel in the vicinity of its singular point is of the same degree as that of the original spherical Poisson kernel. Received: 4 June 1996 / Accepted: 7 April 1997     

On the accuracy of modified Stokes's integration in high-frequency gravimetric geoid determination

 Two numerical techniques are used in recent regional high-frequency geoid computations in Canada: discrete numerical integration and fast Fourier transform. These two techniques have been tested for their numerical accuracy using a synthetic gravity field. The synthetic field was generated by artificially extending the EGM96 spherical harmonic coefficients to degree 2160, which is commensurate with the regular 5^′ geographical grid used in Canada. This field was used to generate self-consistent sets of synthetic gravity anomalies and synthetic geoid heights with different degree variance spectra, which were used as control on the numerical geoid computation techniques. Both the discrete integration and the fast Fourier transform were applied within a 6^∘ spherical cap centered at each computation point. The effect of the gravity data outside the spherical cap was computed using the spheroidal Molodenskij approach. Comparisons of these geoid solutions with the synthetic geoid heights over western Canada indicate that the high-frequency geoid can be computed with an accuracy of approximately 1 cm using the modified Stokes technique, with discrete numerical integration giving a slightly, though not significantly, better result than fast Fourier transform. Received: 2 November 1999 / Accepted: 11 July 2000     

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     

On the accurate numerical evaluation of geodetic convolution integrals

In the numerical evaluation of geodetic convolution integrals, whether by quadrature or discrete/fast Fourier transform (D/FFT) techniques, the integration kernel is sometimes computed at the centre of the discretised grid cells. For singular kernels—a common case in physical geodesy—this approximation produces significant errors near the computation point, where the kernel changes rapidly across the cell. Rigorously, mean kernels across each whole cell are required. We present one numerical and one analytical method capable of providing estimates of mean kernels for convolution integrals. The numerical method is based on Gauss-Legendre quadrature (GLQ) as efficient integration technique. The analytical approach is based on kernel weighting factors, computed in planar approximation close to the computation point, and used to convert non-planar kernels from point to mean representation. A numerical study exemplifies the benefits of using mean kernels in Stokes’s integral. The method is validated using closed-loop tests based on the EGM2008 global gravity model, revealing that using mean kernels instead of point kernels reduces numerical integration errors by a factor of ~5 (at a grid-resolution of 10 arc min). Analytical mean kernel solutions are then derived for 14 other commonly used geodetic convolution integrals: Hotine, Eötvös, Green-Molodensky, tidal displacement, ocean tide loading, deflection-geoid, Vening-Meinesz, inverse Vening-Meinesz, inverse Stokes, inverse Hotine, terrain correction, primary indirect effect, Molodensky’s G1 term and the Poisson integral. We recommend that mean kernels be used to accurately evaluate geodetic convolution integrals, and the two methods presented here are effective and easy to implement.     

Local relationships between the disturbing density, the disturbing potential and the disturbing gravity of the Earth's gravity field

A function having some properties of a wavelet and being harmonic around a given point in R ³ is defined, and three models showing the local relationships between the disturbing density, the disturbing potential and the disturbing gravity are established by using the function as the kernel function of the integrals in the models. The local relationship has two meanings. One is that we can evaluate with a high accuracy the integrals in the models by using mainly high-accuracy and high-resolution data in a local area. The other is that we can obtain a stable solution with high resolution when we invert the integrals in the models because of the rapid decrease of the kernel function of the integrals. As a result, with these models we evaluate one quantity with high resolution, in a band limited by the maximum degree of a set of geopotential coefficients or by the resolution (spacing) of the local data, from another quantity (or quantities) in a local area, and the resulting solution is stable. Received: 6 April 1998 / Accepted: 16 June 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     

Methods of harmonic synthesis for global geopotential models and their first-, second- and third-order gradients

Four widely used algorithms for the computation of the Earth’s gravitational potential and its first-, second- and third-order gradients are examined: the traditional increasing degree recursion in associated Legendre functions and its variant based on the Clenshaw summation, plus the methods of Pines and Cunningham–Metris, which are free from the singularities that distinguish the first two methods at the geographic poles. All four methods are reorganized with the lumped coefficients approach, which in the cases of Pines and Cunningham–Metris requires a complete revision of the algorithms. The characteristics of the four methods are studied and described, and numerical tests are performed to assess and compare their precision, accuracy, and efficiency. In general the performance levels of all four codes exhibit large improvements over previously published versions. From the point of view of numerical precision, away from the geographic poles Clenshaw and Legendre offer an overall better quality. Furthermore, Pines and Cunningham–Metris are affected by an intrinsic loss of precision at the equator and suffer from additional deterioration when the gravity gradients components are rotated into the East-North-Up topocentric reference system. Electronic supplementary material  The online version of this article (doi:) contains supplementary material, which is available to authorized users.     

Wavelet evaluation of the Stokes and Vening Meinesz integrals

The wavelet transform is a powerful tool in evaluating some singular geodetic integrals. Due to its localization properties in both the time (space) and frequency (scale) domains, and because the kernels of some geodetic integrals have singular points and decay smoothly and quickly away from the singularities, many wavelet transform coefficients of the kernels become zeros or negligible, and only a small number of wavelet transform coefficients are significant. It is thus possible to significantly compress the kernels of these integrals on a wavelet basis by neglecting the zero coefficients and the small coefficients below a certain threshold. Therefore, wavelets provide a convenient way of efficiently evaluating these integrals in terms of fast computation and savings of computer memory. A modified algorithm for the wavelet evaluation of Stokes' integral is presented. The same modified algorithm is applied to the evaluation of the Vening Meinesz integral, whose kernel has a stronger singularity than does Stokes' kernel. Numerical examples illustrate the efficiency and accuracy of the wavelet methods.Acknowledgments.The author express their sincere thanks to Dr. Salamonwicz for providing his PhD thesis. E-mail correspondence between the authors and Dr. Barthelmes and Dr. Benciolini contributed to the work. R. Benciolini and the other two anonymous reviewers are thanked for their constructive comments. Support for this research was provided by research grants to Dr. Sideris from the Natural Sciences and Engineering Reserch Council of Canada (NSERC) and the Geomatics for Informed Decisions (GEOIDE) Network of Centres of Excellence. The MATLAB Wavelet Toolbox package was used as the platform to develop the software in this project.     

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.     

Topographic effects by the Stokes–Helmert method of geoid and quasi-geoid determinations

The topographic potential and the direct topographic effect on the geoid are presented as surface integrals, and the direct gravity effect is derived as a rigorous surface integral on the unit sphere. By Taylor-expanding the integrals at sea level with respect to topographic elevation (H) the power series of the effects is derived to arbitrary orders. This study is primarily limited to terms of order H ². The limitations of the various effects in the frequently used planar approximations are demonstrated. In contrast, it is shown that the spherical approximation to power H ² leads to a combined topographic effect on the geoid (direct plus indirect effect) proportional to H˜² (where terms of degrees 0 and 1 are missing) of the order of several metres, while the combined topographic effect on the height anomaly vanishes, implying that current frequent efforts to determine the direct effect to this order are not needed. The last result is in total agreement with Bjerhammar's method in physical geodesy. It is shown that the most frequently applied remove–restore technique of topographic masses in the application of Stokes' formula suffers from significant errors both in the terrain correction C (representing the sum of the direct topographic effect on gravity anomaly and the effect of continuing the anomaly to sea level) and in the term t (mainly representing the indirect effect on the geoidal or quasi-geoidal height). Received: 18 August 1998 / Accepted: 4 October 1999     

Study of a wide-angle laser ranging system for relative positioning of ground-based benchmarks with millimeter accuracy

A wide-angle airborne laser ranging system (WA-ALRS) is developed at the Institut Géographique National (IGN), France, with the aim of providing a new geodesy technique devoted to large (100 km²) networks with a high density (1 km⁻²) of benchmarks. The main objective is to achieve a 1-mm accuracy in relative vertical coordinates from aircraft measurements lasting a few hours. This paper reviews the methodology and analyzes the first experimental data achieved from a specific ground-based experiment. The accuracy in relative coordinate estimates is studied with the help of numerical simulations. It is shown that strong accuracy limitations arise with a small laser beam divergence combined with short range measurements when relatively few simultaneous range data are produced. The accuracy is of a few cm in transverse coordinates and a few mm in radial coordinates. The results from ground-based experimental data are fairly compatible with these predictions. The use of a model for systematic errors in the vehicle trajectory is shown to be necessary to achieve such a high accuracy. This work yields the first complete validation of modelization and methodology of this technique. An accuracy better than 1 mm and a few mm in vertical and horizontal coordinates, respectively, is predicted for aircraft experiments. Received: 19 June 1997 / Accepted: 17 February 1998     

Fast numerical solution of the linearized Molodensky problem

 When standard boundary element methods (BEM) are used in order to solve the linearized vector Molodensky problem we are confronted with two problems: (1) the absence of O(|x|⁻²) terms in the decay condition is not taken into account, since the single-layer ansatz, which is commonly used as representation of the disturbing potential, is of the order O(|x|⁻¹) as x→∞. This implies that the standard theory of Galerkin BEM is not applicable since the injectivity of the integral operator fails; (2) the N×N stiffness matrix is dense, with N typically of the order 10⁵. Without fast algorithms, which provide suitable approximations to the stiffness matrix by a sparse one with O(N(logN) ^s ), s≥0, non-zero elements, high-resolution global gravity field recovery is not feasible. Solutions to both problems are proposed. (1) A proper variational formulation taking the decay condition into account is based on some closed subspace of co-dimension 3 of the space of square integrable functions on the boundary surface. Instead of imposing the constraints directly on the boundary element trial space, they are incorporated into a variational formulation by penalization with a Lagrange multiplier. The conforming discretization yields an augmented linear system of equations of dimension N+3×N+3. The penalty term guarantees the well-posedness of the problem, and gives precise information about the incompatibility of the data. (2) Since the upper left submatrix of dimension N×N of the augmented system is the stiffness matrix of the standard BEM, the approach allows all techniques to be used to generate sparse approximations to the stiffness matrix, such as wavelets, fast multipole methods, panel clustering etc., without any modification. A combination of panel clustering and fast multipole method is used in order to solve the augmented linear system of equations in O(N) operations. The method is based on an approximation of the kernel function of the integral operator by a degenerate kernel in the far field, which is provided by a multipole expansion of the kernel function. Numerical experiments show that the fast algorithm is superior to the standard BEM algorithm in terms of CPU time by about three orders of magnitude for N=65 538 unknowns. Similar holds for the storage requirements. About 30 iterations are necessary in order to solve the linear system of equations using the generalized minimum residual method (GMRES). The number of iterations is almost independent of the number of unknowns, which indicates good conditioning of the system matrix. Received: 16 October 1999 / Accepted: 28 February 2001     

On non-linear low–low SST observation equations for the determination of the geopotential based on an analytical solution

In satellite data analysis, one big advantage of analytical orbit integration, which cannot be overestimated, is missed in the numerical integration approach: spectral analysis or the lumped coefficient concept may be used not only to design efficient algorithms but overall for much better insight into the force-field determination problem. The lumped coefficient concept, considered from a practical point of view, consists of the separation of the observation equation matrix A=BT into the product of two matrices. The matrix T is a very sparse matrix separating into small block-diagonal matrices connecting the harmonic coefficients with the lumped coefficients. The lumped coefficients are nothing other than the amplitudes of trigonometric functions depending on three angular orbital variables; therefore, the matrix N=B ^T B will become for a sufficient length of a data set a diagonal dominant matrix, in the case of an unlimited data string length a strictly diagonal one. Using an analytical solution of high order, the non-linear observation equations for low–low SST range data can be transformed into a form to allow the application of the lumped concept. They are presented here for a second-order solution together with an outline of how to proceed with data analysis in the spectral domain in such a case. The dynamic model presented here provides not only a practical algorithm for the parameter determination but also a simple method for an investigation of some fundamental questions, such as the determination of the range of the subset of geopotential coefficients which can be properly determined by means of SST techniques or the definition of an optimal orbital configuration for particular SST missions. Numerical results have already been obtained and will be published elsewhere. Received: 15 January 1999 / Accepted: 30 November 1999     

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     

计算Legendre函数导数的非奇异方法

引力场关于经度和纬度方向的梯度在两极附近会产生奇异性现象,这将会给诸如重力场和静态洋流探索(GOCE,Gravity field and stesdy-state Oceam Circulation Explorer)数据处理等引力场的研究工作带来诸多不便和困难。这里首先分析了该奇异性产生的原因,即目前采用的球坐标系自身在两极处是奇异的;然后利用Legendre函数的性质推导了一组不含任何奇异性的计算引力场梯度的计算公式;最后与常用的迭代方法进行了实例计算比较,结果表明所导出的公式不仅计算精度大大提高,而且计算用时也不会增加。