3D Euler deconvolution methods (fatiando.gravmag.euler)¶

Euler deconvolution methods for potential fields.

• EulerDeconv: The classic 3D solution to Euler’s equation for potential fields (Reid et al., 1990). Runs on the whole dataset.
• EulerDeconvEW: Run Euler deconvolution on an expanding window over the data set and keep the best estimate.
• EulerDeconvMW: Run Euler deconvolution on a moving window over the data set to produce a set of estimates.

References

Reid, A. B., J. M. Allsop, H. Granser, A. J. Millett, and I. W. Somerton (1990), Magnetic interpretation in three dimensions using Euler deconvolution, Geophysics, 55(1), 80-91, doi:10.1190/1.1442774.

class fatiando.gravmag.euler.EulerDeconv(x, y, z, field, xderiv, yderiv, zderiv, structural_index)[source]¶

Bases: fatiando.inversion.misfit.Misfit

Classic 3D Euler deconvolution of potential field data.

Follows the formulation of Reid et al. (1990). Performs the deconvolution on the whole data set. For windowed approaches, use EulerDeconvMW (moving window) and EulerDeconvEW (expanding window).

Works on any potential field that satisfies Euler’s homogeneity equation (both gravity and magnetic, assuming simple sources):

\[(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),\]

in which \(f_i\) is the given potential field observation at point \((x_i, y_i, z_i)\), \(b\) is the base level (a constant shift of the field, like a regional field), \(\eta\) is the structural index, and \((x_0, y_0, z_0)\) are the coordinates of a point on the source (for a sphere, this is the center point).

The Euler deconvolution estimates \((x_0, y_0, z_0)\) and \(b\) given a potential field and its x, y, z derivatives and the structural index. However, this assumes that the sources are ideal (see the table below). We recommend reading Reid and Thurston (2014) for a discussion on what the structural index means and what it does not mean.

Warning

Please read the paper Reid et al. (2014) to avoid doing horrible things with Euler deconvolution. Uieda et al. (2014) offer a practical tutorial using Fatiando code and show some common misinterpretations.

After Reid et al. (2014), values of the structural index (SI) can be:

Source type	SI (Mag)	SI (Grav)
Point, sphere	3	2
Line, cylinder, thin bed fault	2	1
Thin sheet edge, thin sill, thin dyke	1	0

Use the fit method to run the deconvolution. The estimated coordinates \((x_0, y_0, z_0)\) are stored in the estimate_ attribute and the estimated base level \(b\) is stored in baselevel_.

Note

Using structural index of 0 is not supported yet.

Note

The data does not need to be gridded for this! So long as you can calculate the derivatives of non-gridded data (using an Equivalent Layer, for example).

Note

x is North, y is East, and z is down.

Note

Units of the input data (x, y, z, field, derivatives) must be in SI units! Otherwise, the results will be in strange units. Use functions in fatiando.utils to convert between units.

Parameters:

• x, y, z

  : 1d-arrays

  The x, y, and z coordinates of the observation points

• field

  : 1d-array

  The potential field measured at the observation points

• xderiv, yderiv, zderiv

  : 1d-arrays

  The x-, y-, and z-derivatives of the potential field (measured or calculated) at the observation points

• index

  : float

  The structural index of the source

References:

Reid, A. B., J. M. Allsop, H. Granser, A. J. Millett, and I. W. Somerton (1990), Magnetic interpretation in three dimensions using Euler deconvolution, Geophysics, 55(1), 80-91, doi:10.1190/1.1442774.

Reid, A. B., J. Ebbing, and S. J. Webb (2014), Avoidable Euler Errors – the use and abuse of Euler deconvolution applied to potential fields, Geophysical Prospecting, doi:10.1111/1365-2478.12119.

Reid, A., and J. Thurston (2014), The structural index in gravity and magnetic interpretation: Errors, uses, and abuses, GEOPHYSICS, 79(4), J61-J66, doi:10.1190/geo2013-0235.1.

Uieda, L., V. C. Oliveira Jr., and V. C. F. Barbosa (2014), Geophysical tutorial: Euler deconvolution of potential-field data, The Leading Edge, 33(4), 448-450, doi:10.1190/tle33040448.1.

config(method, **kwargs)¶

Configure the optimization method and its parameters.

This sets the method used by fit and the keyword arguments that are passed to it.

Parameters:

• method

  : string

  The optimization method. One of: 'linear', 'newton', 'levmarq', 'steepest', 'acor'

Other keyword arguments that can be passed are the ones allowed by each method.

Some methods have required arguments:

• newton, levmarq and steepest require the initial argument (an initial estimate for the gradient descent)
• acor requires the bounds argument (min/max values for the search space)

See the corresponding docstrings for more information:

• linear
• newton
• levmarq
• steepest
• acor

copy(deep=False)¶

Make a copy of me together with all the cached methods.

estimate_¶

A nicely formatted version of the estimate.

If the class implements a fmt_estimate method, this will its results. This can be used to convert the parameter vector to a more useful form, like a fatiando.mesher object.

fit()¶

Solve for the parameter vector that minimizes this objective function.

Uses the optimization method and parameters defined using the config method.

The estimated parameter vector can be accessed through the p_ attribute. A (possibly) formatted version (converted to a more manageable type) of the estimate can be accessed through the property estimate_.

fmt_estimate(p)[source]¶

Separate the (x, y, z) point coordinates from the baselevel.

Coordinates are stored in estimate_ and a base level is stored in baselevel_.

gradient(p)¶

The gradient vector of the misfit function.

\[\bar{g} = -2\bar{\bar{J}}^T\bar{r}\]

where \(\bar{\bar{J}}\) is the Jacobian matrix and \(\bar{r}\) is the residual vector.

Parameters:

• p

  : 1d-array

  The parameter vector where the gradient is evaluated

Returns:

• gradient

  : 1d-array

  The gradient vector.

hessian(p)¶

The Hessian of the misfit function with respect to the parameters.

Calculated using the Gauss approximation:

\[\bar{\bar{H}} \approx 2\bar{\bar{J}}^T\bar{\bar{J}}\]

where \(\bar{\bar{J}}\) is the Jacobian matrix.

For linear problems, the Hessian matrix is cached in memory, so calling this method again will not trigger a re-calculation.

Parameters:

• p

  : 1d-array

  The parameter vector where the Hessian is evaluated

Returns:

• hessian

  : 2d-array

  The Hessian matrix

regul_param¶

The regularization parameter (scale factor) for the objetive function.

Defaults to 1.

residuals(p=None)¶

Calculate the residuals vector (observed - predicted data).

Parameters:

• p

  : 1d-array or None

  The parameter vector used to calculate the residuals. If None, will use the current estimate stored in estimate_.

Returns:

• residuals

  : 1d-array or list of 1d-arrays

  The residual vector. If this is the sum of 1 or more Misfit instances, will return the residual vector from each of the summed misfits in the order of the sum.

set_weights(weights)¶

Set the data weights.

Using weights for the data, the least-squares data-misfit function becomes:

\[\phi = \bar{r}^T \bar{\bar{W}}\bar{r}\]

Parameters:

• weights

  : 1d-array or 2d-array or None

  Weights for the data vector. If None, will remove any weights that have been set before. If it is a 2d-array, it will be interpreted as the weight matrix \(\bar{\bar{W}}\). If it is a 1d-array, it will be interpreted as the diagonal of the weight matrix (all off-diagonal elements will default to zero). The weight matrix can be a sparse array from scipy.sparse.

value(p)¶

Calculate the value of the misfit for a given parameter vector.

The value is given by:

\[\phi(\bar{p}) = \bar{r}^T\bar{\bar{W}}\bar{r}\]

where \(\bar{r}\) is the residual vector and \(bar{\bar{W}}\) are optional data weights.

Parameters:

• p

  : 1d-array or None

  The parameter vector.

Returns:

• value

  : float

  The value of the misfit function.

class fatiando.gravmag.euler.EulerDeconvEW(x, y, z, field, xderiv, yderiv, zderiv, structural_index, center, sizes)[source]¶

Bases: fatiando.gravmag.euler.EulerDeconv

Euler deconvolution using an expanding window scheme.

Uses data inside a window of growing size to perform the Euler deconvolution. Keeps the best result, judged by the estimated error.

The deconvolution is performed as in EulerDeconv.

Use the fit method to produce an estimate. The estimated point is stored in the attribute estimate_ and the base level in baselevel_.

Parameters:

• x, y, z

  : 1d-arrays

  The x, y, and z coordinates of the observation points

• field

  : 1d-array

  The potential field measured at the observation points

• xderiv, yderiv, zderiv

  : 1d-arrays

  The x-, y-, and z-derivatives of the potential field (measured or calculated) at the observation points

• index

  : float

  The structural index of the source

• center

  : [x, y]

  The x, y coordinates of the center of the expanding windows.

• sizes

  : list or 1d-array

  The sizes of the windows.

config(method, **kwargs)¶

Configure the optimization method and its parameters.

This sets the method used by fit and the keyword arguments that are passed to it.

Parameters:

• method

  : string

  The optimization method. One of: 'linear', 'newton', 'levmarq', 'steepest', 'acor'

Other keyword arguments that can be passed are the ones allowed by each method.

Some methods have required arguments:

• newton, levmarq and steepest require the initial argument (an initial estimate for the gradient descent)
• acor requires the bounds argument (min/max values for the search space)

See the corresponding docstrings for more information:

• linear
• newton
• levmarq
• steepest
• acor

copy(deep=False)¶

Make a copy of me together with all the cached methods.

estimate_¶

A nicely formatted version of the estimate.

If the class implements a fmt_estimate method, this will its results. This can be used to convert the parameter vector to a more useful form, like a fatiando.mesher object.

fit()[source]¶

Perform the Euler deconvolution with expanding windows.

The estimated point is stored in estimate_, the base level in baselevel_.

fmt_estimate(p)¶

Separate the (x, y, z) point coordinates from the baselevel.

Coordinates are stored in estimate_ and a base level is stored in baselevel_.

gradient(p)¶

The gradient vector of the misfit function.

\[\bar{g} = -2\bar{\bar{J}}^T\bar{r}\]

where \(\bar{\bar{J}}\) is the Jacobian matrix and \(\bar{r}\) is the residual vector.

Parameters:

• p

  : 1d-array

  The parameter vector where the gradient is evaluated

Returns:

• gradient

  : 1d-array

  The gradient vector.

hessian(p)¶

The Hessian of the misfit function with respect to the parameters.

Calculated using the Gauss approximation:

\[\bar{\bar{H}} \approx 2\bar{\bar{J}}^T\bar{\bar{J}}\]

where \(\bar{\bar{J}}\) is the Jacobian matrix.

For linear problems, the Hessian matrix is cached in memory, so calling this method again will not trigger a re-calculation.

Parameters:

• p

  : 1d-array

  The parameter vector where the Hessian is evaluated

Returns:

• hessian

  : 2d-array

  The Hessian matrix

regul_param¶

The regularization parameter (scale factor) for the objetive function.

Defaults to 1.

residuals(p=None)¶

Calculate the residuals vector (observed - predicted data).

Parameters:

• p

  : 1d-array or None

  The parameter vector used to calculate the residuals. If None, will use the current estimate stored in estimate_.

Returns:

• residuals

  : 1d-array or list of 1d-arrays

  The residual vector. If this is the sum of 1 or more Misfit instances, will return the residual vector from each of the summed misfits in the order of the sum.

set_weights(weights)¶

Set the data weights.

Using weights for the data, the least-squares data-misfit function becomes:

\[\phi = \bar{r}^T \bar{\bar{W}}\bar{r}\]

Parameters:

• weights

  : 1d-array or 2d-array or None

  Weights for the data vector. If None, will remove any weights that have been set before. If it is a 2d-array, it will be interpreted as the weight matrix \(\bar{\bar{W}}\). If it is a 1d-array, it will be interpreted as the diagonal of the weight matrix (all off-diagonal elements will default to zero). The weight matrix can be a sparse array from scipy.sparse.

value(p)¶

Calculate the value of the misfit for a given parameter vector.

The value is given by:

\[\phi(\bar{p}) = \bar{r}^T\bar{\bar{W}}\bar{r}\]

where \(\bar{r}\) is the residual vector and \(bar{\bar{W}}\) are optional data weights.

Parameters:

• p

  : 1d-array or None

  The parameter vector.

Returns:

• value

  : float

  The value of the misfit function.

class fatiando.gravmag.euler.EulerDeconvMW(x, y, z, field, xderiv, yderiv, zderiv, structural_index, windows, size, keep=0.2)[source]¶

Bases: fatiando.gravmag.euler.EulerDeconv

Solve an Euler deconvolution problem using a moving window scheme.

Uses data inside a window moving to perform the Euler deconvolution. Keeps only a top percentage of the estimates from all windows.

The deconvolution is performed as in EulerDeconv.

Use the fit method to produce an estimate. The estimated points are stored in estimate_ as a 2D numpy array. Each line in the array is an [x, y, z] coordinate of a point. The base levels are stored in baselevel_.

Parameters:

• x, y, z

  : 1d-arrays

  The x, y, and z coordinates of the observation points

• field

  : 1d-array

  The potential field measured at the observation points

• xderiv, yderiv, zderiv

  : 1d-arrays

  The x-, y-, and z-derivatives of the potential field (measured or calculated) at the observation points

• index

  : float

  The structural index of the source

• windows

  : (ny, nx)

  The number of windows in the y and x directions

• size

  : (dy, dx)

  The size of the windows in the y and x directions

• keep

  : float

  Decimal percentage of solutions to keep. Will rank the solutions by increasing error and keep only the first keep percent.

config(method, **kwargs)¶

Configure the optimization method and its parameters.

This sets the method used by fit and the keyword arguments that are passed to it.

Parameters:

• method

  : string

  The optimization method. One of: 'linear', 'newton', 'levmarq', 'steepest', 'acor'

Other keyword arguments that can be passed are the ones allowed by each method.

Some methods have required arguments:

• newton, levmarq and steepest require the initial argument (an initial estimate for the gradient descent)
• acor requires the bounds argument (min/max values for the search space)

See the corresponding docstrings for more information:

• linear
• newton
• levmarq
• steepest
• acor

copy(deep=False)¶

Make a copy of me together with all the cached methods.

estimate_¶

A nicely formatted version of the estimate.

If the class implements a fmt_estimate method, this will its results. This can be used to convert the parameter vector to a more useful form, like a fatiando.mesher object.

fit()[source]¶

Perform the Euler deconvolution on a moving window.

The estimated points are stored in estimate_, the base levels in baselevel_.

fmt_estimate(p)[source]¶

Separate the (x, y, z) point coordinates from the baselevel.

Coordinates are stored in estimate_ and a base level is stored in baselevel_.

gradient(p)¶

The gradient vector of the misfit function.

\[\bar{g} = -2\bar{\bar{J}}^T\bar{r}\]

where \(\bar{\bar{J}}\) is the Jacobian matrix and \(\bar{r}\) is the residual vector.

Parameters:

• p

  : 1d-array

  The parameter vector where the gradient is evaluated

Returns:

• gradient

  : 1d-array

  The gradient vector.

hessian(p)¶

The Hessian of the misfit function with respect to the parameters.

Calculated using the Gauss approximation:

\[\bar{\bar{H}} \approx 2\bar{\bar{J}}^T\bar{\bar{J}}\]

where \(\bar{\bar{J}}\) is the Jacobian matrix.

For linear problems, the Hessian matrix is cached in memory, so calling this method again will not trigger a re-calculation.

Parameters:

• p

  : 1d-array

  The parameter vector where the Hessian is evaluated

Returns:

• hessian

  : 2d-array

  The Hessian matrix

regul_param¶

The regularization parameter (scale factor) for the objetive function.

Defaults to 1.

residuals(p=None)¶

Calculate the residuals vector (observed - predicted data).

Parameters:

• p

  : 1d-array or None

  The parameter vector used to calculate the residuals. If None, will use the current estimate stored in estimate_.

Returns:

• residuals

  : 1d-array or list of 1d-arrays

  The residual vector. If this is the sum of 1 or more Misfit instances, will return the residual vector from each of the summed misfits in the order of the sum.

set_weights(weights)¶

Set the data weights.

Using weights for the data, the least-squares data-misfit function becomes:

\[\phi = \bar{r}^T \bar{\bar{W}}\bar{r}\]

Parameters:

• weights

  : 1d-array or 2d-array or None

  Weights for the data vector. If None, will remove any weights that have been set before. If it is a 2d-array, it will be interpreted as the weight matrix \(\bar{\bar{W}}\). If it is a 1d-array, it will be interpreted as the diagonal of the weight matrix (all off-diagonal elements will default to zero). The weight matrix can be a sparse array from scipy.sparse.

value(p)¶

Calculate the value of the misfit for a given parameter vector.

The value is given by:

\[\phi(\bar{p}) = \bar{r}^T\bar{\bar{W}}\bar{r}\]

where \(\bar{r}\) is the residual vector and \(bar{\bar{W}}\) are optional data weights.

Parameters:

• p

  : 1d-array or None

  The parameter vector.

Returns:

• value

  : float

  The value of the misfit function.