From a227115035465a9196bcc7f716ff1a7013c548cf Mon Sep 17 00:00:00 2001 From: mrava87 Date: Tue, 28 Nov 2023 21:16:09 +0300 Subject: [PATCH] feature: handles https://github.com/PyLops/pylops/issues/554 --- pylops/waveeqprocessing/oneway.py | 42 ++++++++++++++++++++++-------- pytests/test_oneway.py | 43 ++++++++++++++++++------------- 2 files changed, 56 insertions(+), 29 deletions(-) diff --git a/pylops/waveeqprocessing/oneway.py b/pylops/waveeqprocessing/oneway.py index ee4d9474..b41779db 100644 --- a/pylops/waveeqprocessing/oneway.py +++ b/pylops/waveeqprocessing/oneway.py @@ -193,6 +193,7 @@ def Deghosting( dr: Sequence[float], vel: float, zrec: float, + kind: Optional[str] = "p", pd: Optional[NDArray] = None, win: Optional[NDArray] = None, npad: Union[Tuple[int], Tuple[int, int]] = (11, 11), @@ -206,13 +207,15 @@ def Deghosting( ) -> Tuple[NDArray, NDArray]: r"""Wavefield deghosting. - Apply seismic wavefield decomposition from single-component (pressure) - data. This process is also generally referred to as model-based deghosting. + Apply seismic wavefield decomposition from single-component (pressure or + vertical velocity) 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}) + Pressure (or vertical velocity) 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 @@ -231,6 +234,8 @@ def Deghosting( Velocity along the receiver array (must be constant) zrec : :obj:`float` Depth of receiver array + kind : :obj:`str`, optional + Type of data (``p`` or ``vz``) pd : :obj:`np.ndarray`, optional Direct arrival to be subtracted from ``p`` win : :obj:`np.ndarray`, optional @@ -260,14 +265,19 @@ def Deghosting( Returns ------- pup : :obj:`np.ndarray` - Up-going wavefield + Up-going pressure (or particle velocity) wavefield pdown : :obj:`np.ndarray` - Down-going wavefield + Down-going (or particle velocity) wavefield + + Raises + ------ + ValueError + If ``kind`` is not "p" or "vz". 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 + The up- and down-going components of a 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 @@ -284,16 +294,22 @@ def Deghosting( In a matrix form we can thus write the total wavefield as: .. math:: - \mathbf{p} - \mathbf{p_d} = (\mathbf{I} + \Phi) \mathbf{p}^- + \mathbf{p} - \mathbf{p_d} = (\mathbf{I} \pm \Phi) \mathbf{p}^- where :math:`\Phi` is one-way propagator implemented via the - :class:`pylops.waveeqprocessing.PhaseShift` operator. + :class:`pylops.waveeqprocessing.PhaseShift` operator. Note that :math:`+` is + used for the pressure data, whilst :math:`-` is used for the vertical velocity + data. .. [1] Amundsen, L., 1993, Wavenumber-based filtering of marine point-source data: GEOPHYSICS, 58, 1335–1348. - """ + # Check kind + if kind not in ["p", "vz"]: + raise ValueError("kind must be p or vz") + + # Identify dimensions ndims = p.ndim if ndims == 2: dims = (nt, nr) @@ -328,7 +344,11 @@ def Deghosting( ) # Decomposition operator - Dupop = Identity(nt * nrs, dtype=p.dtype) + Pop + if kind == "p": + Dupop = Identity(nt * nrs, dtype=p.dtype) + Pop + else: + Dupop = Identity(nt * nrs, dtype=p.dtype) - Pop + if dottest: Dottest(Dupop, nt * nrs, nt * nrs, verb=True) diff --git a/pytests/test_oneway.py b/pytests/test_oneway.py index 2e8162d4..48f73a9e 100755 --- a/pytests/test_oneway.py +++ b/pytests/test_oneway.py @@ -22,8 +22,10 @@ "f0": 40, } -par1 = {"ny": 8, "nx": 10, "nt": 20, "dtype": "float32"} # even -par2 = {"ny": 9, "nx": 11, "nt": 21, "dtype": "float32"} # odd +par1 = {"ny": 8, "nx": 10, "nt": 20, "kind": "p", "dtype": "float32"} # even, p +par2 = {"ny": 9, "nx": 11, "nt": 21, "kind": "p", "dtype": "float32"} # odd, p +par1v = {"ny": 8, "nx": 10, "nt": 20, "kind": "vz", "dtype": "float32"} # even, vz +par2v = {"ny": 9, "nx": 11, "nt": 21, "kind": "vz", "dtype": "float32"} # odd, vz # deghosting params vel_sep = 1000.0 # velocity at separation level @@ -34,27 +36,31 @@ wav = ricker(t[:41], f0=parmod["f0"])[0] -@pytest.fixture(scope="module") +@pytest.fixture def create_data2D(): """Create 2d dataset""" - t0_plus = np.array([0.02, 0.08]) - t0_minus = t0_plus + 0.04 - vrms = np.array([1400.0, 1800.0]) - amp = np.array([1.0, -0.6]) - p2d_minus = hyperbolic2d(x, t, t0_minus, vrms, amp, wav)[1].T + def core(datakind): + t0_plus = np.array([0.02, 0.08]) + t0_minus = t0_plus + 0.04 + vrms = np.array([1400.0, 1800.0]) + amp = np.array([1.0, -0.6]) - kx = np.fft.ifftshift(np.fft.fftfreq(parmod["nx"], parmod["dx"])) - freq = np.fft.rfftfreq(parmod["nt"], parmod["dt"]) + p2d_minus = hyperbolic2d(x, t, t0_minus, vrms, amp, wav)[1].T - Pop = -PhaseShift(vel_sep, 2 * zrec, parmod["nt"], freq, kx) + kx = np.fft.ifftshift(np.fft.fftfreq(parmod["nx"], parmod["dx"])) + freq = np.fft.rfftfreq(parmod["nt"], parmod["dt"]) - # Decomposition operator - Dupop = Identity(parmod["nt"] * parmod["nx"]) + Pop + Pop = -PhaseShift(vel_sep, 2 * zrec, parmod["nt"], freq, kx) - p2d = Dupop * p2d_minus.ravel() - p2d = p2d.reshape(parmod["nt"], parmod["nx"]) - return p2d, p2d_minus + # Decomposition operator + Dupop = Identity(parmod["nt"] * parmod["nx"]) + datakind * Pop + + p2d = Dupop * p2d_minus.ravel() + p2d = p2d.reshape(parmod["nt"], parmod["nx"]) + return p2d, p2d_minus + + return core @pytest.mark.parametrize("par", [(par1), (par2)]) @@ -87,10 +93,10 @@ def test_PhaseShift_3dsignal(par): ) -@pytest.mark.parametrize("par", [(par1), (par2)]) +@pytest.mark.parametrize("par", [(par1), (par2), (par1v), (par2v)]) def test_Deghosting_2dsignal(par, create_data2D): """Deghosting of 2d data""" - p2d, p2d_minus = create_data2D + p2d, p2d_minus = create_data2D(1 if par["kind"] is "p" else -1) p2d_minus_inv, p2d_plus_inv = Deghosting( p2d, @@ -100,6 +106,7 @@ def test_Deghosting_2dsignal(par, create_data2D): parmod["dx"], vel_sep, zrec, + kind=par["kind"], win=np.ones_like(p2d), npad=0, ntaper=0,