diff --git a/pylops/waveeqprocessing/kirchhoff.py b/pylops/waveeqprocessing/kirchhoff.py index 2a6b53f4..17e8a23f 100644 --- a/pylops/waveeqprocessing/kirchhoff.py +++ b/pylops/waveeqprocessing/kirchhoff.py @@ -1,12 +1,16 @@ __all__ = ["Kirchhoff"] import logging +import os +import warnings import numpy as np from pylops import LinearOperator from pylops.signalprocessing import Convolve1D +from pylops.utils._internal import _value_or_sized_to_array from pylops.utils.decorators import reshaped +from pylops.utils.tapers import taper try: import skfmm @@ -24,6 +28,10 @@ try: 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 except ModuleNotFoundError: jit = None prange = range @@ -65,6 +73,19 @@ class Kirchhoff(LinearOperator): 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 + + Include dynamic weights in computations (``True``) or not (``False``). + Note that when ``mode=byot``, the user is required to provide such weights + in ``amp``. trav : :obj:`numpy.ndarray`, optional Traveltime table of size :math:`\lbrack (n_y) n_x n_z \times n_s n_r \rbrack` (to be provided if @@ -75,9 +96,32 @@ class Kirchhoff(LinearOperator): Amplitude table of size :math:`\lbrack (n_y) n_x n_z \times n_s n_r \rbrack` (to be provided if ``mode='byot'``) - mode : :obj:`str`, optional - Computation mode (``analytic``, ``eikonal`` or ``byot``, see Notes for - more details) + 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 introduced. See ``aperture`` for implementation + details regarding scalar and tuple cases. + + anglerefl : :obj:`np.ndarray`, optional + .. versionadded:: 2.0.0 + + Angle between the normal of the reflectors and the vertical axis in degrees + snell : :obj:`float` or :obj:`tuple`, optional + .. versionadded:: 2.0.0 + + Threshold on Snell's law evaluation. If larger, the source-receiver-image + point is discarded. If ``None``, no check on the validity of the Snell's + law is performed. See ``aperture`` for implementation details regarding + scalar and tuple cases. engine : :obj:`str`, optional Engine used for computations (``numpy`` or ``numba``). dtype : :obj:`str`, optional @@ -104,49 +148,83 @@ class Kirchhoff(LinearOperator): ----- The Kirchhoff demigration operator synthesizes seismic data given a propagation velocity model :math:`v` and a reflectivity model :math:`m`. - In forward mode: + In forward mode [1]_, [2]_: .. math:: d(\mathbf{x_r}, \mathbf{x_s}, t) = \widetilde{w}(t) * \int_V G(\mathbf{x_r}, \mathbf{x}, t) m(\mathbf{x}) G(\mathbf{x}, \mathbf{x_s}, t)\,\mathrm{d}\mathbf{x} - where :math:`m(\mathbf{x})` is the model that represents the reflectivity + where :math:`m(\mathbf{x})` represents the reflectivity at every location in the subsurface, :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 - a filtered version of the the wavelet :math:`w(t)` [1]_. In our current - implementation, the following high-frequency approximation of the Green's - functions is adopted: + a filtered version of the wavelet :math:`w(t)` [3]_ (or the wavelet itself when + ``wavfilter=False``). 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. In our current - implementations with ``mode=analytic`` and ``mode=eikonal``, this amplitude - scaling is computed as + :math:`t(\mathbf{x_r}, \mathbf{x})` is the traveltime. When ``dynamic=False`` the + amplitude is 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 :math:`a(\mathbf{x}, \mathbf{y})=\frac{1}{\|\mathbf{x} - \mathbf{y}\|}`, that is, the reciprocal of the distance between the two points, approximating the geometrical spreading of the wavefront. - In general, this multiplicative factor could contain other corrections (e.g. - obliquity factors, reflection coefficient of the incident wave at the reflector, - aperture tapers, etc.) [2]_. + Moreover an angle scaling is included in the modelling operator + added as follows: - Depending on the choice of ``mode`` the traveltime of the - Green's function will be also computed differently: + .. math:: + d(\mathbf{x_r}, \mathbf{x_s}, t) = + \tilde{w}(t) * \int_V a(\mathbf{x}, \mathbf{x_s}) a(\mathbf{x}, \mathbf{x_r}) + \frac{|cos \theta_s + cos \theta_r|} {v(\mathbf{x})} e^{j \omega (t(\mathbf{x_r}, \mathbf{x}) + + t(\mathbf{x}, \mathbf{x_s}))} m(\mathbf{x}) \,\mathrm{d}\mathbf{x} + + where :math:`\theta_s` and :math:`\theta_r` are the angles between the source-side + and receiver-side rays and the normal to the reflector at the image point (or + the vertical axis at the image point when ``reflslope=None``), respectively. + + 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``: traveltime curves between - source to receiver pairs are computed for every subsurface point and + * ``mode=analytic`` or ``mode=eikonal``: traveltimes, geometrical spreading, and angles + are computed for every source-image point-receiver triplets and the Green's functions are implemented from traveltime look-up tables, placing - the reflectivity values at corresponding source-to-receiver time in the + scaled reflectivity values at corresponding source-to-receiver time in the data. - * ``byot``: bring your own table. Traveltime table provided + * ``byot``: bring your own tables. Traveltime table are provided directly by user using ``trav`` input parameter. Similarly, in this case one - can provide their own amplitude scaling ``amp``. - - The adjoint of the demigration operator is a *migration* operator which + can provide their own amplitude scaling ``amp`` (which should include the angle + scaling too). + + Three 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 (or the normal to the reflector if ``anglerefl`` is provided. + This aperture limitation also avoid including grazing angles whose contributions + can introduce aliasing effects. Note that for a homogenous medium and slowly varying + heterogenous medium the offset and angle aperture limits may work in the same way; + * ``snell``: the maximum allowed snell's angle is expressed as the absolute value of + the sum between incident and emerging angles defined as in the ``angleaperture`` case. + This aperture limitation is introduced to turn a scattering-based Kirchhoff engine into + a reflection-based Kirchhoff engine where each image point is not considered as scatter + but as a local horizontal (or straight) reflector. + + 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. @@ -157,6 +235,9 @@ class Kirchhoff(LinearOperator): .. [2] Santos, L.T., Schleicher, J., Tygel, M., and Hubral, P. "Seismic modeling by demigration", Geophysics, 65(4), pp. 1281-1289, 2000. + .. [3] Safron, L. "Multicomponent least-squares Kirchhoff depth migration", + MSc Thesis, 2018. + """ def __init__( @@ -170,20 +251,27 @@ def __init__( wav, wavcenter, y=None, + mode="eikonal", + wavfilter=False, + dynamic=False, trav=None, amp=None, - mode="eikonal", + aperture=None, + angleaperture=90, + anglerefl=None, + snell=None, engine="numpy", dtype="float64", name="D", ): + # identify geometry ( self.ndims, _, dims, - ny, - nx, - nz, + self.ny, + self.nx, + self.nz, ns, nr, _, @@ -195,33 +283,123 @@ def __init__( dt = t[1] - t[0] self.nt = len(t) + # store ix-iy locations of sources and receivers + 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: 3D normalized distances + pass + + # compute traveltime + self.dynamic = dynamic if mode in ["analytic", "eikonal", "byot"]: if mode in ["analytic", "eikonal"]: # compute traveltime table - self.trav, _, _, dist = Kirchhoff._traveltime_table( - z, x, srcs, recs, vel, y=y, mode=mode - ) - # need to add a scalar in the denominator to avoid division by 0 - # currently set to 1/100 of max distance to avoid having to large - # scaling around the source. This number may change in future or - # left to the user to define - epsdist = 1e-2 - self.amp = 1 / (dist + epsdist * np.max(dist)) + ( + self.trav, + trav_srcs, + trav_recs, + dist, + 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 to avoid division by 0 + # currently set to 1/100 of max distance to avoid having to large + # scaling around the source. This number may change in future or + # left to the user to define + epsdist = 1e-2 + self.amp = 1 / (dist + epsdist * np.max(dist)) + + # compute angles + if self.ndims == 2: + # 2d with vertical + if anglerefl is None: + 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) + self.cosangle = np.cos(self.angle_srcs).reshape( + np.prod(dims), ns, 1 + ) + np.cos(self.angle_recs).reshape(np.prod(dims), 1, nr) + self.cosangle = self.cosangle.reshape( + np.prod(dims), ns * nr + ) + else: + # TODO: 2D with normal + raise NotImplementedError( + "angle scaling with anglerefl currently not available" + ) + + self.amp *= np.abs(self.cosangle) + if mode == "analytic": + self.amp /= vel + else: + self.amp /= vel.reshape(np.prod(dims), 1) + + else: + # TODO: 3D + raise NotImplementedError( + "dynamic=True currently not available in 3D" + ) else: self.trav = trav - self.amp = amp + if self.dynamic: + self.amp = amp + else: - raise NotImplementedError("method must be analytic or eikonal") + raise NotImplementedError("method must be analytic, eikonal or byot") self.itrav = (self.trav / dt).astype("int32") self.travd = self.trav / dt - self.itrav - self.wav = self._wavelet_reshaping(wav, dt) + + # create wavelet operator + if wavfilter: + self.wav = self._wavelet_reshaping( + wav, dt, srcs.shape[0], recs.shape[0], self.ndims + ) + else: + self.wav = wav self.cop = Convolve1D( - (ns * nr, self.nt), h=wav, offset=wavcenter, axis=1, dtype=dtype + (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 is not None: + warnings.warn( + "Aperture is currently defined as ratio of offset over depth, " + "and may be not ideal for highly heterogenous media" + ) + self.aperture = ( + (2 * self.nx / self.nz,) + 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 and snell law + 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]] + ) + self.snell = ( + (np.pi,) if snell is None else np.deg2rad(_value_or_sized_to_array(snell)) + ) + if len(self.snell) == 1: + self.snell = np.array([0.8 * self.snell[0], self.snell[0]]) + + # dimensions 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) @@ -289,12 +467,22 @@ def _traveltime_table(z, x, srcs, recs, vel, y=None, mode="eikonal"): 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 : :obj:`numpy.ndarray` + Total distance table of size + :math:`\lbrack (n_y*) n_x n_z \times n_s \rbrack` (or constant) + 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` (or constant) + 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, @@ -306,6 +494,8 @@ def _traveltime_table(z, x, srcs, recs, vel, y=None, mode="eikonal"): 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") @@ -383,46 +573,275 @@ def _traveltime_table(z, x, srcs, recs, vel, y=None, mode="eikonal"): else: raise NotImplementedError("method must be analytic or eikonal") - return trav, trav_srcs, trav_recs, dist + # compute traveltime gradients at image points + trav_srcs_grad = np.gradient( + trav_srcs.reshape(*dims, ns), axis=np.arange(ndims) + ) + trav_recs_grad = np.gradient( + trav_recs.reshape(*dims, nr), axis=np.arange(ndims) + ) + + return trav, trav_srcs, trav_recs, dist, trav_srcs_grad, trav_recs_grad - def _wavelet_reshaping(self, wav, dt): + def _wavelet_reshaping(self, wav, dt, dimsrc, dimrec, dimv): """Apply wavelet reshaping as from theory in [1]_""" f = np.fft.rfftfreq(len(wav), dt) W = np.fft.rfft(wav, n=len(wav)) - if self.ndims == 2: - Wfilt = W * np.abs(2 * np.pi * f) - else: + if dimsrc == 2 and dimv == 2: + # 2D + Wfilt = W * (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 _kirch_matvec(x, y, nsnr, nt, ni, itrav, travd, amp): + def _trav_kirch_matvec(x, y, nsnr, nt, ni, itrav, travd): for isrcrec in prange(nsnr): itravisrcrec = itrav[:, isrcrec] travdisrcrec = travd[:, isrcrec] - ampisrcrec = amp[:, isrcrec] for ii in range(ni): index = itravisrcrec[ii] dindex = travdisrcrec[ii] - damp = ampisrcrec[ii] if 0 <= index < nt - 1: - y[isrcrec, index] += x[ii] * (1 - dindex) * damp - y[isrcrec, index + 1] += x[ii] * dindex * damp + y[isrcrec, index] += x[ii] * (1 - dindex) + y[isrcrec, index + 1] += x[ii] * dindex return y @staticmethod - def _kirch_rmatvec(x, y, nsnr, nt, ni, itrav, travd, amp): + def _trav_kirch_rmatvec(x, y, nsnr, nt, ni, itrav, travd): for ii in prange(ni): itravii = itrav[ii] travdii = travd[ii] - ampii = amp[ii] for isrcrec in range(nsnr): if 0 <= itravii[isrcrec] < nt - 1: y[ii] += ( x[isrcrec, itravii[isrcrec]] * (1 - travdii[isrcrec]) + x[isrcrec, itravii[isrcrec] + 1] * travdii[isrcrec] - ) * ampii[isrcrec] + ) + return y + + @staticmethod + def _amp_kirch_matvec( + x, + y, + nsnr, + nt, + ni, + itrav, + travd, + amp, + aperturemin, + aperturemax, + aperturetap, + nz, + six, + rix, + angleaperturemin, + angleaperturemax, + angles_srcs, + angles_recs, + snellmin, + snellmax, + ): + nr = angles_recs.shape[-1] + daperture = aperturemax - aperturemin + dangleaperture = angleaperturemax - angleaperturemin + dsnell = snellmax - snellmin + for isrcrec in prange(nsnr): + # extract traveltime, amplitude, src/rec coordinates at given src/pair + itravisrcrec = itrav[:, isrcrec] + travdisrcrec = travd[:, isrcrec] + ampisrcrec = amp[:, isrcrec] + sixisrcrec = six[isrcrec] + rixisrcrec = rix[isrcrec] + # extract source and receiver angles + angles_src = angles_srcs[:, isrcrec // nr] + angles_rec = angles_recs[:, isrcrec % nr] + for ii in range(ni): + # extract traveltime, amplitude at given image point + index = itravisrcrec[ii] + dindex = travdisrcrec[ii] + damp = ampisrcrec[ii] + # extract source and receiver angle + angle_src = angles_src[ii] + angle_rec = angles_rec[ii] + abs_angle_src = abs(angle_src) + abs_angle_rec = abs(angle_rec) + abs_angle_src_rec = abs(angle_src + angle_rec) + aptap = 1.0 + # angle apertures checks + if ( + abs_angle_src < angleaperturemax + and abs_angle_rec < angleaperturemax + and abs_angle_src_rec < snellmax + ): + 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 + ) + ] + ) + if abs_angle_src_rec >= snellmin: + # extract snell taper value + aptap = ( + aptap + * aperturetap[ + int(20 * (abs_angle_src_rec - snellmin) // dsnell) + ] + ) + + # identify x-index of image point + iz = ii % nz + # aperture check + aperture = abs(sixisrcrec - rixisrcrec) / iz + if aperture < aperturemax: + if aperture >= aperturemin: + # extract aperture taper value + aptap = ( + aptap + * aperturetap[ + int(20 * ((aperture - aperturemin) // daperture)) + ] + ) + # time limit check + if 0 <= index < nt - 1: + # assign values + y[isrcrec, index] += x[ii] * (1 - dindex) * damp * aptap + y[isrcrec, index + 1] += x[ii] * dindex * damp * aptap + return y + + @staticmethod + def _amp_kirch_rmatvec( + x, + y, + nsnr, + nt, + ni, + itrav, + travd, + amp, + aperturemin, + aperturemax, + aperturetap, + nz, + six, + rix, + angleaperturemin, + angleaperturemax, + angles_srcs, + angles_recs, + snellmin, + snellmax, + ): + nr = angles_recs.shape[-1] + daperture = aperturemax - aperturemin + dangleaperture = angleaperturemax - angleaperturemin + dsnell = snellmax - snellmin + for ii in prange(ni): + itravii = itrav[ii] + travdii = travd[ii] + ampii = amp[ii] + # extract source and receiver angles + angle_srcs = angles_srcs[ii] + angle_recs = angles_recs[ii] + # identify x-index of image point + iz = ii % nz + for isrcrec in range(nsnr): + index = itravii[isrcrec] + dindex = travdii[isrcrec] + sixisrcrec = six[isrcrec] + rixisrcrec = rix[isrcrec] + # extract source and receiver angle + angle_src = angle_srcs[isrcrec // nr] + angle_rec = angle_recs[isrcrec % nr] + abs_angle_src = abs(angle_src) + abs_angle_rec = abs(angle_rec) + abs_angle_src_rec = abs(angle_src + angle_rec) + aptap = 1.0 + # angle apertures checks + if ( + abs_angle_src < angleaperturemax + and abs_angle_rec < angleaperturemax + and abs_angle_src_rec < snellmax + ): + 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 + ) + ] + ) + if abs_angle_src_rec >= snellmin: + # extract snell taper value + aptap = ( + aptap + * aperturetap[ + int(20 * (abs_angle_src_rec - snellmin) // dsnell) + ] + ) + + # aperture check + aperture = abs(sixisrcrec - rixisrcrec) / iz + if aperture < aperturemax: + if aperture >= aperturemin: + # extract aperture taper value + aptap = ( + aptap + * aperturetap[ + int(20 * ((aperture - aperturemin) // daperture)) + ] + ) + # time limit check + if 0 <= index < nt - 1: + # assign values + y[ii] += ( + ( + x[isrcrec, index] * (1 - dindex) + + x[isrcrec, index + 1] * dindex + ) + * ampii[isrcrec] + * aptap + ) return y def _register_multiplications(self, engine): @@ -431,20 +850,53 @@ def _register_multiplications(self, engine): if engine == "numba" and jit is not None: # numba numba_opts = dict( - nopython=True, nogil=True, parallel=True + nopython=True, nogil=True, parallel=parallel ) # fastmath=True, - self._kirch_matvec = jit(**numba_opts)(self._kirch_matvec) - self._kirch_rmatvec = jit(**numba_opts)(self._kirch_rmatvec) + if self.dynamic: + self._kirch_matvec = jit(**numba_opts)(self._amp_kirch_matvec) + self._kirch_rmatvec = jit(**numba_opts)(self._amp_kirch_rmatvec) + else: + 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 is None: logging.warning(jit_message) + if self.dynamic: + self._kirch_matvec = self._amp_kirch_matvec + self._kirch_rmatvec = self._amp_kirch_rmatvec + else: + self._kirch_matvec = self._trav_kirch_matvec + self._kirch_rmatvec = self._trav_kirch_rmatvec @reshaped def _matvec(self, x): y = np.zeros((self.nsnr, self.nt), dtype=self.dtype) - y = self._kirch_matvec( - x.ravel(), y, self.nsnr, self.nt, self.ni, self.itrav, self.travd, self.amp - ) + if self.dynamic: + inputs = ( + x.ravel(), + y, + self.nsnr, + self.nt, + self.ni, + self.itrav, + self.travd, + self.amp, + 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, + self.snell[0], + self.snell[1], + ) + else: + 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 @@ -453,7 +905,30 @@ def _rmatvec(self, x): x = self.cop._rmatvec(x.ravel()) x = x.reshape(self.nsnr, self.nt) y = np.zeros(self.ni, dtype=self.dtype) - y = self._kirch_rmatvec( - x, y, self.nsnr, self.nt, self.ni, self.itrav, self.travd, self.amp - ) + if self.dynamic: + inputs = ( + x, + y, + self.nsnr, + self.nt, + self.ni, + self.itrav, + self.travd, + self.amp, + 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, + self.snell[0], + self.snell[1], + ) + else: + inputs = (x, y, self.nsnr, self.nt, self.ni, self.itrav, self.travd) + y = self._kirch_rmatvec(*inputs) return y diff --git a/pylops/waveeqprocessing/lsm.py b/pylops/waveeqprocessing/lsm.py index 0df8cb38..7aab2e7c 100644 --- a/pylops/waveeqprocessing/lsm.py +++ b/pylops/waveeqprocessing/lsm.py @@ -54,15 +54,11 @@ class LSM: y : :obj:`numpy.ndarray` Additional spatial axis (for 3-dimensional problems) kind : :str`, optional - Kind of modelling operator (``kirchhoff``, ``twowayac``) - mode : :obj:`str`, optional - Computation mode (``eikonal``, ``analytic`` - only for - constant velocity) - engine : :obj:`str`, optional - Engine used for computations (``numpy`` or ``numba``) when ``kind=kirchhoff`` - is used + Kind of modelling operator (``kirchhoff``, ``twoway``) dottest : :obj:`bool`, optional Apply dot-test + **kwargs_mod : :obj:`int`, optional + Additional arguments to pass to modelling operators Attributes ---------- @@ -124,25 +120,14 @@ def __init__( wavcenter, y=None, kind="kirchhoff", - mode="eikonal", - engine="numba", dottest=False, + **kwargs_mod, ): self.y, self.x, self.z = y, x, z if kind == "kirchhoff": self.Demop = Kirchhoff( - z, - x, - t, - srcs, - recs, - vel, - wav, - wavcenter, - y=y, - mode=mode, - engine=engine, + z, x, t, srcs, recs, vel, wav, wavcenter, y=y, **kwargs_mod ) elif kind == "twowayac": shape = (len(x), len(z)) @@ -159,6 +144,7 @@ def __init__( recs[1], t[0], len(t), + **kwargs_mod, ) else: diff --git a/pytests/test_kirchhoff.py b/pytests/test_kirchhoff.py index ebc7f882..8cc9ccd2 100755 --- a/pytests/test_kirchhoff.py +++ b/pytests/test_kirchhoff.py @@ -7,7 +7,7 @@ from pylops.waveeqprocessing.kirchhoff import Kirchhoff PAR = { - "ny": 10, + "ny": 3, "nx": 12, "nz": 20, "nt": 50, @@ -16,9 +16,9 @@ "dz": 2, "dt": 0.004, "nsy": 4, - "nry": 8, + "nry": 3, "nsx": 6, - "nrx": 4, + "nrx": 2, } # Check if skfmm is available and by-pass tests using it otherwise. This is @@ -51,9 +51,12 @@ wav, _, wavc = ricker(t[:41], f0=40) -par1 = {"mode": "analytic"} -par2 = {"mode": "eikonal"} -par3 = {"mode": "byot"} +par1 = {"mode": "analytic", "dynamic": False} +par2 = {"mode": "eikonal", "dynamic": False} +par3 = {"mode": "byot", "dynamic": False} +par1d = {"mode": "analytic", "dynamic": True} +par2d = {"mode": "eikonal", "dynamic": True} +par3d = {"mode": "byot", "dynamic": True} def test_identify_geometry(): @@ -120,7 +123,7 @@ def test_traveltime_ana(): """ src = np.array([100, 0])[:, np.newaxis] - _, trav_srcs_ana, trav_recs_ana, dist_ana = Kirchhoff._traveltime_table( + _, trav_srcs_ana, trav_recs_ana, dist_ana, _, _ = Kirchhoff._traveltime_table( np.arange(0, 200, 1), np.arange(0, 200, 1), src, src, v0, mode="analytic" ) assert dist_ana[0, 0] == 200 @@ -132,11 +135,23 @@ def test_traveltime_table(): """Compare analytical and eikonal traveltimes in homogenous medium""" if skfmm_enabled: # 2d - trav_ana, trav_srcs_ana, trav_recs_ana, dist_ana = Kirchhoff._traveltime_table( - z, x, s2d, r2d, v0, mode="analytic" - ) + ( + trav_ana, + trav_srcs_ana, + trav_recs_ana, + dist_ana, + _, + _, + ) = Kirchhoff._traveltime_table(z, x, s2d, r2d, v0, mode="analytic") - trav_eik, trav_srcs_eik, trav_recs_eik, dist_eik = Kirchhoff._traveltime_table( + ( + trav_eik, + trav_srcs_eik, + trav_recs_eik, + dist_eik, + _, + _, + ) = Kirchhoff._traveltime_table( z, x, s2d, r2d, v0 * np.ones((PAR["nx"], PAR["nz"])), mode="eikonal" ) @@ -145,11 +160,23 @@ def test_traveltime_table(): assert_array_almost_equal(trav_ana, trav_eik, decimal=2) # 3d - trav_ana, trav_srcs_ana, trav_recs_ana, dist_ana = Kirchhoff._traveltime_table( - z, x, s3d, r3d, v0, y=y, mode="analytic" - ) + ( + trav_ana, + trav_srcs_ana, + trav_recs_ana, + dist_ana, + _, + _, + ) = Kirchhoff._traveltime_table(z, x, s3d, r3d, v0, y=y, mode="analytic") - trav_eik, trav_srcs_eik, trav_recs_eik, dist_eik = Kirchhoff._traveltime_table( + ( + trav_eik, + trav_srcs_eik, + trav_recs_eik, + dist_eik, + _, + _, + ) = Kirchhoff._traveltime_table( z, x, s3d, @@ -164,13 +191,13 @@ def test_traveltime_table(): assert_array_almost_equal(trav_ana, trav_eik, decimal=2) -@pytest.mark.parametrize("par", [(par1), (par2), (par3)]) +@pytest.mark.parametrize("par", [(par1), (par2), (par3), (par1d), (par2d), (par3d)]) def test_kirchhoff2d(par): """Dot-test for Kirchhoff operator""" vel = v0 * np.ones((PAR["nx"], PAR["nz"])) if par["mode"] == "byot": - trav, _, _, dist = Kirchhoff._traveltime_table( + trav, _, _, dist, _, _ = Kirchhoff._traveltime_table( z, x, s2d, r2d, v0, mode="analytic" ) amp = 1 / (dist + 1e-2 * dist.max()) @@ -194,3 +221,39 @@ def test_kirchhoff2d(par): mode=par["mode"], ) assert dottest(Dop, PAR["nsx"] * PAR["nrx"] * PAR["nt"], PAR["nz"] * PAR["nx"]) + + +@pytest.mark.parametrize("par", [(par1), (par2), (par3)]) +def test_kirchhoff3d(par): + """Dot-test for Kirchhoff operator""" + vel = v0 * np.ones((PAR["ny"], PAR["nx"], PAR["nz"])) + + if par["mode"] == "byot": + trav, _, _, dist, _, _ = Kirchhoff._traveltime_table( + z, x, s3d, r3d, v0, y=y, mode="analytic" + ) + amp = 1 / (dist + 1e-2 * dist.max()) + else: + trav = None + amp = None + + if skfmm_enabled or par["mode"] != "eikonal": + Dop = Kirchhoff( + z, + x, + t, + s3d, + r3d, + vel if par["mode"] == "eikonal" else v0, + wav, + wavc, + y=y, + trav=trav, + amp=amp, + mode=par["mode"], + ) + assert dottest( + Dop, + PAR["nsx"] * PAR["nrx"] * PAR["nsy"] * PAR["nry"] * PAR["nt"], + PAR["nz"] * PAR["nx"] * PAR["ny"], + ) diff --git a/pytests/test_lsm.py b/pytests/test_lsm.py index 26656dfe..33024710 100755 --- a/pytests/test_lsm.py +++ b/pytests/test_lsm.py @@ -50,8 +50,10 @@ wav, _, wavc = ricker(t[:41], f0=40) -par1 = {"mode": "analytic"} -par2 = {"mode": "eikonal"} +par1 = {"mode": "analytic", "dynamic": False} +par2 = {"mode": "eikonal", "dynamic": False} +par1d = {"mode": "analytic", "dynamic": True} +par2d = {"mode": "eikonal", "dynamic": True} def test_unknown_mode(): @@ -60,7 +62,7 @@ def test_unknown_mode(): _ = LSM(z, x, t, s2d, r2d, 0, np.ones(3), 1, mode="foo") -@pytest.mark.parametrize("par", [(par1), (par2)]) +@pytest.mark.parametrize("par", [(par1), (par2), (par1d), (par2d)]) def test_lsm2d(par): """Dot-test and inverse for LSM operator""" if skfmm_enabled or par["mode"] != "eikonal": @@ -79,21 +81,18 @@ def test_lsm2d(par): wav, wavc, mode=par["mode"], + dynamic=par["dynamic"], dottest=True, ) - # Try both v1 and v2 versions of amp - for amp in [lsm.Demop.amp, np.ones_like(lsm.Demop.trav)]: - lsm.Demop.amp = amp + d = lsm.Demop * refl.ravel() + d = d.reshape(PAR["nsx"], PAR["nrx"], PAR["nt"]) - d = lsm.Demop * refl.ravel() - d = d.reshape(PAR["nsx"], PAR["nrx"], PAR["nt"]) + minv = lsm.solve(d.ravel(), **dict(iter_lim=100, show=True)) + minv = minv.reshape(PAR["nx"], PAR["nz"]) - minv = lsm.solve(d.ravel(), **dict(iter_lim=100, show=True)) - minv = minv.reshape(PAR["nx"], PAR["nz"]) + dinv = lsm.Demop * minv.ravel() + dinv = dinv.reshape(PAR["nsx"], PAR["nrx"], PAR["nt"]) - dinv = lsm.Demop * minv.ravel() - dinv = dinv.reshape(PAR["nsx"], PAR["nrx"], PAR["nt"]) - - assert_array_almost_equal(d, dinv, decimal=1) - assert_array_almost_equal(refl, minv, decimal=1) + assert_array_almost_equal(d / d.max(), dinv / d.max(), decimal=2) + assert_array_almost_equal(refl / refl.max(), minv / refl.max(), decimal=1) diff --git a/tutorials/lsm.py b/tutorials/lsm.py index 8d5877f8..7f9fb4b2 100755 --- a/tutorials/lsm.py +++ b/tutorials/lsm.py @@ -95,7 +95,7 @@ ############################################################################### # We can now create our LSM object and invert for the reflectivity using two # different solvers: :py:func:`scipy.sparse.linalg.lsqr` (LS solution) and -# :py:func:`pylops.optimization.sparsity.FISTA` (LS solution with sparse model). +# :py:func:`pylops.optimization.sparsity.fista` (LS solution with sparse model). nt = 651 dt = 0.004 t = np.arange(nt) * dt @@ -103,7 +103,16 @@ lsm = pylops.waveeqprocessing.LSM( - z, x, t, sources, recs, v0, wav, wavc, mode="analytic" + z, + x, + t, + sources, + recs, + v0, + wav, + wavc, + mode="analytic", + engine="numba", ) d = lsm.Demop * refl