Source code for pylops.waveeqprocessing.oneway

import logging

import numpy as np
from scipy.sparse.linalg import lsqr

from pylops import Diagonal, Identity, LinearOperator, Pad
from pylops.signalprocessing import FFT
from pylops.utils import dottest as Dottest
from pylops.utils.backend import to_cupy_conditional
from pylops.utils.tapers import taper2d, taper3d

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

def _phase_shift(phiin, freq, kx, vel, dz, ky=0, adj=False):
    """Phase shift extrapolation for single depth level in 2d/3d constant
    velocity medium

    phiin : :obj:`numpy.ndarray`
        Frequency-wavenumber spectrum of input wavefield
        (only positive frequencies)
    freq : :obj:`numpy.ndarray`
        Positive frequency axis (already gridded with kx)
    kx : :obj:`int`, optional
        Horizontal wavenumber first axis (already gridded with freq)
    ky : :obj:`int`, optional
        Horizontal wavenumber second axis (already gridded with freq). If ``0``
        is provided, this reduces to phase shift in a 2d medium.
    vel : :obj:`float`, optional
        Constant propagation velocity.
    dz : :obj:`float`, optional
        Depth step.

    phiin : :obj:`numpy.ndarray`
        Frequency-wavenumber spectrum of depth extrapolated wavefield
        (only positive frequencies)

    # vertical slowness
    kz = (freq / vel) ** 2 - kx ** 2 - ky ** 2
    kz = np.sqrt(kz.astype(phiin.dtype))
    # ensure evanescent region is complex positive
    kz = np.real(kz) - 1j * np.sign(dz) * np.abs(np.imag(kz))
    # create and apply propagator
    gazx = np.exp(-1j * 2 * np.pi * dz * kz)
    if adj:
        gazx = np.conj(gazx)
    phiout = phiin * gazx
    return phiout

class _PhaseShift(LinearOperator):
    """Phase shift operator in frequency-wavenumber domain

    Apply positive phase shift directly in frequency-wavenumber domain.
    See :class:`pylops.waveeqprocessingPhaseShift` for more details on the
    input parameters.


    def __init__(self, vel, dz, freq, kx, ky=None, dtype="complex64"):
        self.vel = vel = dz
        # define frequency and horizontal wavenumber axes
        if ky is None:
            ky = 0
            [freq, kx] = np.meshgrid(freq, kx, indexing="ij")
            [freq, kx, ky] = np.meshgrid(freq, kx, ky, indexing="ij")
        # define vertical wavenumber axis
        kz = (freq / vel) ** 2 - kx ** 2 - ky ** 2
        kz = np.sqrt(kz.astype(dtype))
        # ensure evanescent region is complex positive
        kz = np.real(kz) - 1j * np.sign(dz) * np.abs(np.imag(kz))
        # create propagator
        self.gazx = np.exp(-1j * 2 * np.pi * dz * kz)

        self.dims = freq.shape
        self.shape = (,
        self.dtype = np.dtype(dtype)
        self.explicit = False

    def _matvec(self, x):
        if not isinstance(self.gazx, type(x)):
            self.gazx = to_cupy_conditional(x, self.gazx)
        y = x.reshape(self.dims) * self.gazx
        return y.ravel()

    def _rmatvec(self, x):
        if not isinstance(self.gazx, type(x)):
            self.gazx = to_cupy_conditional(x, self.gazx)
        y = x.reshape(self.dims) * np.conj(self.gazx)
        return y.ravel()

[docs]def PhaseShift(vel, dz, nt, freq, kx, ky=None, dtype="float64"): r"""Phase shift operator Apply positive (forward) phase shift with constant velocity in forward mode, and negative (backward) phase shift with constant velocity in adjoint mode. Input model and data should be 2- or 3-dimensional arrays in time-space domain of size :math:`[n_t \times n_x \;(\times n_y)]`. Parameters ---------- vel : :obj:`float`, optional Constant propagation velocity dz : :obj:`float`, optional Depth step nt : :obj:`int`, optional Number of time samples of model and data freq : :obj:`numpy.ndarray` Positive frequency axis kx : :obj:`int`, optional Horizontal wavenumber axis (centered around 0) of size :math:`[n_x \times 1]`. ky : :obj:`int`, optional Second horizontal wavenumber axis for 3d phase shift (centered around 0) of size :math:`[n_y \times 1]`. dtype : :obj:`str`, optional Type of elements in input array Returns ------- Pop : :obj:`pylops.LinearOperator` Phase shift operator Notes ----- The phase shift operator implements a one-way wave equation forward propagation in frequency-wavenumber domain by applying the following transformation to the input model: .. math:: d(f, k_x, k_y) = m(f, k_x, k_y) e^{-j \Delta z \sqrt{\omega^2/v^2 - k_x^2 - k_y^2}} where :math:`v` is the constant propagation velocity and :math:`\Delta z` is the propagation depth. In adjoint mode, the data is propagated backward using the following transformation: .. math:: m(f, k_x, k_y) = d(f, k_x, k_y) e^{j \Delta z \sqrt{\omega^2/v^2 - k_x^2 - k_y^2}} Effectively, the input model and data are assumed to be in time-space domain and forward Fourier transform is applied to both dimensions, leading to the following operator: .. math:: \mathbf{d} = \mathbf{F}^H_t \mathbf{F}^H_x \mathbf{P} \mathbf{F}_x \mathbf{F}_t \mathbf{m} where :math:`\mathbf{P}` perfoms the phase-shift as discussed above. """ dtypefft = (np.ones(1, dtype=dtype) + 1j * np.ones(1, dtype=dtype)).dtype if ky is None: dims = (nt, kx.size) dimsfft = (freq.size, kx.size) else: dims = (nt, kx.size, ky.size) dimsfft = (freq.size, kx.size, ky.size) Fop = FFT(dims, dir=0, nfft=nt, real=True, dtype=dtype) Kxop = FFT( dimsfft, dir=1, nfft=kx.size, real=False, fftshift_after=True, dtype=dtypefft ) if ky is not None: Kyop = FFT( dimsfft, dir=2, nfft=ky.size, real=False, fftshift_after=True, dtype=dtypefft, ) Pop = _PhaseShift(vel, dz, freq, kx, ky, dtypefft) if ky is None: Pop = Fop.H * Kxop * Pop * Kxop.H * Fop else: Pop = Fop.H * Kxop * Kyop * Pop * Kyop.H * Kxop.H * Fop # Recasting of type is required to avoid FFT operators to cast to complex. # We know this is correct because forward and inverse FFTs are applied at # the beginning and end of this combined operator Pop.dtype = dtype return LinearOperator(Pop)
[docs]def Deghosting( p, nt, nr, dt, dr, vel, zrec, pd=None, win=None, npad=(11, 11), ntaper=(11, 11), restriction=None, sptransf=None, solver=lsqr, dottest=False, dtype="complex128", **kwargs_solver ): r"""Wavefield deghosting. Apply seismic wavefield decomposition from single-component (pressure) data. This process is also generally referred to as model-based deghosting. Parameters ---------- p : :obj:`np.ndarray` Pressure data of of size :math:`\lbrack n_{r_x}\,(\times n_{r_y}) \times n_t \rbrack` (or :math:`\lbrack n_{r_{x,\text{sub}}}\, (\times n_{r_{y,\text{sub}}}) \times n_t \rbrack` in case a ``restriction`` operator is provided. Note that :math:`n_{r_{x,\text{sub}}}` (and :math:`n_{r_{y,\text{sub}}}`) must agree with the size of the output of this operator) nt : :obj:`int` Number of samples along the time axis nr : :obj:`int` or :obj:`tuple` Number of samples along the receiver axis (or axes) dt : :obj:`float` Sampling along the time axis dr : :obj:`float` or :obj:`tuple` Sampling along the receiver array of the separated pressure consituents vel : :obj:`float` Velocity along the receiver array (must be constant) zrec : :obj:`float` Depth of receiver array pd : :obj:`np.ndarray`, optional Direct arrival to be subtracted from ``p`` win : :obj:`np.ndarray`, optional Time window to be applied to ``p`` to remove the direct arrival (if ``pd=None``) ntaper : :obj:`float` or :obj:`tuple`, optional Number of samples of taper applied to propagator to avoid edge effects npad : :obj:`float` or :obj:`tuple`, optional Number of samples of padding applied to propagator to avoid edge effects angle restriction : :obj:`pylops.LinearOperator`, optional Restriction operator sptransf : :obj:`pylops.LinearOperator`, optional Sparsifying operator solver : :obj:`float`, optional Function handle of solver to be used if ``kind='inverse'`` dottest : :obj:`bool`, optional Apply dot-test dtype : :obj:`str`, optional Type of elements in input array. If ``None``, directly inferred from ``p`` **kwargs_solver Arbitrary keyword arguments for chosen ``solver`` Returns ------- pup : :obj:`np.ndarray` Up-going wavefield pdown : :obj:`np.ndarray` Down-going wavefield Notes ----- Up- and down-going components of seismic data :math:`p^-(x, t)` and :math:`p^+(x, t)` can be estimated from single-component data :math:`p(x, t)` using a ghost model. The basic idea [1]_ is that of using a one-way propagator in the f-k domain (also referred to as ghost model) to predict the down-going field from the up-going one (excluded the direct arrival and its source ghost referred here to as :math:`p_d(x, t)`): .. math:: p^+ - p_d = e^{-j k_z 2 z_\text{rec}} p^- where :math:`k_z` is the vertical wavenumber and :math:`z_\text{rec}` is the depth of the array of receivers In a matrix form we can thus write the total wavefield as: .. math:: \mathbf{p} - \mathbf{p_d} = (\mathbf{I} + \Phi) \mathbf{p}^- where :math:`\Phi` is one-way propagator implemented via the :class:`pylops.waveeqprocessing.PhaseShift` operator. .. [1] Amundsen, L., 1993, Wavenumber-based filtering of marine point-source data: GEOPHYSICS, 58, 1335–1348. """ ndims = p.ndim if ndims == 2: dims = (nt, nr) nrs = nr nkx = nr + 2 * npad kx = np.fft.ifftshift(np.fft.fftfreq(nkx, dr)) ky = None else: dims = (nt, nr[0], nr[1]) nrs = nr[0] * nr[1] nkx = nr[0] + 2 * npad[0] kx = np.fft.ifftshift(np.fft.fftfreq(nkx, dr[0])) nky = nr[1] + 2 * npad[1] ky = np.fft.ifftshift(np.fft.fftfreq(nky, dr)) nf = nt freq = np.fft.rfftfreq(nf, dt) # Phase shift operator zprop = 2 * zrec if ndims == 2: taper = taper2d(nt, nr, ntaper).T Padop = Pad(dims, ((0, 0), (npad, npad))) else: taper = taper3d(nt, nr, ntaper).transpose(2, 0, 1) Padop = Pad(dims, ((0, 0), (npad[0], npad[0]), (npad[1], npad[1]))) Pop = ( -Padop.H * PhaseShift(vel, zprop, nt, freq, kx, ky) * Padop * Diagonal(taper.ravel(), dtype=dtype) ) # Decomposition operator Dupop = Identity(nt * nrs, dtype=p.dtype) + Pop if dottest: Dottest(Dupop, nt * nrs, nt * nrs, verb=True) # Add restriction if restriction is not None: Dupop_norestr = Dupop Dupop = restriction * Dupop # Add sparsify transform if sptransf is not None: Dupop_norestr = Dupop_norestr * sptransf Dupop = Dupop * sptransf # Define data if pd is not None: d = p - pd else: d = win * p # Inversion pup = solver(Dupop, d.ravel(), **kwargs_solver)[0] # Apply sparse transform if sptransf is not None: p = Dupop_norestr * pup # reconstruct p at finely sampled spatial axes pup = sptransf * pup p = np.real(p).reshape(dims) # Finalize estimates pup = np.real(pup).reshape(dims) pdown = p - pup return pup, pdown