Euler deconvolution of potential field data

In this tutorial (second article in the TLE series) we'll talk about a widely used method of interpretation for potential-field data called Euler deconvolution. Our goal is to demonstrate its usefulness and, most importantly, call attention to some pitfalls encountered in the interpretation of the results. The code and synthetic data required to reproduce our results and figures can be found in the accompanying IPython notebooks at figshare or GitHub. The notebooks also expand the analysis presented in this article. We encourage you to download the data and try it on your software of choice. For this tutorial we'll use the implementation in the open-source Python package Fatiando a Terra.

Brief overview

Euler deconvolution was first developed by Thompson (1982) and later extended by Reid et al. (1990). Since then, it has been adapted and improved upon by Keating (1998), Mushayandebvu et al. (2004), and many others. This popularity is largely due to its great simplicity of implementation and use, making it the tool of choice for a quick initial interpretation.

In many cases, maps of gravity and magnetic data (and transformations thereof) provide good constraints on the horizontal location of an anomaly source. Euler deconvolution adds an extra dimension to the interpretation. It estimates a set of (x, y, z) points that, ideally, fall inside the source of the anomaly. Euler deconvolution requires the x-, y-, and z-derivatives of the data and a parameter called the structural index (SI). The SI is an integer number that is related to the homogeneity of the potential field and varies for different fields and source types (Stavrev and Reid, 2007). For example, in the case of total field magnetic anomaly data, a sphere is represented by an SI of 3, whereas a dike is represented by an SI of 1. There are methods that can estimate the SI and we refer the reader to Barbosa et al. (1999) and Melo et al. (2013) for more information.

An example with synthetic data

The best way to fully grasp the usefulness and limitations of any geophysical method is to apply it in a controlled setting using synthetic data. To make things more interesting, we'll apply the Euler deconvolution on data generated by a model that differs slightly from ideal sources (Figure 1).

Figure 1: Our model (a) and synthetic total field anomaly data (b). The model simulates a dike (blue), a sill (green), and an intrusive body (red). The data is corrupted with 5 nT pseudo-random Gaussian noise.

First, we need to calculate the 3 directional derivatives of the data. For gridded data, this can be done with the Fourier transform using module fourier from Fatiando a Terra. Assuming that the data has been loaded into Numpy arrays x, y, z with the respective x-, y-, z-coordinates of each data point and the array data with our magnetic data (see the accompanying notebook):

   from fatiando.gravmag import fourier
   from fatiando import utils
   # Need data in SI units
   data_si = utils.nt2si(data)
   # shape = (ny, nx) number of points in y and x
   dx = fourier.derivx(x, y, data_si, shape)
   dy = fourier.derivy(x, y, data_si, shape)
   dz = fourier.derivz(x, y, data_si, shape)

Now that we have calculated the derivatives, we have to create a solver to run the deconvolution. Fatiando provides the Classic solver which uses the whole dataset to estimate a single point (source location). The usual strategy for multiple sources, introduced by Reid et al. (1990), is to use a moving window over the data. A point is estimated for each window using the data that fall inside it, producing a large number of solutions. The challenge then becomes separating the good solutions from the spurious ones. FitzGerald et al. (2004) provide an overview of different selection criteria. Fatiando provides the MovingWindow solver that can run a given Euler solver on a moving window (see the accompanying notebook for more details on the selection criterion used):

   from fatiando.gravmag.euler import Classic, MovingWindow
   classic = Classic(x, y, z, data_si, dx, dy, dz, struct_index)
   solver = MovingWindow(classic, windows=(50, 50), size=(1000, 1000))
   solver.fit()

Calling solver.fit() runs the deconvolution using the specified 50 x 50 windows of size 1 km. The estimated points are stored in the solver.estimate_ variable.

For comparison, we ran the deconvolution using various different combinations of structural index and window size. Figure 2 shows scatter plots of the solutions (the typical way Euler solutions are presented). Notice how changing the window size influences the spread of the solutions and can affect our interpretation. For example, for a 1 km window and structural index 3, we cannot tell apart the two model bodies in the left. Furthermore, increasing the structural index will increase the depths of the solutions.

Figure 2: Euler deconvolution solutions for varying structural index (SI) and moving window size.

We can view these solutions in 3D using the Fatiando module myv, a wrapper for the Mayavi visualization software (http://code.enthought.com/projects/mayavi):

   from fatiando.vis import myv
   myv.figure()
   myv.points(solver.estimate_)
   myv.axes(myv.outline(model_bounds))
   myv.show()

This will create something similar to Figure 3. Notice that the Euler solutions do fall approximately inside the sources. However, it is curious how the solutions make arch-like structures that have no relation to the actual forms of the sources. In fact, Silva et al. (2001) have shown that Euler solutions are biased and that the bias depends on the moving-window scheme and selection criterion used.

Figure 3: 3D view of the Euler deconvolution solutions (black dots) for a moving window of 3 km and structural index 3.

Last few words of caution

We repeat here for emphasis: Euler deconvolution solutions do not represent the three-dimensional outline of the sources. They indicate, at most, an approximate source location. It is crucial to remember that Euler solutions are subject to non-uniqueness and bias just like any other geophysical inverse problem. Tests using synthetic data are indispensable to assert the plausibility of our interpretations.

External links

• IPython notebook that generates the synthetic data - Uses the forward modeling in Fatiando a Terra to create the synthetic data.
• IPython notebook for Euler deconvolution - Notebook that runs the Euler deconvolution on the synthetic data and produces the figures. Has a more detailed explanation of the method and generates a few more results.
• Tutorials repository - The GitHub repository for the Geophysical tutorials section of TLE has the notebooks from all tutorials published so far.
• GitHub repository - Has the text of the tutorial, IPython notebook (source code), synthetic data, and figures.
• figshare repository - A persistent archive for the synthetic data. Use this to download the data and for citation purposes.
• Fatiando a Terra - The open-source Python package used to produce the data and results in this tutorial.
• Corresponding author: Leonardo Uieda, leouiedagmail.com

find literature about Euler deconvolution of potential field data