Discussion on some FFT problems to determine the geoid

The aim of this investigation is to study some FFT problems related to the application of FFT to gravity field convolution integrals. And the others, such as the effect of spectral leakage, edge effects, cyclic convolution and effect of padding, are also discussed. A numerical test for these problems is made. A large area of Western China selected for the test is located between 30°N～36°N and 96°E～102°E and includes 1 858 gravity observations on land. The results show that the removal of the bias in the residual gravity anomalies is important to avoid spectral leakage. One hundred percent zero padding is highly recommended for further research of the geoid to remove cyclic convolution errors and edge effects. 1-D FFT is recommended for precise local geoid determination because it does not use kernel approximation.     

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     

Local geoid determination combining gravity disturbances and GPS/levelling: a case study in the Lake Nasser area, Aswan, Egypt

 The use of GPS for height control in an area with existing levelling data requires the determination of a local geoid and the bias between the local levelling datum and the one implicitly defined when computing the local geoid. If only scarse gravity data are available, the heights of new data may be collected rapidly by determining the ellipsoidal height by GPS and not using orthometric heights. Hence the geoid determination has to be based on gravity disturbances contingently combined with gravity anomalies. Furthermore, existing GPS/levelling data may also be used in the geoid determination if a suitable general gravity field modelling method (such as least-squares collocation, LSC) is applied. A comparison has been made in the Aswan Dam area between geoids determined using fast Fourier transform (FFT) with gravity disturbances exclusively and LSC using only the gravity disturbances and the disturbances combined with GPS/levelling data. The EGM96 spherical harmonic model was in all cases used in a remove–restore mode. A total of 198 gravity disturbances spaced approximately 3 km apart were used, as well as 35 GPS/levelling points in the vicinity and on the Aswan Dam. No data on the Nasser Lake were available. This gave difficulties when using FFT, which requires the use of gridded data. When using exclusively the gravity disturbances, the agreement between the GPS/levelling data were 0.71 ± 0.17 m for FFT and 0.63 ± 0.15 for LSC. When combining gravity disturbances and GPS/levelling, the LSC error estimate was ±0.10 m. In the latter case two bias parameters had to be introduced to account for a possible levelling datum difference between the levelling on the dam and that on the adjacent roads. Received: 14 August 2000 / Accepted: 28 February 2001     

Residual geoid and free-air gravity over the indian offshore from ers-1 high resolution altimeter data

In this study, ERS-1 altimeter data over the Indian offshore have been processed for deriving marine geoid and gravity. Processing of altimeter data involves corrections for various atmospheric and oceanographic effects, stacking and averaging of repeat passes, cross-over correction, removal of deeper earth and bathymetric effects, spectral analyses and conversion of geoid into free-air gravity anomaly. Methods for generation of residual geoid and free-air gravity anomaly using high resolution ERS-1 168 day repeat altimeter data were developed. High resolution detailed geoid maps, gravity anomaly and their spectral components have been generated over the Indian offshore using ERS-I altimeter data and ARCGIS system. A number of known megastructures over the study area have been successfully interpreted e.g. Bombay High, Saurastra platform, 90° east ridge etc. from these maps.     

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     

Improving Orthometric Heights Using Precise GPS Measurements

The vertical component obtained from the Global Positioning System (GPS) observations is from the ellipsoid (a mathematical surface), and therefore needs to be converted to the orthometric height, which is from the geoid (represented by the mean sea level). The common practice is to use existing bench marks (around the four corners of a project area and interpolate for the rest of the area), but in many areas bench marks may not be available, in which case an existing geoid undulation is used. Present available global geoid undulation values are not generally as detailed as needed, and in many areas they are not known better than ±1 to ±5 m, because of many limitations. This article explains the difficulties encountered in obtaining precise geoid undulation with some example computations, and proposes a technique of applying corrections to the best available global geoid undulations using detailed free-air gravity anomalies (within a 2° × 2° area) to get relative centimeter accuracy. Several test computations have been performed to decide the optimal block sizes and the effective spherical distances to compute the regional and the local effects of gravity anomalies on geoid undulations by using the Stokes integral. In one test computation a 2° × 2° area was subdivided into smaller surface elements. A difference of 37.34 ± 1.6 cm in geoid undulation was obtained over the same 2° × 2° area when 1° × 1° block sizes were replaced by a combination of 5' × 5' and 1' × 1' subdivision integration elements (block sizes).     

Local geoid computation from gravity using the fast fourier transform technique

Model computations were performed for the study of numerical errors which are interjected into local geoid computations byFFT. The gravity field model was generated through the attractions of granitic prisms derived from actual geology. Changes in sampling interval introduced only0.3 cm variation in geoid heights. Although zero padding alone provided an improvement of more than5 cm in theFFT generated geoid, the combination of spectral windowing (tapering) and padding further reduced numerical errors. For theGPS survey of Franklin County, Ohio, the parameters selected as a result of model computations, allow large reduction in local data requirements while still retaining the centimeter accuracy when tapering and padding is applied.     

Gravity anomalies from satellite altimetry: comparison between computation via geoid heights and via deflections of the vertical

The accumulation of good quality satellite altimetry missions allows us to have a precise geoid with fair resolution and to compute free air gravity anomalies easily by fast Fourier transform (FFT) techniques.In this study we are comparing two methods to get gravity anomalies. The first one is to establish a geoid grid and transform it into anomalies using inverse Stokes formula in the spectral domain via FFT. The second one computes deflection of the vertical grids and transforms them into anomalies.The comparison is made using different data sets: Geosat, ERS-1 and Topex-Poseidon exact repeat misions (ERMs) north of 30°S and Geosat geodetic mission (GM) south of 30°S. The second method which transforms the geoid gradients converted into deflection of the vertical values is much better and the results have been favourably evaluated by comparison with marine gravity data.     

Improvements on existing geopotential coefficients models for precise determination of 1° global geoid and mean gravity anomalies

The concept of an idealised earth having 1° averaged heights over its land surface is introduced as a means to improve upon the existing geopotential coefficient solutions without the use of additional observed data, in order to provide more precise knowledge of the earth’s gravity field in the form of 1° global geoid and 1° mean free-air gravity anomalies especially over the mountainous regions with the visible topography condensed into the actual geoid, first by referring them to the idealised earth and then by reducing the same to the actual earth on applying appropriate corrections for the differences between the two earths.     

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

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

Detailed gravimetric geoid for the North Atlantic

A detailed gravimetric geoid in the North Atlantic Ocean, named DGGNA-77, has been computed, based on a satellite and gravimetry derived earth potential model (consisting in spherical harmonic coefficients up to degree and order 30) and mean free air surface gravity anomalies (35180 1°×1° mean values and 245000 4′×4′ mean values). The long wavelength undulations were computed from the spherical harmonics of the reference potential model and the details were obtained by integrating the residual gravity anomalies through the Stokes formula: from 0 to 5° with the 4′×4′ data, and from 5° to 20° with the 1°×1° data. For computer time reasons the final grid was computed with half a degree spacing only. This grid extends from the Gulf of Mexico to the European and African coasts. Comparisons have been made with Geos 3 altimetry derived geoid heights and with the 5′×5′ gravimetric geoid derived byMarsh andChang [8] in the northwestern part of the Atlantic Ocean, which show a good agreement in most places apart from some tilts which porbably come from the satellite orbit recovery.     

On the inversion of geodetic integrals defined over the sphere using 1-D FFT

An iterative method is presented which performs inversion of integrals defined over the sphere. The method is based on one-dimensional fast Fourier transform (1-D FFT) inversion and is implemented with the projected Landweber technique, which is used to solve constrained least-squares problems reducing the associated 1-D cyclic-convolution error. The results obtained are as precise as the direct matrix inversion approach, but with better computational efficiency. A case study uses the inversion of Hotine’s integral to obtain gravity disturbances from geoid undulations. Numerical convergence is also analyzed and comparisons with respect to the direct matrix inversion method using conjugate gradient (CG) iteration are presented. Like the CG method, the number of iterations needed to get the optimum (i.e., small) error decreases as the measurement noise increases. Nevertheless, for discrete data given over a whole parallel band, the method can be applied directly without implementing the projected Landweber method, since no cyclic convolution error exists.     

Precise geoid determination over Sweden using the Stokes–Helmert method and improved topographic corrections

 Four different implementations of Stokes' formula are employed for the estimation of geoid heights over Sweden: the Vincent and Marsh (1974) model with the high-degree reference gravity field but no kernel modifications; modified Wong and Gore (1969) and Molodenskii et al. (1962) models, which use a high-degree reference gravity field and modification of Stokes' kernel; and a least-squares (LS) spectral weighting proposed by Sj?berg (1991). Classical topographic correction formulae are improved to consider long-wavelength contributions. The effect of a Bouguer shell is also included in the formulae, which is neglected in classical formulae due to planar approximation. The gravimetric geoid is compared with global positioning system (GPS)-levelling-derived geoid heights at 23 Swedish Permanent GPS Network SWEPOS stations distributed over Sweden. The LS method is in best agreement, with a 10.1-cm mean and ±5.5-cm standard deviation in the differences between gravimetric and GPS geoid heights. The gravimetric geoid was also fitted to the GPS-levelling-derived geoid using a four-parameter transformation model. The results after fitting also show the best consistency for the LS method, with the standard deviation of differences reduced to ±1.1 cm. For comparison, the NKG96 geoid yields a 17-cm mean and ±8-cm standard deviation of agreement with the same SWEPOS stations. After four-parameter fitting to the GPS stations, the standard deviation reduces to ±6.1 cm for the NKG96 geoid. It is concluded that the new corrections in this study improve the accuracy of the geoid. The final geoid heights range from 17.22 to 43.62 m with a mean value of 29.01 m. The standard errors of the computed geoid heights, through a simple error propagation of standard errors of mean anomalies, are also computed. They range from ±7.02 to ±13.05 cm. The global root-mean-square error of the LS model is the other estimation of the accuracy of the final geoid, and is computed to be ±28.6 cm. Received: 15 September 1999 / Accepted: 6 November 2000     

Gravity field convolutions without windowing and edge effects

A new set of formulas has been developed for the computation of geoid undulations and terrain corrections by FFT when the input gravity anomalies and heights are mean gridded values. The effects of the analytical and the discrete spectra of kernel functions and that of zero-padding on the computation of geoid undulations and terrain corrections are studied in detail.Numerical examples show that the discrete spectrum is superior to the analytically-defined one. By using the discrete spectrum and 100% zero-padding, the RMS differences are 0.000 m for the FFT geoid undulations and 0.200 to 0.000 mGal for the FFT terrain corrections compared with results obtained by numerical integration.     

全球及局部海洋扰动重力反演的快速解析方法

从经典边值问题理论及球谐函数理论出发,在空域推导获得了由大地水准面高以及垂线偏差计算扰动重力的解析计算公式,为利用卫星测高数据反演海洋扰动重力提供了理论基础。针对全球海洋区域和局部海洋区域的扰动重力反演,在前人已有工作基础上,提出了改进的基于一维FFT的精确快速算法,保证了计算结果与原解析方法完全一致,且计算速度提高约20倍。该算法在提高计算效率的同时避免了由于引入FFT而产生的混叠、边缘效应问题,而且对观测数据的序列长度没有硬性要求,使得应用更加灵活。利用EGM2008地球重力场模型分别生成了2.5'分辨率大地水准面高数据和垂线偏差数据,按照本文提出的改进方法(采用全球积分计算)分别反演获得了全球及局部海洋区域的扰动重力。经比较分析,由大地水准面和垂线偏差分别反演获得的扰动重力其差异在0.8×10^-5 m/s²以内,这说明两种反演方法是基本一致的,但在数据包含系统误差的情况下,由垂线偏差反演扰动重力具有一定优势。     

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

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

A bias-free geodetic boundary value problem approach to height datum unification

A geodetic boundary value problem (GBVP) approach has been formulated which can be used for solving the problem of height datum unification. The developed technique is applied to a test area in Southwest Finland with approximate size of 1.5° × 3° and the bias of the corresponding local height datum (local geoid) with respect to the geoid is computed. For this purpose the bias-free potential difference and gravity difference observations of the test area are used and the offset (bias) of the height datum, i.e., Finnish Height Datum 2000 (N2000) fixed to Normaal Amsterdams Peil (NAP) as origin point, with respect to the geoid is computed. The results of this computation show that potential of the origin point of N2000, i.e., NAP, is (62636857.68 ± 0.5) (m²/s²) and as such is (0.191 ± 0.003) (m) under the geoid defined by W ₀ = 62636855.8 (m²/s²). As the validity test of our methodology, the test area is divided into two parts and the corresponding potential difference and gravity difference observations are introduced into our GBVP separately and the bias of height datums of the two parts are computed with respect to the geoid. Obtaining approximately the same bias values for the height datums of the two parts being part of one height datum with one origin point proves the validity of our approach. Besides, the latter test shows the capability of our methodology for patch-wise application.     

The northern European geoid: a case study on long-wavelength geoid errors

 The long-wavelength geoid errors on large-scale geoid solutions, and the use of modified kernels to mitigate these effects, are studied. The geoid around the Nordic area, from Greenland to the Ural mountains, is considered. The effect of including additional gravity data around the Nordic/Baltic land area, originating from both marine, satellite and ground-based measurements, is studied. It is found that additional data appear to increase the noise level in computations, indicating the presence of systematic errors. Therefore, the Wong–Gore modification to the Stokes kernel is applied. This method of removing lower-order terms in the Stokes kernel appears to improve the geoid. The best fit to the global positioning system (GPS) leveling points is obtained with a degree of modification of approximately 30. In addition to the study of modification errors, the results of different methods of combining satellite altimetry gravity and other gravimetry are presented. They all gave comparable results, at the 6-cm level, when evaluated for the Nordic GPS networks. One dimensional (1-D) and 2-D fast Fourier transform (FFT) methods are also compared. It is shown that even though methods differ by up to 6 cm, the fit to the GPS is essentially the same. A surprising conclusion is that the addition of more data does not always produce a better geoid, illustrating the danger of systematic errors in data. Received: 4 July 2001 / Accepted: 21 February 2002     

ON THE EVALUATION OF DEFLECTIONS OF THE VERTICAL USING FFT TECHNIQUE

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

Updating and computing the geoid using two dimensional fast Hartley transform and fast T transform

In the analyses of 2D real arrays, fast Hartley (FHT), fast T (FTT) and real-valued fast Fourier transforms are generally preferred in lieu of a complex fast Fourier transform due to the advantages of the former with respect to disk storage and computation time. Although the FHT and the FTT in one dimension are identical, they are different in two or more dimensions. Therefore, first, definitions and some properties of both transforms and the related 2D FHT and FTT algorithms are stated. After reviewing the 2D FHT and FTT solutions of Stokes' formula in planar approximation, 2D FHT and FTT methods are developed for geoid updating to incorporate additional gravity anomalies. The methods are applied for a test area which includes a 64×64 grid of 3^′×3^′ point gravity anomalies and geoid heights calculated from point masses. The geoids computed by 2D FHT and FTT are found to be identical. However, the RMS value of the differences between the computed and test geoid is ±15 mm. The numerical simulations indicate that the new methods of geoid updating are practical and accurate with considerable savings on storage requirements. Received: 15 February 1996; Accepted: 22 January 1997