On the evaluation of deflections of the vertical using FFT technique

There exist three types of convolution formulae for the efficient evaluation of gravity field convolution integrals, i.e., the planar 2D convolution, the spherical 2D convolution and the spherical 1D convolution. The largest drawback of both the planar and the spherical 2D FFT methods is that, due to the approximations in the kernel function, only inexact results can be achieved. Apparently, the reason is the meridian convergence at higher latitudes. As the meridians converge, the ??,?λ blocks do not form a rectangular grid, as is assumed in 2D FFT methods. It should be pointed out that the meridian convergence not only leads to an approximation error in the kernel function, but also causes an approximation error during the implementation of 2D FFT in computer. In order to meet the increasing need for precise determination of the vertica deflections, this paper derives a more precise planar 2D FFT formula for the computation of the vertical deflections. After having made a detailed comparison between the planar and the spherical 2D FFT formulae, we find out the main source of errors causing the loss in accuracy by applying the conventional spherical 2D FFT method. And then, a modified spherical 2D FFT formula for the computation of the vertical deflections is developed in this paper. A series of numerical tests have been carried out to illustrate the improvement made upon the old spherical 2D FFT. The second part of this paper is to discuss the influences of the spherical harmonic reference field, the limited capsize, and the singular integral on the computation of the vertical deflections. The results of the vertical deflections over China by applying the spherical 1D FFT formula with different integration radii have been compared to the astro-observed vertical deflections in the South China Sea to obtain a set of optimum deflection computation parameters.     

A formula for computing the gravity disturbance from the second radial derivative of the disturbing potential

 A formula for computing the gravity disturbance and gravity anomaly from the second radial derivative of the disturbing potential is derived in detail using the basic differential equation with spherical approximation in physical geodesy and the modified Poisson integral formula. The derived integral in the space domain, expressed by a spherical geometric quantity, is then converted to a convolution form in the local planar rectangular coordinate system tangent to the geoid at the computing point, and the corresponding spectral formulae of 1-D FFT and 2-D FFT are presented for numerical computation. Received: 27 December 2000 / Accepted: 3 September 2001     

Evaluation of deflections of the vertical on the sphere and the plane: a comparison of FFT techniques

This paper presents a set of efficient formulas to evaluate the deflections of the vertical on the sphere using gridded data. The Vening-Meinesz formula, the topographic indirect effect on the deflections of the vertical as well as the terrain corrections are expressed as both 2D and 1D convolutions on the sphere, and consequently can be evaluated by the 2D and the 1D fast Fourier transform (FFT). When compared with the results obtained from pointwise integration, the use of the 1D FFT gives identical results, and therefore these results were used as control values in this paper. The use of the spherical 2D FFT improves significantly the computational efficiency with little sacrifice of accuracy (0.6^″ rms difference from the 1D FFT results). The planar 2D FFT, which is as efficient as the spherical 2D FFT, gives worse results (1.2^″ rms difference from the 1D FFT results) because of the extra approximations. Received: 27 February 1996; Accepted: 24 January 1997     

关于Stokes公式的球面卷积和平面卷积的注记

讨论了Stokes公式球面卷积和平面卷积形式的近似性和严密性问题,分析了Stokes函数球面卷积形式和平面卷积形式的关系,推导了其间的差值表达式,估算了最大差值及其对计算大地水准面差距的误差影响。同时指出,将顾及Stokes函数全项的平面卷积公式称为严密公式的提法,仅仅是相对仅顾及Stokes函数首项的简单平面卷积公式而言,认为更合理的提法应该是“高精度Stokes平面近似卷积公式”。理论分析表明,球面卷积不可能严格转化为等效的平面卷积。     

Inverse Vening Meinesz formula and deflection-geoid formula: applications to the predictions of gravity and geoid over the South China Sea

Using the spherical harmonic representations of the earth's disturbing potential and its functionals, we derive the inverse Vening Meinesz formula, which converts deflection of the vertical to gravity anomaly using the gradient of the H function. The deflection-geoid formula is also derived that converts deflection to geoidal undulation using the gradient of the C function. The two formulae are implemented by the 1D FFT and the 2D FFT methods. The innermost zone effect is derived. The inverse Vening Meinesz formula is employed to compute gravity anomalies and geoidal undulations over the South China Sea using deflections from Seasat, Geosat, ERS-1 and TOPEX//POSEIDON satellite altimetry. The 1D FFT yields the best result of 9.9-mgal rms difference with the shipborne gravity anomalies. Using the simulated deflections from EGM96, the deflection-geoid formula yields a 4-cm rms difference with the EGM96-generated geoid. The predicted gravity anomalies and geoidal undulations can be used to study the tectonic structure and the ocean circulations of the South China Sea. Received: 7 April 1997 / Accepted: 7 January 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.     

Evaluation of deflections of the vertical using topographic-isostatic and astrogeodetic data for the area of Greece

The evaluation of deflections of the vertical for the area of Greece is attempted using a combination of topographic and astrogeodetic data. Tests carried out in the area bounded by 35°≤ϕ≤42°, 19°≤λ≤27° indicate that an accuracy of ±3″.3 can be obtained in this area for the meridian and prime vertical deflection components when high resolution topographic data in the immediate vicinity of computation points are used, combined with high degree spherical harmonic expansions of the geopotential and isostatic reduction potential. This accuracy is about 25% better than the corresponding topographic-Moho deflection components which are evaluated using topographic and Moho data up to 120 km around each station, without any combination with the spherical harmonic expansion of the geopotential or isostatic reduction potential. The accuracy in both cases is increased to about 2″.6 when the astrogeodetic data available in the area mentioned above are used for the prediction of remaining values. Furthermore the estimation of datum-shift parameters is attempted using least squares collocation.     

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     

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     

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.     

Comparisons of spectral techniques for geoid computations over large regions

The Stokes formula is efficiently evaluated by the one-and two- dimensional (1D, 2D) fast Fourier transform (FFT) technique in the plane and on the sphere in order to obtain precise geoid determinatiover a large area such as Europe. Using a high-pass filtered spherical harmonic reference model (OSU91A truncated to different degrees), gridded gravity anomalies and geoid heights were produced and the anomalies were used as input in the FFT software. Various tests were performed with respect to the different kernel functions used, to the spherical computations in bands, as well as to windowing, edge effects and extent of the area. It is thus demonstrated that, in geoid computations over large regions, the 1D spherical FFT and the 2D multiband spherical FFT in combination with discrete spectra for the kernel functions and 100% zero-padding give better results than those obtained by the other transform techniques. Additionally, numerical tests were carried out at the same test area using the planar fast Hartley transform (FHT) instead of the FFT and the results obtained by the two attractive alternatives were compared regarding the requirements in both computer time and computer memory needed in geoid height computations.A slightly modified version of the paper has been presented at the XX EGS General Assembly, Hamburg, 3–7 April, 1995     

利用垂线偏差计算大地水准面中央区效应的改进方法

为提高利用Molodensky公式反演测高大地水准面中央区效应的精度,视中央区为矩形域,将垂线偏差分量表示成双二次多项式插值形式,引入非奇异变换,推导出了大地水准面的计算公式。垂线偏差理论模型下的分析表明本文导出公式误差为零,而传统公式的误差与纬度以及垂线偏差子午分量与卯酉分量之间的比值有关;以中纬度区域分辨率为2'*2'的垂线偏差数据为背景场进行了实际计算,结果表明在反演计算点本身所在的1个网格对大地水准面的贡献时,传统公式与本文导出公式计算结果差值的最大值达数厘米。本文导出公式可为测高大地水准面的高精度反演提供理论依据。     

Discussion on some FFT problems to determine the geoid

1　ＩｎｔｒｏｄｕｃｔｉｏｎＩｎｔｈｅｍｉｄ＿1 980ｓ,ｔｈｅｆａｓｔＦｏｕｒｉｅｒｔｒａｎｓｆｏｒｍａｔｉｏｎ(ＦＦＴ)ｂｅｇａｎｔｏｆｉｎｄｗｉｄｅｓｐｒｅａｄｕｓｅｉｎｇｅｏｉｄｄｅｔｅｒ ｍｉｎａｔｉｏｎｂｅｃａｕｓｅｏｆｉｔｓｅｆｆｉｃｉｅｎｔｅｖａｌｕａｔｉｏｎｏｆｃｏｎ ｖｏｌｕｔｉｏｎｉｎｔｅｇｒａｌｓ,ｗｈｅｎｃｏｍｐａｒｅｄｔｏｃｌａｓｓｉｃａｌｎｕ ｍｅｒｉｃａｌｉｎｔｅｇｒａｔｉｏｎ .Ｆｏｒｍａｎｙｙｅａｒｓ,ｔｈｅｐｌａｎａｒ,2＿ＤＦＦＴｈａｄｂｅｅｎｕｓｅｄ (Ｓｃｈｗａｒｚ ,1 …     

不同扰动位泛函间的积分变换广义核函数

基于物理大地测量边值问题的解,利用一阶边界算子定义,推导重力异常Δg、单层密度μ、大地水准面高N,垂线偏差ε、扰动重力δg等扰动场元的解。利用球谐函数的正交特性,通过对核函数的算子运算,可以得到上述扰动场元的有关逆变换公式。相对经典物理大地测量公式应用的边界面条件,笔者将含有因子r的对应扰动场元反演关系的公式称为广义积分公式。针对常用的重力异常Δg、大地水准面高N,垂线偏差ε、扰动重力δg计算,重点分析它们之间的变换关系,给出利用某个选定扰动场元计算其他扰动场元的广义积分公式。同时,通过对积分边界面的讨论,分析经典公式与广义积分公式的差异和联系。最后,给出所有外部扰动场元与核函数映射的关系表。     

APPLICATION OF THE HARTLEY TRANSFORM IN THE COMPUTATION OF THE GEOID IN CHINA

This paper presents a method for the computation of the Stokes for-mula using the Fast Hartley Transform（FHT）techniques.The algorithm is mostsuitable for the computation of real sequence transform,while the Fast FourierTransform（FFT）techniques are more suitable for the computaton of complex se-quence transform.A method of spherical coordinate transformation is presented inthis paper.By this method the errors,which are due to the approximate term inthe convolution of Stokes formula,can be effectively eliminated.Some numericaltests are given.By a comparison with both FFT techniques and numerical integra-tion method,the results show that the resulting values of geoidal undulations byFHT techniques are almost the same as by FFT techniques,and the computation-al speed of FHT techniques is about two times faster than that of FFT techniques.     

Spherical panel clustering and its numerical aspects

In modern approximation methods, linear combinations in terms of (space localizing) radial basis functions play an essential role. Areas of application are numerical integration formulae on the unit sphere Ω corresponding to prescribed nodes, spherical spline interpolation and spherical wavelet approximation. The evaluation of such a linear combination is a time-consuming task, since a certain number of summations, multiplications and the calculation of scalar products are required. A generalization of the panel clustering method in a spherical setup is presented. The economy and efficiency of panel clustering are demonstrated for three fields of interest, namely upward continuation of the Earth's gravitational potential, geoid computation by spherical splines and wavelet reconstruction of the gravitational potential. Received: 1 October 1997 / Accepted: 1 April 1998     

Relation between geoidal undulation, deflection of the vertical and vertical gravity gradient revisited

The vertical gradients of gravity anomaly and gravity disturbance can be related to horizontal first derivatives of deflection of the vertical or second derivatives of geoidal undulations. These are simplified relations of which different variations have found application in satellite altimetry with the implicit assumption that the neglected terms—using remove-restore—are sufficiently small. In this paper, the different simplified relations are rigorously connected and the neglected terms are made explicit. The main neglected terms are a curvilinear term that accounts for the difference between second derivatives in a Cartesian system and on a spherical surface, and a small circle term that stems from the difference between second derivatives on a great and small circle. The neglected terms were compared with the dynamic ocean topography (DOT) and the requirements on the GOCE gravity gradients. In addition, the signal root-mean-square (RMS) of the neglected terms and vertical gravity gradient were compared, and the effect of a remove-restore procedure was studied. These analyses show that both neglected terms have the same order of magnitude as the DOT gradient signal and may be above the GOCE requirements, and should be accounted for when combining altimetry derived and GOCE measured gradients. The signal RMS of both neglected terms is in general small when compared with the signal RMS of the vertical gravity gradient, but they may introduce gradient errors above the spherical approximation error. Remove-restore with gravity field models reduces the errors in the vertical gravity gradient, but it appears that errors above the spherical approximation error cannot be avoided at individual locations. When computing the vertical gradient of gravity anomaly from satellite altimeter data using deflections of the vertical, the small circle term is readily available and can be included. The direct computation of the vertical gradient of gravity disturbance from satellite altimeter data is more difficult than the computation of the vertical gradient of gravity anomaly because in the former case the curvilinear term is needed, which is not readily available.     

The solution of the general geodetic boundary value problem by least squares

 A general scheme is given for the solution in a least-squares sense of the geodetic boundary value problem in a spherical, constant-radius approximation, both uniquely and overdetermined, for a large class of observations. The only conditions are that the relation of the observations to the disturbing potential is such that a diagonalization in the spectrum can be found and that the error-covariance function of the observations is isotropic and homogeneous. Most types of observations used in physical geodesy can be adjusted to fit into this approach. Examples are gravity anomalies, deflections of the vertical and the second derivatives of the gravity potential. Received: 3 November 1999 / Accepted: 25 September 2000     

Prediction of vertical deflections from high-degree spherical harmonic synthesis and residual terrain model data

This study demonstrates that in mountainous areas the use of residual terrain model (RTM) data significantly improves the accuracy of vertical deflections obtained from high-degree spherical harmonic synthesis. The new Earth gravitational model EGM2008 is used to compute vertical deflections up to a spherical harmonic degree of 2,160. RTM data can be constructed as difference between high-resolution Shuttle Radar Topography Mission (SRTM) elevation data and the terrain model DTM2006.0 (a spherical harmonic terrain model that complements EGM2008) providing the long-wavelength reference surface. Because these RTM elevations imply most of the gravity field signal beyond spherical harmonic degree of 2,160, they can be used to augment EGM2008 vertical deflection predictions in the very high spherical harmonic degrees. In two mountainous test areas—the German and the Swiss Alps—the combined use of EGM2008 and RTM data was successfully tested at 223 stations with high-precision astrogeodetic vertical deflections from recent zenith camera observations (accuracy of about 0.1 arc seconds) available. The comparison of EGM2008 vertical deflections with the ground-truth astrogeodetic observations shows root mean square (RMS) values (from differences) of 3.5 arc seconds for ξ and 3.2 arc seconds for η, respectively. Using a combination of EGM2008 and RTM data for the prediction of vertical deflections considerably reduces the RMS values to the level of 0.8 arc seconds for both vertical deflection components, which is a significant improvement of about 75%. Density anomalies of the real topography with respect to the residual model topography are one factor limiting the accuracy of the approach. The proposed technique for vertical deflection predictions is based on three publicly available data sets: (1) EGM2008, (2) DTM2006.0 and (3) SRTM elevation data. This allows replication of the approach for improving the accuracy of EGM2008 vertical deflection predictions in regions with a rough topography or for improved validation of EGM2008 and future high-degree spherical harmonic models by means of independent ground truth data.     

On the indirect effect in the Stokes–Helmert method of geoid determination

The classical integral formula for determining the indirect effect in connection with the Stokes–Helmert method is related to a planar approximation of the sea level. A strict integral formula, as well as some approximations to it, are derived. It is concluded that the cap- size truncated integral formulas will suffer from the omission of some long-wavelength contributions, of the order of 50 cm in high mountains for the classical formula. This long-wavelength information can be represented by a set of spherical harmonic coefficients of the topography to, say, degree and order 360. Hence, for practical use, a combination of the classical formula and a set of spherical harmonics is recommended. Received: 10 March 1998 / Accepted: 16 November 1998