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')[source]#
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
- z
numpy.ndarray
Depth axis
- x
numpy.ndarray
Spatial axis
- t
numpy.ndarray
Time axis for data
- srcs
numpy.ndarray
Sources in array of size \(\lbrack 2 (3) \times n_s \rbrack\) The first axis should be ordered as (
y
,)x
,z
.- recs
numpy.ndarray
Receivers in array of size \(\lbrack 2 (3) \times n_r \rbrack\) The first axis should be ordered as (
y
,)x
,z
.- vel
numpy.ndarray
orfloat
Velocity model of size \(\lbrack (n_y\,\times)\; n_x \times n_z \rbrack\) (or constant)
- wav
numpy.ndarray
Wavelet.
- wavcenter
int
Index of wavelet center
- y
numpy.ndarray
Additional spatial axis (for 3-dimensional problems)
- mode
str
, optional Computation mode (
analytic
,eikonal
orbyot
, see Notes for more details)- wavfilter
bool
, optional New in version 2.0.0.
Apply wavelet filter (
True
) or not (False
)- dynamic
bool
, 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).- trav
numpy.ndarray
ortuple
, 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, onlymode='dynamic'
is only possible when traveltimes are provided in the latter form.- amp
numpy.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 whenmode='dynamic'
is chosen.- aperture
float
ortuple
, 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.- angleaperture
float
ortuple
, 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. Seeaperture
for implementation details regarding scalar and tuple cases.- snell
float
ortuple
, 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.
- engine
str
, optional Engine used for computations (
numpy
ornumba
).- dtype
str
, optional Type of elements in input array.
- name
str
, optional New in version 2.0.0.
Name of operator (to be used by
pylops.utils.describe.describe
)
- z
- Raises
- NotImplementedError
If
mode
is neitheranalytic
,eikonal
, orbyot
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 as2D
: \(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 formode=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
ormode=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 usingtrav
input parameter. Similarly, in this case one can also provide their own amplitude scalingamp
.
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
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
#
19. Image Domain Least-squares migration