pylops.waveeqprocessing.Kirchhoff

class pylops.waveeqprocessing.Kirchhoff(z, x, t, srcs, recs, vel, wav, wavcenter, y=None, mode='eikonal', wavfilter=False, dynamic=False, trav=None, amp=None, aperture=None, angleaperture=90.0, snell=None, engine='numpy', dtype='float64', name='K')

Kirchhoff demigration operator.

Kirchhoff-based demigration/migration operator. Uses a high-frequency approximation of the Green’s function propagators based on traveltimes and amplitudes that are either computed internally by solving the Eikonal equation, or passed directly by the user (which can use any other propagation engine of choice).

Parameters

znumpy.ndarray

Depth axis

xnumpy.ndarray

Spatial axis

tnumpy.ndarray

Time axis for data

srcsnumpy.ndarray

Sources in array of size \(\lbrack 2 (3) \times n_s \rbrack\) The first axis should be ordered as (y,) x, z.

recsnumpy.ndarray

Receivers in array of size \(\lbrack 2 (3) \times n_r \rbrack\) The first axis should be ordered as (y,) x, z.

velnumpy.ndarray or float

Velocity model of size \(\lbrack (n_y\,\times)\; n_x \times n_z \rbrack\) (or constant)

wavnumpy.ndarray

Wavelet.

wavcenterint

Index of wavelet center

ynumpy.ndarray

Additional spatial axis (for 3-dimensional problems)

modestr, optional

Computation mode (analytic, eikonal or byot, see Notes for more details)

wavfilterbool, optional

New in version 2.0.0.

Apply wavelet filter (True) or not (False)

dynamicbool, optional

New in version 2.0.0.

Apply dynamic weights in computations (True) or not (False). This includes both the amplitude terms of the Green’s function and the reflectivity-related scaling term (see equations below).

travnumpy.ndarray or tuple, optional

Traveltime table of size \(\lbrack (n_y) n_x n_z \times n_s n_r \rbrack\) or pair of traveltime tables of size \(\lbrack (n_y) n_x n_z \times n_s \rbrack\) and \(\lbrack (n_y) n_x n_z \times n_r \rbrack\) (to be provided if mode='byot'). Note that the latter approach is recommended as less memory demanding than the former. Moreover, only mode='dynamic' is only possible when traveltimes are provided in the latter form.

ampnumpy.ndarray, optional

New in version 2.0.0.

Pair of amplitude tables of size \(\lbrack (n_y) n_x n_z \times n_s \rbrack\) and \(\lbrack (n_y) n_x n_z \times n_r \rbrack\) (to be provided if mode='byot'). Note that this parameter is only required when mode='dynamic' is chosen.

aperturefloat or tuple, optional

New in version 2.0.0.

Maximum allowed aperture expressed as the ratio of offset over depth. If None, no aperture limitations are introduced. If scalar, a taper from 80% to 100% of aperture is applied. If tuple, apertures below the first value are accepted and those after the second value are rejected. A tapering is implemented for those between such values.

angleaperturefloat or tuple, optional

New in version 2.0.0.

Maximum allowed angle (either source or receiver side) in degrees. If None, angle aperture limitations are not introduced. See aperture for implementation details regarding scalar and tuple cases.

snellfloat or tuple, optional

Deprecated, will be removed in v3.0.0. Simply kept for back-compatibility with previous implementation, but effectively not affecting the behaviour of the operator.

enginestr, optional

Engine used for computations (numpy or numba).

dtypestr, optional

Type of elements in input array.

namestr, optional

New in version 2.0.0.

Name of operator (to be used by pylops.utils.describe.describe)

Raises

NotImplementedError

If mode is neither analytic, eikonal, or byot

Notes

The Kirchhoff demigration operator synthesizes seismic data given a propagation velocity model \(v(\mathbf{x})\) and a reflectivity model \(m(\mathbf{x})\). In forward mode [1], [2], [3]:

\[d(\mathbf{x_r}, \mathbf{x_s}, t) = \widetilde{w}(t) * \int_V \frac{2 \cos\theta} {v(\mathbf{x})} G(\mathbf{x_r}, \mathbf{x}, t) G(\mathbf{x}, \mathbf{x_s}, t) m(\mathbf{x}) \,\mathrm{d}\mathbf{x}\]

where \(G(\mathbf{x}, \mathbf{x_s}, t)\) and \(G(\mathbf{x_r}, \mathbf{x}, t)\) are the Green’s functions from source-to-subsurface-to-receiver and finally \(\widetilde{w}(t)\) is either a filtered version of the wavelet \(w(t)\) as explained below (wavfilter=True) or the wavelet itself when (wavfilter=False). Moreover, an angle scaling is included in the modelling operator, where the reflection angle \(\theta=(\theta_s-\theta_r)/2\) is half of the opening angle, with \(\theta_s\) and \(\theta_r\) representing the angles between the source-side and receiver-side rays and the vertical at the image point, respectively.

In our implementation, the following high-frequency approximation of the Green’s functions is adopted:

\[G(\mathbf{x_r}, \mathbf{x}, \omega) = a(\mathbf{x_r}, \mathbf{x}) e^{j \omega t(\mathbf{x_r}, \mathbf{x})}\]

where \(a(\mathbf{x_r}, \mathbf{x})\) is the amplitude and \(t(\mathbf{x_r}, \mathbf{x})\) is the traveltime. When dynamic=False, the amplitude correction terms are disregarded leading to a kinematic-only Kirchhoff operator.

\[d(\mathbf{x_r}, \mathbf{x_s}, t) = \tilde{w}(t) * \int_V e^{j \omega (t(\mathbf{x_r}, \mathbf{x}) + t(\mathbf{x}, \mathbf{x_s}))} m(\mathbf{x}) \,\mathrm{d}\mathbf{x}\]

On the other hand, when dynamic=True, the amplitude scaling is defined as

• 2D: \(a(\mathbf{x}, \mathbf{y})=\frac{1}{\sqrt{\text{dist}(\mathbf{x}, \mathbf{y})}}\)

• 3D: \(a(\mathbf{x}, \mathbf{y})=\frac{1}{\text{dist}(\mathbf{x}, \mathbf{y})}\)

approximating the geometrical spreading of the wavefront. For mode=analytic, \(\text{dist}(\mathbf{x}, \mathbf{y})=\|\mathbf{x} - \mathbf{y}\|\), whilst for mode=eikonal, this is computed internally by the Eikonal solver.

The wavelet filtering is applied as follows [4]:

• 2D: \(\tilde{W}(f)=\sqrt{j\omega} \cdot W(f)\)

• 3D: \(\tilde{W}(f)=-j\omega \cdot W(f)\)

Depending on the choice of mode the traveltime and amplitude of the Green’s function will be also computed differently:

• mode=analytic or mode=eikonal: traveltimes, amplitudes, and angles are computed for every source-image point-receiver triplets upfront and the Green’s functions are implemented from traveltime and amplitude look-up tables, placing scaled reflectivity values at corresponding source-to-receiver time in the data.

• mode=byot: bring your own tables. Traveltime table are provided directly by user using trav input parameter. Similarly, in this case one can also provide their own amplitude scaling amp.

Two aperture limitations have been also implemented as defined by:

• aperture: the maximum allowed aperture is expressed as the ratio of offset over depth. This aperture limitation avoid including grazing angles whose contributions can introduce aliasing effects. A taper is added at the edges of the aperture;

• angleaperture: the maximum allowed angle aperture is expressed as the difference between the incident or emerging angle at every image point and the vertical axis. This aperture limitation also avoid including grazing angles whose contributions can introduce aliasing effects. Note that for a homogenous medium and slowly varying heterogeneous medium the offset and angle aperture limits may work in the same way.

Finally, the adjoint of the demigration operator is a migration operator which projects data in the model domain creating an image of the subsurface reflectivity.

1

Bleistein, N., Cohen, J.K., and Stockwell, J.W. “Mathematics of Multidimensional Seismic Imaging, Migration and Inversion”, 2000.

2

Santos, L.T., Schleicher, J., Tygel, M., and Hubral, P. “Seismic modeling by demigration”, Geophysics, 65(4), pp. 1281-1289, 2000.

3

Yang, K., and Zhang, J. “Comparison between Born and Kirchhoff operators for least-squares reverse time migration and the constraint of the propagation of the background wavefield”, Geophysics, 84(5), pp. R725-R739, 2019.

4

Safron, L. “Multicomponent least-squares Kirchhoff depth migration”, MSc Thesis, 2018.

Attributes

shapetuple

Operator shape

explicitbool

Operator contains a matrix that can be solved explicitly (True) or not (False)

Methods

__init__(z, x, t, srcs, recs, vel, wav, ...)
adjoint()
apply_columns(cols)	Apply subset of columns of operator
cond([uselobpcg])	Condition number of linear operator.
conj()	Complex conjugate operator
div(y[, niter, densesolver])	Solve the linear problem \(\mathbf{y}=\mathbf{A}\mathbf{x}\).
dot(x)	Matrix-matrix or matrix-vector multiplication.
eigs([neigs, symmetric, niter, uselobpcg])	Most significant eigenvalues of linear operator.
matmat(X)	Matrix-matrix multiplication.
matvec(x)	Matrix-vector multiplication.
reset_count()	Reset counters
rmatmat(X)	Matrix-matrix multiplication.
rmatvec(x)	Adjoint matrix-vector multiplication.
todense([backend])	Return dense matrix.
toimag([forw, adj])	Imag operator
toreal([forw, adj])	Real operator
tosparse()	Return sparse matrix.
trace([neval, method, backend])	Trace of linear operator.
transpose()

Examples using pylops.waveeqprocessing.Kirchhoff

15. Least-squares migration

15. Least-squares migration

20. Image Domain Least-squares migration

20. Image Domain Least-squares migration