Notebook

Euler deconvolution of synthetic magnetic data¶

Euler deconvolution of potential fields is a widely used interpretation technique. This is, in no small part, due to its simplicity and ease of implementation and use. This IPython notebook will apply standard Euler deconvolution on synthetic total field magnetic anomaly data in an attempt to highlight some common misconceptions in the interpretation of the results.

We use the Euler deconvolution implemented in the open-source Python library Fatiando a Terra. The library also provides routines for plotting, forward modeling, calculating derivatives through FFT (required for Euler deconvolution) and more.

Euler deconvolution is based on Euler's homogeneity equation,

$$ (x_i - x_0)\dfrac{\partial f_i}{\partial x} + (y_i - y_0)\dfrac{\partial f_i}{\partial y} + (z_i - z_0)\dfrac{\partial f_i}{\partial z} = \eta (b - f_i), $$

where $(x_i, y_i, z_i)$ are the coordinates of observation $f_i$ (the total field anomaly), $(x_0, y_0, z_0)$ is a point on the source (where that point is depends on the nature of the source), $b$ is a constant base level (shift in the data), and $\eta$ is called the structural index (SI) and is related to the nature of the source.

Given the observed total field anomaly $f_i$, its derivatives, and the structural index, we can solve the above equation for the unknown point $(x_0, y_0, z_0)$ and base level $b$. This is implemented as the class Classic from the euler module of Fatiando a Terra.

However, this is only valid if there is only one source. When there are multiple sources, the most common approach is to select a small window over the data and solve using the observations that fall inside the window. The window is then moved and another solution is found. Thus, we end up with a set of solutions. This is called a moving window scheme and is implemented in class MovingWindow in the euler module.

A drawback of this approach is the sheer number of solutions generated, many of which are spurious. Thus, the solutions have to be filtered so that only consistent solutions are kept. How to best perform this filtering is an active field of research and will not be covered here. The MovingWindow class implements a simple filter that keeps only a percentage of the solutions, ranked by the average uncertainty of the estimate. This is an approach similar to that adopted by Beiki and Pedersen (2010).

Returning to the subject of the structural index, it is a number that characterizes the decay of the potential field with distance. The structural indexes of some ideal sources (sphere, cylinder, vertical dike, etc) have been calculated (e.g., Reid et al., 1990). However, most, if not all, geologic sources cannot be properly approximated by such ideal sources. There are methods that have been developed to estimate the structural index during Euler deconvolution. However, that is outside the scope of our study and will not be covered here.

We will show an application of Euler deconvolution using the moving window scheme on synthetic total field magnetic anomaly data. This application will illustrate how the deconvolution results behave for non-ideal sources. We will use the results obtained to highlight how the output of Euler deconvolution should and, most importantly, should not be interpreted.

Loading the data¶

The first step is to load (import) the required libraries.

In [1]:

%matplotlib inline
import numpy as np
import json
import cPickle as pickle
from IPython.display import Image
from fatiando.gravmag.euler import Classic, ExpandingWindow, MovingWindow
from fatiando.gravmag import fourier, sphere, polyprism
from fatiando.mesher import PolygonalPrism, Sphere
from fatiando.vis import mpl, myv
from fatiando import utils
import fatiando
print("Using Fatiando a Terra version {0}".format((fatiando.version)))
mpl.rc('font', size=9)

Using Fatiando a Terra version 0.2

Now, we load the synthetic data and accompanying metadata using modules json and numpy.

In [2]:

x, y, z, tf = np.loadtxt('data/synthetic_data.txt', unpack=True)
with open('data/metadata.json') as f:
    metadata = json.load(f)
    area = metadata['area']
    bounds = metadata['bounds']
    shape = metadata['shape']
    inc = metadata['inc']
    dec = metadata['dec']

In [3]:

metadata

Out[3]:

{u'area': [5000, 25000, 5000, 25000],
 u'bounds': [0, 30000, 0, 30000, 0, 5000],
 u'dec': 30,
 u'inc': -15,
 u'shape': [100, 100]}

Module mpl of Fatiando a Terra interfaces with the 2D plotting library matplotlib. We can use it to generate a filled contour plot of the data. Note that the x coordinates (array x) are North-South and y (array y) are East-West.

In [4]:

mpl.figure(figsize=(3.33, 2.7))
mpl.axis('scaled')
mpl.contourf(y, x, tf, shape, 60, cmap=mpl.cm.RdBu_r)
mpl.colorbar().set_label('nT')
mpl.m2km()
mpl.xlabel('East (km)')
mpl.ylabel('North (km)')
mpl.savefig('fig/data.png', dpi=600)

The polygonal prism model used to generate the data is stored in a Python pickle file.

In [5]:

with open('data/model.pickle') as f:
    model = pickle.load(f)
    dike, sill, batholith = model

We can use module myv to interface with 3D plotting library Mayavi and render the model. Notice that the model is composed of 3 different structures: a slightly crooked vertical dike, an irregular sill and a structure resembling an intrusive body, like a batholith.

In [6]:

scene = myv.figure(size=(1000, 600))
myv.polyprisms([dike], linewidth=1, color=(0, 0, 1))
myv.polyprisms([sill], linewidth=1, color=(0, 1, 0))
myv.polyprisms([batholith], linewidth=1, color=(1, 0, 0))
ax = myv.axes(myv.outline(bounds), ranges=[b*0.001 for b in bounds], nlabels=3, fmt='%.0f')
ax.axes.x_label, ax.axes.y_label, ax.axes.z_label = 'x', 'y', 'z'
ax.axes.font_factor = 2
myv.wall_north(bounds)
myv.wall_bottom(bounds)
scene.scene.camera.position = [-37282.427190894407, -35382.676458226633, -36760.719551179718]
scene.scene.camera.focal_point = [13286.06353429469, 14199.525065081265, 5611.7451127001987]
scene.scene.camera.view_angle = 19.2
scene.scene.camera.view_up = [0.37050322249789591, 0.35546774432369316, -0.85812006436401433]
scene.scene.camera.clipping_range = [42859.55137481264, 132512.66650687021]
scene.scene.camera.compute_view_plane_normal()
scene.scene.render()
myv.savefig('fig/model.png')
myv.show()
Image(filename='fig/model.png')

Out[6]:

Lets have a look at some of the properties of the model:

In [7]:

print('Dike:')
print(' top (m): {}'.format(dike.z1))
print(' bottom (m): {}'.format(dike.z2))
print(' magnetization (SI): {}'.format(dike.props['magnetization']))

print('Sill:')
print(' top (m): {}'.format(sill.z1))
print(' bottom (m): {}'.format(sill.z2))
print(' magnetization (SI): {}'.format(sill.props['magnetization']))

print('Batholith:')
print(' top (m): {}'.format(batholith.z1))
print(' bottom (m): {}'.format(batholith.z2))
print(' magnetization (SI): {}'.format(batholith.props['magnetization']))

Dike:
 top (m): 0.0
 bottom (m): 5000.0
 magnetization (SI): 10
Sill:
 top (m): 1000.0
 bottom (m): 1500.0
 magnetization (SI): 10
Batholith:
 top (m): 500.0
 bottom (m): 4000.0
 magnetization (SI): 2

Structural index¶

The structural index (SI) characterizes the nature of the source of the anomaly. It will vary depending on the potential field, hence, the SI for a type of source will not be same for gravity and magnetic fields. SI values have been calculated for ideal sources. After Reid et al. (1990), the SIs for the magnetic case are

Source	SI
Dike	1
Sill	1
Sphere	3
Vertical pipe	2

We'll assume that the SI values that best represent our model are: 1 for the dike and sill, 3 for the batholith.

Calculating derivatives¶

Euler deconvolution requires the derivatives of the total-field anomaly. They could be directly measured (if possible) or calculated numerically. The important thing is that the derivatives need to be calculated/measured at the same observation points $(x_i, y_i, z_i)$ as the data. We will use module fourier of Fatiando a Terra to calculate derivatives using the FFT. Note that this requires the total field anomaly to be in the Système International units and on a regular grid. Output derivatives will be in units of T/m.

In [8]:

tf_si = utils.nt2si(tf)
xderiv = fourier.derivx(x, y, tf_si, shape)
yderiv = fourier.derivy(x, y, tf_si, shape)
zderiv = fourier.derivz(x, y, tf_si, shape)

In [9]:

mpl.figure(figsize=(10, 4))
mpl.subplots_adjust()
levels = 30
cbargs = dict(orientation='horizontal')

ax = mpl.subplot(1, 3, 1)
mpl.axis('scaled')
mpl.title('x derivative')
mpl.contourf(y, x, xderiv, shape, levels, cmap=mpl.cm.RdBu_r)
mpl.colorbar(**cbargs).set_label('T/m')
mpl.m2km()
mpl.xlabel('East (km)')
mpl.ylabel('North (km)')

ax = mpl.subplot(1, 3, 2)
mpl.axis('scaled')
mpl.title('y derivative')
mpl.contourf(y, x, yderiv, shape, levels, cmap=mpl.cm.RdBu_r)
mpl.colorbar(**cbargs).set_label('T/m')
mpl.m2km()
mpl.xlabel('East (km)')
mpl.ylabel('North (km)')

ax = mpl.subplot(1, 3, 3)
mpl.axis('scaled')
mpl.title('z derivative')
mpl.contourf(y, x, zderiv, shape, levels, cmap=mpl.cm.RdBu_r)
cb = mpl.colorbar(**cbargs).set_label('T/m')
mpl.m2km()
mpl.xlabel('East (km)')
mpl.ylabel('North (km)')

Out[9]:

<matplotlib.text.Text at 0x7711ed0>

Euler deconvolution¶

Now we have all that we need to perform an Euler deconvolution with a moving-window scheme. We will use varying structural indexes (SI) and moving window sizes to illustrate the dependence of the results on these parameters.

First, we'll use a 1 km window and vary the structural index from 1 to 3. We'll store the results in a Python list for later use. The parameter keep=0.3 tells the moving window scheme to keep only the best 30% of solutions using the criterion of Beiki and Pedersen (2010). Note that we have to pass the total field anomaly and the derivatives in Système International units.

In [10]:

euler_s1 = [MovingWindow(Classic(x, y, z, tf_si, xderiv, yderiv, zderiv, i), 
                         windows=(50, 50), size=(1000, 1000), keep=0.3).fit()
            for i in [1, 2, 3]]

Now we do the same for a 3 km window.

In [11]:

euler_s3 = [MovingWindow(Classic(x, y, z, tf_si, xderiv, yderiv, zderiv, i), 
                         windows=(50, 50), size=(3000, 3000), keep=0.3).fit()
            for i in [1, 2, 3]]

For the sake of completeness, let's run on a 5 km window as well.

In [12]:

euler_s5 = [MovingWindow(Classic(x, y, z, tf_si, xderiv, yderiv, zderiv, i), 
                         windows=(50, 50), size=(5000, 5000), keep=0.3).fit()
            for i in [1, 2, 3]]

The typical way that deconvolution results are show is using a scatter map with color (or point size) representing the depth (z coordinate). We'll make 6 plots, one for each combination of window size and structural index. While we plot, we'll also filter out unrealistic results that lie above the data.

In [13]:

# Get the maximum value of z that we got from all runs.
# We'll use this to set the limits of the color bars
maxz = max(p[2] for euler in [euler_s1, euler_s3] for e in euler for p in e.estimate_)
labels = ['(a)', '(b)', '(c)', '(d)', '(e)', '(f)']
fig = mpl.figure(figsize=(7, 4.8))
mpl.subplots_adjust(left=0.06, top=0.98, bottom=0.05, wspace=0, hspace=0.02)
for j, euler in enumerate([euler_s1, euler_s3]):
    for i, e in enumerate(euler):
        xc, yc, zc = e.estimate_.T
        xc, yc, zc = xc[zc > 0], yc[zc > 0], zc[zc > 0]
        ax = mpl.subplot(2, 3, 3*j + i + 1)
        mpl.axis('scaled')
        mpl.text(6000, 23000, labels[3*j + i], fontsize=12, fontweight='bold')
        mpl.title('SI={0} | Window={1:.0f} km'.format(e.euler.model['structural_index'], e.size[0]*0.001))
        #for body in model:
        #    mpl.polygon(body, xy2ne=True, style='-b', linewidth=1.5)
        plot = mpl.scatter(yc, xc, s=30, c=zc*0.001, cmap=mpl.cm.RdYlGn_r, vmin=0, vmax=maxz*0.001)
        mpl.set_area(area)
        mpl.m2km()
        if i in [1, 2, 4, 5]:
            ax.set_yticklabels([])
        if j == 1:
            mpl.xlabel('East (km)')
            ax.set_xticklabels(['{:.0f}'.format(t*0.001) for t in ax.get_xticks()[:-1]])
        if j in [0]:
            ax.set_xticklabels([])
        if i in [0, 3]:
            mpl.ylabel('North (km)')
fig.colorbar(plot, cax=mpl.axes([0.91, 0.09, 0.03, 0.85]), orientation='vertical').set_label('Depth (km)')
mpl.savefig('fig/euler-solutions-low.png')
mpl.savefig('fig/euler-solutions.png', dpi=600)

In [14]:

# Make the figure for the 5km window separately because it won't go in the tutorial
maxz = max(p[2] for euler in [euler_s5] for e in euler for p in e.estimate_)
fig = mpl.figure(figsize=(12, 3.5))
for i, e in enumerate(euler_s5):
        xc, yc, zc = e.estimate_.T
        xc, yc, zc = xc[zc > 0], yc[zc > 0], zc[zc > 0]
        ax = mpl.subplot(1, 3, i + 1)
        mpl.axis('scaled')
        mpl.title('SI={0} | Window={1:.0f} km'.format(e.euler.model['structural_index'], e.size[0]*0.001))
        plot = mpl.scatter(yc, xc, s=30, c=zc*0.001, cmap=mpl.cm.RdYlGn_r, vmin=0, vmax=maxz*0.001)
        mpl.set_area(area)
        mpl.m2km()
        mpl.xlabel('East (km)')
        mpl.ylabel('North (km)')
fig.colorbar(plot, cax=mpl.axes([0.91, 0.09, 0.03, 0.85]), orientation='vertical').set_label('Depth (km)')

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 sill from the intrusive body. The solutions get even more tangled when using a 5 km window.

Varying the structural index will give different depths for the solutions. A general trend is that, the larger the structural index, the deeper the solutions.

Most importantly, we call attention to the fact that the Euler deconvolution solutions do not represent the three-dimensional shape of the sources. This is even more evident when we view the solutions in 3D along with the original model. For convenience, we'll plot only the solutions for a 3 km window and structural index 3.

In [15]:

scene = myv.figure(size=(1000, 600))
myv.polyprisms(model, opacity=0.3, linewidth=2)
myv.points([e for e in euler_s3[2].estimate_ if e[2] > 0])
ax = myv.axes(myv.outline(bounds), ranges=[b*0.001 for b in bounds], nlabels=3, fmt='%.0f')
ax.axes.x_label, ax.axes.y_label, ax.axes.z_label = 'x', 'y', 'z'
ax.axes.font_factor = 2
myv.wall_north(bounds)
myv.wall_bottom(bounds)
scene.scene.camera.position = [-10753.852296021096, -59325.007533753123, -22266.66194804628]
scene.scene.camera.focal_point = [16722.015964135739, 13784.075401721391, 4529.0411289607291]
scene.scene.camera.view_angle = 19.2
scene.scene.camera.view_up = [0.12859903242950413, 0.29830909042254711, -0.94576634293543571]
scene.scene.camera.clipping_range = [43526.915602369831, 131500.63847030781]
scene.scene.camera.compute_view_plane_normal()
scene.scene.render()
myv.savefig('fig/euler-solutions-3d.png')
scene.scene.camera.position = [61783.513706029873, -51534.690503972124, -18750.237814442207]
scene.scene.camera.focal_point = [14324.659515610345, 12669.434351699714, 2305.8434486142896]
scene.scene.camera.view_angle = 5.033164800000001
scene.scene.camera.view_up = [-0.14301336238889156, 0.21138117409951956, -0.96688426267807892]
scene.scene.camera.clipping_range = [41498.787075197979, 137750.00692720726]
scene.scene.camera.compute_view_plane_normal()
scene.scene.render()
myv.savefig('fig/euler-solutions-3d-zoom.png')
myv.show()
Image(filename='fig/euler-solutions-3d.png')

Out[15]: