Source code for pylops.waveeqprocessing.kirchhoff

__all__ = ["Kirchhoff"]


import logging
import os
import warnings
from typing import Optional, Tuple, Union

import numpy as np

from pylops import LinearOperator
from pylops.signalprocessing import Convolve1D
from pylops.utils import deps
from pylops.utils._internal import _value_or_sized_to_array
from pylops.utils.decorators import reshaped
from pylops.utils.tapers import taper
from pylops.utils.typing import DTypeLike, NDArray

skfmm_message = deps.skfmm_import("the kirchhoff module")
jit_message = deps.numba_import("the kirchhoff module")

if skfmm_message is None:
    import skfmm

if jit_message is None:
    from numba import jit, prange

    # detect whether to use parallel or not
    numba_threads = int(os.getenv("NUMBA_NUM_THREADS", "1"))
    parallel = True if numba_threads != 1 else False
else:
    prange = range

logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.WARNING)


[docs]class Kirchhoff(LinearOperator): r"""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 : :obj:`numpy.ndarray` Depth axis x : :obj:`numpy.ndarray` Spatial axis t : :obj:`numpy.ndarray` Time axis for data srcs : :obj:`numpy.ndarray` Sources in array of size :math:`\lbrack 2 (3) \times n_s \rbrack` The first axis should be ordered as (``y``,) ``x``, ``z``. recs : :obj:`numpy.ndarray` Receivers in array of size :math:`\lbrack 2 (3) \times n_r \rbrack` The first axis should be ordered as (``y``,) ``x``, ``z``. vel : :obj:`numpy.ndarray` or :obj:`float` Velocity model of size :math:`\lbrack (n_y\,\times)\; n_x \times n_z \rbrack` (or constant) wav : :obj:`numpy.ndarray` Wavelet. wavcenter : :obj:`int` Index of wavelet center y : :obj:`numpy.ndarray` Additional spatial axis (for 3-dimensional problems) mode : :obj:`str`, optional Computation mode (``analytic``, ``eikonal`` or ``byot``, see Notes for more details) wavfilter : :obj:`bool`, optional .. versionadded:: 2.0.0 Apply wavelet filter (``True``) or not (``False``) dynamic : :obj:`bool`, optional .. versionadded:: 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 : :obj:`numpy.ndarray` or :obj:`tuple`, optional Traveltime table of size :math:`\lbrack (n_y) n_x n_z \times n_s n_r \rbrack` or pair of traveltime tables of size :math:`\lbrack (n_y) n_x n_z \times n_s \rbrack` and :math:`\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. amp : :obj:`numpy.ndarray`, optional .. versionadded:: 2.0.0 Pair of amplitude tables of size :math:`\lbrack (n_y) n_x n_z \times n_s \rbrack` and :math:`\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. aperture : :obj:`float` or :obj:`tuple`, optional .. versionadded:: 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 : :obj:`float` or :obj:`tuple`, optional .. versionadded:: 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. snell : :obj:`float` or :obj:`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. engine : :obj:`str`, optional Engine used for computations (``numpy`` or ``numba``). dtype : :obj:`str`, optional Type of elements in input array. name : :obj:`str`, optional .. versionadded:: 2.0.0 Name of operator (to be used by :func:`pylops.utils.describe.describe`) Attributes ---------- shape : :obj:`tuple` Operator shape explicit : :obj:`bool` Operator contains a matrix that can be solved explicitly (``True``) or not (``False``) Raises ------ NotImplementedError If ``mode`` is neither ``analytic``, ``eikonal``, or ``byot`` Notes ----- The Kirchhoff demigration operator synthesizes seismic data given a propagation velocity model :math:`v(\mathbf{x})` and a reflectivity model :math:`m(\mathbf{x})`. In forward mode [1]_, [2]_, [3]_: .. math:: 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 :math:`G(\mathbf{x}, \mathbf{x_s}, t)` and :math:`G(\mathbf{x_r}, \mathbf{x}, t)` are the Green's functions from source-to-subsurface-to-receiver and finally :math:`\widetilde{w}(t)` is either a filtered version of the wavelet :math:`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 :math:`\theta=(\theta_s-\theta_r)/2` is half of the opening angle, with :math:`\theta_s` and :math:`\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: .. math:: G(\mathbf{x_r}, \mathbf{x}, \omega) = a(\mathbf{x_r}, \mathbf{x}) e^{j \omega t(\mathbf{x_r}, \mathbf{x})} where :math:`a(\mathbf{x_r}, \mathbf{x})` is the amplitude and :math:`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. .. math:: 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``: :math:`a(\mathbf{x}, \mathbf{y})=\frac{1}{\sqrt{\text{dist}(\mathbf{x}, \mathbf{y})}}` * ``3D``: :math:`a(\mathbf{x}, \mathbf{y})=\frac{1}{\text{dist}(\mathbf{x}, \mathbf{y})}` approximating the geometrical spreading of the wavefront. For ``mode=analytic``, :math:`\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``: :math:`\tilde{W}(f)=\sqrt{j\omega} \cdot W(f)` * ``3D``: :math:`\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. """ def __init__( self, z: NDArray, x: NDArray, t: NDArray, srcs: NDArray, recs: NDArray, vel: NDArray, wav: NDArray, wavcenter: int, y: Optional[NDArray] = None, mode: str = "eikonal", wavfilter: bool = False, dynamic: bool = False, trav: Optional[NDArray] = None, amp: Optional[NDArray] = None, aperture: Optional[Tuple[float, float]] = None, angleaperture: Union[float, Tuple[float, float]] = 90.0, snell: Optional[Tuple[float, float]] = None, engine: str = "numpy", dtype: DTypeLike = "float64", name: str = "K", ) -> None: warnings.warn( "A new implementation of Kirchhoff is provided in v2.1.0. " "This currently affects only the inner working of the " "operator, end-users can continue using the operator in " "the same way. Nevertheless, it is now recommended to provide" "the variables trav (and amp) as a tuples containing the " "traveltime (and amplitude) tables for sources and receivers " "separately. This behaviour will eventually become default in " "version v3.0.0.", FutureWarning, ) # identify geometry ( self.ndims, _, dims, self.ny, self.nx, self.nz, ns, nr, _, _, _, _, _, ) = Kirchhoff._identify_geometry(z, x, srcs, recs, y=y) self.dt = t[1] - t[0] self.nt = len(t) # store ix-iy locations of sources and receivers for aperture filter self.dynamic = dynamic if self.dynamic: dx = x[1] - x[0] if self.ndims == 2: self.six = ( np.tile((srcs[0] - x[0]) // dx, (nr, 1)).T.astype(int).ravel() ) self.rix = np.tile((recs[0] - x[0]) // dx, (ns, 1)).astype(int).ravel() elif self.ndims == 3: # TODO: compute 3D indices for aperture filter # currently no aperture filter in 3D... just make indices 0 # so check if always passed self.six = np.zeros(nr * ns) self.rix = np.zeros(nr * ns) # compute traveltime and distances self.travsrcrec = True # use separate tables for src and rec traveltimes if mode in ["analytic", "eikonal", "byot"]: if mode in ["analytic", "eikonal"]: # compute traveltime table ( self.trav_srcs, self.trav_recs, dist_srcs, dist_recs, trav_srcs_grad, trav_recs_grad, ) = Kirchhoff._traveltime_table(z, x, srcs, recs, vel, y=y, mode=mode) if self.dynamic: # need to add a scalar in the denominator of amplitude term to avoid # division by 0, currently set to 1e-2 of max distance to avoid having # too large scaling around the source. This number may change in future # or left to the user to define epsdist = 1e-2 self.maxdist = epsdist * (np.max(dist_srcs) + np.max(dist_recs)) if self.ndims == 2: self.amp_srcs, self.amp_recs = 1.0 / np.sqrt( dist_srcs + self.maxdist ), 1.0 / np.sqrt(dist_recs + self.maxdist) else: self.amp_srcs, self.amp_recs = 1.0 / ( dist_srcs + self.maxdist ), 1.0 / (dist_recs + self.maxdist) else: if isinstance(trav, tuple): self.trav_srcs, self.trav_recs = trav else: self.travsrcrec = False self.trav = trav if self.dynamic and not self.travsrcrec: raise NotImplementedError( "separate traveltime tables must be provided " "when selecting mode=dynamic" ) if self.dynamic: if isinstance(amp, tuple): self.amp_srcs, self.amp_recs = amp else: raise NotImplementedError( "separate amplitude tables must be provided " ) if self.travsrcrec: # compute traveltime gradients at image points trav_srcs_grad = np.gradient( self.trav_srcs.reshape(*dims, ns), axis=np.arange(self.ndims), ) trav_recs_grad = np.gradient( self.trav_recs.reshape(*dims, nr), axis=np.arange(self.ndims), ) else: raise NotImplementedError("method must be analytic, eikonal or byot") # compute angles with vertical if self.dynamic: if self.ndims == 2: self.angle_srcs = np.arctan2( trav_srcs_grad[0], trav_srcs_grad[1] ).reshape(np.prod(dims), ns) self.angle_recs = np.arctan2( trav_recs_grad[0], trav_recs_grad[1] ).reshape(np.prod(dims), nr) else: trav_srcs_grad = np.concatenate( [trav_srcs_grad[i][np.newaxis] for i in range(3)] ) trav_recs_grad = np.concatenate( [trav_recs_grad[i][np.newaxis] for i in range(3)] ) self.angle_srcs = ( np.sign(trav_srcs_grad[1]) * np.arccos( trav_srcs_grad[-1] / np.sqrt(np.sum(trav_srcs_grad**2, axis=0)) ) ).reshape(np.prod(dims), ns) self.angle_recs = ( np.sign(trav_srcs_grad[1]) * np.arccos( trav_recs_grad[-1] / np.sqrt(np.sum(trav_recs_grad**2, axis=0)) ) ).reshape(np.prod(dims), nr) # pre-compute traveltime indices if total traveltime is used if not self.travsrcrec: self.itrav = (self.trav / self.dt).astype("int32") self.travd = self.trav / self.dt - self.itrav # create wavelet operator if wavfilter: self.wav = self._wavelet_reshaping( wav, self.dt, srcs.shape[0], recs.shape[0], self.ndims ) else: self.wav = wav self.cop = Convolve1D( (ns * nr, self.nt), h=self.wav, offset=wavcenter, axis=1, dtype=dtype ) # create fixed-size aperture taper for all apertures self.aperturetap = taper(41, 20, "hanning")[20:] # define aperture # if aperture=None, we want to ensure the check is always matched (no aperture limits...) # if aperture!=None in 3d, force to None as aperture checks are not yet implemented if aperture is not None and self.ndims == 3: aperture = None warnings.warn( "Aperture is forced to None as currently not implemented in 3D" ) if aperture is not None: warnings.warn( "Aperture is currently defined as ratio of offset over depth, " "and may be not ideal for highly heterogeneous media" ) self.aperture = ( (self.nx + 1,) if aperture is None else _value_or_sized_to_array(aperture) ) if len(self.aperture) == 1: self.aperture = np.array([0.8 * self.aperture[0], self.aperture[0]]) # define angle aperture angleaperture = [0.0, 1000.0] if angleaperture is None else angleaperture self.angleaperture = np.deg2rad(_value_or_sized_to_array(angleaperture)) if len(self.angleaperture) == 1: self.angleaperture = np.array( [0.8 * self.angleaperture[0], self.angleaperture[0]] ) # dimensions self.ns, self.nr = ns, nr self.nsnr = ns * nr self.ni = np.prod(dims) dims = tuple(dims) if self.ndims == 2 else (dims[0] * dims[1], dims[2]) dimsd = (ns, nr, self.nt) super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, name=name) # save velocity if using dynamic to compute amplitudes if self.dynamic: self.vel = ( vel.flatten() if not isinstance(vel, (float, int)) else vel * np.ones(np.prod(dims)) ) self._register_multiplications(engine) @staticmethod def _identify_geometry( z: NDArray, x: NDArray, srcs: NDArray, recs: NDArray, y: Optional[NDArray] = None, ) -> Tuple[ int, int, NDArray, int, int, int, int, int, float, float, float, NDArray, NDArray, ]: """Identify geometry and acquisition size and sampling""" ns, nr = srcs.shape[1], recs.shape[1] nz, nx = len(z), len(x) dz = np.abs(z[1] - z[0]) dx = np.abs(x[1] - x[0]) if y is None: ndims = 2 shiftdim = 0 ny = 1 dy = None dims = np.array([nx, nz]) dsamp = np.array([dx, dz]) origin = np.array([x[0], z[0]]) else: ndims = 3 shiftdim = 1 ny = len(y) dy = np.abs(y[1] - y[0]) dims = np.array([ny, nx, nz]) dsamp = np.array([dy, dx, dz]) origin = np.array([y[0], x[0], z[0]]) return ndims, shiftdim, dims, ny, nx, nz, ns, nr, dy, dx, dz, dsamp, origin @staticmethod def _traveltime_table( z: NDArray, x: NDArray, srcs: NDArray, recs: NDArray, vel: Union[float, NDArray], y: Optional[NDArray] = None, mode: str = "eikonal", ) -> Tuple[NDArray, NDArray, NDArray, NDArray, NDArray, NDArray]: r"""Traveltime table Compute traveltimes along the source-subsurface-receivers triplet in 2- or 3-dimensional media given a constant or depth- and space variable velocity. Parameters ---------- z : :obj:`numpy.ndarray` Depth axis x : :obj:`numpy.ndarray` Spatial axis srcs : :obj:`numpy.ndarray` Sources in array of size :math:`\lbrack 2 (3) \times n_s \rbrack` recs : :obj:`numpy.ndarray` Receivers in array of size :math:`\lbrack 2 (3) \times n_r \rbrack` vel : :obj:`numpy.ndarray` or :obj:`float` Velocity model of size :math:`\lbrack (n_y \times)\, n_x \times n_z \rbrack` (or constant) y : :obj:`numpy.ndarray` Additional spatial axis (for 3-dimensional problems) mode : :obj:`numpy.ndarray`, optional Computation mode (``eikonal``, ``analytic`` - only for constant velocity) Returns ------- trav_srcs : :obj:`numpy.ndarray` Source-to-subsurface traveltime table of size :math:`\lbrack (n_y*) n_x n_z \times n_s \rbrack` trav_recs : :obj:`numpy.ndarray` Receiver-to-subsurface traveltime table of size :math:`\lbrack (n_y) n_x n_z \times n_r \rbrack` dist_srcs : :obj:`numpy.ndarray` Source-to-subsurface distance table of size :math:`\lbrack (n_y*) n_x n_z \times n_s \rbrack` dist_recs : :obj:`numpy.ndarray` Receiver-to-subsurface distance table of size :math:`\lbrack (n_y) n_x n_z \times n_r \rbrack` trav_srcs_gradient : :obj:`numpy.ndarray` Source-to-subsurface traveltime gradient table of size :math:`\lbrack (n_y*) n_x n_z \times n_s \rbrack` trav_recs_gradient : :obj:`numpy.ndarray` Receiver-to-subsurface traveltime gradient table of size :math:`\lbrack (n_y) n_x n_z \times n_r \rbrack` """ # define geometry ( ndims, shiftdim, dims, ny, nx, nz, ns, nr, _, _, _, dsamp, origin, ) = Kirchhoff._identify_geometry(z, x, srcs, recs, y=y) # compute traveltimes if mode == "analytic": if not isinstance(vel, (float, int)): raise ValueError("vel must be scalar for mode=analytical") # compute grid if ndims == 2: X, Z = np.meshgrid(x, z, indexing="ij") X, Z = X.ravel(), Z.ravel() else: Y, X, Z = np.meshgrid(y, x, z, indexing="ij") Y, X, Z = Y.ravel(), X.ravel(), Z.ravel() dist_srcs2 = np.zeros((ny * nx * nz, ns)) dist_recs2 = np.zeros((ny * nx * nz, nr)) for isrc, src in enumerate(srcs.T): dist_srcs2[:, isrc] = (X - src[0 + shiftdim]) ** 2 + ( Z - src[1 + shiftdim] ) ** 2 if ndims == 3: dist_srcs2[:, isrc] += (Y - src[0]) ** 2 for irec, rec in enumerate(recs.T): dist_recs2[:, irec] = (X - rec[0 + shiftdim]) ** 2 + ( Z - rec[1 + shiftdim] ) ** 2 if ndims == 3: dist_recs2[:, irec] += (Y - rec[0]) ** 2 trav_srcs = np.sqrt(dist_srcs2) / vel trav_recs = np.sqrt(dist_recs2) / vel dist_srcs = trav_srcs * vel dist_recs = trav_recs * vel elif mode == "eikonal": if skfmm_message is None: dist_srcs = np.zeros((ny * nx * nz, ns)) dist_recs = np.zeros((ny * nx * nz, nr)) trav_srcs = np.zeros((ny * nx * nz, ns)) trav_recs = np.zeros((ny * nx * nz, nr)) for isrc, src in enumerate(srcs.T): src = np.round((src - origin) / dsamp).astype(np.int32) phi = np.ones_like(vel) if ndims == 2: phi[src[0], src[1]] = -1 else: phi[src[0], src[1], src[2]] = -1 dist_srcs[:, isrc] = (skfmm.distance(phi=phi, dx=dsamp)).ravel() trav_srcs[:, isrc] = ( skfmm.travel_time(phi=phi, speed=vel, dx=dsamp) ).ravel() for irec, rec in enumerate(recs.T): rec = np.round((rec - origin) / dsamp).astype(np.int32) phi = np.ones_like(vel) if ndims == 2: phi[rec[0], rec[1]] = -1 else: phi[rec[0], rec[1], rec[2]] = -1 dist_recs[:, irec] = (skfmm.distance(phi=phi, dx=dsamp)).ravel() trav_recs[:, irec] = ( skfmm.travel_time(phi=phi, speed=vel, dx=dsamp) ).ravel() else: raise NotImplementedError(skfmm_message) else: raise NotImplementedError("method must be analytic or eikonal") # compute traveltime gradients at image points trav_srcs_grad = np.gradient( trav_srcs.reshape(*dims, ns), *dsamp, axis=np.arange(ndims) ) trav_recs_grad = np.gradient( trav_recs.reshape(*dims, nr), *dsamp, axis=np.arange(ndims) ) return ( trav_srcs, trav_recs, dist_srcs, dist_recs, trav_srcs_grad, trav_recs_grad, ) def _wavelet_reshaping( self, wav: NDArray, dt: float, dimsrc: int, dimrec: int, dimv: int, ) -> NDArray: """Apply wavelet reshaping to account for omega scaling factor originating from the wave equation""" f = np.fft.rfftfreq(len(wav), dt) W = np.fft.rfft(wav, n=len(wav)) if dimsrc == 2 and dimv == 2: # 2D Wfilt = W * np.sqrt(1j * 2 * np.pi * f) elif (dimsrc == 2 or dimrec == 2) and dimv == 3: # 2.5D raise NotImplementedError("2.D wavelet currently not available") elif dimsrc == 3 and dimrec == 3 and dimv == 3: # 3D Wfilt = W * (-1j * 2 * np.pi * f) wavfilt = np.fft.irfft(Wfilt, n=len(wav)) return wavfilt @staticmethod def _trav_kirch_matvec( x: NDArray, y: NDArray, nsnr: int, nt: int, ni: int, itrav: NDArray, travd: NDArray, ) -> NDArray: for isrcrec in prange(nsnr): itravisrcrec = itrav[:, isrcrec] travdisrcrec = travd[:, isrcrec] for ii in range(ni): itravii = itravisrcrec[ii] travdii = travdisrcrec[ii] if 0 <= itravii < nt - 1: y[isrcrec, itravii] += x[ii] * (1 - travdii) y[isrcrec, itravii + 1] += x[ii] * travdii return y @staticmethod def _trav_kirch_rmatvec( x: NDArray, y: NDArray, nsnr: int, nt: int, ni: int, itrav: NDArray, travd: NDArray, ) -> NDArray: for ii in prange(ni): itravii = itrav[ii] travdii = travd[ii] for isrcrec in range(nsnr): itravisrcrecii = itravii[isrcrec] travdisrcrecii = travdii[isrcrec] if 0 <= itravisrcrecii < nt - 1: y[ii] += ( x[isrcrec, itravisrcrecii] * (1 - travdisrcrecii) + x[isrcrec, itravisrcrecii + 1] * travdisrcrecii ) return y @staticmethod def _travsrcrec_kirch_matvec( x: NDArray, y: NDArray, ns: int, nr: int, nt: int, ni: int, dt: float, trav_srcs: NDArray, trav_recs: NDArray, ) -> NDArray: for isrc in prange(ns): travisrc = trav_srcs[:, isrc] for irec in range(nr): travirec = trav_recs[:, irec] trav = travisrc + travirec itrav = (trav / dt).astype("int32") travd = trav / dt - itrav for ii in range(ni): itravii = itrav[ii] travdii = travd[ii] if 0 <= itravii < nt - 1: y[isrc * nr + irec, itravii] += x[ii] * (1 - travdii) y[isrc * nr + irec, itravii + 1] += x[ii] * travdii return y @staticmethod def _travsrcrec_kirch_rmatvec( x: NDArray, y: NDArray, ns: int, nr: int, nt: int, ni: int, dt: float, trav_srcs: NDArray, trav_recs: NDArray, ) -> NDArray: for ii in prange(ni): trav_srcsii = trav_srcs[ii] trav_recsii = trav_recs[ii] for isrc in prange(ns): trav_srcii = trav_srcsii[isrc] for irec in range(nr): trav_recii = trav_recsii[irec] travii = trav_srcii + trav_recii itravii = int(travii / dt) travdii = travii / dt - itravii if 0 <= itravii < nt - 1: y[ii] += ( x[isrc * nr + irec, itravii] * (1 - travdii) + x[isrc * nr + irec, itravii + 1] * travdii ) return y @staticmethod def _ampsrcrec_kirch_matvec( x: NDArray, y: NDArray, ns: int, nr: int, nt: int, ni: int, dt: float, vel: NDArray, trav_srcs: NDArray, trav_recs: NDArray, amp_srcs: NDArray, amp_recs: NDArray, aperturemin: float, aperturemax: float, aperturetap: NDArray, nz: int, six: NDArray, rix: NDArray, angleaperturemin: float, angleaperturemax: float, angles_srcs: NDArray, angles_recs: NDArray, ) -> NDArray: daperture = aperturemax - aperturemin dangleaperture = angleaperturemax - angleaperturemin for isrc in prange(ns): travisrc = trav_srcs[:, isrc] ampisrc = amp_srcs[:, isrc] angleisrc = angles_srcs[:, isrc] for irec in range(nr): travirec = trav_recs[:, irec] trav = travisrc + travirec itrav = (trav / dt).astype("int32") travd = trav / dt - itrav ampirec = amp_recs[:, irec] angleirec = angles_recs[:, irec] sixisrcrec = six[isrc * nr + irec] rixisrcrec = rix[isrc * nr + irec] # compute cosine of half opening angle and total amplitude scaling cosangle = np.cos((angleisrc - angleirec) / 2.0) amp = 2.0 * cosangle * ampisrc * ampirec / vel for ii in range(ni): itravii = itrav[ii] travdii = travd[ii] damp = amp[ii] # extract source and receiver angle at given image point angle_src = angleisrc[ii] angle_rec = angleirec[ii] abs_angle_src = abs(angle_src) abs_angle_rec = abs(angle_rec) # angle apertures checks aptap = 1.0 if ( abs_angle_src < angleaperturemax and abs_angle_rec < angleaperturemax ): if abs_angle_src >= angleaperturemin: # extract source angle aperture taper value aptap = ( aptap * aperturetap[ int( 20 * (abs_angle_src - angleaperturemin) // dangleaperture ) ] ) if abs_angle_rec >= angleaperturemin: # extract receiver angle aperture taper value aptap = ( aptap * aperturetap[ int( 20 * (abs_angle_rec - angleaperturemin) // dangleaperture ) ] ) # identify x-index of image point iz = ii % nz # aperture check aperture = abs(sixisrcrec - rixisrcrec) / (iz + 1) if aperture < aperturemax: if aperture >= aperturemin: # extract aperture taper value aptap = ( aptap * aperturetap[ int( 20 * ((aperture - aperturemin) // daperture) ) ] ) # time limit check if 0 <= itravii < nt - 1: y[isrc * nr + irec, itravii] += ( x[ii] * (1 - travdii) * damp * aptap ) y[isrc * nr + irec, itravii + 1] += ( x[ii] * travdii * damp * aptap ) return y @staticmethod def _ampsrcrec_kirch_rmatvec( x: NDArray, y: NDArray, ns: int, nr: int, nt: int, ni: int, dt: float, vel: NDArray, trav_srcs: NDArray, trav_recs: NDArray, amp_srcs: NDArray, amp_recs: NDArray, aperturemin: float, aperturemax: float, aperturetap: NDArray, nz: int, six: NDArray, rix: NDArray, angleaperturemin: float, angleaperturemax: float, angles_srcs: NDArray, angles_recs: NDArray, ) -> NDArray: daperture = aperturemax - aperturemin dangleaperture = angleaperturemax - angleaperturemin for ii in prange(ni): trav_srcsii = trav_srcs[ii] trav_recsii = trav_recs[ii] amp_srcsii = amp_srcs[ii] amp_recsii = amp_recs[ii] velii = vel[ii] angle_srcsii = angles_srcs[ii] angle_recsii = angles_recs[ii] # identify x-index of image point iz = ii % nz for isrc in range(ns): trav_srcii = trav_srcsii[isrc] amp_srcii = amp_srcsii[isrc] angle_src = angle_srcsii[isrc] for irec in range(nr): trav_recii = trav_recsii[irec] travii = trav_srcii + trav_recii itravii = int(travii / dt) travdii = travii / dt - itravii amp_recii = amp_recsii[irec] angle_rec = angle_recsii[irec] sixisrcrec = six[isrc * nr + irec] rixisrcrec = rix[isrc * nr + irec] abs_angle_src = abs(angle_src) abs_angle_rec = abs(angle_rec) # compute cosine of half opening angle and total amplitude scaling cosangle = np.cos((angle_src - angle_rec) / 2.0) damp = 2.0 * cosangle * amp_srcii * amp_recii / velii # angle apertures checks aptap = 1.0 if ( abs_angle_src < angleaperturemax and abs_angle_rec < angleaperturemax ): if abs_angle_src >= angleaperturemin: # extract source angle aperture taper value aptap = ( aptap * aperturetap[ int( 20 * (abs_angle_src - angleaperturemin) // dangleaperture ) ] ) if abs_angle_rec >= angleaperturemin: # extract receiver angle aperture taper value aptap = ( aptap * aperturetap[ int( 20 * (abs_angle_rec - angleaperturemin) // dangleaperture ) ] ) # aperture check aperture = abs(sixisrcrec - rixisrcrec) / (iz + 1) if aperture < aperturemax: if aperture >= aperturemin: # extract aperture taper value aptap = ( aptap * aperturetap[ int( 20 * ((aperture - aperturemin) // daperture) ) ] ) # time limit check if 0 <= itravii < nt - 1: # assign values y[ii] += ( ( x[isrc * nr + irec, itravii] * (1 - travdii) + x[isrc * nr + irec, itravii + 1] * travdii ) * damp * aptap ) return y def _register_multiplications(self, engine: str) -> None: if engine not in ["numpy", "numba"]: raise KeyError("engine must be numpy or numba") if engine == "numba" and jit_message is None: # numba numba_opts = dict( nopython=True, nogil=True, parallel=parallel ) # fastmath=True, if self.dynamic and self.travsrcrec: self._kirch_matvec = jit(**numba_opts)(self._ampsrcrec_kirch_matvec) self._kirch_rmatvec = jit(**numba_opts)(self._ampsrcrec_kirch_rmatvec) elif self.travsrcrec: self._kirch_matvec = jit(**numba_opts)(self._travsrcrec_kirch_matvec) self._kirch_rmatvec = jit(**numba_opts)(self._travsrcrec_kirch_rmatvec) elif not self.travsrcrec: self._kirch_matvec = jit(**numba_opts)(self._trav_kirch_matvec) self._kirch_rmatvec = jit(**numba_opts)(self._trav_kirch_rmatvec) else: if engine == "numba" and jit_message is not None: logging.warning(jit_message) if self.dynamic and self.travsrcrec: self._kirch_matvec = self._ampsrcrec_kirch_matvec self._kirch_rmatvec = self._ampsrcrec_kirch_rmatvec elif self.travsrcrec: self._kirch_matvec = self._travsrcrec_kirch_matvec self._kirch_rmatvec = self._travsrcrec_kirch_rmatvec elif not self.travsrcrec: self._kirch_matvec = self._trav_kirch_matvec self._kirch_rmatvec = self._trav_kirch_rmatvec @reshaped def _matvec(self, x: NDArray) -> NDArray: y = np.zeros((self.nsnr, self.nt), dtype=self.dtype) if self.dynamic and self.travsrcrec: inputs = ( x.ravel(), y, self.ns, self.nr, self.nt, self.ni, self.dt, self.vel, self.trav_srcs, self.trav_recs, self.amp_srcs, self.amp_recs, self.aperture[0], self.aperture[1], self.aperturetap, self.nz, self.six, self.rix, self.angleaperture[0], self.angleaperture[1], self.angle_srcs, self.angle_recs, ) elif self.travsrcrec: inputs = ( x.ravel(), y, self.ns, self.nr, self.nt, self.ni, self.dt, self.trav_srcs, self.trav_recs, ) elif not self.travsrcrec: inputs = (x.ravel(), y, self.nsnr, self.nt, self.ni, self.itrav, self.travd) y = self._kirch_matvec(*inputs) y = self.cop._matvec(y.ravel()) return y @reshaped def _rmatvec(self, x: NDArray) -> NDArray: x = self.cop._rmatvec(x.ravel()) x = x.reshape(self.nsnr, self.nt) y = np.zeros(self.ni, dtype=self.dtype) if self.dynamic and self.travsrcrec: inputs = ( x, y, self.ns, self.nr, self.nt, self.ni, self.dt, self.vel, self.trav_srcs, self.trav_recs, self.amp_srcs, self.amp_recs, self.aperture[0], self.aperture[1], self.aperturetap, self.nz, self.six, self.rix, self.angleaperture[0], self.angleaperture[1], self.angle_srcs, self.angle_recs, ) elif self.travsrcrec: inputs = ( x, y, self.ns, self.nr, self.nt, self.ni, self.dt, self.trav_srcs, self.trav_recs, ) elif not self.travsrcrec: inputs = (x, y, self.nsnr, self.nt, self.ni, self.itrav, self.travd) y = self._kirch_rmatvec(*inputs) return y