Total magnetization estimation methods (fatiando.gravmag.magdir)¶

Estimation of the total magnetization vector of homogeneous bodies.

It estimates parameters related to the magnetization vector of homogeneous bodies.

Algorithms

• DipoleMagDir: This class estimates the Cartesian components of the magnetization vector of homogeneous dipolar bodies with known center. The estimated magnetization vector is converted to dipole moment, inclination (positive down) and declination (with respect to x, North).

class fatiando.gravmag.magdir.DipoleMagDir(x, y, z, data, inc, dec, points)[source]¶

Bases: fatiando.inversion.misfit.Misfit

Estimate the magnetization vector of a set of dipoles from magnetic total field anomaly using the method of Oliveira Jr. et al. (2015).

By using the well-known first-order approximation of the total field anomaly (Blakely, 1996, p. 179) produced by a set of dipoles, the estimation of the Cartesian components of the magnetization vectors is formulated as linear inverse problem. After estimating the magnetization vectors, they are converted to dipole moment, inclination (positive down) and declination (with respect to x, North).

After solving, use the estimate_ attribute to get the estimated magnetization vectors in dipole moment, inclination and declination. The estimated magnetization vectors in Cartesian coordinates can be accessed through the p_ attribute.

Note

Assumes x = North, y = East, z = Down.

Parameters:

• x, y, z

  : 1d-arrays

  The x, y, z coordinates of each data point.

• data

  : 1d-array

  The total field magnetic anomaly data at each point.

• inc, dec

  : floats

  The inclination and declination of the inducing field

• points

  : list of points [x, y, z]

  Each point [x, y, z] is the center of a dipole. Will invert for the Cartesian components of the magnetization vector of each dipole. Subsequently, the estimated magnetization vectors are converted to dipole moment, inclination and declination.

Note

Inclination is positive down and declination is measured with respect to x (North).

References:

Blakely, R. (1996), Potential theory in gravity and magnetic applications: CUP

Oliveira Jr., V. C., D. P. Sales, V. C. F. Barbosa, and L. Uieda (2015), Estimation of the total magnetization direction of approximately spherical bodies, Nonlin. Processes Geophys., 22(2), 215-232, doi:10.5194/npg-22-215-2015.

Examples:

Estimation of the total magnetization vector of dipoles with known centers

>>> import numpy as np
>>> from fatiando import gridder, utils
>>> from fatiando.gravmag import sphere
>>> from fatiando.mesher import Sphere, Prism
>>> # Produce some synthetic data
>>> area = (0, 10000, 0, 10000)
>>> x, y, z = gridder.scatter(area, 500, z=-150, seed=0)
>>> model = [Sphere(3000, 3000, 1000, 1000,
...              {'magnetization': utils.ang2vec(6.0, -20.0, -10.0)}),
...          Sphere(7000, 7000, 1000, 1000,
...              {'magnetization': utils.ang2vec(6.0, 30.0, -40.0)})]
>>> inc, dec = -9.5, -13
>>> tf = sphere.tf(x, y, z, model, inc, dec)
>>> # Give the coordinates of the dipoles
>>> points = [[3000.0, 3000.0, 1000.0], [7000.0, 7000.0, 1000.0]]
>>> p_true = np.hstack((ang2vec(CM*(4.*np.pi/3.)*6.0*1000**3,
...                                             -20.0, -10.0),
...                        ang2vec(CM*(4.*np.pi/3.)*6.0*1000**3,
...                                              30.0, -40.0)))
>>> estimate_true = [utils.vec2ang(p_true[3*i : 3*i + 3]) for i
...                                in range(len(points))]
>>> # Make a solver and fit it to the data
>>> solver = DipoleMagDir(x, y, z, tf, inc, dec, points).fit()
>>> # Check the fit
>>> np.allclose(tf, solver.predicted(), rtol=0.001, atol=0.001)
True
>>> # solver.p_ returns the Cartesian components of the
>>> # estimated magnetization vectors
>>> for p in solver.p_: print "%.10f" % p
2325.8255393651
-410.1057950109
-859.5903757213
1667.3411086852
-1399.0653093445
1256.6370614359
>>> # Check the estimated parameter vector
>>> np.allclose(p_true, solver.p_, rtol=0.001, atol=0.001)
True
>>> # The parameter vector is not that useful so use solver.estimate_
>>> # to convert the estimated magnetization vectors in dipole moment,
>>> # inclination and declination.
>>> for e in solver.estimate_:
...    print "%.10f %.10f %.10f" % (e[0], e[1], e[2])
2513.2741228718 -20.0000000000 -10.0000000000
2513.2741228718 30.0000000000 -40.0000000000
>>> # Check the converted estimate
>>> np.allclose(estimate_true, solver.estimate_, rtol=0.001,
...                                                 atol=0.001)
True

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]¶

Convert the estimate parameters from Cartesian to inclination, declication, and intensity.

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.