The fatiando package has been deprecated. Please check out the new tools in the Fatiando a Terra website: www.fatiando.org

fatiando 0.5

Seismic profiling (fatiando.seismic.profile)

Forward modeling and inversion of vertical seismic profiling (VSP) data.

In this kind of profiling, the wave source is located at the surface on top of the well. The seismic waves are then measured at different depths along the well.

Forward modeling

Inversion

  • LayeredStraight: Inverts for the slownesses of a layered model assuming straight ray paths.
class fatiando.seismic.profile.LayeredStraight(traveltimes, zp, thickness)[source]

Bases: fatiando.inversion.misfit.Misfit

Inversion of straight-ray travel-times for the velocity of a layered medium

Assumes that the source is at the top of the well and that rays follow a straight path (no reflection or refraction). Also assumes known thicknesses (may be a fine discretization if real thickness is not known).

Actually solves for the slowness (1/velocity) so that the problem becomes linear and more manageable. Use the estimate_ attribute to get the estimated velocities. Slowness with stored in the estimated parameter vector p_.

Uses fatiando.seismic.ttime2d.straight for forward modeling.

Note

In most cases requires regularization. The recommended types are Damping and Smoothness1D.

Parameters:

  • traveltimes
    : list

    The first-arrival travel-times calculated at the measurement stations

  • zp
    : list

    The depths of the measurement stations (seismometers)

  • thickness
    : list

    The thickness of each layer in order of increasing depth

Notes:

The ith travel-time \(t_i\) measured at depth \(z_i\) is a function of the wave velocity \(v_j\) and distance \(d_{ij}\) that it traveled in each layer

\[t_i(z_i) = \sum\limits_{j=1}^M \frac{d_{ij}}{v_j}\]

The distance \(d_{ij}\) is smaller or equal to the thickness of the layer \(s_j\). Notice that \(d_{ij} = 0\) if the jth layer is below \(z_i\), \(d_{ij} = s_j\) if the jth layer is above \(z_i\), and \(d_{ij} < s_j\) if \(z_i\) is inside the jth layer.

To make \(t_i\) linear with respect to \(v_j\), we can use the slowness \(w_j = 1/v_j\) instead of velocity

\[t_i(z_i) = \sum\limits_{j=1}^M d_{ij} w_j\]

Thus, the parameters we want to estimate in this inversion are the slownesses of each layer.

From the above equation, we can see that the element \(G_{ij}\) of the Jacobian (sensitivity) matrix is given by

\[G_{ij} = d_{ij}\]

Examples:

Using some synthetic data produced by layered_straight_ray and assuming that the thickness of the layers is known:

>>> import numpy as np
>>> # Make a 4 layer model
>>> thicks = [10, 20, 10, 30]
>>> vels = [2, 4, 10, 8]
>>> # Set the recording stations
>>> zp = range(1, sum(thicks), 5)
>>> # Calculate the travel-times
>>> tts = layered_straight_ray(thicks, vels, zp)
>>> # Solve for the slowness assuming known thicknesses
>>> solver = LayeredStraight(tts, zp, thicks).fit()
>>> # The estimated velocities
>>> solver.estimate_
array([ 2.,  4., 10.,  8.])
>>> # and the corresponding slownesses
>>> solver.p_
array([ 0.5  ,  0.25 ,  0.1  ,  0.125])
>>> # Check the fit
>>> np.all(np.abs(solver.residuals()) < 10**-10)
True

See the Cookbook for more complex examples that use regularization and unknown thicknesses.

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:

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 estimated slowness to velocity.

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.

fatiando.seismic.profile.layered_straight_ray(thickness, velocity, zp)[source]

Calculates straight-ray (no refraction) travel-times in a layered model.

The source is assumed to be at z = 0 and on the well. The z-axis is positive downward.

Parameters:

  • thickness
    : list

    The thickness of each layer in order of increasing depth

  • velocity
    : list

    The velocity of each layer in order of increasing depth

  • zp
    : list

    The depths of the measurement stations (seismometers)

Returns:

  • travel_times
    : array

    The first-arrival travel-times calculated at the measurement stations.

Examples:

>>> import numpy as np
>>> # Make a 4 layer model
>>> thicks = [10, 20, 10, 30]
>>> vels = [2, 4, 10, 5]
>>> # Set the recording depths
>>> zs = [10, 30, 40, 70]
>>> # Calculate the travel-times from a surface source
>>> tt = layered_straight_ray(thicks, vels, zs)
>>> tt
array([  5.,  10.,  11.,  17.])