Nonlinear Processes in Geophysics Grid preparation for magnetic and gravity data using fractal fields

Most interpretive methods for potential field (magnetic and gravity) measurements require data in a gridded format. Many are also based on using fast Fourier transforms to improve their computational efficiency. As such, grids need to be full (no undefined values), rectangular and periodic. Since potential field surveys do not usually provide data sets in this form, grids must first be prepared to satisfy these three requirements before any interpretive method can be used. Here, we use a method for grid preparation based on a fractal model for predicting field values where necessary. Using fractal field values ensures that the statistical and spectral character of the measured data is preserved, and that unwanted discontinuities at survey boundaries are minimized. The fractal method compares well with standard extrapolation methods using gridding and maximum entropy filtering. The procedure is demonstrated on a portion of a recently flown aeromagnetic survey over a volcanic terrane in southern British Columbia, Canada.


Introduction
Magnetic and gravity data are often subject to a number of processing or enhancement techniques designed to improve their interpretive value.These procedures, such as calculating derivatives, downward continuation, reduction to the pole, can all be efficiently carried out in the frequency domain using fast Fourier transforms (FFTs) (e.g., Kanasewich, 1981).FFT algorithms require that data exist at all points on the grid (i.e., there are no undefined values) and that the data are periodic with a period equal to the grid dimensions.Since magnetic and gravity data sets do not usually satisfy either of these requirements, they must be appropriately prepared before FFT-based grid processing algorithms are used (Cordell and Grauch, 1982;Ricard and Blakely, 1988).FFT algo-rithms can be based on powers of two or arbitrary factors.In the following we assume that a power of two FFT is used.
The shape of a survey is driven by a number of factors including the geology of the region being investigated, acquisition cost, flightline orientation and the configuration of the targeted geological feature, if known.As an example, Fig. 1 shows a portion of an aeromagnetic survey carried out in British Columbia.Standard gridding of this set of flightline data produces a complex shaped grid (blue line, Fig. 1) that must then be expanded to a rectangular form (black rectangle, Fig. 1) so that there are no undefined values in the grid.The undefined parts of the grid must be populated by extrapolated values usually generated from the original gridded data.This expanded and now fully defined rectangular grid is, however, not yet in a form ready for applying FFT algorithms.The assumption of grid periodicity assumed in FFT algorithms must first be satisfied.In the simple case of a survey covering a rectangular area, one can remove the mean of the data, apply a tapering window to the data along the edges and pad the grid with zero values.The width of the tapering window should be about one tenth of the size of the grid.Use of a tapering window results in loss of data along the edges.Moreover, real surveys are seldom rectangular.
For non-rectangular data grids, periodicity can be achieved in a variety of ways.Initially, some form of extrapolation is used to add data along survey edges and avoid the loss of information that a tapering window applied inside the survey limits would cause.Care must be taken with the extrapolation in order to avoid problems with discontinuities in the values at or near to the grid edges.Edge discontinuities will cause ringing (Gibb's phenomenon), when grids are filtered or transformed in the course of standard potential field data processing, such as continuation or calculating derivatives.Because of the assumed periodicity of the data grid, ringing will often "bleed" from one edge of the grid to the opposite edge (wrap-around).To avoid this, extrapolated values must  be tapered within a specified zone at the grid edges down to some average grid value, which then exists at all points around the grid perimeter.In practice the mean of the grid and in some cases a first-order trend should be removed before processing.Linear trends are often present in data sets and cause large differences in values at opposite edges of the resulting grids; thus, the most expedient way to treat these effects is to simply fit a low-order surface to the input grid and remove the trend.If necessary, this can be added back to the data after further processing.The width of the tapered zone must be large enough to smooth out any large magnitude differences between values at opposite grid edges.Padding or extrapolated regions at grid edges must have large widths, especially when three-dimensional (3-D) inversions of the data are done since fitting sources at depth will affect data at large horizontal distances.An alternative to tapering is to introduce extrapolated values that are already guaranteed to match and be periodic at the grid edges and do not contain large magnitude differences near to the grid perimeter.This is the technique used in the fractal approach below.
Where real and synthetic data meet (blue line, Fig. 1), extrapolated values need to fit seamlessly to the measured data and should preserve the statistical character of the input data as closely as possible.This ensures that discontinuities caused by a sudden change in wavelength and amplitude content are avoided and that spurious wavelength components are not introduced into the extrapolated grid, which after filtering could affect values within the survey limits.
Common methods that are used for the extrapolation include minimum curvature gridding (Briggs, 1974), mirroring the input data (Baranov, 1975, p. 42) and linear prediction filtering (Gibert and Galdeano, 1985).The two latter approaches are 1-D in that they act on the grid rows and columns independently.Mirroring the data will preserve the frequency content of the original grid values, but when the input grid edges are irregular, the extrapolated values may require some smoothing to avoid large offsets between adjacent rows or columns in areas where extrapolated data are mirrored.Predicting data using maximum entropy preserves the original frequency content of the data along the columns and rows of the grid only up to how well the underlying autoregressive data model fits the observations.Operating on columns followed by rows or vice versa presents the same problems as the mirroring approach.Extrapolation using gridding has the advantage of being a 2-D approach.Nonetheless, the extrapolated values become much smoother than the input data as the distance from the survey grid edges increases.This smoothing changes the frequency content of the input data, which can have serious consequences on the final processed data grid.
The aim of this paper is to introduce an approach to grid preparation based on a fractal description of the data, whether gravitational or magnetic.The parameters of the fractal model are determined from the observed data and the synthesized fractal field used to extrapolate grid values into undefined regions where required.

Fractal description of potential fields
An alternative to existing approaches to grid preparation discussed in the previous section is to use knowledge of the character of the field so that any predicted values have the same statistical behaviour and frequency content as the known gridded values.Magnetic and gravitational fields have been shown to be well-described as fractal (Gregotski et al., 1991;Pilkington and Todoeschuck, 1993;2004;Maus and Dimri, 1994;Lovejoy and Schertzer, 2007).Specifically, the power spectrum of both fields is proportional to f β , where f is the spatial frequency and β is the scaling exponent (the slope of the spectrum in log-log space).For gravity data, published values of β range from −4.5 to −5 based on regional and continent-wide data compilations (Maus and Dimri, 1996;Maus et al., 1998;Pilkington and Todoeschuck, 2004).Magnetic data show a wide range of β values between −1 and −4.8, based on sample areas with scales from <10 km up to >1000 km (Gregotski et al., 1991;Pilkington and Todoeschuck, 1993;Maus and Dimri, 1995;1996;Maus et al., 1997;Bouligand et al., 2009).Although not a necessary requirement for the fractal (or scaling) potential field description (Fedi et al., 1997;Quarta et al., 2000), the spatial distribution of rock properties (density and magnetic susceptibility) that cause these fields are also fractal with behaviours summarized in Bouligand et al. (2009).These particular rock properties have also been described in multifractal terms (Lovejoy et al., 2001;Fedi, 2003;Gettings, 2005).Regardless of whether a simple fractal or multi-fractal model is the better description of the potential field being processed, the main outcome in practical terms from these studies is that a scaling or correlated type of description is more realistic than the earlier, commonly used assumption of an uncorrelated density/susceptibility distribution within the Earth's crust.Fedi et al. (1997) and Quarta et al. (2000) showed that magnetisation distributions other than a scaling one can also produce a scaling magnetic field power spectrum.Their models of uncorrelated distributions of blocks with constant magnetisation are spectrally equivalent to a scaling medium.A "blocky" distribution for magnetisation is not unrealistic since the basic approach of mapping magnetic units from magnetic field data involves assuming regions have similar susceptibilities and are associated with a single lithology.Nevertheless, susceptibility measurements from drill holes and rock sample suites over a wide range of scales support truly scaling properties of the underlying physical properties.Furthermore, evidence of scaling behaviour has been provided through analyses using non-spectral methods.Susceptibility logs were analysed with the rescaled range method and shown to be scaling by Leonardi and Kumpel (1996) The fractal model of magnetic and gravity field data has already been exploited in a variety of uses including kriging of aeromagnetic data using a fractal covariance model (Pilkington et al., 1994), inversion for fractally magnetised source distributions (Maus and Dimri, 1995), Curie depth determination (Maus et al., 1997;Bouligand et al., 2009;Bansal et al., 2011), deriving accurate covariance models for satellite gravity data (Bansal and Dimri, 2005) and synthetic model-making to determine filtering parameters (Pilkington and Cowan, 2006).In the following, we outline a further use of fractals in the preparation of gravity and magnetic data grids prior to FFT-based processing and enhancement algorithms.

Method
Since we have a reliable model for predicting the character of crustal magnetic fields, we can use this knowledge to more realistically "fill in" and extrapolate grids before they are processed.The method we use is that of conditional simulation (Journel and Huijbregts, 1978;Tubman and Crane, 1995).This approach aims to produce synthesized values that have a specified statistical character and that match real values where known.For example, petroleum reservoir simulation to predict fluid flow requires a synthetic porosity distribution that matches those values measured in well-logs (Tubman and Crane, 1995).For the magnetic and gravity cases, the known values occur at all points within the survey area, but the important or conditioning data only occur at the survey edges, external or internal.The simulation approach needs to satisfy two requirements: one is that the simulated field has the same or very similar character to the measured field  and the other is that the simulated values match those at the original survey grid boundaries.
To create a synthetic fractal field, we generate Gaussian white noise (which has a flat power spectrum) with a specified mean (M) and standard deviation (σ ).These values are Fourier transformed and multiplied by f β/2 where β is the required scaling exponent (Pilkington et al., 1994).To account for the distance between the top of the source distribution (usually assumed to be ground level) and the measurement altitude, the synthetic field values are upward continued in the frequency domain.The grid is also projected to the same geomagnetic latitude as the observed data (Blakely, 1996).Finally, inverse Fourier transformation gives the desired field.This synthetic grid is then scaled to the same standard deviation as the measured grid.

Data example
Figure 1 shows a map of aeromagnetic data from the Bonaparte Lake area of British Columbia, Canada.This area was flown by helicopter in 2006 with a line spacing of 420 m and a nominal mean terrain clearance of 125 m (Thomas and Pilkington, 2008).The magnetic data were gridded with an interval of 100 m, resulting in a grid with 200×200 cells.The geology in this region is dominated by the Chilcotin Group basaltic volcanics and the Kamloops Group calc-alkaline volcanic rocks.Minor amounts of Nicola Group volcanic rocks (andesitic to basaltic) occur in the western half of the  area.The major northwest-southeast trending positive magnetic anomaly is interpreted to be caused by Nicola volcanics buried beneath thin (<25 m), Chilcotin Group cover.Elsewhere in the region, Nicola Group volcanics are commonly associated with high-amplitude (commonly >1000 nT and rarely >4000 nT) anomalies.Regions with Kamloops or Chilcotin volcanic rocks generally exhibit a spotty magnetic fabric, lacking any dominant coherent trends.
Figure 2 shows the power spectrum of the field in Fig. 1.A fractal field with similar character to the measured field was then determined based on the observed power spectrum.The parameter β can be estimated from the slope of the long wavelength part of the measured spectrum.This slope will be modified slightly when the fractal field is upward continued to the average terrain clearance (125 m).Matching the power level is easily done once a trial spectrum is plotted.A scaling parameter β = −2.5 was estimated from the observed data spectrum and used for the generated fractal field shown in Fig. 3.This grid has a 256 × 256 cell dimension and is the size of the final extrapolated grid.How much to extend the measured data will be constrained by the shape of the original grid, but should be large enough to allow for a smooth transition from original data out to the edge of the extended grid.In our experience, an extrapolation zone of width equal to 10 % of the original grid dimensions appears to be sufficient (cf., Fig. 1).The power spectrum of the synthetic data grid (Fig. 3) is shown in Fig. 2, where a good match with the observed spectrum is apparent.The observed spectrum shows a fall-off in power at the longest wavelengths (> 20 km) compared to the synthetic values.This difference is often seen and is interpreted to be caused by the limited depth extent of magnetic sources (Maus et al., 1997;Bouligand et al., 2009).Using an FFT approach to compute the fractal field ensures the periodicity requirement at the extended grid edges is met.In order to fill in the region between the original grid (Fig. 1) and the final grid edges with fractal values, there must be no discontinuities along the original grid perimeter.Therefore, using only those fractal field values coincident with the observed field grid (Fig. 1), values are extrapolated outwards using minimum curvature.To ensure the edges of the expanded grid match opposite edges (guaranteeing periodicity), zero values are assigned to the boundary of the expanded grid (black line in Fig. 1).This extrapolated grid, shown in Fig. 4, is then subtracted from the synthetic fractal field grid, resulting in a conditioned field grid that is now zero at the observed field grid edges (Fig. 5).This conditioned field has the desired statistical character within the extrapolated regions.In order for this conditioned field to fit the observed grid at the latter's edges, the observed field is extrapolated outwards using minimum curvature, in the same way as the fractal grid (Fig. 6).Now the conditioned field and expanded observed field are added to give the final extrapolated grid (Fig. 7).The results show that (1) there are no discontinuities at the original grid edges, (2) extrapolated values have the same character as the initial gridded data and (3) this character persists even at large distances from the original grid edges.The original grid in Fig. 1 contains 33162 defined values, while the extrapolated grid in Fig. 7 is 256 × 256, so just less than 50 % of the grid consists of extrapolated values.This level of extrapolation is not extreme, but still shows that fractal extrapolation is an efficient way to preserve the frequency content of the measured data while ensuring periodicity.A flow chart of the complete procedure for fractal grid preparation is given in Fig. 8.

Conclusions
Fractal extrapolation is an alternative to commonly used grid extrapolation techniques.In these techniques, periodicity is obtained by padding the grid with zeros after extrapolation, based on maximum entropy prediction or extrapolation of the grid using minimum curvature or inverse distance gridding.In the case of non-rectangular grids with irregular edges, this can lead to complex algorithms.Mirroring rectangular grids is rather simple, but can be very complex for irregularly shaped grids.The proposed technique is obtained from a series of simple steps that are easy to implement.The technique is independent of the shape of the survey.As seen in the real data example, the extrapolated grid does not show      any discontinuities along the edges of the survey; for other extrapolation techniques this can only be obtained by significant programming effort.The main advantage of the fractal method is that the extrapolated grid has the same frequency content as the original grid.

Fig. 1
Fig. 1 Portion of the Bonaparte Lake aeromagnetic survey in British Columbia.Blue line is edge of gridded survey data.Black rectangle line is the outline of expanded rectangular grid.Coordinates in metres are UTM zone 21.

Fig. 1 .
Fig. 1.Portion of the Bonaparte Lake aeromagnetic survey in British Columbia.Blue line is edge of gridded survey data.Black rectangle line is the outline of expanded rectangular grid.Coordinates in metres are UTM zone 21.
, while Dolan et al. (1998) demonstrated broad-band fractal scaling for several petrophysical logs based on four different calculation methods.

Fig. 2 .
Fig. 2. Spectra of the fractal grid upward continued to 125 m (continuous line) and the measured magnetic data (dotted line).

Fig. 2 .
Fig. 2. Spectra of the fractal grid upward continued to 125 m (continuous line) and the measured magnetic data (dotted line).

Fig. 3 .
Fig. 3. Synthetic fractal magnetic grid upward continued to a height of 125 m.Note the periodicity of the grid.Top and bottom edges are identical as well as left and right edges.Coordinates in metres are UTM zone 21.

Fig. 3 .
Fig. 3. Synthetic fractal magnetic grid upward continued to a height of 125 m.Note the periodicity of the grid.Top and bottom edges are identical as well as left and right edges.Coordinates in metres are UTM zone 21. 19

Fig. 4 .
Fig. 4. Synthetic fractal grid extrapolated to rectangular edges after removal of fractal data located outside the survey area outline in blue in Fig. 1.Edges are zeros.Coordinates in metres are UTM zone 21.

Fig. 4 .
Fig. 4. Synthetic fractal grid extrapolated to rectangular edges after removal of fractal data located outside the survey area outline in blue in Fig. 1.Edges are zeros.Coordinates in metres are UTM zone 21.

20Fig. 5 .
Fig. 5. Conditioned field given by the field in Fig. 4 subtracted from the field in Fig. 3.This conditioned grid is zero at the observed field grid edges (blue line, Fig. 1) and within the survey area.Coordinates in metres are UTM zone 21.

Fig. 5 .
Fig.5.Conditioned field given by the field in Fig.4subtracted from the field in Fig.3.This conditioned grid is zero at the observed field grid edges (blue line, Fig.1) and within the survey area.Coordinates in metres are UTM zone 21.

Fig. 6 .
Fig. 6.Magnetic data after extrapolation to a rectangular area.Measured data remain unchanged.Edges are zero.Coordinates in metres are UTM zone 21. 22

Fig. 7 .
Fig. 7. Final grid after fractal extrapolation.The grid is now periodic.Coordinates in metres are UTM zone 21.

Fig. 7 .
Fig. 7. Final grid after fractal extrapolation.The grid is now periodic.Coordinates in metres are UTM zone 21. 23

Fig. 8 .
Fig. 8. Flow chart summarizing the steps required for fractal grid extrapolation.Where quantities have been plotted, the appropriate figure number is indicated.The expressions used for different quantities (e.g., Fobs) are not used in the text but used here for brevity.

Fig. 8 .
Fig. 8. Flow chart summarizing the steps required for fractal grid extrapolation.Where quantities have been plotted, the appropriate figure number is indicated.The expressions used for different quantities (e.g., F obs ) are not used in the text but used here for brevity.