The scaling behaviour of the edges of gravity data

Measurements of the earth's gravity field are widely used in geophysical exploration programs. The geological interpretation process often involves the identification of the boundaries, or edges, of different regions. This can be achieved through a variety of techniques. This paper examines the statistical distribution of the size of the edges produced by a synthetic gravity model, and compares the results with those obtained from a gravity dataset from South Africa.


Introduction
The identification of the boundaries between different geological units is often performed using potential field data, which can be obtained quickly and relatively inexpensively by a variety of airborne platforms.Techniques such as directional derivatives or sun shading (Horn, 1982) can be used to emphasise edges with particular orientations.Alternatively edges with any orientation can be enhanced using the gradient amplitude (or total horizontal derivative, TDX) of the field, f (Jahne, 2005, p. 339): The regions of the data with the largest gradients will obviously produce the greatest values of the TDX.The zero contours of the second horizontal derivative can also be used to locate edges within geophysical datasets, and they have the additional property of being able to track edges with any value of TDX.For map datasets the second horizontal derivative can be computed in different ways, such as the Laplacian L (Jahne, 2005, p. 345): (2)

The statistical distribution of edges produced by an ensemble of point sources
The gravity response of a buried sphere is given by; where G is the gravitational constant and m is the mass of the sphere.Its second horizontal derivative is given by where r 2 = (x 2 + y 2 ).Hence the zero contour of the Laplacian L of the anomaly from a sphere will have a radius of z/2. Figure 1a shows a synthetic gravity dataset produced by a set of spheres.Their masses and depths are identical, and their locations are random.Their edges, as determined by the zero contour of the Laplacian, are plotted in Fig. 1b, and Fig. 1g shows the histogram of the edge contour lengths.Figure 1 c, d, e, f show the plots obtained when different numbers of spheres were used.The contour locations were determined by linearly interpolating the Laplacian values of the data.When the spheres are few and relatively shallow then the overlap of anomalies is infrequent and the histogram has a prominent peak at a contour length of π.z, where z is the depth of the spheres.However, as the number of spheres increases, their anomalies will interact more often and the histogram shape becomes a power law.The presence of a power law in the behaviour of a system is a measure of its selfsimilarity, i.e. portions of an object look similar to the object itself, when suitably magnified.In this case, small contours have a similar appearance to large contours.Many objects in nature, such as clouds and trees, display fractal characteristics (Mandelbrot, 1982).This type of scale invariant behaviour is exhibited by many geoscientific phemonema, such Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.as seismicity (Turcotte, 1997, p. 56), coastlines (Richardson, 1961), topography (Gagnon et al., 2006), faulting (Bohnenstiehl and Kleinrock, 1999), and non-linear least-squares inverse theory (Cooper, 2000).Fractal analysis has also been used to determine an optimum grid size in the interpolation of geophysical data (Dimri et al., 2005).
The fractal behaviour of a system is characterised by its fractal dimension, which is determined from the exponent of the power law that describes its behaviour in some manner.Customarily this is performed by taking logarithms of some measure of the data (in this case, the number of edge contours within a range of lengths is plotted against the mean value of those lengths) and then a straight line is fitted to the result.The (negative of the) gradient of the line then gives the frac- tal dimension of the object or system.Figure 1i shows a plot of the fractal dimension obtained as a function of the number of spheres used.As more spheres are added their responses overlap more and hence the degree of spatial averaging of the gravity field gradually increases, the result eventually being a small number of relatively large contours.Boundary effects then become important as most contours intersect the spatial limits of the dataset.The statistics on the lengths of the contours then become unreliable and the power-law relationship breaks down.The shallower the sources are, the smaller their interaction and the greater the number of spheres required before the breakdown occurs; conversely, the deeper they are, the lower the fractal dimension obtained and the smaller the number of spheres required to produce a breakdown in the power-law relationship.
All calculations were performed in Matlab (under Windows 7 32 bit) on an Intel Core i7 processor running at 3.33 GHz .It was found initially that computing the gravity of 100 000 sources over a 1000 × 1000 point grid would take approximately 75 min, so the process had to be speeded up.Firstly, since each computation is essentially the same, one prior computation of the anomaly of a single source centred in a grid 2000 × 2000 points in size was made.In all, 1000 × 1000 point portions of this grid were then used when its centre was located randomly within the desired 1000 × 1000 point grid.Hence no repeated computations of the anomaly were made.Secondly, the computations were moved from the CPU to the GPU (nVidia GeForce GTX295) using the CUDA architecture.The time of computation was then reduced to a more reasonable 90 s.
Figure 1j shows the result of allowing the sphere masses to vary (randomly) as well as their location.The result is similar to that shown in Fig. 1i, although a higher fractal dimension is obtained.If the sphere depths are varied instead of their mass a similar behaviour was also observed.Fractal dimension obtained as a function of the number of randomly located buried sphere sources, at a depth of 5 m.The mass of each sphere was generated randomly between 0 and 10 11 kg.Because the fractal dimension depended on the locations of the spheres and these were calculated randomly, the experiment was repeated five times and the mean value and the standard deviation determined.
3 The statistical distribution of the edges of a gravity dataset from South Africa Figure 2a shows a Bouguer gravity dataset from South Africa.The image is 650 km ×700 km in size, and the original approximately 90 000 gravity measurements (made by the Council for Geoscience, Pretoria) were gridded to a cell size of 1 km.The prominent circular anomaly on the middle right side is the Trompsburg high (Buchmann, 1960), and the Vredefort dome impact structure is just visible in the upper right corner of the figure.This portion of the South African gravity dataset was chosen because it is the largest rectangular area that does not intersect South Africa's borders, beyond which no data were available.Figure 2b overlays the zero contours of the Laplacian on the gravity data, and Fig. 2c and d show the histogram of the edge contour lengths and a power-law fit.As suggested by the theoretical framework in the previous section, the scaling behaviour of the edges of the gravity data followed a power law, which had a fractal dimension of 1.87 for this particular area.

Conclusions
The scaling behaviour of the edges of ensembles of point gravity sources were studied, and the distribution of the lengths was found to depend upon the number of sources in the map.A power law could be reliably fit to the length distribution once sufficient sources were present (the number was dependant on the source depth).This suggested that the scaling behaviour of the edges of a real gravity dataset would show similar power-law behaviour, and this was found to be the case.

Fig. 1a .
Fig. 1a.The gravity response from 100 000 point sources, each with a mass of 10 11 kg and a depth of 5 m.

Fig. 1b .
Fig. 1b.Edge contours of the data shown in Fig. 1a.The contours are coloured according to their length.

Fig
Fig. 1c, d.Histogram and log-log histogram of the lengths of the edge contours obtained from the gravity response from 10 000 sources at a depth of 5 m.A fitted power-law distribution is shown as a solid line.

Fig
Fig. 1e, f.Histogram and log-log histogram of the lengths of the edge contours obtained from the gravity response from 50 000 sources at a depth of 5 m.A fitted power-law distribution is shown as a solid line.

Fig
Fig. 1g, h.Histogram and log-log histogram of the lengths of the edge contours in Fig. 1b.A fitted power-law distribution is shown as a solid line.

Fig. 1i .
Fig. 1i.Fractal dimension obtained as a function of the number of randomly located buried sphere sources of mass 10 11 kg, at a depth of 5 m.Because the fractal dimension depended on the locations of the spheres and these were determined randomly, the calculations were repeated five times and the mean value and the standard deviation determined.

Fig. 2a .
Fig. 2a.Gravity data over a portion of South Africa.The grid spacing is 1 km.

Fig. 2b .
Fig. 2b.The gravity dataset from Fig. 2a is overlain by the zero contour of the Laplacian function.The contours are colour coded by their length.

Fig
Fig. 2c, d.Histogram and log-log histogram of the lengths of the edge contours in Fig. 2b.A fitted power-law distribution is shown as a solid line.