Nested anisotropic geostatistical gridding of airborne geophysical data Open Access

Practicalities of geophysical surveys dictate that data are acquired at a dense sample spacing along the direction of the traverse, and cost limitations of surveying dictate that traverse lines can be separated by distances of two or more orders of magnitude greater than the along-line or inline spacing. This anisotropic sampling results in aliasing of geophysical features with wavelengths shorter than the line spacing when producing maps or gridded products. This problem can be exacerbated when data are acquired from airborne platforms, a commonly cost-effective choice when surveying large regions. As a result, despite the orthodox advantages of close sampling, maps produced from airborne geophysical surveys are limited to offering spatial resolution at the scale of approximately one half to one third of the line spacing of the survey with the risk of aliasing in the final maps. The visual noise of aliasing is sometimes referred to as boudinage or string-of-beads patterns (Smith and O’Connell, 2005). Furthermore, the geophysical parameters that we wish to measure with a compromised sampling strategy may exhibit azimuthally anisotropic features (e.g., a thin feature not entirely perpendicular to the survey lines). This complicates production of the final products of geophysical surveys (such as total magnetic intensity [TMI] maps) because many gridding algorithms rely upon isotropic gridding techniques. This paper attempts to address both problems with a novel approach to interpolation.

There are several techniques available for gridding geophysical survey data. Inverted layers of electrical conductivity from airborne electromagnetic (AEM) data are often gridded using splines. Akima (1970) proposes an extension of the technique by interpolating irregularly spaced data without the instability inherent in cubic splines. One of the most popular methods for aeromagnetic data is minimum-curvature interpolation, which attempts to produce a smooth surface through the measured data points by balancing data reproduction against a tension parameter (Briggs, 1974). The technique is conceptually simple, but, unlike a cubic spline, it does not always honor the data at measurement locations, and it suffers from boudinage artifacts due to anisotropic sampling (Smith and O’Connell, 2005; Naprstek and Smith, 2019). More generally, Sambridge et al. (1995) propose the natural neighbor method for interpolating earth property measurements. This technique, which relies on Delaunay triangulation, produces a locally referenced grid that reproduces data exactly at the reference points, and it is continuous at all locations other than the reference points. However, the natural neighbor method also suffers from aliasing and boudinage effects.

Outside of the geophysics community, kriging is perhaps one of the most commonly used geostatistical interpolation algorithms (Webster and Oliver, 2007). Kriging relies on spatial covariance models to estimate values at sites remote from the measurement locations and, in the correct conditions, produces the best linear unbiased predictor of the unknown values at those locations.

Some geophysical applications using kriging include Pryet et al. (2011), who use kriging to create 3D conductivity maps of a portion of the Galapagos Islands from stitched 1D layered-earth inversions of an AEM survey; Steinmetz et al. (2015) combine kriged AEM results with other data sets to improve their interpretation of a buried valley near Cuxhaven, Germany; Beucher et al. (2020) use handheld electromagnetic (EM) sensors to estimate peat bog thickness; and Christensen et al. (2015) combine AEM data with ancillary borehole data to update depth-to-bedrock maps for road constructions in Norway. In all cases, the ostensible method used was ordinary kriging of the logarithm of the electrical conductivity (also known as lognormal kriging) using an isotropic covariance function. However, isotropic kriging assumes that the spatial variation of the estimated parameter is the same in every direction. Furthermore, using a single estimate of the covariance of the data over the entire region assumes that the spatial rate of change of the observed data is uniform. In my experience, these assumptions necessarily lead to boudinage structures in the gridded data, even when the conductivity is log-transformed. Christensen et al. (2015) specifically mention using an isotropic variogram for their work, primarily because the data sets studied were directed along narrow corridors. Hansen (1993) compares minimum-curvature gridding of aeromagnetic data to isotropic and anisotropic kriging over the Juan de Fuca Ridge in the Pacific Northwest region and shows that more reasonably interpretable grids can be produced considering anisotropy in the data.

The motivation for this research is due to recent experience with a regional AEM survey designed to map paleochannels in Mid West, Western Australia (WA) (Davis, 2016; Davis et al., 2016a). The logistics of the survey dictated that the flight lines intersected the paleochannels at glancing angles throughout the survey region. This resulted, as expected, in unacceptable boudinage patterns in the gridded results. As a first attempt at addressing the unwanted artifacts, I have used anisotropic variograms at every interpolation point to create a final grid that appeared to ameliorate the problem (Davis, 2016). In a similar time frame, Fouedjio et al. (2016) present a comprehensive approach to nonstationary kriging that included anisotropy. They provide the mathematical framework for the use of anchor points that contain the necessary anisotropy information spread throughout the survey region using a smoothing kernel. The adaptation of the Fouedjio et al. (2016) approach to airborne geophysical surveys is the starting point of this paper.

Due to the volume of data acquired during an airborne geophysical survey, kriging can be computationally expensive. This problem can be alleviated by preconditioning the survey data into blocks that are of the same size as the desired output grid, typically approximately one quarter of the line separation. This has the added advantage of providing a regular grid for the variogram calculations. Distances between observed data are converted to indexes, and the calculation of differences between data samples is reduced to a set of lookup values, conceptually similar to block kriging. The estimation of variograms at anchor points becomes a repetitive task, well suited to parallelization.

This paper revisits the concepts of kriging and the estimation of the covariance weighting functions applied to interpolation. Anisotropy is briefly discussed, as well as the concept of adding additional variogram functions (known as nesting). Building from Fouedjio et al. (2016), I introduce the Laplace approximation method for model selection. This determines which variogram is used at every location in the interpolation region: An isotropic, anisotropic, or nested anisotropic variogram is selected by a using method described by Gregory (2010).

A synthetic example, involving reconstruction of satellite imagery data, is shown to demonstrate the application of the method. This is followed by two case histories, both situated in WA, that relied heavily on airborne geophysical data. The first example is in Mid West WA, and it uses nested anisotropic kriging to determine the depths of paleochannels interpreted from AEM data. The second case history examines AEM and aeromagnetic data in the Pilbara region north of Newman. Both examples show how nested anisotropic kriging improves the final products and eliminates or significantly reduces boudinage string-of-bead effects from airborne geophysical survey grids.

A discussion of the constitutive equations of gridding is required. In this section, I will introduce the system of equations for ordinary kriging, explore the concept of the variogram, and introduce the concepts of stationarity and nonstationarity in geostatistics.

The empirical variogram is defined for the set of observations ${y(x1)$⁠, $y(x2)$⁠, ... $y(xN)}$ as a means to provide estimates of the random function $Y(x)$ at location $x$⁠. This is achieved by providing a substitute variogram function $γ(h)$ that approximately matches the stationary variogram $γ^(h)$ over some range. The choice of function $γ(h)$ depends on the variability seen in the empirical variogram with respect to the lag distance. Implicit in this development are some assumptions about the nature of the random function $Y(x)$⁠. Normally, it is assumed to be locally stationary or quasistationary if the mean value of $Y(x)$ and the covariance between observations vary slowly in a neighborhood around location $x$ (e.g., Webster and Oliver, 2007). Fouedjio et al. (2016) describe a generalized framework that relaxes the assumption of stationarity in $Y(x)$ and allows nonstationary variogram structures to be used for predicting $Y(x)$ at unobserved locations. The essential difference is that covariance between observations can become spatially dependent, even though the same covariance function is used throughout the region of interest.

As in the stationary case, the covariance between variables and their mean are assumed to be slowly varying in the areas of interest for local prediction. Local empirical variograms can be calculated from the observations, and some function $γ(h)$ can be calculated for that local region. However, the parameters of $γ(h)$ that best fit the empirical variograms are allowed to slowly vary. The set of parameters for the estimated variogram functions can be propagated over the survey area by a smoothing function such as a Gaussian kernel.

In these equations, $Ψx0$ and $Λx0$ are the eigenvector matrix and eigenvalues of a geometric parameterization of the survey space surrounding $x0$⁠, where $Ψx0$ operates as a matrix that rotates an ellipse with principle axes $λ1$ and $λ2$ (contained in $Λx0$⁠) by angle $ψx0$⁠. Distances between points located at $x$ and $y$ near $x0$ are calculated according to $λ(λ1(x0)λ2(x0),ψ(x0))=(x−y)Σx0−(x−y)T$⁠, where $λ$ is the transformed distance parameter and $Σx0−$ is the inverse of the deformation matrix in equation 5. As before, it is typical to group the calculated distances into lagged vector distances (as in equation 4).

Geophysical surveys are generally conducted along straight lines separated by a nominal distance $d$ to maximize the coverage and minimize the survey cost. Typically, data collected in the direction along a survey line are sampled much more densely than data across the lines, leading to increased resolution of along-line geophysical parameters compared to any other direction. Interpolating, or gridding, parameters over the entire survey region leads to a necessary loss of resolution over the study area meaning that along-line sections are still an important part of the interpretation process. Typical values of grid cell size are from one fifth to one third of the line spacing distance (Billings and Richards, 2000).

The significance of the preceding argument is that we approach the problem of gridding by starting with a set of observations that is reduced in size and complexity from the original survey data. Grid cell locations are chosen so that the cell centers coincide as much as possible with survey lines. Data are reduced by computing the mean value of survey data that falls within a given grid cell. The cell mean may be arithmetic or geometric, depending on the survey data acquired (e.g., for electrical conductivity values, the appropriate mean is geometric). This has the added benefit that distances between input data for variogram computations (equation 4) are simplified, and lag vectors $h$ for a variogram calculated at position $x$ are indexed based on their position in the regional grid.

The mostly linear alignment of geophysical surveys, and preconditioning of data into a support grid, also gives us an idea not only of where to place anchor points, but also to determine the maximum radius of investigation for the variogram at each anchor point. For a line spacing of $d$⁠, I use grid cell sizes of $d/4$ meaning that there are grid centers at {0, d/4, d/2, 3d/4, d}. Anchor points are chosen to be at distance $d$ apart (in the along- and crossline directions), centered at $d/2$ between and along the lines (resulting in a grid of anchor points). Anchor points are therefore also interpolation centers, and the parameters are easily smoothed across lines. The maximum search radius for computation of the empirical variogram can be found through multiple trials in the ordinary kriging stage; however, good results can be obtained using a radius between $2.5d$ and $3.0d$⁠.

These simplifications greatly reduce the computational burden of computing empirical variograms and estimating the best-fitting parameters for the modeled parameter covariances. Because variogram computation and parameter estimation at each anchor point are considered independently, the computation becomes embarrassingly parallel.

To demonstrate the arguments developed here, an ersatz model replicating data sampling characteristics similar to an airborne geophysical survey is constructed from satellite imagery data of the earth (Google, 2021). The image (Figure 1a) is of a fluvial system near Meekatharra, WA. It is downsampled to a 20.1 × 20.1 km grid with 100 m resolution, and the red, green, blue (RGB) color bands are flattened to a single band, which is rescaled from 0 to 500.

A synthetic data set is constructed from the image data in Figure 1a. Flight lines, trending north–south, are spaced 400 m apart to replicate a moderately spaced survey. The survey samples pixels at 100 m intervals along the flight direction. These are shown with thin black lines in the figure. It is not necessary, for this example, to precondition the data into a regular grid. Figure 1b shows the resulting survey data that is used to construct grids.

Figure 1c shows a detailed portion of box 1 in Figure 1b. The survey data are shown, and the black dots indicate the positions of the anchor points, set in a 400 × 400 m grid. Similarly, Figure 1d shows a detail of box 2; only this time the search radii for a given anchor point are shown. The dark circle shows a search radius of 1600 m, whereas the gold circle has a radius of 800 m.

The method is straightforward. First, a global isotropic variogram is calculated for the entire region, the best-fitting sill $s$ and range $λ$ are estimated for the single variogram (equations 4, 8, and 9), and a map is created for the survey area (equation 1). The empirical variogram is calculated using all survey locations and a maximum search radius of 1600 m. The resulting empirical and modeled variogram are shown with black circles and a black line, respectively, in Figure 2.

Next, isotropic empirical variograms are calculated for every anchor point. The best-fitting regional isotropic parameters are used as starting values to estimate the isotropic best-fitting parameters $s(x)$ and $λ(x)$ at each anchor point in the survey, the results of which are smoothed using equation 10 for every interpolation location. The variogram and best-fitting model for anchor point 1 located at (13,200 m, 15,200 m) are shown in blue, whereas the variogram and best-fitting model for anchor point 2 (5600 m, 5600 m) are shown gold in Figure 2. The gridded maps for these examples are shown in Figure 3a and 3b. Both of these maps exhibit the string-of-beads artifacts for narrow features that are not strictly parallel to the flight lines (Smith and O’Connell, 2005).

Using the best-fitting isotropic parameters as the starting values, the anisotropic parameters ${s(x),λ1(x),λ2(x),ψ(x)}$ are estimated using equations 10 and 12. The semimajor radius starting values are copied from $λ(x)$⁠, and the semiminor starting values are some factor (e.g., 0.8) of $λ(x)$⁠. Angles $ψ(x)$ are proposed by examining the previous maps or by setting an expected angle of anisotropy. The best-fitting parameters are smoothed throughout the survey region using equations 10 and 12. The first nonstationary anisotropic map is then produced from the smoothed parameters (Figure 3c). We see that many of the string-of-beads features have been eliminated in this figure as a result of accounting for anisotropy.

Nested anisotropic parameters are estimated at each anchor point by building on the results from the anisotropy step. Starting sills $s1(x)$ and $s2(x)$ are given the value of $2s(x)/5$⁠, distance parameters $λ1,1(x)$ and $λ2,1(x)$ are recycled from $λ1(x)$ and $λ2(x)$⁠, whereas distance parameters $λ1,2(x)$ and $λ2,2(x)$ are proposed with values of one half and one third the maximum search radius, respectively. The starting angles are $ψ1(x)=ψ(x)$ and $ψ2(x)=ψ(x)+π/2$⁠. An example of the best-fitting nested anisotropic ellipses, and the predicted variogram, for the anchor point at (13,200 m, 15,200 m) are shown in Figure 4b. The empirical variogram for this point is displayed in Figure 4a, showing distinct anisotropy. The estimated parameters are smoothed using equations 10 and 12 as before, and the nested anisotropic map is produced for the survey area (Figure 3d).

Recalling that equation 13 is used to calculate the odds ratio of preferring one model over another, the final step is to apply the equation throughout the entire survey area to select the preferred variogram models at each estimation point. That is, for point $x$⁠, equation 13 is used to determine if the nested anisotropic model is preferable to the single anisotropic model and, if not, if the single anisotropic model is preferable to the simple isotropic model. Once the preferred model is determined for each point, a composite map is produced that only accepts an increased model complexity where it is warranted by the data. The result of the analysis is shown in Figure 3e. The original image is included for reference in Figure 3f.

One advantage of using geostatistical interpolation methods such as kriging is that we can also obtain the variance of the estimated geophysical parameters. Figure 5a and 5b shows the kriging values and error estimates for the ordinary isotropic kriging estimates, whereas Figure 5c and 5d presents the globally selected local anisotropic kriging estimates. As we expect from ordinary kriging, the error in the estimates is lowest along the flight lines and increases between them. The error in the estimates between lines is quite uniform because the covariance function is the same throughout, and all interpolation points between lines have approximately the same number of points contributing to the estimate. There is much more variation in the error of the interpolated parameters when we consider anisotropy (Figure 5c and 5d). The spatial distribution of the error is due to anisotropy. The spatially dependent, anisotropic covariance functions decrease the estimation errors in some locations while increasing it in others. Figure 6 shows the distribution of errors for the isotropic ordinary kriging (blue) and the local anisotropic interpolation (gold). The mean and median values of error for the anisotropic kriging are lower than for ordinary kriging indicating that anisotropic kriging generally better resolves the original image.

It is also useful to examine the total time taken to calculate each step. Keeping in mind that each increase in complexity relies upon the previous step, the total time for each phase of the estimation process is shown in Table 1. All calculations were conducted in parallel on a laptop computer with six processors. As expected, the processing time of each step increases for diminishing reduction in error.

Having demonstrated the interpolation procedure and results with synthetic data, it is appropriate to further illustrate with case histories. In this section, I will use two case histories to show how anisotropic gridding is applied to airborne magnetic and electromagnetic data, and also to interpret products derived from geophysical survey data.

The Murchison Paleochannels Project, funded by the government of WA’s Royalties for Regions initiative, used AEM geophysics to map paleovalleys at a regional scale to assess and provide advice on water information, water availability, and supply options (Government of Western Australia, Department of Water, 2015). The survey comprised more than 14,000 line-km of data flown at 4000 m spacing in the Mid West region of WA. The survey data were inverted for subsurface electrical conductivity and interpreted for paleochannel depth, an example of which is shown for line 134,001 in Figure 7. (Davis et al., 2016b, 2020). The cross section illustrates the depth of the paleochannel clearly as the interface between the high and low conductors. The paleochannel depths picked from the section are shown in white and the pregridded (1000 m) depths are shown in red.

Because the survey lines were oriented to minimize ferry and survey costs, many lines were flown oblique to the strike of the paleochannel thalwegs, which is suboptimal. The gridding produced with isotropic kriging resulted in a blobby paleochannel depth map with clear gridding artifacts between flight lines. The artifacts are obvious in Figure 8a, particularly in the east of the region. The white line in the figure indicates the location and orientation of line 134,001. Flight lines are shown with narrow black lines, illustrating that they are not perpendicular to the paleochannel thalweg.

The problem of interpolating for continuous palaeochannel depth inspired the idea of applying anisotropic kriging to geophysical survey data. Figure 8b shows the gridded paleochannel depth map following all the processing steps of the previous section. The result is a depth map with continuous linear trends connecting the paleochannel networks together and draining toward the northwest. All data used for the Murchison Palechannels Project may be found in the CSIRO Data Access Portal (Davis, 2020).

The second case study involves AEM and TMI data for groundwater exploration and managed aquifer recharge development in the Pilbara region of WA. One of the goals of the Transforming Agriculture in the Pilbara (TAP) Project, sponsored by a Royalties for Regions and managed by the WA Department of Primary Industries and Regional Development, is to map near-surface electrical conductivity for aquifer detection, characterization, and groundwater resource development (Government of Western Australia, 2018). Specific targets include the delineation of paleochannel flow paths and deposition patterns, alluvial fan transported materials, fault detection, and cover depth mapping. The TAP airborne survey was conducted along east–west-trending flight lines separated by 400 m chosen to provide good regional coverage while still providing hectare-scale resolution of near-surface features (Davis et al., 2021a, 2021b).

The postprocessed AEM data were inverted using a 30-layer spatially constrained inversion package (Auken et al., 2015). The logarithm of the inversion model conductivities was then pregridded into 100 × 100 m cells and run through the procedure described previously (details in Davis et al., 2021b). Figure 9a and 9c shows the results of the isotropic kriging for a shallow (6.6–9.2 m) layer (Figure 9a) and a deep (35.4–40.6 m) layer (Figure 9c). Figure 9b and 9d shows the same layers using the best-fitting results of the isotropic, anisotropic, and nested anisotropic kriging methods. Between-line boudinage artifacts are reduced using the advanced methods. The shallow conductivity map of Figure 9b shows smoothly varying conductivity indicating near-surface features of an alluvial fan that transported regolith materials from the Fortescue paleoriver system across the Wittenoom Formation. The fan trends northward toward Fortescue Marsh, but it spreads outward to the east across the Formation. Higher conductivities are attributed to clay-rich sediments, although presently there is little drilling in the area to support this. Improvements in crossline connections between features can be seen by comparing the areas in rectangles in Figure 9a. Figure 9d shows a deeper cover sequence across the Wittenoom Formation. The abrupt change in conductivity in the east marks the boundary of the Formation from the Pinjian Chert Breccia and is governed by the Fortescue River Fault. There is evidence of the Fortescue Palaeoriver emptying to the north from approximately 182,000 m east at the bottom of the figure. The north-trending linear feature near the center of the figure is from powerlines along the Marble Bar Road. The round, low-conductivity feature at (180,000 m; 7,445,000 m; marked with a circle in Figure 9b) has been identified as a possible groundwater resource. As with the shallow conductivity map, features in Figure 9d are generally smoothly varying even across lines. Linear features are enhanced (see the arrow in Figure 9b), and boundaries with low-conductivity features are sharper (see the rectangle in Figure 9b). Details of the interpretation of this AEM survey can be found in Davis et al. (2021a).

The final example of anisotropic nested kriging relates to TMI data acquired as part of the TAP Project AEM survey. The TMI data were mean-stripped and gridded to 100 × 100 m cells, resulting in an anomaly map. Figure 10a and 10b shows the results of isotropic gridding and probability-selected nested anisotropic kriging, respectively. Figure 10–10f shows magnified regions, marked with the black boxes in Figure 10a. Figure 10c and 10e shows the magnified regions 1 and 2 for the isotropic kriging, whereas Figure 10d and 10f shows the magnified regions for the anisotropic gridding. The regional bedrock geology of the area is also shown in Figure 10a (State of Western Australia [Department of Mines and Petroleum], 2020).

As observed in the conductivity maps from the AEM data, the anisotropic kriged maps of the TMI anomaly are much smoother than the isotropic maps and the between-line blobs are minimized or eliminated. Figure 10c and 10d shows a magnified region to the west of the survey, which highlights the interface between Boolgeeda Iron Formations to the south and the Wittenoom Formation to the north. The northwest–southeast-striking linear feature is the Poonda Fault, a major break in the geology. Although both gridding techniques show the fault, the break in Figure 10d is slightly smoother. Furthermore, stippling of the grids is reduced in Figure 10d compared to Figure 10c. Perhaps more interesting is the presence of another fault in the center of these images, striking just east of north from the center of the figures (marked with an arrow in Figure 10c. In Figure 10d, the feature is slightly enhanced by the anisotropic gridding. Figure 10e and 10f illustrates another interesting feature in the southwest of the survey, north of Jimblebar. These images show the intersection of two dikes: one striking to the northeast and the other just south of east. Although the dikes may be interpreted in the isotropic gridding, boudinage effects make interpretation more difficult. In contrast, the anisotropic kriging significantly reduces boudinage effects and reduces the ambiguity of interpretation. There are some gridding artifacts still present in Figure 10f, but they are reduced markedly from the ordinary kriging results. All grids and data used to make the maps are available from the CSIRO Data Access Portal (Davis and Department of Primary Industries and Regional Development, 2020a, 2020b). Figures showing the global isotropic variogram and local variograms for the anchor point shown with a black dot in Figure 10a, 10b, 10e, and 10f can be found in the supplemental material, which can be accessed through the following links: S1, S2, S3, S4, S5, S6, S7, and S8).

Airborne geophysical surveys are noteworthy in that they offer greater resolution along the direction of flight than across it. For this reason, it is best to plan surveys at right angles to the expected geologic strike of the features to be detected. This is not always possible or practical during the planning stages of the survey, with the repercussion that grids and maps can produce boudinage artifacts in the final products. This paper focuses on the development of nested anisotropic kriging to sharpen features and improve between-line continuity in geophysical survey data.

The anisotropic gridding algorithm proposed in this paper is based on a method of nonstationary kriging that has been extended to airborne geophysical data by providing a simple method for systematically refining grid products while also taking advantage of the highly regular spacing of the surveys. I have also developed the concept of including higher order nesting of anisotropic variograms, and I provide an objective method to select which variogram structure to use at each location of the area to be interpolated. This extension comes at the cost of increased computing time, but the results are an improvement on traditional gridding.

In the examples shown, the method cleans up crossline structures remarkably well and often reduces or removes the boudinage effects from the grids. This has been demonstrated with the paleochannel depth map in the Mid West region of WA and in the AEM and magnetics maps in the TAP project. The gridded results can be used for visual interpretation and also for structure mapping. In the TAP EM data, the near-surface and deeper features have greater lateral continuity. There are still artifacts present in the results, although they have been reduced. The results show that greater model complexity can reduce the MS error in the predicted results (at least for the satellite imagery example shown).

Although it has not been demonstrated here, it would be possible to use the techniques developed on other airborne geophysical data sets such as radiometric or gravity surveys. In principle, the methods would be the same: Transform the data to reflect the expected change in parameters with respect to distance (such as taking the log of the data), reduce the data set to the resolution desired in the final products, and compute the grids using the described procedure. Along-line density and regular spacing can be used to our advantage because the data can be indexed for the variogram calculations. The method is straightforward and easy to apply in simple consecutive steps. By pregridding and indexing the data, another potential application of anisotropic variogram analysis becomes apparent, and that is in the field of 3D inversion. If, say, AEM survey data were to be converted to conductivity-depth information through the use of one of the widely available fast transform techniques, the anisotropic analysis could be conducted to improve the spatial regularization required in 3D modeling. This also could be done for 3D inversion of airborne magnetic data. It is expected that the anisotropic analysis would improve the resulting 3D models.

There are parameters that need to be adjusted and tweaked during gridding, which will require some experimentation. For example, there are plenty of different functions available to model the variograms. Although I have found that the Whittle function offers reasonable results for satellite imagery, AEM data, and magnetics, other functions may perform better in certain circumstances. Selection of variogram functions can be accomplished through the model selection equations given here with minimal changes. If this is desired, the model selection can be calculated starting at the isotropic fitting of the anchor points. The Whittle function offers low variance at short lag distances and an approximately exponential increase in variance at larger lags, which appear to be appropriate for EM (which was inverted using horizontal and vertical spatial regularization) and magnetic data.

There are other parameters that must be adjusted depending on the line spacing, the geophysical parameters of interest, and the expected geologic variation; for example, anchor point spacing, maximum search radius for the variograms, and characteristic distance for the Gaussian kernels all require at least a preliminary investigation for each data set considered. However, good results can be achieved with anchor point spacing equivalent to the line separation, a maximum search distance of three-to-four-line spacings, and a characteristic length of one half the line spacing gives reasonable results. Estimation of anisotropy angles can be improved by seeding the starting parameters $ϕ(x)$⁠. In each example shown, however, a common starting angle was used for every anchor point. As described previously in this paper, the selection of starting parameters is straightforward. Although each new step in the process does rely upon the previous one, so that some experimentation must be made with the initial parameters, it is easily envisaged that the whole process can be automated for electromagnetic and aeromagnetic surveys. The resulting enhancement of across-line features and the sharpening of edges and regions of increased variation will aid geologists and geophysicists in their interpretation of geophysical survey data, despite the increased computational cost.

I thank CSIRO Mineral Resources, CSIRO Land and Water, and the Deep Earth Imaging Future Science Platform for support in publishing this article. F. Fouedjio, T. Munday, and A. King are acknowledged for interesting and stimulating discussions about gridding. All software was developed in the MATLAB programming environment from first principles, with the exceptions of the least-squares parameter searching functions.

Data is available in CSIRO Data Access Portal, which is referenced in this paper.

A biography and photograph of the author are not available.