Evaluation of terrain corrections through FFT and classical integration in two selected areas of the Andes and their impact on geoidal heights

Avaliação das correções do terreno através da FFT e da integração clássica em duas áreas selecionadas dos Andes e do seu impacto sobre a altura geoidal

ARTIGOS

Maria Eugenia Gomez^{I, II}; Diego Ruben Bagu^I; Daniel Del Cogliano^{I, II}; Raul A. Perdomo^{I, II}

^IDto. de Astrometría. Facultad de Ciencias Astronómicas y Geofísicas. Universidad Nacional de La Plata (UNLP) - La Plata - Argentina E-mails: {megomez; dbagu; daniel}@fcaglp.unlp.edu.ar; [e-mail]

^IIConsejo Nacional de Investigaciones Científicas y Técnicas - Buenos Aires - Argentina E-mails: {megomez; dbagu; daniel}@fcaglp.unlp.edu.ar; [e-mail]

ABSTRACT

As part of the regional geoid modeling, terrain corrections were computed in Tierra del Fuego island and in the west side of the province of Mendoza. The first place is located in the southernmost region of Argentina and Mendoza is in the center-west of this country. Considering both study areas, elevations range from 0 m to 6500 m. The classical integration of prism contribution and the 2-D FFT technique were used to estimate terrain corrections. This study aims at discussing the results obtained by both approaches and their applicability considering their advantages and disadvantages according to the regions under investigation. The analysis allowed us to conclude that classical integration has a better performance than FFT methods, especially in the highest regions where terrain corrections can be overestimated in more than 20 mGals by FFT. Both techniques described show similar results in flat areas. Finally, the effect that the error of terrain corrections computation has on geoidal heights is also discussed and numerically tested. It is proved that an error in gravity anomalies of 20 mGals may cause up to 2 m geoid height error.

Keywords: Terrain Corrections; Classical Integration; Tierra del Fuego; Mendoza.

RESUMO

Como parte da modelagem regional do geoide, as correções do terreno foram calculados na ilha da Tierra del Fuego e no lado ocidental da província de Mendoza. O primeiro lugar está localizado na região sul da Argentina e Mendoza está na parte centro-oeste do país. Considerando as duas áreas de estudo, as elevações variam de 0 m a 6500 m. A integração clássica de contribuição por prismas e a técnica 2-D FFT foram utilizadas pra estimar as correções do terreno. Este trabalho tem como objetivo discutir os resultados obtidos pelas duas abordagens e sua aplicabilidade, considerando suas vantagens e suas desvantagens de acordo com as regiões de estudo. A análise permitiu concluir que a integração clássica tem um desempenho melhor do que os métodos da FFT, especialmente nas regiões mais elevadas onde as correções do terreno podem ser superestimadas em mais de 20 mGals pela técnica do FFT. Ambos os métodos mostram resultados semelhantes em áreas planas. Finalmente, o efeito que o erro do cálculo nas correções do terreno tem nas alturas geoidais também é discutido e numericamente testado. Mostra-se que um erro de 20 mGals em anomalias da gravidade pode causar um erro de até 2 m nas alturas geoidais.

Palavras-chave: Correções do Terreno; Integração Clássica; Tierra del Fuego; Mendoza.

1. INTRODUCTION

Terrain corrections (TCR) constitute an important part in the reduction scheme of gravimetric data, particularly in mountainous regions. They are a useful tool for smoothing gravity anomalies and their application allows the use of interpolation techniques in order to calculate a geoid model.

There is no single way to determine TCR. Broadly speaking, it is possible to distinguish two cases: methods which involve the use of prisms and estimate TCR through a classical integration and those employing FFT (Fast Fourier Transform) techniques.

Classical integration has the advantage of being straightforward to interpret while FFT algorithms are faster. The estimation of TCR can take hours or even days of processing in the case of the classical integration due to how the gravitational contribution of each prism is integrated.

The speed of computation depends mainly on the resolution of the Digital Elevation Model (DEM) employed, its extension and the capacity of the CPU processor.

In practice, the FFT is often given preference to the classical approach due to its higher integration speed. Regarding the computed values, the differences between both techniques are negligible in flat regions. However, in high mountainous areas they can exceed tenths of mGals. The reason for these discrepancies relies mainly on the steep topographic gradients which sometimes exceed 45° and cause the Fourier series to fail to converge (KIRBY and FEATHERSTONE, 2001). It must be clearly stated that this kind of topography can cause any method to give wrong results due to the rough terrain or the inadequate representation of the surface relief by the DEM.

In order to apply FFT despite the convergence problem, usually the resolution of the DEM is decreased leading to smaller topographic gradients. This smoothing can imply a resolution of 1' or less and can be considered a safe strategy in regions of moderate gradients. In mountainous areas, the denser the DEM, the better the representation of the highest frequencies. It is then necessary to reach a compromise between the best representation of the topography and the applicability of the TCR computation methods.

Two GRAVSOFT routines, namely TC and TCFOUR (TSCHERNING et al., 1992) were applied to quantify the effect of using classical integration or FFT for terrain correction estimations.

One focus lies on Tierra del Fuego island. In this region the elevation ranges from 0 to 2600 m. The second test area is situated in the province of Mendoza. Here, the elevations range between 200 and 6500 m. Both regions represent parts of the Andean Cordilleras, one of the highest mountain ranges on Earth (Figure 1).

The following sections are devoted to the explanation of each method and their results. The practical relevance of the results obtained is demonstrated by a numerical estimation of the effect of terrain correction errors on geoid determinations for the area of Tierra del Fuego.

1.1 Terrain Corrections Estimated by Classical Integration (TC)

TC estimates the effect of topography with respect to the Bouguer plate in a conventional manner. It integrates the gravitational effect produced by each prism constructed from the DEM (FORSBERG, 1984).

For a prism bounded by x= xp-x, y=yp-y, z=hp-z, the gravitational effect of terrain at a point P is given by equation (1):

where is the density, G is the gravitational constant and r is the spherical geocentric radius.

(xp,yp,hp) and (x,y,z) represent the Cartesian coordinates of the computation point of height hp and the corresponding coordinates of mass elements respectively.

The total TCR is given by the sum of the terrain effect from all prisms.

Due to the instabilities of these formulas over large distances, McMillan (1958) formulas are used in such case (equation 2).

and depend on the distance between the mass element and its distance to point P. For more details, see McMillan (1958).

One way to speed up the calculation is to use two DEM grids of different resolution: one high-resolution DEM covering the region close to the point of calculation and a less dense grid to take into account the gravitational effect caused by distant structures. For this purpose, two different radii must be defined.

When the proximity to the point of computation, P, is small enough (3x3 grid points), the calculation is performed using bilinear interpolation (TSCHERNING et al., 1992).

The described expressions are the formalism applied by the TC program to compute TCR.

1.2 Terrain Corrections Determined by FFT (TCFOUR)

Under a planar approximation, TCR at point P is written as equation (3).

where

In the FFT approach, the solution is based on splitting equation (3) into three convolutions.

Distance l can be written as equation (5):

with

If l is written as (5), integration of equation (1) along the z component, leads to:

A series development of the term in square brackets shows that it converges only if Δh/s <1. This implies that gradients must not exceed 45°.

This development can be transformed into three convolutions in terms of h and h², such as TCR is calculated by equation (8):

with f and f₀, given by equations (9) and (10).

f₀ is a singular integral which presents no problem at zero distance when discrete data is used (FORSBERG, 1985).

Equation (8) can be solved applying a two-dimensional FFT, as is the case of the TCFOUR routine used in this work.

As in the calculation with prisms, two different DEM grids may be introduced here for the same reasons stated above.

TCFOUR is a suitable routine and FFT methods are widely in use due to the speed of FFT computation (FORSBERG, 1985; SIDERIS, 1985; SCHWARZ et al., 1990).

2. ANALYISIS OF THE ESTIMATED TERRAIN CORRECTIONS

2.1 Topographic Analysis

Two topographic maps of the corresponding regions are shown in Figures 2a and 2b. Both were derived from SRTM3 data (Shuttle Radar Topography Mission, FARR et al., 2007).

When computing TCR in Tierra del Fuego, a 20 km radius was used for the interior zone to take into consideration the gravitational effect of the closest topography. A 160 km radius according to the Hayford 167 km zone was applied for the external part (FORSBERG, 1984). Both radii were obtained after experimentations regarding the relationship between the minimum and maximum radius to achieve a required TCR accuracy of 0.1 mGal. For the inner grid a 3" resolution was chosen while for the external grid a 30" resolution was applied.

The radii for the inner and outer zones of Mendoza were the same as those for Tierra del Fuego. However, the resolution was 6" for the inner part and 12" for the external. As can be noticed, this resolution was less than that of the SRTM3. Although there was no significant difference when using 3" resolution for the inner grid, a 6" grid provided the best result. These values were adopted while applying the techniques described above.

The adopted resolutions for the inner and outer grids differ in both regions because it is a compromise between the most detailed representation of the topography and the minimization of the effect of 45° gradients. These gradients are larger in the area of Mendoza because of the steep topography in the whole region under investigation. Under this situation, a lower resolution for the inner grid is mandatory. On the contrary, the outer grid requires higher resolution than that applied to the corresponding grid of Tierra del Fuego.

The heights do not exceed 2600 m in Tierra del Fuego, which makes possible to use a resolution of 3" for the inner grid, even when this region can be elevated. Due to the fact that the highest mountains are located in the south-west part and the remaining heights are of 500 m or less, it is reasonable to use a 30" resolution for the outer grid. This number was determined through a few experiments concerning the variation of both radii.

The density value used in both regions for the TCR determination was 2.67 g/cm³ according to Lodolo et al. (2007) and Blitzcow et al. (2009).

In Figures 3a and 3b terrain correction maps for both study regions are shown. In both cases, TCR were computed using the TC program; the classical approach.

Terrain corrections in Tierra del Fuego reach 60 mGals and 85 % are under 35 mGals.

In the case of Mendoza's province it is observed that terrain corrections exceed 120 mGals and 85 % are under 28 mGals for the selected region.

Kirby and Featherstone (2001) concluded that TCR should not exceed 75 mGals. However, their results were obtained in regions like Hawaii, but not in the Andes. As it is seen here and in Blitzcow et al. (2009), TCR can exceed 100 mGals in the highest parts of the Andes.

Figures 4a and 4b show the results of TCR determination for both regions in analogy to Figures 3a and 3b, but based on the application of the 2D-FFT. A remarkable increase in the maximum values is observed in both areas: the maximums values of TCR in Tierra del Fuego and Mendoza are 80 mGals and 175 mGals, respectively.

In the maps of topographic gradients (Figures 5a and 5b), it is possible to appreciate that the steeper slopes, often exceeding 45°, coincide with the highest mountains. The correlation coefficients between TCR and the topographic gradients are of 60 % in both cases.

Figure 6 confirms that the greatest differences (more than 50 mGals) between both methods correspond to the highest zones. These differences also show a correlation with the largest gradients. This condition may imply that FFT solution is not the most satisfactory one in high mountainous areas even if the resolution of the DEM is lowered. FFT overestimates TCR values in the steepest areas due to the convergence problem.

In flat regions the difference between both methods are negligible like in the center-east part of Tierra del Fuego. This means that if the computation area is extended and flat, FFT can be advantageous because of its speed.

3. THE EFFECT ON GEOID DETERMINATION

Terrain corrections, as stated before, are part of the procedure of geoid determination when gravity data is used. Especially in mountainous regions, they may have a significant impact on the gravity anomalies.

Assuming that the observed gravity anomalies have been accurately estimated (better than 0.1 mGals) there are still two error sources related to the DEM accuracy and the error in the estimation of topographic effects. Then, the focus should be the effect that terrain correction errors have on geoidal height determinations.

It was demonstrated that the differences between both methods for terrain correction determinations exceed 15 mGals in the southwest part of Tierra del Fuego island (Figure 6a). In the case of Mendoza, they are also larger than 20 mGals in the highest region (Figure 6b). In both places the largest differences between both techniques cover extended high mountain regions. Those discrepancies can be regarded as a systematic error of nearly 20 mGal in those areas.

In order to estimate the numerical effect that a 20 mGals error in mountainous zones has on geoid estimation, a test employing the Equivalent Source Technique (EST) was done. The EST method attempts to solve the inverse gravimetric problem: the distribution of masses below a given surface (geoid or topography) which are responsible for the observed external gravitational potential is non-unique. The EST technique consists in the determination of the location and mass of point masses, the configuration of which is consistent to the observed gravity anomalies, geoidal undulations and other observables of the Earth gravity field. For further details, see Dampney (1969), Del Cogliano (2006) and Guspí et al. (2004).

In this case, 654 Helmert (Faye) anomalies distributed in the province of Tierra del Fuego together with the remove-compute-restore method were used. These anomalies resulted from gravimetric reduction and interpolation of 512 observed gravimetric anomalies distributed in the province (Figure 7).

20 mGals were added to produce a bias of that magnitude in part of these 654 gravity anomalies. As it was shown in the previous section, the differences between methods were around 20 mGals in the hightest places of Tierra del Fuego. In this way, the impact of erroneous TCR derived by the FFT approach is simulated. The affected gravity anomalies were located on the south-west side of the province, below 54°S and to the west of 67.6°S. This procedure provided a solution which was compared to an undisturbed gravimetric geoid, without applying the 20 mGals bias. The differences in geoid undulation are shown in Figure 8. The maximum error reaches 1 m in the highest places. Since the difference between both methods is about 0.5 mGals in flat areas (Figure 6a), their contribution would be up to 0.05 m in geoidal heights.

A different result was obtained when 20 mGals were systematically added to all the 654 gravity anomaly values. In this case, the maximum effect on geoidal heights was 2 m (Figure 9).

The results shown in Figure 9 confirm the results published by Vaniček and Martinec (1994). They indicated that a 0.01 mGals systematic error in gravity anomalies produces 1 mm effect in geoidal height.

It would be of great interest to carry out in Mendoza the same experiences as for Tierra del Fuego but the coverage of gravimetric observed data is not sufficient to such purpose.

In Mendoza, where large gradients are distributed all over the region, the error introduced by FFT based TCR can be regarded as systematic. This may imply an error of 2 m or more in geoid determination, depending on the error of the resulting gravity anomalies

4. DISCUSSION AND CONCLUSION

TCR were computed in two different regions of Argentina with two different techniques.

In both regions the differences between both methods correlate with the topography.

The fact that FFT overestimates TCR by tenths of mGals in regions of rough topography, makes classical integration the method which provides the most reliable results in those places. It should be mentioned that the difference between both methods can be neglected in flat areas.

In the case of Mendoza's province, TCR computed by FFT are in excess by more than 50 mGals, which can produce up to 5 m error in geoid determination and may affect geologic interpretations.

A numerical test done in Tierra del Fuego determined up to 1 m geoidal undulation effect caused by a 20 mGals error in gravity anomalies located in the highest places. This error also affected geoidal height determinations in lower areas, in a sensible way.

When the 20 mGals constitute a systematic error in all gravity anomalies, there is a systematic effect of 1m to 2m in geoidal heights too.

Considering that TCR participate in the smoothing procedure of gravity anomalies, when TCR are needed for high-accuracy geoid modelling in high mountains, the use of classical integration instead of FFT is recommended.

ACKNOWLEDGEMENTS

This research was partially financed by a CONICET scholarship and UNLP funds. The authors would like to thank Dr. Renee Forsberg for providing the TC and TCFOUR routines for terrain correction computations.

All figures were generated using Generic Mapping Tools v. 4.5.8 (WESSEL and SMITH, 2012).

BIBLIOGRAPHICAL REFERENCES

BLITZKOW, D.; OLIVEIRA CANCORO DE MATOS, A. C.; CINTRA, J. P. Digital terrain model evaluation and computation of terrain correction and indirect effect in South America, GEOACTA, v. 34, p. 59-74, 2009.

DAMPNEY, C. N. G. The equivalent source technique. Geophysics, v. 34- 1, p. 39-53, 1969.

DEL COGLIANO, D. Modelado del Geoide con GPS y Gravimetría. Caracterización de la Estructura Geológica de Tandil. Tesis Doctoral. UNR, Argentina. 2006.

FARR, T. G.; ROSEN, P. A.; CARO, E.; CRIPPEN, R.; DUREN, R.; HENSLEY, S.; KOBRICK, M.; PALLER, M.; RODRIGUEZ, E.; ROTH, L.; SEAL, D.; SHAFFER, S.; SHIMADA, J.; UMLAND, J.; WERNER, M.; OSKIN, M.; BURBANK, D.; ALSDORF, D. The Shuttle Radar Topography Mission. 2007. Accesed http://www2.jpl.nasa.gov/srtm/SRTM_paper.pdf

FORSBERG, R. A Study of Terrain Reductions, density anomalies and geophysical inversion methods in gravity field modelling. Report 355, Department of Geodetic Science and Surveying, The Ohio State University, Columbus, 1984.

FORSBERG, R. Gravity Field Terrain Effect Computations by FFT. Bulletin Geodesique, v. 59, p. 342-360, 1985.

GUSPI, F.; Introcaso, A., Introcaso B. Gravity-enhanced representation of measured geoid undulations using equivalent sources. Geophysical Journal International, v. 158, p. 1-8, 2004.

KIRBY, J. F.; FEATHERSTONE, W. E. Anomalously large gradients in the "GEODATA 9 SECOND" Digital Elevation Model of Australia, and their effects on gravimetric terrain corrections. Cartography, v. 29-1, 2001.

LODOLO, E.; LIPPAI, H.; TASSONE, A.; ZANOLLA, C.; MENICHETTI, M.; HORMAECHEA, J. L. Gravity map of the Isla Grande de Tierra Del Fuego, and morphology of Lago Fagnano. Geologica Acta v. 5, p. 307-314, 2007.

MCMILLAN, W. D. The theory of potential. Theoretical mechanics, v. 2, 1958.

SCHWARZ, K. P.; SIDERIS, M. G.; FORSBERG, R. The use of FFT techniques in physical geodesy. Geophys. J. Int., 100, 485-514, 1990.

SIDERIS, M. A Fast Fourier Transform method for computing terrain corrections. Man. Geod., v. 10, no. 1, p. 66- 71, 1985.

TSCHERNING, C.C.; FORSBERG, R.; KNUDSEN, P. The GRAVSOFT package for geoid determination. Proc. 1st Continental Workshop on the Geoid in Europe, Prague, p. 327-334, 1992.

VANIČEK, P.; MARTINEC, Z. The Stokes-Helmert scheme for the evaluation of a precise geoid. Man. Geod., v.19, p. 119 - 128, 1994.

WESSEL, P.; SMITH, W. H. F. The Generic Mapping Tools Technical Reference and Cookbook, Version 4.5.8, 2012.

Recebido em abril de 2013.

Aceito em julho de 2013.

Evaluation of terrain corrections through FFT and classical integration in two selected areas of the Andes and their impact on geoidal heights

Avaliação das correções do terreno através da FFT e da integração clássica em duas áreas selecionadas dos Andes e do seu impacto sobre a altura geoidal

Publication Dates

• Publication in this collection
  27 Sept 2013
• Date of issue
  Sept 2013

History

• Received
  Apr 2012
• Accepted
  July 2013