Inversion for the relief of 2D basins (fatiando.gravmag.basin2d)¶

2D inversion for the basement relief of sedimentary basins.

There are different parametrizations available.

The simplest are meant more as an exercise and initiation in inverse problems:

• Triangular: assumes a basin with a triangular cross-section (think “foreland”).
• Trapezoidal: assumes a basin with a trapezoidal cross-section (think “grabben”).

More complex parametrizations are:

• PolygonalBasinGravity: approximate the basin by a polygon.

class fatiando.gravmag.basin2d.PolygonalBasinGravity(x, z, data, npoints, props, top=0, xlim=None)[source]¶

Bases: fatiando.inversion.misfit.Misfit

Estimate the relief of a sedimentary basin approximating by a polygon.

Currently only works with gravity data.

The top of the basin is straight and fixed at a given height. Polygon vertices are distributed evenly in the x-direction. The inversion estimates the depths of each vertex.

This is a non-linear inversion. Therefore you must configure it before running to choose a solver method and set the initial estimate. Use the config method for this.

Recommended configuration: Levemberg-Marquardt algorithm ('levmarq') with initial estimate to the average expected depth of the basin.

Typical regularization to use with this class are: Damping, Smoothness1D, TotalVariation1D.

The forward modeling is done using talwani. Derivatives are calculated using a 2-point finite difference approximation.

Tip

Use the estimate_ attribute to get a Polygon version of the estimated parameters (attribute p_).

Parameters:

• x, z

  : 1d-arrays

  The x and z coordinates of the observations. In meters.

• data

  : 1d-array

  The observed data.

• npoints

  : int

  Number of points to use

• props

  : dict

  The physical properties dictionary that will be assigned to the basin Polygon. Ex: to give the basin a density contrast of 500 kg/m3 props={'density': 500}.

• top

  : float

  The value of the z-coordinate where the top of the basin will be fixed. In meters. Default: 0.

• xlim

  : None or list = [xmin, xmax]

  The horizontal limits of the model. If not given, will use the limits of the data (i.e., [x.min(), x.max()]).

Examples:

Lets run an inversion on synthetic data from a simple model of a trapezoid basin (a polygon with 4 vertices). We’ll assume that the horizontal limits of the basin are the same as the limits of the data:

>>> from fatiando.mesher import Polygon
>>> from fatiando.gravmag import talwani
>>> import numpy as np
>>> # Make some synthetic data from a simple basin
>>> props = {'density': -500}
>>> model = [Polygon([[3000, 0], [2000, 800], [1000, 500], [0, 0]],
...                  props)]
>>> x = np.linspace(0, 3000, 50)
>>> z = -np.ones_like(x)  # Put data at 1m height
>>> data = talwani.gz(x, z, model)
>>> # Make the solver, configure, and invert.
>>> # Will use only 2 points because the two in the corners are
>>> # considered fixed in the inversion (at 'top').
>>> misfit = PolygonalBasinGravity(x, z, data, 2, props, top=0)
>>> _ = misfit.config('levmarq', initial=100*np.ones(misfit.nparams)).fit()
>>> misfit.p_
array([ 800.,  500.])
>>> type(misfit.estimate_)
<class 'fatiando.mesher.Polygon'>
>>> misfit.estimate_.vertices
array([[ 3000.,     0.],
       [ 2000.,   800.],
       [ 1000.,   500.],
       [    0.,     0.]])

If the x range of the data points is larger than the basin, you can specify a horizontal range for the basin model. When this is not specified, it is deduced from the data:

>>> x = np.linspace(-500, 3500, 80)
>>> z = -np.ones_like(x)
>>> data = talwani.gz(x, z, model)
>>> # Specify that the model used for inversion should be within
>>> # x => [0, 3000]
>>> misfit = PolygonalBasinGravity(x, z, data, 2, props, top=0,
...                                xlim=[0, 3000])
>>> _ = misfit.config('levmarq', initial=100*np.ones(misfit.nparams)).fit()
>>> misfit.p_
array([ 800.,  500.])
>>> misfit.estimate_.vertices
array([[ 3000.,     0.],
       [ 2000.,   800.],
       [ 1000.,   500.],
       [    0.,     0.]])

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 parameter vector to a fatiando.mesher.Polygon so that it can be used for plotting and forward modeling.

Examples:

>>> import numpy as np
>>> # Make some arrays to create the estimator clas
>>> x = np.linspace(-100, 300, 50)
>>> z = np.zeros_like(x)
>>> data = z
>>> misfit = PolygonalBasinGravity(x, z, data, 3, {}, top=-100)
>>> poly = misfit.fmt_estimate([1, 2, 3])
>>> type(poly)
<class 'fatiando.mesher.Polygon'>
>>> poly.vertices
array([[ 300., -100.],
       [ 200.,    1.],
       [ 100.,    2.],
       [   0.,    3.],
       [-100., -100.]])

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

jacobian(p)[source]¶

Calculate the Jacobian (sensitivity) matrix for a parameter vector.

p2vertices(p)[source]¶

Convert a parameter vector into vertices a Polygon.

Parameters:

• p

  : 1d-array

  The parameter vector with the depth of the polygon vertices

Returns:

• vertices

  : 2d-array

  Like a list of [x, z] coordinates of each vertex

Examples:

>>> import numpy as np
>>> # Make some arrays to create the estimator clas
>>> x = np.linspace(-100, 300, 50)
>>> z = np.zeros_like(x)
>>> data = z
>>> misfit = PolygonalBasinGravity(x, z, data, 3, {}, top=-100)
>>> misfit.p2vertices([1, 2, 3])
array([[ 300., -100.],
       [ 200.,    1.],
       [ 100.,    2.],
       [   0.,    3.],
       [-100., -100.]])

predicted(p)[source]¶

Calculate the predicted data for a parameter vector.

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.basin2d.Trapezoidal(x, z, gz, verts, density)[source]¶

Bases: fatiando.inversion.misfit.Misfit

Estimate the relief of a trapezoidal basin.

Use when the basin can be approximated by a 2D body with trapezoidal vertical cross-section, like in rifts.

The trapezoid is assumed to have 2 known vertices at the surface (the edges of the basin) and two unknown vertices in the subsurface. We assume that the x coordinates of the unknown vertices are the same as the x coordinates of the known vertices (i.e., the unknown vertices are directly under the known vertices). The inversion will then estimate the z coordinates of the unknown vertices.

The forward modeling is done using talwani. Derivatives are calculated using a 2-point finite difference approximation.

Tip

Use the estimate_ attribute to produce a polygon from the estimated parameter vector (p_).

Parameters:

• x, z

  : array

  Arrays with the x and z coordinates of the profile data points

• gz

  : array

  The profile gravity anomaly data

• verts

  : list of lists

  [[x1, z1], [x2, z2]] List of the [x, z] coordinates of the left and right know vertices, respectively.

  Warning

  Very important that the vertices in the list be ordered from left to right! Otherwise the forward model will give results with an inverted sign and terrible things may happen!

• density

  : float

  Density contrast of the basin

• delta

  : float

  Interval used to calculate the approximate derivatives

Note

The recommended solver for this inverse problem is the Levemberg-Marquardt method. Since this is a non-linear problem, set the desired method and initial solution using the config method. See the example bellow.

Example with synthetic data:

>>> import numpy as np
>>> from fatiando.mesher import Polygon
>>> from fatiando.gravmag import talwani
>>> # Make a trapezoidal basin model (will estimate the z coordinates
>>> # of the last two points)
>>> verts = [[10000, 1], [90000, 1], [90000, 5000], [10000, 3000]]
>>> model = Polygon(verts, {'density':500})
>>> # Generate the synthetic gz profile
>>> x = np.linspace(0, 100000, 50)
>>> z = np.zeros_like(x)
>>> gz = talwani.gz(x, z, [model])
>>> # Make a solver and fit it to the data
>>> solver = Trapezoidal(x, z, gz, verts[0:2], 500).config(
...     'levmarq', initial=[1000, 500]).fit()
>>> # p_ is the estimated parameter vector (z1 and z2 in this case)
>>> z1, z2 = solver.p_
>>> print('{:.1f}, {:.1f}'.format(z1, z2))
5000.0, 3000.0
>>> # The parameter vector is not that useful so use estimate_ to get a
>>> # Polygon object
>>> poly = solver.estimate_
>>> poly.vertices
array([[  1.00000000e+04,   1.00000000e+00],
       [  9.00000000e+04,   1.00000000e+00],
       [  9.00000000e+04,   5.00000000e+03],
       [  1.00000000e+04,   3.00000000e+03]])
>>> poly.props
{'density': 500}
>>> # Check is the residuals are all small
>>> np.all(np.abs(solver.residuals()) < 10**-10)
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 parameter vector to a fatiando.mesher.Polygon so that it can be used for plotting and forward modeling.

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.basin2d.Triangular(x, z, gz, verts, density)[source]¶

Bases: fatiando.inversion.misfit.Misfit

Estimate the relief of a triangular basin.

Use when the basin can be approximated by a 2D body with triangular vertical cross-section, like foreland basins.

The triangle is assumed to have 2 known vertices at the surface (the edges of the basin) and one unknown vertex in the subsurface. The inversion will estimate the (x, z) coordinates of the unknown vertex.

The forward modeling is done using talwani. Derivatives are calculated using a 2-point finite difference approximation.

Tip

Use the estimate_ attribute to produce a polygon from the estimated parameter vector (p_).

Parameters:

• x, z

  : array

  Arrays with the x and z coordinates of the profile data points

• gz

  : array

  The profile gravity anomaly data

• verts

  : list of lists

  [[x1, z1], [x2, z2]] List of the [x, z] coordinates of the left and right know vertices, respectively.

  Warning

  Very important that the vertices in the list be ordered from left to right! Otherwise the forward model will give results with an inverted sign and terrible things may happen!

• density

  : float

  Density contrast of the basin

• delta

  : float

  Interval used to calculate the approximate derivatives

Note

The recommended solver for this inverse problem is the Levemberg-Marquardt method. Since this is a non-linear problem, set the desired method and initial solution using the config method of this class. See the example bellow.

Example using synthetic data:

>>> import numpy as np
>>> from fatiando.mesher import Polygon
>>> from fatiando.gravmag import talwani
>>> # Make a triangular basin model (will estimate the last point)
>>> verts = [(10000, 1), (90000, 1), (50000, 5000)]
>>> left, middle, right = verts
>>> model = Polygon(verts, {'density':500})
>>> # Generate the synthetic gz profile
>>> x = np.linspace(0, 100000, 50)
>>> z = np.zeros_like(x)
>>> gz = talwani.gz(x, z, [model])
>>> # Make a solver and fit it to the data
>>> solver = Triangular(x, z, gz, [left, middle], 500).config(
...     'levmarq', initial=[10000, 1000]).fit()
>>> # p_ is the estimated parameter vector (x and z in this case)
>>> x, z = solver.p_
>>> print('{:.1f}, {:.1f}'.format(x, z))
50000.0, 5000.0
>>> # The parameter vector is not that useful so use estimate_ to get a
>>> # Polygon object
>>> poly = solver.estimate_
>>> poly.vertices
array([[  1.00000000e+04,   1.00000000e+00],
       [  9.00000000e+04,   1.00000000e+00],
       [  5.00000000e+04,   5.00000000e+03]])
>>> poly.props
{'density': 500}
>>> # Check is the residuals are all small
>>> np.all(np.abs(solver.residuals()) < 10**-10)
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 parameter vector to a Polygon so that it can be used for plotting and forward modeling.

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

jacobian(p)[source]¶

Calculate the Jacobian (sensitivity) matrix for a given parameter vector.

predicted(p)[source]¶

Calculate predicted data for a given parameter vector.

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.