The current terrain correction procedure corrects gravity measurements for the effect of terrain from a distance Rmin to a distance Rmax. Each gravity station is processed independently, and therefore corrections calculated for a particular station do not depend upon possible location errors of other stations. The correction procedure is designed to provide corrections which can be added to the "outer-zone" corrections calculated by existing terrain correction programs which are in common use in the United States, such as that of Plouff (1966). The distance Rmax is specified by the user during execution, and should be chosen to be large enough that the corrections calculated by the outer-zone correction program are accurate. The value selected for Rmax. therefore depends upon the smallest compartment size available to the outer-zone terrain correction program. Typically, that size is 15", 30", or 1' of arc, or roughly 450 to 1800 m. The distance Rmin is also supplied by the user; in areas of relatively low relief, values for Rmin of 5 m or so seem appropriate. In areas of very high local relief, the 30 m digitization cannot accurately represent the local terrain, and thus Rmin should be selected accordingly. A value for Rmin of 45 m [about the radius of Hammer's (1939) Zone C] seems to provide good results where local relief is large.
The correction procedure divides the circular area enclosed by
Rmax into an inner zone and an outer zone; the
radius separating the inner and outer zones is denoted Rmed.
A surface is fit to the elevations between Rmin and
Rmed; this surface is numerically integrated to
calculate that portion of the terrain correction that is due to the
terrain located in the interval Rmin <= R
<= Rmed. The figure below illustrates the
subdivision of areas used in the calculation.
The current version of the program uses the triangulation and interpolation procedures of Renka (1984) to interpolate elevations in this interval; earlier versions of the programs fit a multiquadric surface to the elevation data within the inner zone. The multiquadric surface approach, while quite accurate, was exceedingly slow and required substantial memory. Between Rmed and Rmax, the terrain effect is calculated by assuming that each elevation sample represents the elevation of a rectangular compartment; a line element formula is used to calculate the effect of each such compartment. Terrain compartments lying between Rmed and Rmax that partially intersect the circular radius Rmed are treated in such a way that only the effect of that portion lying outside Rmed is added to the overall terrain effect. A similar procedure is applied to those compartments that partially intersect the outer radius Rmax
The use of multiquadric surfaces to fit the terrain data, while not the most efficient method from a computational standpoint, is a rather accurate method of representing topography (Hardy, 1971), and has been used for calculating terrain corrections for gravity stations by others (e.g, Krohn, 1976; Gettings, 1981). A comparison of several scattered data interpolation methods by Franke (1982) concludes that the multiquadric quation method, from the viewpoint of smoothness and accuracy, was one of the best scattered data interpolation methods of those that were tested. A disadvantage of using the multiquadric equation surface technique is that both time and memory requirements for fitting the multiquadric surface are proportional to N², where N represents the number of elevation samples used to construct the multiquadric surface; for the case under consideration, N is the number of DTM samples falling within or close to the innermost radius Rmin. Moreover, the evaluation of the surface at a single location requires summing the N terms defining the surface; as such a sum must be constructed many times for each innermost terrain correction, the evaluation phase of the calculation can take as long as the surface fitting phase. Finally, all the computations related to the construction and evaluation of the multiquadric equations should be performed in double precision on 32-bit machines.
Note that the numerical integration procedure used to calculate the effect of terrain lying between Rmin and Rmed essentially integrates the effect of a line element between the two radii. This integration is repeated every 5° to obtain the effect of all the terrain lying within the circular region. The radial portion of the integration is performed using an adaptive technique called QUAN8; this routine is described by Forsythe and others (1977, page 97). The technique is really the same as that of Gettings (1981); however, Gettings derived his method in an extremely complicated manner. An very simple derivation can show that the method used here is exactly equivalent.
A terrain surface is used close to the station location because the elevations provided by the DTM are samples and do not actually represent mean elevations of rectangular compartments. The use of a surface provides a terrain representation that should be much closer to reality than the use of compartments of constant mean elevation. At a certain distance from the station, the procedure of considering the elevation samples as representing the mean elevation of a rectangular compartment should yield numerical results than are not distinguishable from the results that would be obtained by actually using mean elevations of compartments whose size would necessarily be larger than 30 m on a side. In fact, the method used is equivalent to numerically integrating the terrain effect (at distances greater than Rmed) using a rectangular rule. Because the compartment size is relatively small, use of the rectangular rule should be rather accurate. Typical values for Rmed and Rmax are 200 m and 4000 m. From a practical standpoint, the value of Rmed is restricted to relatively small values (usually less than 400 m), due to the increased memory and time required to fit and integrate the multiquadric surface to the elevation data. The value of Rmax can rather easily be increased to greater values, though, as no surface fitting is involved. However, little increase in accuracy of the calculation of the overall terrain effect is expected by increasing Rmax to values much greater than 2500 m.
The elevation of the gravity station is not directly used during the computation of the effect of terrain effects. Instead, the elevation at the horizontal location of the gravity station is calculated from the multiquadric representation of the terrain surface, and this calculated elevation is used in the computation of the effect of terrain. The actual elevation of the gravity station is not used at all.
The reason for not using the actual elevation of the gravity station for the computation of the terrain effect is simply that the DTM elevations often differ from the actual elevations by amounts of up to 7 m. The use of the actual elevation, rather than one consistent with the DTM terrain, effectively leads to gravity stations being located in deep holes or on very steep hills, whereas in fact such terrain features probably do not exist in the immediate vicinity of the gravity station. The departure of the DTM elevations from the actual elevations is a relatively local phenomenon, as the mean difference between the DTM and the actual terrain surface should be much closer to 0 than to 7 m.
Numerically, the procedures used for the calculation of the terrain effect are extremely accurate, especially when compared to terrain corrections calculated using template methods. However, the corrections calculated are no better than the terrain data that are used to represent the terrain about each gravity station. These data are typically at least as accurate as the elevation data present on standard USGS topographic maps, though, and thus it seems that the resulting terrain corrections will always be more accurate than those calculated using template methods, when the basic source of terrain data consists of USGS topographic maps.
Up to 5 DTMs may be provided. If more than one DTM is provided, the grids should be provided in increasing coarseness, with sufficient overlap that continuous terrain coverage is obtained for all gravity stations to be processed. For example, one might provide two terrain grids, one of discretization 20 m and another of discretization 50 m: the 20-m grid could be used for terrain corrections to a distance of 2 km, and the 50-m grid could be used to a distance of 10 km. It is certainly possible to use a single DTM having a relatively fine discretization, but computation times can become lengthy (and grid file sizes quite large) if finely-discretized grids are used to calculate terrain corrections to great distances. Note that the terrain grid must have the same horizontal resolution in both the east and north directions; that is, the grids must have square compartments.
The user has a number of choices to make prior to calculating the inner-zone terrain corrections. The use of the default radii for the terrain correction calcuations is often adequate, but the user may want to experiment with the radii used. If a DTM having a resolution different from 30 m is used, the user will almost surely wish to use different values for the correction radii than the default values. Another important factor is to determine if a bias exists between the DTM elevations and the elevations provided with the gravity measurement.
The default values for Rmin and Rmed are 30 m and 250 m, respectively. These values are appropriate for USGS DTM's that are normally distributed with a resolution of 30 m. However, the grid file supplied by the user can have an arbitrary grid increment, and it is likely that it will differ from 30 m. In such cases, the value of Rmin, especially, should be changed to a much smaller value. I suggest a value 1-3 times the spatial resolution of the DTM, depending upon the overall quality of the DTM.
The value chosen for Rmed should be selected in conjunction with the value used for the integration angular increment and the spatial resolution of the first DTM used (which will always have the highest spatial resolution). At a radius of Rmed from the gravity station, the numerical sampling increment will never be closer than Rmed×sin(angle), where angle is the angular increment used. The horizontal distance implied by this relation should not be much greater than the spatial resolution of the DTM. Normally, values of angle and Rmed should be selected such that the horizontal distance implied by this relation is roughly the same as the spatial resolution of the DTM. For example, the default values for Rmed and angle imply a horizontal distance (at a radius of Rmed from the gravity station) of 26 m, slightly less than 30 m.
It is very important for the user to attempt to determine what bias, if any, exists between the elevations in the DTM and the elevations accompanying the gravity measurements. The author's experience has been that such biases are relatively common. To assist in the identification and estimation of such biases, RasterTC notes the actual and the interpolated value for all elevations in the gravity data set. Additionally, any gravity measurement that is "close" (within 2.5 m) to a DTM elevation is noted and the two elevations are compared. Such information appears in the listing file generated by the program. If a systematic elevation bias exists between the gravity data elevations (as well as the auxiliary elevation data, if used) and the DTM elevations, the user should either (a) always use interpolated elevations rather than actual elevations in the terrain correction calculation or (b) apply a constant shift to the gravity data elevations in order to minimize errors associated with the elevation bias. Note that if auxiliary elevation data are used, it is probably best to apply an elevation shift (which is always applied to both the gravity station data as well as the auxiliary elevation data) rather than simply use interpolated elevations, as the interpolated elevation data could vary rapidly horizontally if the DTM data and the auxiliary elevation data, which collectively form the basis for the interpolation, are used when a bias exists. Note that one can always apply a bias correction and use interpolated elevations.