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

Added in version 2.0.0.

Apply wavelet filter (True) or not (False)

dynamicbool, optional

Added 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, mode='dynamic' and engine='cuda' are only possible when traveltimes are provided in the latter form.

ampnumpy.ndarray, optional

Added 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

Added 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

Added 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, numba or cuda).

dtypestr, optional

Type of elements in input array.

namestr, optional

Added in version 2.0.0.

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

Attributes:

ndimsint

Number of spatial dimensions (2 or 3).

nyint

Number of samples in y axis (1 if ndims=2).

nxint

Number of samples in x axis.

nzint

Number of samples in z axis.

dtfloat

Sampling step in time axis.

ntint

Number of samples in time axis.

nsnrint

Number of source-receiver pairs.

niint

Number of image points (ni=nx*nz if ndims=2 and ni=ny*nx*nz if ndims=3).

sixnumpy.ndarray

ix locations of sources (1d array of size ns*nr).

rixnumpy.ndarray

ix locations of receivers (1d array of size ns*nr).

travsrcrecbool

Whether separate traveltime tables for sources and receivers are used (True) or a single traveltime table for source-receiver pairs (False).

trav_srcsnumpy.ndarray

Traveltime table from sources to image points of size \(\lbrack (n_y) n_x n_z \times n_s \rbrack\).

trav_recsnumpy.ndarray

Traveltime table from image points to receivers of size \(\lbrack (n_y) n_x n_z \times n_r \rbrack\).

travnumpy.ndarray

Traveltime table from sources to receivers of size \(\lbrack (n_y) n_x n_z \times n_s n_r \rbrack\).

itravnumpy.ndarray

Integer traveltime table (in samples) from sources to receivers of size \(\lbrack (n_y) n_x n_z \times n_s n_r \rbrack\) (only if travsrcrec=False).

travdnumpy.ndarray

Fractional part of the traveltime table from sources to receivers of size \(\lbrack (n_y) n_x n_z \times n_s n_r \rbrack\) (only if travsrcrec=False).

maxdistfloat

Maximum distance between sources/receivers and image points.

amp_srcsnumpy.ndarray

Amplitude table from sources to image points of size \(\lbrack (n_y) n_x n_z \times n_s \rbrack\).

amp_recsnumpy.ndarray

Amplitude table from image points to receivers of size \(\lbrack (n_y) n_x n_z \times n_r \rbrack\).

angle_srcsnumpy.ndarray

Incident angle table from sources to image points of size \(\lbrack (n_y) n_x n_z \times n_s \rbrack\) ( only if dynamic=True).

angle_recsnumpy.ndarray

Emerging angle table from image points to receivers of size \(\lbrack (n_y) n_x n_z \times n_r \rbrack\) ( only if dynamic=True).

coppylops.Convolve1D

Wavelet convolution operator of size \(\lbrack n_t \times n_t \rbrack\).

aperturetapnumpy.ndarray

Aperture taper

aperturetuple or numpy.ndarray

Aperture limits (in terms of offset over depth)

apertureanglenumpy.ndarray

Angle aperture limits (in degrees)

velnumpy.ndarray

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

dimstuple

Shape of the array after the adjoint, but before flattening.

For example, x_reshaped = (Op.H * y.ravel()).reshape(Op.dims).

dimsdtuple

Shape of the array after the forward, but before flattening.

For example, y_reshaped = (Op * x.ravel()).reshape(Op.dimsd).

shapetuple

Operator shape.

Raises:

NotImplementedError

If engine="cuda" and trav is provided as a single table

ValueError

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.

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

19. Image Domain Least-squares migration

19. Image Domain Least-squares migration