Application of terrain corrections in Australia (2001)

N.G. Direen (with contributions from T. Luyendyk)

• Abstract
• What is terrain correction? Why is it desirable?
• Some philosophies for calculating terrain corrections
• Intrepid^TM algorithm
• Implementation
• Elevation determination
• Algorithm testing
• General conclusions
• Acknowledgments
• References
• Related links

Abstract

Recent improvements in quality of position and elevation determinations due to Differential GPS have led to substantial improvements in the accuracy and precision of observed gravity data, and a general renaissance in the use of gravimetry in Australia. Increasingly, users of the data are specifying levels of precision in derived measurements such as the Bouguer Anomaly, that reflect only these positioning improvements without considering topographic effects on precision. The consequence has been that most gravity data routinely submitted to the National Gravity Database are Simple Bouguer Anomalies only, and contain a range of unquantifiable acquisition noise that degrades the quality of the overall dataset.

This paper briefly describes the principles behind the terrain correction of observed gravity data, and outlines one method of rapidly calculating the value of a "first pass" terrain correction for gravity values in Australia. This method makes use of the widely available 9 second Digital Elevation Model of Australia and the Intrepid^TM geophysical software package developed by Desmond Fitzgerald & Associates. The assumptions and pitfalls of this particular method are also described.

A comparison between the Intrepid method and traditional Hammer-chart manual methods is described using a test dataset from northern Tasmania. The comparison shows that the Intrepid ^TM method is suitable for applying coarse corrections using the generally available DEM. However, the limitations of the DEM itself prevent use in extreme topographic situations, where a more detailed DEM is needed. Overall, in Australian conditions, the method should prove practical.

[back to top]

What is terrain correction? Why is it desirable?

Gravimetric methods measure variations in the gravitational field in order to deduce changes in the structure and composition of the Earth. The gravity field measured at a given location is proportional to the size of attracting masses and inversely proportional to the square of distance from those masses. Variations of gravity due to these two effects may mask the field effects from sources of interest.

For example, when the source of part of the gravity field lies at a constant position in the Earth, the meter may be nearer or closer to the source as it is moved over topography. In addition, other masses above or around the meter may attract the small mass inside the meter. This would cause both the wavelength and amplitude of the field of the source of interest to be distorted, or obscured altogether. For example, removal of a significant terrain correction to reveal potentially economic mineralisation targets in western Tasmania is discussed by Roberts & Mudge (1997).

To avoid these types of spurious effects, the Bouguer Anomaly is derived by a series of corrections to the observed data (Figure 1). These corrections produce a gravity value for the observation point free from the effects of observing above or below a given datum (Figure 1a), and any masses interposed between the meter and the datum (Figure 1b). The complete Bouguer correction also accounts for the effects of any masses that are higher than the observation point, or voids lower than the observation point (Figure 1c). These effects can be imagined by considering the effect on a pendulum at the observation point; the bob is attracted by excess mass (m) above, and further deflected by mass deficit (-m) below (see Figure 1c inset). Both effects produce a lower value for the vertical component of the field. The correction for these masses is called the terrain correction.

In principle, the complete Bouguer correction should also contain a further correction for the curvature of the Earth, which causes the terrain correction for a flat "Bouguer Slab" to be overestimated at distances far from the meter (Figure 1d). In practice, the subdued nature of most Australian topography results in negligible terrain corrections (ie less than data acquisition noise) at distances of >22 km from most stations, where earth curvature corrections are required (Leaman, 1998). Both terrain corrections past 22 km and their earth curvature recorrections are therefore ignored for the purpose of this discussion.

Past practice in Australia has focussed on the method of Hammer (1939) for calculating terrain corrections; for example all Tasmanian gravity stations have a correction calculated using this method. This method involves the use of a graticule of segmented annuli (Figure 2), representing sections of cylinders of different radii centred on the station. The graticule is overlain on topographic maps of the terrain about the station, and the average height in each annular segment estimated with the corresponding terrain effect calculated from the gravimetric attraction of a cylindrical segment. As noted by Sandberg (1958) and Leaman (1998), Hammer's method is manual, tedious and prone to estimation error.

Whichever method is chosen, the main effect of the terrain correction is to raise the value of the Bouguer Anomaly in the presence of significant terrain, where the measured value will be too low. But what constitutes significant terrain? The answer to this question is determined by the size of the target sources the gravity method is being used to define:

• a significant piece of terrain is anything big enough and close enough to the meter to mask the signal of the target.

This relationship is one of diminishing returns – terrain of low relief close to the meter may be significant, but not when it is further away. Significant features close to the meter may include rocks, small cliffs, slope breaks, and ditches with relief of up to a metre or two (Leaman, 1998). Significant features further away from the meter would usually include hills, mountains and valleys of increasing size with distance.

The current noise envelope due to errors in height in Differential GPS gravity surveys positioned to 10 cm, equates to an error in the Simple Bouguer Anomaly of approximately 0.3 µms^-2. To maintain this precision requires correcting for topographic features of heights that exceed the limits given in Figure 3, based on Hammer's (1939) tables. Calculations have been made for features that:

1. occupy one sector of the hollow cylinder annulus eg a small hill;
2. occupy one complete annulus eg a range forming a ring around a basin
3. occupy all annuli eg applicable when the observation is made on top of a hill on an otherwise featureless plain.

Most real topography would usually fall somewhere between these end-member cases.

In Australia, the terrain correction is rarely applied to gravity data, whatever the scale of survey. This is not because there is no "significant terrain", but because users of the data are prepared to accept an unspecified amount of random, topographically-induced noise in their data rather than embark on the previously time-consuming task of computing the relevant terrain corrections.

[back to top]

Some philosophies for calculating terrain corrections

All terrain corrections work by calculating the gravitational effect of excess or deficit mass based on the topography around the station. Typically this is achieved by dividing the topography around any station into a series of sectors, having an assumed average height (eg Hammer, 1939), although other methods and approximations may be applied (eg Parker, 1995, 1996)

Four main philosophies apply to the calculation of terrain effects (Figure 4):

• Splitting the terrain into dipping planes (Figure 4a), and calculating the gravitational effects of these surfaces (Sandberg, 1958; Zhou et al., 1990). An analytic solution can be calculated only where a high-resolution topographic model (5–10 m cell size) is available for the inner zones (out to hundreds of metres). For this reason, these corrections are not widely applicable in Australian conditions, as topographic models of sufficient precision are rarely available.
• Splitting the terrain into a series of vertical prisms, such as right rectangular prisms (Figure 4b), cylindrical segments (Figure 4c), rods, or conical elements, and summing the effects of the analytic solutions for these elements (Hammer, 1939; Bott, 1959; Kane, 1962; Olivier & Simard, 1981) or combinations of these (Roberts & Mudge, 1997). These methods have less accuracy in the inner zones, because the tops of the prisms are usually considered to be flat with an average height. However, greater accuracy has been obtained using sloping topped prisms (Blais & Ferland, 1984) and conical element methods, but both are more computationally intensive, and additionally require accurate merging with any flat-topped solutions further from the station.
• Splitting the terrain into a series of vertical "mass-lines" (Figure 4d), and integrating these for an area around the station, but not in the area nearest the station, where assumptive errors may be very large (St.John & Green, 1967).
• Using the mathematical properties of any of the above solutions and approximating them using a Fast Fourier Transform on a gridded digital terrain model, in order to speed up the process (Sideris, 1985; Cogbill, 1990; Parker, 1995,1996; Kirby & Featherstone, 1999). The limitations of this method depend upon the properties of the analytic solution transformed into the frequency domain.

All of these variants are amenable to calculation by computer, based on digital elevation models. However, as Leaman (1998) has noted, the greatest effects in the gravity field come from terrain near the meter, which is usually not resolved by most terrain models in use. This results in an undercorrection, but any correction is better than none, especially where large terrain effects occur outside this inner zone.

The pixel resolution of the Australia-wide DEM, available from Geoscience Australia, is approximately 250 m. The standard deviation of vertical cells in the model from known spot heights is around 35 m ( https://www.ga.gov.au/products/servlet/controller?event=DEFINE_PRODUCTS ; viewed 2/5/00). The horizontal resolution approximately covers Hammer zones A to E, in which height differences of up to 36 m are significant (Figure 3). Thus any correction developed from this DEM will be considered a minimum approximation to corrections for Hammer zones F to M, and further estimates must be made for the important inner zones (A–E). Note that even using 25 m cell DEM's, available in some parts of Australia, height differences of up to 4 m would still be significant in these inner zones. These effects necessitate careful station positioning away from, and/or detailed descriptions of, near-station relief, as outlined by Leaman (1998).

[back to top]

Intrepid^TM algorithm

The method trialled for applying automated terrain corrections using the Australia-wide DEM and incorporated into the Intrepid^TM processing system is outlined below.

The method divides the region surrounding a gravity station into cells in the x, y plane. A mean elevation is assigned to each cell by interpolation from the DEM grid regridded with available gravity heights to improve accuracy. Prisms are then formed by projecting the cells up or down to the plane of the station elevation (Figure 5). Each prism is assigned a standard density (2.67 t/m³) and the terrain correction is then calculated at the station as the sum of the gravitational effects due to each prism.

To improve calculation times, cell area is increased in proportion to distance from the station. This is effected by considering concentric circles centred on the station. In the circle closest to the station, each cell has a user-defined area, allowing the user to take advantage of any high-resolution DEM available. Moving out, the area of the cells in each subsequent circle is doubled.

Two methods are used to calculate the gravitational effect of a prism. For prisms falling within the innermost circle, calculations are based on the analytic solution of a right rectangular prism (Coggon, 1976). In all other cases, a thin rod solution (Coggon, 1976) is used to speed up the calculation.

This method provides maximum precision in the zones nearest to the station, while allowing more efficient calculation further away. In addition, the algorithm is flexible, allowing the cell size to be chosen according to the resolution of the DEM, and also allowing the gravity station heights to be merged with the DEM, thus reducing a major source of error (Leaman, 1998; Cogbill, 1990).

[back to top]

Implementation

The user specifies the number of circles to use and a minimum cell size for the inner prisms, from which the size of the concentric circles surrounding each station are also calculated. The minimum cell size chosen should be a function of the resolution of the DEM. For example, the 9'' DEM has a basic cell size of approximately 250 m; any cell size specified smaller than this will require interpolating beyond the accuracy of the data. A maximum of 5 circles can be used. The factors to convert from cell size to circle radii are currently 16, 32, 64, 256 and 1024. Note that using all five rings with the 9'' DEM would provide a terrain correction out to a distance of 256 km, necessitating a further correction for the curvature of the earth beyond 22 km.

The DEM grid is sub-divided into contiguous sub-blocks, each sub-block having a size equal to that of the largest cell size — 16 times the minimum cell size. Since cell sizes are multiples of 2, a sub-block will always be an exact multiple of that of any allowed cell. Instead of calculating the terrain correction for each station in turn, the algorithm proceeds by sequentially considering each sub-block of the DEM grid and assigning its effect to each station. The procedure for each station follows.

The smallest distance from the current sub-block of the grid to the station is calculated and is used to determine the circle about the station in which the sub-block lies. This in turn determines the cell size of the prisms to be used for this sub-block/station pairing. If the sub-block falls beyond the largest circle, it is ignored; if not, the effect of each prism comprising the sub-block is computed and added as an on-going total for the station.

[back to top]

Elevation determination

Two methods are available to interpolate DEM grid values to an elevation at the (x, y) point at the centre of the base of a prism:

1. A value is linearly interpolated directly from the four grid points that surround the (x, y) point.
2. Inverse distance averaging is used. A search radius about the (x, y) point is nominated and all grid points within this radius are used to determine a weighted average. At least 8 grid points are required before an average is accepted, with the search radius being increased until this condition is met. The initial search radius is three times the grid cell size. Failure occurs if the search radius exceeds either the prism cell size or 9 times the grid cell size, resulting in the prism being ignored. If the search radius is R and the distance from the (x, y) point to the grid point is D, then the weight assigned to that grid point is (R-D)/D.

If inverse distance averaging is used, the actual height of the station can be included.

In evaluating computer generated terrain corrections with those computed manually, significant differences were found in those cases where the elevations calculated from the DEM grid disagreed with those measured with the gravity readings. This occurs where the gridding algorithm is unable to exactly honour the original gravity station heights due to the grid sampling interval and curvature parameters. Good agreement can be obtained by simply replacing the gravity station elevation with that calculated from the DEM grid. Consequently, the option also exists to replace all station elevations by those interpolated from the DEM grid for calculating the terrain correction.

[back to top]

Algorithm testing

Test data were chosen from gravity data covering part of northern Tasmania, around the area of the Longford Basin (Direen & Leaman, 1997). The area chosen has three characteristics useful for testing a terrain correction algorithm:

1. it contains a range of both rugged and gentle topography in a relatively small area;
2. it has a reasonably homogeneous areal spread of gravity stations, including on the Western Tiers escarpment area, which rises from 150 masl to around 1300 masl;
3. it has a station density better than 1 station /km² over most of the area.

The area of study has also been the subject of at least three detailed gravity surveys and has been gravity modelled by Hinch (1965), Longman & Leaman (1971), Direen (1995) and Direen & Leaman (1997). The area chosen was also the first region of Tasmania where the need for the application of a Terrain Correction to the Bouguer Anomaly was recognised in 1970 (D.E. Leaman, Pers. Comm., 1995).

The dataset used in testing was taken from that provided to AGSO by Mineral Resources Tasmania. It consists of 3200 stations in the area delineated by 146.5E–147.3E, 41.16S  41.83S. For each station, values for the Simple Bouguer Anomaly (un-terrain corrected), and Total Bouguer Anomaly (ie. with the Terrain Correction applied) at 2.67 t/m³ have been calculated. Terrain corrections were calculated by hand, using the Hammer (1939) graticule method; this assumed a density of 2.67 t/m³ and was carried out to a radius of 22 km (Hammer Zone M). The range of terrain corrections manually applied to the data was 0 to 234 µms^-2.

The digital elevation model (DEM) used for testing Intrepid was the AUSLIG 9'' (~ 250 m) DEM. The gravity measurements, spatially sampled down to 400 m, support an elevation model cell size of approximately 100 m, so the gravity station heights were gridded into a new DEM with a cell size of 100 m, or more than twice the horizontal resolution of the original DEM (Figure 6). An indication of the ruggedness of topography in the area is shown in the profile from the western Tiers to Mt Barrow (Figure 6b).

The SBA and the TBA (with manual TC's) for the gravity data were then gridded in Intrepid^TM (Figures 7, 8) to provide the benchmark data against which the automatic terrain correction module was tested. There are significant differences in these original grids (Figure 6b), highlighting the importance of the terrain correction in the test area, particularly around the Western Tiers, where the SBA is too low by the order of 100 µm/s^-2.

The terrain correction generated using Intrepid^TM on the SBA data was then calculated and gridded, implementing modifications to the code to force the station heights onto the DEM surface (Figure 9). The initial cell size chosen was 125 m, with four rings chosen. This produces a terrain correction out to 256 × 100 m or 25.6 km, which is approximately equivalent to the Hammer Zone M corrections applied manually.

To good approximation, these constraints should replicate the details included in manual corrections, with the exception of the innermost Zones (out to 125 m, or Hammer zones A to ~D). In the original manual calculations, Zones from C to G were calculated using 1:25 000 scale maps, whereas G  to M were calculated using 1:100 000 scale maps with less precision. Zones A to C were manually calculated from detailed terrain observations made during acquisition, and are not able to be modelled at the resolution of the DEM — thus any automated estimates of the terrain correction will be minima.

Differencing of the original TBA data and the automated approximation was carried out and the results gridded (Figure 10). Good results were achieved by the program in areas of low relief or of small topographic change (Figure 10b). The algorithm performed poorly on the Tiers (overcorrection), and at the Tamar mouth and on the Central Plateau (undercorrection). The Tamar Valley and environs (Mt Barrow, Mt Arnon) also appear to disrupt the routine (usually overcorrecting).

The large undercorrection at the Tamar mouth gives some clues to the reasons for this poor performance. In this area, the AUSLIG DEM cuts out, due to Bass Strait. The true terrain correction value is affected by the water and bathymetry in this area, but the algorithm is unable to consider this data. Thus, the failure of the algorithm in this area is really due to the lack of data in the DEM rather than the assumptions the program is forced to adopt. It does, however, highlight a limitation on the general use of the algorithm in near-coastal areas if the AUSLIG 9'' DEM is being used. In the case of the Western Tiers and Tamar Valley regions, the errors may also be related to the inability of the DEM to fully model the precipitous changes in slope.

[back to top]

General conclusions

The Intrepid^TM algorithm performs to reasonable expectations, and is able to broadly replicate results of manual terrain correction. However, the extreme nature of the topography included in the test data, and the inability of the DEM to fully model this effect led to significant failure (errors of up to 140 µms^-2 in a dataset with a dynamic range of 120 to –400 µms^-2). This is a general constraint on automated systems, as has been noted by Leaman (1998), and Leaman (1991). The lack of a definite tendency to either under- or over-correct is more problematic; in theory all automated corrections should be minimum estimates (ie under-corrected).

Tests are still to be conducted on a continent-scale dataset to determine the efficiency of the algorithm when dealing with larger blocks of stations.

[back to top]

Acknowledgments

The author would like to thank Tony Luyendyk for his dedication in testing various parts of the Intrepid algorithm. Alice Murray, Ian Hone and Michael Morse are thanked for reviews that substantially improved the manuscript. Anna Dyke is thanked for providing a copy of her M.Sc thesis and for engaging discussions on the subject. This paper is published with the permission of the Chief Executive Officer, AGSO.

[back to top]

References

Blais J. A. R. & Ferland R. 1984. Optimization in gravimetric terrain corrections. Canadian Journal of Earth Sciences 21, 505-515.

Bott M. H. P. 1959. The use of electronic digital computers for the evaluation of gravimetric terrain corrections. Geophysical Prospecting 7, 46-54.

Cogbill A. H. 1990. Gravity terrain corrections calculated using digital elevation models. Geophysics 55, 102-106.

Coggon J. H. 1976. Magnetic and gravity anomalies of polyhedra. Geoexploration 14, 93-105.

Direen N. G. 1995. Geophysical modelling of the Longford Basin, northern Tasmania. Unpublished B.Sc Hons thesis, University of Tasmania

Direen N. G. & Leaman D. E. 1997. Geophysical modelling of structure and tectonostratigraphic history of the Longford Basin, northern Tasmania. Exploration Geophysics 28, 29-33

Dyke A. 1997. Assessment and application of gravity terrain correction methods. Unpublished M.Sc Thesis, University of Leeds.

Hammer S. 1939. Terrain corrections for gravimeter stations. Geophysics 4, 184-194.

Hinch A. J. 1965. A gravity survey of the Cressy area, northern Tasmania. Unpublished B.Sc Hons thesis, University of Tasmania.

Kane M. F. 1962. A comprehensive system of terrain corrections using a digital computer Geophysics 27, 455-462.

Kirby J. F. & Featherstone W. E. 1999. Terrain correcting Australian gravity observations using the national digital elevation model and the fast Fourier transform. Australian Journal of Earth Sciences 46, 555-562.

Leaman D. E. 1998. The gravity terrain correction –practical considerations. Exploration Geophysics 29, 467-471.

Leaman D. E. 1991. Exploration geophysics, remote sensing and image processing. CODES M.Econ Geol Short Course Manual 13.

Longman M. J. & Leaman D. E. 1971. Gravity survey of the Tertiary basins in northern Tasmania. Tasmanian Geological Survey Bulletin 59, Tasmanian Department of Mines, Hobart.

Olivier R. J. & Simard R. G. 1981. Improvement of the conic prism model for terrain correction in rugged topography. Geophysics 46, 1054-1056.

Parker R. L. 1995. Improved Fourier terrain correction, Part I. Geophysics 60, 1007-1017.

Parker R. L. 1996. Improved Fourier terrain correction, Part II. Geophysics 61, 365-372.

Peng M., Li Y. C. & Sideris M. G. 1995. First results on the computation of terrain corrections by the 3D FFT method. manuscripta geodaetica 20, 475-488.

Roberts S. S. J. & Mudge S. T. 1997. Magnetic and gravity modelling of the Renison Tin Mine, Tasmania. Exploration Geophysics 28, 292-295.

Sandberg C. H. 1958. Terrain corrections for an inclined plane in gravity computations. Geophysics 23, 701-711.

Sideris M. G. 1985. A fast fourier transform method for computing terrain corrections. manuscripta geodaetica 10, 66-73.

St John V. P. & Green R. 1967. Topographic and isostatic corrections to gravity surveys in mountainous areas. Geophysical Prospecting 15, 151-162.

Zhou X., Zhong B. & Li X. 1990. Gravimetric terrain corrections by triangular-element method. Geophysics 55, 232-238.

[back to top]

Related links

• Geophysics home page
• Continental Geophysics
• Australian Fundamental Gravity Network (AFGN)
• Seismic Acquisition and Processing
• AEM Acquisition and Interpretation
• Magnetic Anomaly Grid Database of Australia (MAGDA)
• Airborne Geophysics and Gravity Data
• Gravimetry
• A transect through central Australia