Merge pull request #555 from mrava87/feat-vzdegh

feature: vz deghosting

Author: mrava87

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)
 

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,