Source code for pylops.waveeqprocessing.seismicinterpolation

__all__ = ["SeismicInterpolation"]

import logging
from typing import List, Optional, Sequence, Tuple, Union

import numpy as np

from pylops import Laplacian, Restriction, SecondDerivative
from pylops.optimization.leastsquares import regularized_inversion
from pylops.optimization.sparsity import fista
from pylops.signalprocessing import (
    FFT2D,
    FFTND,
    ChirpRadon2D,
    ChirpRadon3D,
    Interp,
    Radon2D,
    Radon3D,
    Sliding2D,
    Sliding3D,
)
from pylops.utils.dottest import dottest as Dottest
from pylops.utils.typing import InputDimsLike, NDArray

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


[docs]def SeismicInterpolation( data: NDArray, nrec: Union[int, Tuple[int, int]], iava: Union[List[Union[int, float]], NDArray], iava1: Optional[Union[List[Union[int, float]], NDArray]] = None, kind: str = "fk", nffts: Optional[Union[int, InputDimsLike]] = None, sampling: Optional[Sequence[float]] = None, spataxis: Optional[NDArray] = None, spat1axis: Optional[NDArray] = None, taxis: Optional[NDArray] = None, paxis: Optional[NDArray] = None, p1axis: Optional[NDArray] = None, centeredh: bool = True, nwins: InputDimsLike = None, nwin: InputDimsLike = None, nover: InputDimsLike = None, engine: str = "numba", dottest: bool = False, **kwargs_solver, ) -> Tuple[NDArray, NDArray, NDArray]: r"""Seismic interpolation (or regularization). Interpolate seismic data from irregular to regular spatial grid. Depending on the size of the input ``data``, interpolation is either 2- or 3-dimensional. In case of 3-dimensional interpolation, data can be irregularly sampled in either one or both spatial directions. Parameters ---------- data : :obj:`np.ndarray` Irregularly sampled seismic data of size :math:`[n_{r_y} \,(\times n_{r_x} \times n_t)]` nrec : :obj:`int` or :obj:`tuple` Number of elements in the regularly sampled (reconstructed) spatial array, :math:`n_{R_y}` for 2-dimensional data and :math:`(n_{R_y}, n_{R_x})` for 3-dimensional data iava : :obj:`list` or :obj:`numpy.ndarray` Integer (or floating) indices of locations of available samples in first dimension of regularly sampled spatial grid of interpolated signal. The :class:`pylops.basicoperators.Restriction` operator is used in case of integer indices, while the :class:`pylops.signalprocessing.Iterp` operator is used in case of floating indices. iava1 : :obj:`list` or :obj:`numpy.ndarray`, optional Integer (or floating) indices of locations of available samples in second dimension of regularly sampled spatial grid of interpolated signal. Can be used only in case of 3-dimensional data. kind : :obj:`str`, optional Type of inversion: ``fk`` (default), ``spatial``, ``radon-linear``, ``chirpradon-linear``, ``radon-parabolic`` , ``radon-hyperbolic``, ``sliding``, or ``chirp-sliding`` nffts : :obj:`int` or :obj:`tuple`, optional nffts : :obj:`tuple`, optional Number of samples in Fourier Transform for each direction. Required if ``kind='fk'`` sampling : :obj:`tuple`, optional Sampling steps ``dy`` (, ``dx``) and ``dt``. Required if ``kind='fk'`` or ``kind='radon-linear'`` spataxis : :obj:`np.ndarray`, optional First spatial axis. Required for ``kind='radon-linear'``, ``kind='chirpradon-linear'``, ``kind='radon-parabolic'``, ``kind='radon-hyperbolic'``, can also be provided instead of ``sampling`` for ``kind='fk'`` spat1axis : :obj:`np.ndarray`, optional Second spatial axis. Required for ``kind='radon-linear'``, ``kind='chirpradon-linear'``, ``kind='radon-parabolic'``, ``kind='radon-hyperbolic'``, can also be provided instead of ``sampling`` for ``kind='fk'`` taxis : :obj:`np.ndarray`, optional Time axis. Required for ``kind='radon-linear'``, ``kind='chirpradon-linear'``, ``kind='radon-parabolic'``, ``kind='radon-hyperbolic'``, can also be provided instead of ``sampling`` for ``kind='fk'`` paxis : :obj:`np.ndarray`, optional First Radon axis. Required for ``kind='radon-linear'``, ``kind='chirpradon-linear'``, ``kind='radon-parabolic'``, ``kind='radon-hyperbolic'``, ``kind='sliding'``, and ``kind='chirp-sliding'`` p1axis : :obj:`np.ndarray`, optional Second Radon axis. Required for ``kind='radon-linear'``, ``kind='chirpradon-linear'``, ``kind='radon-parabolic'``, ``kind='radon-hyperbolic'``, ``kind='sliding'``, and ``kind='chirp-sliding'`` centeredh : :obj:`bool`, optional Assume centered spatial axis (``True``) or not (``False``). Required for ``kind='radon-linear'``, ``kind='radon-parabolic'`` and ``kind='radon-hyperbolic'`` nwins : :obj:`int` or :obj:`tuple`, optional Number of windows. Required for ``kind='sliding'`` and ``kind='chirp-sliding'`` nwin : :obj:`int` or :obj:`tuple`, optional Number of samples of window. Required for ``kind='sliding'`` and ``kind='chirp-sliding'`` nover : :obj:`int` or :obj:`tuple`, optional Number of samples of overlapping part of window. Required for ``kind='sliding'`` and ``kind='chirp-sliding'`` engine : :obj:`str`, optional Engine used for Radon computations (``numpy/numba`` for ``Radon2D`` and ``Radon3D`` or ``numpy/fftw`` for ``ChirpRadon2D`` and ``ChirpRadon3D``) dottest : :obj:`bool`, optional Apply dot-test **kwargs_solver Arbitrary keyword arguments for :py:func:`pylops.optimization.leastsquares.regularized_inversion` solver if ``kind='spatial'`` or :py:func:`pylops.optimization.sparsity.FISTA` solver otherwise Returns ------- recdata : :obj:`np.ndarray` Reconstructed data of size :math:`[n_{R_y}\,(\times n_{R_x} \times n_t)]` recprec : :obj:`np.ndarray` Reconstructed data in the sparse or preconditioned domain in case of ``kind='fk'``, ``kind='radon-linear'``, ``kind='radon-parabolic'``, ``kind='radon-hyperbolic'`` and ``kind='sliding'`` cost : :obj:`np.ndarray` Cost function norm Raises ------ KeyError If ``kind`` is neither ``spatial``, ``fl``, ``radon-linear``, ``radon-parabolic``, ``radon-hyperbolic`` nor ``sliding`` Notes ----- The problem of seismic data interpolation (or regularization) can be formally written as .. math:: \mathbf{y} = \mathbf{R} \mathbf{x} where a restriction or interpolation operator is applied along the spatial direction(s). Here :math:`\mathbf{y} = [\mathbf{y}_{R1}^T, \mathbf{y}_{R2}^T,\ldots, \mathbf{y}_{RN^T}]^T` where each vector :math:`\mathbf{y}_{Ri}` contains all time samples recorded in the seismic data at the specific receiver :math:`R_i`. Similarly, :math:`\mathbf{x} = [\mathbf{x}_{r1}^T, \mathbf{x}_{r2}^T,\ldots, \mathbf{x}_{rM}^T]`, contains all traces at the regularly and finely sampled receiver locations :math:`r_i`. Several alternative approaches can be taken to solve such a problem. They mostly differ in the choice of the regularization (or preconditining) used to mitigate the ill-posedness of the problem: * ``spatial``: least-squares inversion in the original time-space domain with an additional spatial smoothing regularization term, corresponding to the cost function :math:`J = ||\mathbf{y} - \mathbf{R} \mathbf{x}||_2 + \epsilon_\nabla \nabla ||\mathbf{x}||_2` where :math:`\nabla` is a second order space derivative implemented via :class:`pylops.basicoperators.SecondDerivative` in 2-dimensional case and :class:`pylops.basicoperators.Laplacian` in 3-dimensional case * ``fk``: L1 inversion in frequency-wavenumber preconditioned domain corresponding to the cost function :math:`J = ||\mathbf{y} - \mathbf{R} \mathbf{F} \mathbf{x}||_2` where :math:`\mathbf{F}` is frequency-wavenumber transform implemented via :class:`pylops.signalprocessing.FFT2D` in 2-dimensional case and :class:`pylops.signalprocessing.FFTND` in 3-dimensional case * ``radon-linear``: L1 inversion in linear Radon preconditioned domain using the same cost function as ``fk`` but with :math:`\mathbf{F}` being a Radon transform implemented via :class:`pylops.signalprocessing.Radon2D` in 2-dimensional case and :class:`pylops.signalprocessing.Radon3D` in 3-dimensional case * ``radon-parabolic``: L1 inversion in parabolic Radon preconditioned domain * ``radon-hyperbolic``: L1 inversion in hyperbolic Radon preconditioned domain * ``sliding``: L1 inversion in sliding-linear Radon preconditioned domain using the same cost function as ``fk`` but with :math:`\mathbf{F}` being a sliding Radon transform implemented via :class:`pylops.signalprocessing.Sliding2D` in 2-dimensional case and :class:`pylops.signalprocessing.Sliding3D` in 3-dimensional case """ dtype = data.dtype ndims = data.ndim if ndims == 1 or ndims > 3: raise ValueError("data must have 2 or 3 dimensions") if ndims == 2: dimsd = data.shape dims = (nrec, dimsd[1]) else: dimsd = data.shape dims = (nrec[0], nrec[1], dimsd[2]) # sampling if taxis is not None: dt = taxis[1] - taxis[0] if spataxis is not None: dspat = np.abs(spataxis[1] - spataxis[0]) if spat1axis is not None: dspat1 = np.abs(spat1axis[1] - spat1axis[0]) # create restriction/interpolation operator if iava.dtype == float: Rop = Interp(dims, iava, axis=0, kind="linear", dtype=dtype) if ndims == 3 and iava1 is not None: dims1 = (len(iava), nrec[1], dimsd[2]) Rop1 = Interp(dims1, iava1, axis=1, kind="linear", dtype=dtype) Rop = Rop1 * Rop else: Rop = Restriction(dims, iava, axis=0, dtype=dtype) if ndims == 3 and iava1 is not None: dims1 = (len(iava), nrec[1], dimsd[2]) Rop1 = Restriction(dims1, iava1, axis=1, dtype=dtype) Rop = Rop1 * Rop # create other operators for inversion if kind == "spatial": prec = False dotcflag = 0 if ndims == 3 and iava1 is not None: Regop = Laplacian(dims=dims, axes=(0, 1), dtype=dtype) else: Regop = SecondDerivative(dims, axis=0, dtype=dtype) SIop = Rop elif kind == "fk": prec = True dimsp = nffts dotcflag = 1 if ndims == 3: if sampling is None: if spataxis is None or spat1axis is None or taxis is None: raise ValueError( "Provide either sampling or spataxis, " f"spat1axis and taxis for kind=%{kind}" ) else: sampling = ( np.abs(spataxis[1] - spataxis[1]), np.abs(spat1axis[1] - spat1axis[1]), np.abs(taxis[1] - taxis[1]), ) Pop = FFTND(dims=dims, nffts=nffts, sampling=sampling) Pop = Pop.H else: if sampling is None: if spataxis is None or taxis is None: raise ValueError( "Provide either sampling or spataxis, " f"and taxis for kind={kind}" ) else: sampling = ( np.abs(spataxis[1] - spataxis[1]), np.abs(taxis[1] - taxis[1]), ) Pop = FFT2D(dims=dims, nffts=nffts, sampling=sampling) Pop = Pop.H SIop = Rop * Pop elif "chirpradon" in kind: prec = True dotcflag = 0 if ndims == 3: Pop = ChirpRadon3D( taxis, spataxis, spat1axis, (np.max(paxis) * dspat / dt, np.max(p1axis) * dspat1 / dt), ).H dimsp = (spataxis.size, spat1axis.size, taxis.size) else: Pop = ChirpRadon2D(taxis, spataxis, np.max(paxis) * dspat / dt).H dimsp = (spataxis.size, taxis.size) SIop = Rop * Pop elif "radon" in kind: prec = True dotcflag = 0 kindradon = kind.split("-")[-1] if ndims == 3: Pop = Radon3D( taxis, spataxis, spat1axis, paxis, p1axis, centeredh=centeredh, kind=kindradon, engine=engine, ) dimsp = (paxis.size, p1axis.size, taxis.size) else: Pop = Radon2D( taxis, spataxis, paxis, centeredh=centeredh, kind=kindradon, engine=engine, ) dimsp = (paxis.size, taxis.size) SIop = Rop * Pop elif kind in ("sliding", "chirp-sliding"): prec = True dotcflag = 0 if ndims == 3: nspat, nspat1 = spataxis.size, spat1axis.size spataxis_local = np.linspace( -dspat * nwin[0] // 2, dspat * nwin[0] // 2, nwin[0] ) spat1axis_local = np.linspace( -dspat1 * nwin[1] // 2, dspat1 * nwin[1] // 2, nwin[1] ) dimsslid = (nspat, nspat1, taxis.size) if kind == "sliding": npaxis, np1axis = paxis.size, p1axis.size Op = Radon3D( taxis, spataxis_local, spat1axis_local, paxis, p1axis, centeredh=True, kind="linear", engine=engine, ) else: npaxis, np1axis = nwin[0], nwin[1] Op = ChirpRadon3D( taxis, spataxis_local, spat1axis_local, (np.max(paxis) * dspat / dt, np.max(p1axis) * dspat1 / dt), ).H dimsp = (nwins[0] * npaxis, nwins[1] * np1axis, dimsslid[2]) Pop = Sliding3D( Op, dimsp, dimsslid, nwin, nover, (npaxis, np1axis), tapertype="cosine" ) # to be able to reshape correctly the preconditioned model dimsp = (nwins[0], nwins[1], npaxis, np1axis, dimsslid[2]) else: nspat = spataxis.size spataxis_local = np.linspace(-dspat * nwin // 2, dspat * nwin // 2, nwin) dimsslid = (nspat, taxis.size) if kind == "sliding": npaxis = paxis.size Op = Radon2D( taxis, spataxis_local, paxis, centeredh=True, kind="linear", engine=engine, ) else: npaxis = nwin Op = ChirpRadon2D(taxis, spataxis_local, np.max(paxis) * dspat / dt).H dimsp = (nwins * npaxis, dimsslid[1]) Pop = Sliding2D(Op, dimsp, dimsslid, nwin, nover, tapertype="cosine") SIop = Rop * Pop else: raise KeyError( "kind must be spatial, fk, radon-linear, " "radon-parabolic, radon-hyperbolic, sliding or chirp-sliding" ) # dot-test if dottest: Dottest( SIop, np.prod(dimsd), np.prod(dimsp) if prec else np.prod(dims), complexflag=dotcflag, raiseerror=True, verb=True, ) # inversion if kind == "spatial": recdata = regularized_inversion(SIop, data.ravel(), [Regop], **kwargs_solver) if isinstance(recdata, tuple): recdata = recdata[0] recdata = recdata.reshape(dims) recprec = None cost = None else: recprec = fista(SIop, data.ravel(), **kwargs_solver) if len(recprec) == 3: cost = recprec[2] else: cost = None recprec = recprec[0] recdata = np.real(Pop * recprec) recprec = recprec.reshape(dimsp) recdata = recdata.reshape(dims) return recdata, recprec, cost