Spatial wavefield gradient-based seismic wavefield separation

Spatial wavefield gradient-based seismic wavefield separation Summary Measurements of the horizontal and vertical components of particle motion combined with estimates of the spatial gradients of the seismic wavefield enable seismic data to be acquired and processed using single dedicated multicomponent stations (e.g. rotational sensors) and/or small receiver groups instead of large receiver arrays. Here, we present seismic wavefield decomposition techniques that use spatial wavefield gradient data to separate land and ocean bottom data into their upgoing/downgoing and P/S constituents. Our method is based on the elastodynamic representation theorem with the derived filters requiring local measurements of the wavefield and its spatial gradients only. We demonstrate with synthetic data and a land seismic field data example that combining translational measurements with spatial wavefield gradient estimates allows separating seismic data recorded either at the Earth’s free-surface or at the sea bottom into upgoing/downgoing and P/S wavefield constituents for typical incidence angle ranges of body waves. A key finding is that the filter application only requires knowledge of the elastic properties exactly at the recording locations and is valid for a wide elastic property range. Numerical modelling, Body waves, Rotational seismology, Wave propagation 1 INTRODUCTION The decomposition of multicomponent land and ocean bottom seismic (OBS) recordings into their upgoing/downgoing and P/S wavefield constituents is a key processing step for multicomponent data analysis. Important applications of wavefield separation techniques are water layer multiple removal and true amplitude imaging. However, wavefield separation of seismic data remains one of the main challenges in seismic data processing to date. The simplest conventional approach to multicomponent data wavefield separation is based on the assumption that the vertical component measurements correspond to P-wave recordings, whereas S-waves are recorded solely on the horizontal components. However, this assumption breaks down at non-normal incidence angles for which P- and S-waves are recorded both on the vertical and horizontal components. More sophisticated approaches to the separation of wave modes on multicomponent recordings have been proposed using array-based filtering techniques (e.g. Dankbaar 1985; Wapenaar et al. 1990; Wapenaar & Haimé 1990; Amundsen et al. 2000). Dankbaar (1985) derived P/S separation filters as a function of the geophone receiving characteristics. Wapenaar et al. (1990) formulated a wave equation-based approach to decompose multicomponent surface seismic data into their primary P- and S-wave components. Wapenaar & Haimé (1990) used the elastodynamic representation theorem to derive upgoing/downgoing and P/S wavefield decomposition filters. Amundsen et al. (2000) extended this approach to multicomponent OBS recordings for receiver-side water layer multiple removal and P/S separation. Array-based wavefield separation filters typically require long arrays of densely spaced multicomponent stations, which significantly increases the field effort and cost and is not feasible in areas with limited access (e.g. the ocean bottom). Additionally, the application of the filters requires knowledge of the near-surface elastic properties. Lateral variations in elastic properties and statics are likely to take place along such extended receiver lines and can result in incorrect separation results. These limitations of array-based decomposition techniques generally render them less practical for field applications. The use of spatial wavefield gradient data as a new observable in seismic exploration has gained increasing attention in recent years and could open up new ways of acquiring and processing multicomponent seismic data (Igel et al. 2012; Schmelzbach et al. 2016b). If multicomponent receivers are closely spaced, not only the vector wavefield but also the spatial gradients of the seismic wavefield can be estimated from the array data (e.g. Spudich et al. 1995; Langston 2007a,b; Liu & Holt 2015; Porter et al. 2016; Sollberger et al. 2016a). Wavefield separation techniques based on spatial wavefield gradient data have been proposed by Robertsson & Muyzert (1999) and Robertsson & Curtis (2002), using spatial wavefield gradient estimates from spatially compact groups of three-component (3C) receivers. Robertsson & Muyzert (1999) proposed a compact volume acquisition design consisting of only four stations to separate the curl- and divergence-free components of the wavefield, but their method requires the presence of a receiver at depth to compute the vertical spatial derivatives. A major step forward was the realization by Robertsson & Curtis (2002) that the free-surface boundary condition allows expressing the vertical derivatives of the particle velocity components in terms of horizontal derivatives. Hence, recording stations at depth to estimate vertical gradients are not required to isolate the upgoing P- and S-waves of land surface seismic data. In multicomponent OBS exploration, spatial wavefield gradient data have mainly been used for deghosting of pressure data (e.g. Osen et al. 2002; Amundsen et al. 2013). Barak et al. (2013, 2015) recently discussed the potential benefit of using rotational data acquired on the sea bottom for enhanced wavefield characterization. In this paper, we present spatial wavefield gradient-based techniques to separate land and OBS recordings into their (1) upgoing/downgoing and (2) P/S wavefield constituents. Up/down separation of OBS data is important for water layer multiple suppression, whereas in the land seismic case we discuss up/down decomposition as an important step towards P/S separation for imaging. In contrast to Robertsson & Curtis (2002), we present filters to isolate the horizontal and vertical particle velocity components of the upgoing S- and P-waves, respectively, which represent more straightforward quantities for interpretation instead of the originally proposed P- and S-wave potentials. For OBS data, we extend the seabed wavefield decomposition filters by Amundsen et al. (2000) to become applicable to local measurements of the wavefield and its spatial gradients only. Furthermore, we investigate the effect of strong near-surface velocity gradients on the wavefield separation results and make the important observation that the P- and S-wave velocities required for our technique are only those exactly at the receiver locations. Finally, we apply our method to field data acquired on land, which, to the best of our knowledge, is the first application of spatial wavefield gradient-based filters to land seismic data. 2 METHODOLOGY Our wavefield separation filters are derived from the elastodynamic representation theorem, which allows to calculate the seismic wavefield inside a volume from recordings along the surface enclosing the volume (Wapenaar & Haimé 1990; Amundsen et al. 2000). Adjusting the surface geometry with corresponding boundary conditions to either the free-surface (zero vertical traction) or the fluid/solid boundary at the ocean bottom (zero tangential stress) and introducing analytical expressions for the Green’s tensor in a homogeneous medium allow for the derivation of formulations for the upgoing wavefield at an infinitesimal distance below the recording interface. 2.1 Wavefield separation of land seismic data using spatial wavefield gradients 2.1.1 Up/down separation Land seismic data acquired at the Earth’s free-surface comprise of superimposed recordings of upgoing waves (P and S) and their downgoing reflections and mode-conversions (Aki & Richards 2002). As a consequence, the recorded composite ground motion differs in polarization from the incident waves dependent on the incidence angle and the medium properties. Up/down separation of land seismic data allows to extract the true amplitude upgoing wavefield and can therefore significantly improve, for example, amplitude-versus-offset analysis and full-waveform inversion. Expressions for the upgoing horizontal ($$v_x^{U}$$ and $$v_y^{U}$$) and vertical ($$v_z^{U}$$) particle velocity components that are exact for homogeneous and isotropic media at an infinitesimal distance below the Earth’s surface are given by (Wapenaar & Haimé 1990):   \begin{equation} v_x^{U} = \frac{1}{2}\left(v_x + F^{v_x}_{v_z} \ast v_z\right), \end{equation} (1)  \begin{equation} v_y^{U} = \frac{1}{2}\left(v_y + F^{v_y}_{v_z} \ast v_z\right) , \end{equation} (2)  \begin{equation} v_z^{U} = \frac{1}{2}\left(v_z - F^{v_z}_{v_x} \ast v_x - F^{v_z}_{v_y} \ast v_y\right) , \end{equation} (3)where vx and vy, and vz are the total horizontal and vertical particle velocities, respectively, recorded at the free-surface. The superscripts of the filter F denote the wavefield component to be extracted, whereas the subscripts denote the wavefield quantity on which the filter acts. Spatiotemporal convolutions are denoted by an asterisk. In the wavenumber–frequency domain, the filters can be expressed as (Wapenaar & Haimé 1990):   \begin{equation} \widetilde{F}^{v_x}_{v_z} = \frac{k_x}{k_z^{(\alpha )}} \left(1-\frac{2}{k_\beta ^2}\left(k_x^2 + k_y^2 + k_z^{(\alpha )}k_z^{(\beta )}\right)\right), \end{equation} (4)  \begin{equation} \widetilde{F}^{v_y}_{v_z} = \frac{k_y}{k_z^{(\alpha )}} \left(1-\frac{2}{k_\beta ^2}\left(k_x^2 + k_y^2 + k_z^{(\alpha )}k_z^{(\beta )}\right)\right) , \end{equation} (5)  \begin{equation} \widetilde{F}^{v_z}_{v_x} = \frac{k_x}{k_z^{(\beta )}} \left(1-\frac{2}{k_\beta ^2}\left(k_x^2 + k_y^2 + k_z^{(\alpha )}k_z^{(\beta )}\right)\right), \end{equation} (6)  \begin{equation} \widetilde{F}^{v_z}_{v_y} = \frac{k_y}{k_z^{(\beta )}} \left(1-\frac{2}{k_\beta ^2}\left(k_x^2 + k_y^2 + k_z^{(\alpha )}k_z^{(\beta )}\right)\right), \end{equation} (7)where $$k_z^{(\alpha )} = k_\alpha \sqrt{1-\frac{k_x^2 + k_y^2}{k_\alpha ^2}}$$ and $$k_z^{(\beta )} = k_\beta \sqrt{1-\frac{k_x^2 + k_y^2}{k_\beta ^2}}$$ are the unsigned vertical wavenumbers for P- and S-waves, respectively, and kx and ky are the horizontal wavenumbers. The P- and S-wavenumbers are denoted by kα = ω/α and kβ = ω/β, respectively, where α is the P-wave velocity, β the S-wave velocity and ω the angular frequency. The implementation of eqs (4)–(7) would require horizontal wavenumber estimates conventionally derived from long receiver arrays. Instead, we seek formulations that only contain first-order horizontal spatial derivatives of the wavefield (corresponding to kx and ky) and that can be estimated using spatially compact receiver arrays. This is achieved by keeping the first term in the Taylor approximations around kx = 0 and ky = 0 to eqs (4)–(7):   \begin{equation} \widetilde{F}^{v_x}_{v_z} \approx \frac{k_x}{k_\alpha }\left(1 - \frac{2\beta }{\alpha }\right) , \end{equation} (8)  \begin{equation} \widetilde{F}^{v_y}_{v_z} \approx \frac{k_y}{k_\alpha }\left(1 - \frac{2\beta }{\alpha }\right) , \end{equation} (9)  \begin{equation} \widetilde{F}^{v_z}_{v_x} \approx \frac{k_x}{k_\beta }\left(1 - \frac{2\beta }{\alpha }\right), \end{equation} (10)  \begin{equation} \widetilde{F}^{v_z}_{v_y} \approx \frac{k_y}{k_\beta }\left(1 - \frac{2\beta }{\alpha }\right). \end{equation} (11) Inserting eqs (8)–(11) into eqs (1)–(3) and transforming the expressions back to the space–frequency domain yields the following expressions for the upgoing particle velocity components:   \begin{equation} v_x^U \approx \frac{1}{2}\left(v_x + \frac{1}{i\omega }(\alpha -2\beta )\frac{\partial v_z}{\partial x}\right), \end{equation} (12)  \begin{equation} v_y^U \approx \frac{1}{2}\left(v_y + \frac{1}{i\omega }(\alpha -2\beta )\frac{\partial v_z}{\partial y}\right), \end{equation} (13)  \begin{equation} v_z^U \approx \frac{1}{2}\left(v_z - \frac{1}{i\omega }\left(\beta -\frac{2\beta ^2}{\alpha }\right)\left(\frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y}\right)\right). \end{equation} (14) Eqs (12)–(14) combine the total recorded wavefield with spatial gradient estimates scaled by the seismic velocities to suppress the downgoing reflections and mode conversions generated at the Earth’s surface. It is also worthwhile to note that the first terms in eqs (12)–(14), which are equal to half of the total recorded wavefield, correspond to the conventional interpretation of multicomponent data assuming that vertical component measurements correspond to P-waves, whereas the horizontal components correspond to S-waves. This only holds for vertically incident waves. In contrast, eqs (12)–(14) better approximate the upgoing wavefield for non-normal incidence angles. We test the performance of our up/down separation filters on synthetic data generated using a 2-D finite-difference (FD) code (Virieux 1986). A force source emitting a 50 Hz Ricker wavelet was placed at 100 m depth in a homogeneous medium with P- and S-wave velocities of 1800 and 600 m s−1, respectively. The incident P- and S-wave arrivals were recorded at the free-surface. The free-surface condition was implemented by imaging the stresses and particle velocities on several grid points above the free-surface (Levander 1988; Crase 1990), which we found to be in good agreement with an analytical solution (Sherwood 1958). A reference solution of the upgoing wavefield at the free-surface was derived using the FD-injection method presented by Robertsson et al. (2015) that allows decomposing FD-computed seismograms into their upgoing and downgoing constituents to within machine precision. Fig. 1 displays the recorded horizontal and vertical particle velocity components and the up/down separated results using eqs (12) and (14), respectively. To assess the performance of the spatial wavefield gradient filters we calculated the relative root-mean-square error (RMSE) against incidence angle between the estimated upgoing wavefield and the FD-injection reference solution for the isolated P- or S-wave arrivals at the free-surface (Fig. 2). The dashed lines correspond to taking half of the total recorded wavefield (denoted here as the conventional method), whereas the solid lines correspond to the error from using our spatial gradient-based eqs (12) and (14). Fig. 2(a) shows that our spatial gradient filters significantly improve the approximation to the horizontal component of the upgoing P-wave compared to the conventional approach, with errors below 10 per cent up to incidence angles of about 30°. In contrast, simply taking half of the total wavefield to approximate the vertical component of the upgoing P-wavefield already results in errors below 10 per cent up to 50° incidence (Fig. 2b) and our spatial gradient filters do not significantly reduce the error. In case of an incident P-wave recorded on the vertical component, the amplitudes of the upgoing P-wave and its downgoing reflections have similar magnitudes, while the contributions of P-to-S conversions to the vertical component recordings are very small and predominantly project onto the horizontal component recording (Aki & Richards 2002). This observation explains why taking half of the total vertical component recording is already a good approximation to the vertical component of the upgoing P-wave, whereas our spatial wavefield gradient-based filters are required to better estimate the horizontal component of the upgoing P-wave. Figure 1. View largeDownload slide Up/down separation of a synthetic land seismic data example. (a) Total horizontal component recording. (b) Total vertical component recording. (c) Estimated upgoing horizontal component. (d) Estimated upgoing vertical component. Figure 1. View largeDownload slide Up/down separation of a synthetic land seismic data example. (a) Total horizontal component recording. (b) Total vertical component recording. (c) Estimated upgoing horizontal component. (d) Estimated upgoing vertical component. Figure 2. View largeDownload slide Relative RMSE between the FD-injection reference solution and the estimated particle velocity components of isolated upgoing P- and S-waves plotted against incidence angle for the synthetic land seismic data example in Fig. 1. The dashed line corresponds to the error from using the conventional method (taking half of the total wavefield), while the solid line represents the error from using spatial wavefield gradients (eqs 12 and 14). (a) Error for the horizontal component of the upgoing P-wave. (b) Error for the vertical component of the upgoing P-wave. (c) Error for the horizontal component of the upgoing S-wave. (d) Error for the vertical component of the upgoing S-wave. Figure 2. View largeDownload slide Relative RMSE between the FD-injection reference solution and the estimated particle velocity components of isolated upgoing P- and S-waves plotted against incidence angle for the synthetic land seismic data example in Fig. 1. The dashed line corresponds to the error from using the conventional method (taking half of the total wavefield), while the solid line represents the error from using spatial wavefield gradients (eqs 12 and 14). (a) Error for the horizontal component of the upgoing P-wave. (b) Error for the vertical component of the upgoing P-wave. (c) Error for the horizontal component of the upgoing S-wave. (d) Error for the vertical component of the upgoing S-wave. In case of an incident S-wave, the spatial wavefield gradient filters do not provide a better approximation to the upgoing horizontal component compared to taking half of the total wavefield (Fig. 2c), whereas our spatial gradient filters enable to extract more accurate amplitude information on the vertical component of the upgoing S-wave in comparison to the conventional solution (Fig. 2d) with errors below 10 per cent up to about 20° incidence angle. The reason why the incidence angle range over which the upgoing S-wave can be reliably estimated is much smaller than for a P-wave is due to the free-surface effect. In particular, S-wave arrivals at the free-surface generate horizontally propagating S-to-P mode conversions for incidence angles beyond the critical angle and these waveforms cannot be removed by our approximated filter expressions. The critical angle is a function of the near-surface P- and S-wave velocities and is about 20° in our example. However, we anticipate that in case of field data, due to strong vertical velocity gradients in the shallow subsurface, body waves will mostly arrive at steep incidence angles close to vertical and well below the critical angle for S-waves. 2.1.2 P/S separation The up/down separation of land seismic data removes the effect of the free-surface, but for imaging purposes an additional P/S separation step is desirable. Robertsson & Curtis (2002) derived expressions to estimate the upgoing P- and S-wave potentials using spatially compact receiver groups (see Appendix A). In contrast to Robertsson & Curtis (2002), we prefer to work with particle velocities instead of wavefield potentials to obtain the horizontal component of the upgoing S-wave constituents and the vertical component of the upgoing P-wave constituents. Similar to the up/down separation technique, we seek expressions that remove undesired wavefield constituents from the total recorded particle velocity components. This can be achieved by normalizing eqs (A14)–(A16) as follows:   \begin{equation} v^{S^U}_x \approx \frac{\beta }{i\omega } (\nabla \times v^U)_y \approx \frac{1}{2} \left(v_x - \frac{1}{i \omega } 2 \beta \frac{\partial v_z}{\partial x}\right), \end{equation} (15)  \begin{equation} v^{S^U}_y \approx \frac{\beta }{i\omega } (\nabla \times v^U)_x \approx \frac{1}{2} \left(v_y - \frac{1}{i \omega } 2 \beta \frac{\partial v_z}{\partial y}\right), \end{equation} (16)  \begin{equation} v^{P^U}_z \approx \frac{\alpha }{i\omega } (\nabla \cdot v^U) \approx \frac{1}{2} \left(v_z + \frac{1}{i\omega } \frac{2\beta ^2}{\alpha } \left(\frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} \right)\right). \end{equation} (17) Eqs (15)–(17) separate the recorded particle velocity components into their P- and S-wave constituents and additionally remove the free-surface effect. For example, the spatial wavefield gradient term ∂vx/∂x + ∂vy/∂y in eq. (17) is proportional to the divergence at the free-surface (Robertsson & Curtis 2002). The divergence measured at the free-surface is sensitive to S-to-P conversions and properly scaled by the medium parameters it effectively removes the total vertical component S-wave constituent. In case of an incident P-wave, this divergence term does not only contain the upgoing P-waves but also the downgoing reflected P-waves and is close to zero (destructive interference) when the amplitudes of both constituents have similar magnitudes. Because the vertical components of the upgoing and downgoing P-waves are about equal in amplitude for a considerable range of incidence angles (Aki & Richards 2002), the vertical component of the upgoing P-wave is retained by scaling the total vertical component P-wave arrival by a factor 1/2 (see Fig. 2b). To demonstrate the P/S separation technique on synthetic land seismic data we use the same 2-D model as for the up/down separation. Fig. 3 displays the total horizontal and vertical particle velocity recordings and the P/S separated results using eqs (15) and (17). Fig. 4 shows the relative RMSE versus incidence angle between the P/S separated results using spatial wavefield gradients (eqs 15 and 17) and the FD-injection reference solution for an upgoing P-wave arrival recorded on the vertical component and an upgoing S-wave arrival recorded on the horizontal component (solid lines). For comparison, we also plotted the error from using the conventional method (dashed lines), which assumes P-waves to be recorded on the vertical component and S-waves to be recorded on the horizontal component only. Our spatial wavefield gradient method significantly improves the P/S separated results for the upgoing horizontal and vertical components, reducing the errors below 10 per cent up to about 25° and 20° incidence, respectively. Although, the P-wave arrival recorded on the horizontal component is fully removed up to 50° incidence, the main error source arises from extracting the horizontal component of the upgoing S-wave beyond the critical angle. Similarly, the vertical component of the incident S-wave is only suppressed up to the critical angle. Nevertheless, the improved P/S separated results demonstrate that our spatial wavefield gradient filters result in good quality recordings for P/S imaging and are particularly attractive to enhance S-wave analyses on horizontal component recordings. Figure 3. View largeDownload slide P/S separation of a synthetic land seismic data example. (a) Total horizontal component recording. (b) Total vertical component recording. (c) Estimated upgoing horizontal component with the P-wave arrival suppressed. (d) Estimated upgoing vertical component with the S-wave arrival suppressed. Figure 3. View largeDownload slide P/S separation of a synthetic land seismic data example. (a) Total horizontal component recording. (b) Total vertical component recording. (c) Estimated upgoing horizontal component with the P-wave arrival suppressed. (d) Estimated upgoing vertical component with the S-wave arrival suppressed. Figure 4. View largeDownload slide Relative RMSE between the FD-injection reference solution and the estimated P/S separation results for the synthetic land seismic data example in Fig. 3. The dashed line corresponds to the error from using the conventional method (P- and S-waves only recorded on the vertical and horizontal components, respectively), while the solid line represents the error from using spatial wavefield gradients (eqs 15 and 17). (a) Error for the P-wave suppression on the upgoing horizontal component. (b) Error for the S-wave suppression on the upgoing vertical component. Figure 4. View largeDownload slide Relative RMSE between the FD-injection reference solution and the estimated P/S separation results for the synthetic land seismic data example in Fig. 3. The dashed line corresponds to the error from using the conventional method (P- and S-waves only recorded on the vertical and horizontal components, respectively), while the solid line represents the error from using spatial wavefield gradients (eqs 15 and 17). (a) Error for the P-wave suppression on the upgoing horizontal component. (b) Error for the S-wave suppression on the upgoing vertical component. 2.2 Wavefield separation of OBS data using spatial wavefield gradients 2.2.1 Up/down separation Up/down separation based on the elastodynamic representation theorem and using spatial wavefield gradients can be extended to OBS data for receiver-side water layer multiple removal. Amundsen et al. (2000) derived analytical expressions for the upgoing horizontal and vertical wavefields that are exact for homogeneous and isotropic media at an infinitesimal distance below the seabed:   \begin{equation} v_x^{U} = \frac{1}{2}\left(v_x + F^{v_x}_{v_z} \ast v_z\right), \end{equation} (18)  \begin{equation} v_y^{U} = \frac{1}{2}\left(v_y + F^{v_y}_{v_z} \ast v_z\right), \end{equation} (19)  \begin{equation} v_z^{U} = \frac{1}{2}\left(v_z - F^{v_z}_{v_x} \ast v_x - F^{v_z}_{v_y} \ast v_y - F^{v_z}_{p} \ast p\right). \end{equation} (20) The filters operating on the particle velocity components are identical to those derived for the land seismic case (eqs 4–7). Eq. (20) contains an additional spatial filter term $$F^{v_z}_{p}$$ that operates on the pressure recordings p. In the wavenumber–frequency domain this filter is given by (Amundsen et al. 2000):   \begin{equation} \widetilde{F}^{v_z}_{p} = \frac{1}{\rho \omega k^{(\beta )}_z}\left(k^2_x + k^2_y + k^{(\alpha )}_z k^{(\beta )}_z \right), \end{equation} (21) where ρ is the density of the seafloor. In order to adapt the up/down separation technique of OBS data to local measurements of the wavefield only, we follow the same approach as for the land seismic case by expanding the filters up to first-order spatial gradients. For $$\widetilde{F}^{v_z}_{p}$$ this results in   \begin{equation} \widetilde{F}^{v_z}_{p} \approx \frac{k_\alpha }{\rho \omega }. \end{equation} (22) Inserting the truncated filter expressions (8)–(11) and (22) into eqs (18)–(20) and transforming back to the space–frequency domain results in the following expressions for the upgoing particle velocity components at an infinitesimal distance below the seabed:   \begin{equation} v_x^U \approx \frac{1}{2}\left(v_x + \frac{1}{i\omega }(\alpha -2\beta )\frac{\partial v_z}{\partial x}\right), \end{equation} (23)  \begin{equation} v_y^U \approx \frac{1}{2}\left(v_y + \frac{1}{i\omega }(\alpha -2\beta )\frac{\partial v_z}{\partial y}\right), \end{equation} (24)  \begin{equation} v_z^U \approx \frac{1}{2}\left(v_z -\frac{1}{i\omega } \left(\beta -\frac{2\beta ^2}{\alpha }\right)\left(\frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y}\right) - \frac{1}{\alpha \rho } p\right). \end{equation} (25) Eqs (23)–(25) suppress all downgoing events, including the direct wave, receiver ghosts and accompanying water layer reverberations (receiver-side multiples) in addition to remaining downgoing reflections and mode conversions generated at the seabed. We demonstrate the spatial gradient-based multiple attenuation technique on synthetic FD modelled vertical particle velocity recordings at the seabed. Our 2-D model consists of a 100 m thick water layer bounded above by a free-surface. Below the water layer is a homogeneous medium with a P-wave velocity of 1800 m s−1, an S-wave velocity of 600 m s−1 and a density of 1600 m kg−3. An explosive point source emitting a 50 Hz Ricker wavelet was placed at 5 m depth below the sea surface. In order to generate overlapping downgoing and upgoing arrivals recorded by a multicomponent OBS sensor, a second explosive point source was located at 150 m below the seafloor. Figs 5(a) and (b) show the total as well as the estimated upgoing vertical particle velocity recordings, obtained by evaluating eq. (25). The downgoing direct wave (black arrow) is effectively suppressed, while the upgoing arrivals (green arrows) are preserved, even in case of full overlap of the upgoing and downgoing events as observed for example at about 100 m offset. The later arriving water layer reverberations (grey arrows) are also effectively suppressed. Figure 5. View largeDownload slide Receiver-side water layer multiple suppression of a synthetic OBS data example. (a) Total vertical component recording. (b) Estimated upgoing vertical component. The downgoing direct wave is denoted by the black arrow, whereas the upgoing waves are indicated by the green arrows. Downgoing water layer reverberations are highlighted by the grey arrows. Figure 5. View largeDownload slide Receiver-side water layer multiple suppression of a synthetic OBS data example. (a) Total vertical component recording. (b) Estimated upgoing vertical component. The downgoing direct wave is denoted by the black arrow, whereas the upgoing waves are indicated by the green arrows. Downgoing water layer reverberations are highlighted by the grey arrows. The performance of eq. (25) to suppress the downgoing direct wave is shown in Fig. 6, which displays the relative RMSE between the filtered data and the FD-injection reference solution for the upgoing wavefield (equal to zero) (solid line). The dashed line represents the error from using the conventional multiple suppression approach, which corresponds to combining the particle velocity and pressure recordings only (PZ summation, Amundsen 1993) and ignoring the spatial gradient term in eq. (25). Both approaches effectively suppress vertical component water layer multiples with errors below 10 per cent up to incidence angles of about 40°, although our spatial wavefield gradient-based technique does not significantly improve the results compared to the conventional PZ summation. Figure 6. View largeDownload slide Relative RMSE between the FD-injection reference solution and the estimated upgoing vertical component plotted against incidence angle for the synthetic OBS data example in Fig. 5. The dashed line corresponds to the error from using the conventional method (PZ summation), while the solid line represents the error from using spatial wavefield gradients (eq. 25). Figure 6. View largeDownload slide Relative RMSE between the FD-injection reference solution and the estimated upgoing vertical component plotted against incidence angle for the synthetic OBS data example in Fig. 5. The dashed line corresponds to the error from using the conventional method (PZ summation), while the solid line represents the error from using spatial wavefield gradients (eq. 25). 2.2.2 P/S separation P/S separation based on the elastodynamic representation theorem using spatial wavefield gradients can be extended to OBS data and the corresponding derivations are given in Appendix B. Similar to the land seismic case, we prefer to interpret the horizontal particle velocity component of the upgoing S-wave constituents and the vertical particle velocity component of the upgoing P-wave constituents, which can be achieved by normalizing the P- and S-wave potential equations (B22)–(B24) as follows:   \begin{equation} v_x^{S^U} \approx \frac{1}{\rho \beta } (\nabla \times v^U)_y \approx \frac{1}{2} \left(v_x - \frac{1}{i\omega }2 \beta \frac{\partial v_z}{\partial x} + \frac{1}{i \omega } \frac{1}{\rho } \frac{\partial p}{\partial x}\right), \end{equation} (26)  \begin{equation} v_y^{S^U} \approx \frac{1}{\rho \beta } (\nabla \times v^U)_x \approx \frac{1}{2} \left(v_y - \frac{1}{i \omega } 2 \beta \frac{\partial v_z}{\partial y} + \frac{1}{i \omega } \frac{1}{\rho } \frac{\partial p}{\partial y}\right), \end{equation} (27)  \begin{eqnarray} v_z^{P^U} \approx \frac{1}{\alpha \rho } (\nabla \cdot v) \approx \frac{1}{2} \left(v_z + \frac{1}{i\omega } \frac{2 \beta ^2}{\alpha } \left(\frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y}\right) - \frac{1}{\alpha \rho }p \right).\nonumber\\ \end{eqnarray} (28) In addition to extracting the upgoing P- or S-waves, eqs (26)–(28) also suppress all downgoing events such as the receiver-side water layer multiples. To illustrate the P/S separation method, OBS data were modelled using the same 2-D subsurface model as for the up/down separation but now using a force source emitting a 50 Hz Ricker wavelet at 100 m below the seafloor to generate both P- and S-waves. Fig. 7 shows the total horizontal and vertical particle velocity recordings at the seafloor and the P/S separated results using eqs (26) and (28). Our technique effectively suppresses the horizontal component of the P-wave and the vertical component of the S-wave up to moderate incidence angles. The sea surface multiple (receiver ghost), which is observed on the vertical component at about 0.21 s, has also been removed by the combined up/down and P/S separation. Figure 7. View largeDownload slide P/S separation of a synthetic OBS data example. (a) Total horizontal component recording. (b) Total vertical component recording. (c) Estimated upgoing horizontal component with the P-wave arrival suppressed. (d) Estimated upgoing vertical component with the S-wave arrival suppressed. Figure 7. View largeDownload slide P/S separation of a synthetic OBS data example. (a) Total horizontal component recording. (b) Total vertical component recording. (c) Estimated upgoing horizontal component with the P-wave arrival suppressed. (d) Estimated upgoing vertical component with the S-wave arrival suppressed. The performance of the seabed P/S separation filters (solid line, eqs 26 and 28) as well as the conventional method (dashed line) against FD-injection as a function of incidence angle is displayed in Fig. 8 and shows very similar behaviour as for the land seismic case (Fig. 4). Figure 8. View largeDownload slide Relative RMSE between the FD-injection reference solution and the estimated P/S separation results for the synthetic OBS data example in Fig. 7. The dashed line corresponds to the error from using the conventional method (P- and S-waves only recorded on the vertical and horizontal components, respectively), while the solid line represents the error from using spatial wavefield gradients (eqs 26 and 28). (a) Error for the P-wave suppression on the upgoing horizontal component. (b) Error for the S-wave suppression on the upgoing vertical component. Figure 8. View largeDownload slide Relative RMSE between the FD-injection reference solution and the estimated P/S separation results for the synthetic OBS data example in Fig. 7. The dashed line corresponds to the error from using the conventional method (P- and S-waves only recorded on the vertical and horizontal components, respectively), while the solid line represents the error from using spatial wavefield gradients (eqs 26 and 28). (a) Error for the P-wave suppression on the upgoing horizontal component. (b) Error for the S-wave suppression on the upgoing vertical component. 3 VELOCITIES FOR WAVEFIELD SEPARATION The application of wavefield separation techniques using spatial filters requires the knowledge of the ‘near-surface’ P- and S-wave velocities (e.g. Dankbaar 1985; Robertsson & Curtis 2002). It has been unclear whether these ‘near-surface’ velocities are the velocities at or below the free-surface, possibly averaged over some volume, and whether they are frequency dependent or not. Robertsson & Curtis (2002) speculated that the required P- and S-wave velocities likely are those experienced by the wavefield, whereby longer wavelengths are influenced by the elastic properties of a larger volume. In order to experimentally investigate the velocity dependence of our wavefield separation filters, we apply our techniques to a synthetic data example with a strong vertical velocity gradient in the shallow subsurface. The P- and S-wave velocities in our synthetic model linearly increase from 1800 to 4800 m s−1 and from 600 to 1800 m s−1, respectively, between the receivers at the free-surface and a depth of 100 m. An explosive source emitting a 50 Hz Ricker wavelet was excited 100 m below the surface. We used an iterative frequency-domain Gauss–Newton inversion scheme (Pratt et al. 1998) to determine the velocities that minimize the RMSE between the estimated upgoing horizontal and vertical wavefield obtained using eqs (12)–(14) and the reference upgoing wavefield obtained by means of FD-injection. The inversion approach consistently yields the P- and S-wave velocities exactly at the free-surface for frequencies above 30 Hz (Fig. 9). Furthermore, the accuracy of the wavefield separation results in case of a strong vertical velocity gradient is the same as in case of a homogeneous subsurface (see Figs 2a and b). Below 30 Hz, the optimum velocity seems to be frequency-dependent with the inverted P-wave velocities being lower than the minimum P-wave velocity in our subsurface model, whereas the inverted S-wave velocities are higher. The frequency dependence can be attributed to near-field effects: for the lowest frequencies in our example, the plane wave assumption inherent to the up/down separation algorithm is not met. In the far field, the velocities required to separate seismic data into their desired wavefield constituents using spatial wavefield gradients are exactly those at the receiver locations. Our synthetic data example demonstrates that highly accurate results can be achieved in the presence of a vertical velocity gradient. Bearing in mind that our wavefield decomposition filters were derived under a homogeneous space Green’s function assumption, this is perhaps not an entirely intuitive result. Figure 9. View largeDownload slide Inverted P- and S-wave velocities as a function of frequency in order to minimize the error of the up/down wavefield separation technique compared to the FD-injection reference solution in case of a strong vertical velocity gradient in the velocity model. Figure 9. View largeDownload slide Inverted P- and S-wave velocities as a function of frequency in order to minimize the error of the up/down wavefield separation technique compared to the FD-injection reference solution in case of a strong vertical velocity gradient in the velocity model. Our results are consistent with observations recently made by Robertsson et al. (2015). Robertsson et al. (2015) found that in order to isolate the upgoing/downgoing constituents of seismic recordings using the FD-injection method, only knowledge of the material properties exactly at the recording location is needed. This observation has so far not been fully appreciated, but has significant implications for field data processing. The shallow subsurface is typically characterized by strong velocity variations and the fact that only the knowledge of the P- and S-wave velocities at the recording locations are required significantly simplifies wavefield separation in practice. 4 FIELD DATA EXAMPLE We applied our spatial wavefield gradient-based technique to a multicomponent land seismic data set acquired in northern Switzerland to characterize the geometry and infill of a Quaternary valley (Schmelzbach et al. 2014; Sollberger et al. 2014). Near-surface studies benefit from the enhanced resolution and improved lithological characterization achieved by S-wave and converted-wave exploration (Helbig & Mesdag 1982). We therefore test our P/S separation filters to suppress P-waves on the horizontal particle velocity component recordings for S-wave imaging. The survey consisted of a 250 m long 2-D line of 3C geophones spaced at 1.5 m intervals. Seismic data were generated by striking a prismatic wedge (3C vector source, Schmelzbach et al. 2016a) with a sledgehammer at 3 m intervals along the recording line. Fig. 10(a) displays a horizontal component recording of a vertically directed source. Based on hodogram analyses and arrival identifications previously carried out by Sollberger et al. (2016b), the P-wave first arrival is marked by the red arrow, whereas the S-wave first arrival is highlighted by the blue arrow. At later times, the shot gather is mainly dominated by Rayleigh waves (green arrow). Figure 10. View largeDownload slide P/S separation of field data. The P- and S-wave first arrivals are indicated by the red and blue arrows, respectively, whereas a Rayleigh wave is indicated by the green arrow. (a) Total horizontal component recording. (b) Estimated upgoing horizontal component with the P-wave arrival suppressed. Figure 10. View largeDownload slide P/S separation of field data. The P- and S-wave first arrivals are indicated by the red and blue arrows, respectively, whereas a Rayleigh wave is indicated by the green arrow. (a) Total horizontal component recording. (b) Estimated upgoing horizontal component with the P-wave arrival suppressed. The P/S separation of the recorded horizontal particle velocity was carried out using eq. (15). The horizontal spatial wavefield gradients of the vertical particle velocities were calculated by finite-differencing the recordings from three neighbouring stations and centring the estimate at the middle receiver. The S-wave velocity was obtained from the traveltime tomography results in Schmelzbach et al. (2014) and is about 250 m s−1 at the surface. Fig. 10(b) shows that our spatial wavefield gradient-based P/S separation filters effectively suppress the undesired P-wave arrival to obtain a horizontal component recording that only contains S-waves. To verify the wavefield separation results we carried out the P/S separation on synthetic data that were generated from the traveltime tomography results in Schmelzbach et al. (2014). Fig. 11 displays the total and P/S separated horizontal component synthetic data and shows an overall very good agreement to both the raw and separated field data (Fig. 10). Hodogram analyses and snapshots of the FD-simulated wavefield confirm that the isolated body wave arrival in Fig. 11(b) is indeed an S-wave arrival. Figure 11. View largeDownload slide P/S separation of synthetic FD data generated from the traveltime tomography results in Schmelzbach et al. (2014). The P- and S-wave first arrivals are indicated by the red and blue arrows, respectively, whereas a Rayleigh wave is indicated by the green arrow. (a) Total horizontal component recording. (b) Estimated upgoing horizontal component with the P-wave arrival suppressed. Figure 11. View largeDownload slide P/S separation of synthetic FD data generated from the traveltime tomography results in Schmelzbach et al. (2014). The P- and S-wave first arrivals are indicated by the red and blue arrows, respectively, whereas a Rayleigh wave is indicated by the green arrow. (a) Total horizontal component recording. (b) Estimated upgoing horizontal component with the P-wave arrival suppressed. 5 DISCUSSION The techniques presented in this paper were tested on synthetic data examples with a typical near-surface P-to-S-wave velocity ratio equal to 3. It is important to stress that additional synthetic data tests showed that our observations generally hold in case of realistic near-surface Poisson’s ratios (between 0.2 and 0.45), as also demonstrated by the field data example. The spatial gradient data in our synthetic and field data examples were estimated by finite-differencing the particle velocity and pressure recordings from neighbouring stations, which in practice require spatially compact receiver groups. Recently, reliable and portable rotational sensors have been made commercially available for seismic purposes (Pierson et al. 2016; Bernauer et al. 2017). Such sensors allow to directly measure the spatial wavefield gradient quantities ∂vz/∂x and ∂vz/∂y, which correspond to the rotational motion around the two horizontal axes at the free-surface. These spatial gradients of the vertical particle velocity component are required for the up/down and P/S separation of land and seabed horizontal particle velocity component data. Collocated measurements of particle velocities and spatial wavefield gradients (e.g. 6C sensors; Sollberger et al. 2017) allow for single station wavefield separation using the expressions derived in this study, which is very attractive for field applications. 6 CONCLUSIONS We present a wavefield decomposition technique combining local measurements of particle velocity and spatial wavefield gradient data to separate multicomponent land and seabed seismic data into their upgoing/downgoing and P/S constituents. Generally, our wavefield separation filters involving first-order derivatives are more effective in suppressing P-waves than S-waves and, hence, are of particular interest for S-wave and P-to-S-converted wave analysis for steep to moderate incidence angles. The filter performance as a function of incidence angle could be improved by including higher-order derivatives, which would require larger arrays and the gradient estimation may be more prone to errors compared to estimating first-order derivatives. Spatial gradients may be derived from small sensor arrays or, given that the horizontal gradient of the vertical particle velocity corresponds at the free-surface to a rotational motion, some spatial gradients can directly be measured with rotational sensors. Therefore, our filters are particularly interesting for the application to single station multicomponent measurements with dedicated sensors at sites where one or only a few stations can be installed such as at the sea bottom or in extraterrestrial seismology. Synthetic data tests revealed that the filter performance only shows a weak dependence on the subsurface elastic parameter variations. Hence, the derived wavefield separation approach has a wide potential application range in terms of environment. An interesting observation is that the filters require as input parameters the elastic parameters found exactly at the receivers. This observation may have implications for other wave equation-based techniques that either depend on or are used to infer on material properties. As only the receiver elastic properties need to be known, the estimation of these material parameters using standard geophysical techniques (e.g. traveltime tomography, surface wave tomography) is straightforward and facilitates the application of our wavefield separation approaches to field data. In principle, the presented wavefield separation technique can easily be adopted to estimate the material properties from identified isolated arrivals. ACKNOWLEDGEMENTS This work was supported by Swiss National Science Foundation grant 200021_156996 and the sponsors of the CARNEVAL consortium (Nagra, OMV and Schlumberger Gould Research). We wish to thank editor Herve Chauris and two anonymous reviewers for constructive reviews that improved the paper. REFERENCES Aki K., Richards P.G., 2002. Quantitative Seimology , University Science Books. Amundsen L., 1993. Wavenumber-based filtering of marine point-source data, Geophysics , 58( 9), 1335– 1348. https://doi.org/10.1190/1.1443516 Google Scholar CrossRef Search ADS   Amundsen L., Ikelle L., Martin J., 2000. Multiple attenuation and P/S splitting of multicomponent OBC data at a heterogeneous sea floor, Wave Motion , 32, 67– 78. https://doi.org/10.1016/S0165-2125(99)00047-5 Google Scholar CrossRef Search ADS   Amundsen L., Ikelle L., Robertsson J.O.A., Westerdahl H., 2013. Wavefield decomposition of seabed node marine recordings, in 83rd Annual International Meeting, SEG, Expanded Abstracts , pp. 146– 150. Barak O., Brune R., Herkenhoff F., Dash R., Rector J., Ronen S., 2013. Seven-component seismic data, in 83rd Annual International Meeting, SEG, Expanded Abstracts , pp. 5151– 5155. Barak O., Key K., Constable S., Milligan P., Ronen S., 2015. Acquiring rotation data on the ocean bottom without rotation sensors, in 85th Annual International Meeting, SEG, Expanded Abstracts , pp. 2148– 2152. Bernauer F. et al.  , 2017. BlueSeis3A - full characterization of a 3C broadband rotational ground motion sensor for seismology, in 19th EGU General Assembly Conference Abstracts , p. 15512. Crase E., 1990. High-order (space and time) finite-difference modeling of the elastic wave equation, in 60th Annual International Meeting, SEG, Expanded Abstracts , pp. 987– 991. Dankbaar J., 1985. Separation of P- and S-waves, Geophys. Prospect. , 33, 970– 986. https://doi.org/10.1111/j.1365-2478.1985.tb00792.x Google Scholar CrossRef Search ADS   Helbig K., Mesdag C.S., 1982. The potential of shear wave observations, Geophys. Prospect. , 30, 413– 431. https://doi.org/10.1111/j.1365-2478.1982.tb01314.x Google Scholar CrossRef Search ADS   Igel H., Brokesova J., Evans J.R., Zembaty Z., 2012. Preface to the special issue on “Advances in rotational seismology: instrumentation, theory, observations, and engineering”, J. Seismol. , 16( 4), 571– 572. https://doi.org/10.1007/s10950-012-9307-6 Google Scholar CrossRef Search ADS   Langston C., 2007a. Spatial Gradient Analysis for Linear Seismic Arrays, Bull. seism. Soc. Am. , 97( 1), 265– 280. https://doi.org/10.1785/0120060100 Google Scholar CrossRef Search ADS   Langston C., 2007b. Wave gradiometry in two dimensions, Bull. seism. Soc. Am. , 97( 2), 401– 416. https://doi.org/10.1785/0120060138 Google Scholar CrossRef Search ADS   Levander A., 1988. Fourth-order finite-difference P–SV seismograms, Geophysics , 53, 1425– 1436. https://doi.org/10.1190/1.1442422 Google Scholar CrossRef Search ADS   Liu Y., Holt W.E., 2015. Wave gradiometry and its link with Helmholtz equation solutions applied to USArray in the eastern U.S., J. geophys. Res. , 120( 8), 5717– 5746. https://doi.org/10.1002/2015JB011982 Google Scholar CrossRef Search ADS   Osen A., Amundsen L., Reitan A., 2002. Toward optimal spatial filters for demultiple and wavefield splitting of ocean-bottom seismic data, Geophysics , 67, 1983– 1990. https://doi.org/10.1190/1.1527098 Google Scholar CrossRef Search ADS   Pierson R., Laughlin D., Brune R., 2016. Advances in rotational seismic measurements, in 86th Annual International Meeting, SEG, Expanded Abstracts , pp. 2263– 2267. Porter R., Liu Y., Holt W.E., 2016. Lithospheric records of orogeny within the continental U.S., Geophys. Res. Lett. , 43( 1), 144– 153. https://doi.org/10.1002/2015GL066950 Google Scholar CrossRef Search ADS   Pratt R., Shin C., Hicks G., 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion, Geophys. J. Int. , 133, 341– 362. https://doi.org/10.1046/j.1365-246X.1998.00498.x Google Scholar CrossRef Search ADS   Robertsson J., Curtis A., 2002. Wavefield separation using densely deployed three-component single-sensor groups in land surface-seismic recordings, Geophysics , 67( 5), 1624– 1633. https://doi.org/10.1190/1.1512809 Google Scholar CrossRef Search ADS   Robertsson J., Muyzert E., 1999. Wavefield separation using a volume distribution of three component recordings, Geophys. Res. Lett. , 26( 18), 2821– 2824. https://doi.org/10.1029/1999GL010472 Google Scholar CrossRef Search ADS   Robertsson J., van Manen D.-J., Schmelzbach C., Van Renterghem C., Amundsen L., 2015. Finite-difference modelling of wavefield constituents, Geophys. J. Int. , 203( 2), 1334– 1342. https://doi.org/10.1093/gji/ggv379 Google Scholar CrossRef Search ADS   Schmelzbach C. et al.  , 2014. Multi-method geophysical imaging of a Quaternary valley in northern Switzerland, in 84th Annual International Meeting, SEG, Expanded Abstracts , pp. 2083– 2087. Schmelzbach C., Sollberger D., Greenhalgh S.A., Horstmeyer H., Maurer H., Robertsson J.O.A., 2016a. 9C seismic data acquisition for near-surface applications: recording, waveform reciprocity and 4C rotation, in 78th Conference and Exhibition , EAGE, Workshop Programme, WS04 B03. Schmelzbach C., Sollberger D., Van Renterghem C., Häusler M., Robertsson J., Greenhalgh S., 2016b. Seismic spatial wavefield gradient and rotational rate measurements as new observables in land seismic exploration, in 18th EGU General Assembly Conference Abstracts , p. 3787. Sherwood J., 1958. Elastic wave propagation in a semi-infinite solid medium, Proc. Phys. Soc. , 71, 207– 219. https://doi.org/10.1088/0370-1328/71/2/308 Google Scholar CrossRef Search ADS   Sollberger D., Schmelzbach C., Horstmeyer H., Reiser F., Rabenstein L., Maurer H., Robertsson J.O.A., Greenhalgh S., 2014. Mapping near surface-elastic parameters of Quaternary sediments using multicomponent seismic techniques, in 16th EGU General Assembly Conference Abstracts , p. 4014. Sollberger D., Schmelzbach C., Robertsson J., Greenhalgh S., Nakamura Y., Khan A., 2016a. The shallow elastic structure of the lunar crusts: new insights from seismic wavefield gradient analysis, Geophys. Res. Lett. , 43, 1– 10. https://doi.org/10.1002/2016GL070883 Google Scholar CrossRef Search ADS   Sollberger D., Schmelzbach C., Van Renterghem C., Robertsson J., Greenhalgh S., 2016b. Single-component elastic wavefield separation at the free surface using source- and receiver-side gradients, in 86th Annual International Meeting, SEG, Expanded Abstracts , pp. 2268– 2273. Sollberger D., Schmelzbach C., Van Renterghem C., Robertsson J., Greenhalgh S., 2017. Automated, six-component, single-station ground-roll identification and suppression by combined processing of translational and rotational ground motion, in 87th Annual International Meeting, SEG, Expanded Abstracts , pp. 5064– 5068. Spudich P., Steck L.K., Hellweg M., Fletcher J.B., Baker L.M., 1995. Transient stresses at Parkfield, California, produced by the M 7.4 Landers earthquake of June 28, 1992: observations from the UPSAR dense seismograph array, J. geophys. Res. , 100, 675– 690. https://doi.org/10.1029/94JB02477 Google Scholar CrossRef Search ADS   Virieux J., 1986. P–SV wave propagation in heterogeneous media: velocity-stress finite-difference method, Geophysics , 51( 4), 889– 901. https://doi.org/10.1190/1.1442147 Google Scholar CrossRef Search ADS   Wapenaar C., Haimé G., 1990. Elastic extrapolation of primary seismic P- and S-waves, Geophys. Prospect. , 38, 23– 60. https://doi.org/10.1111/j.1365-2478.1990.tb01833.x Google Scholar CrossRef Search ADS   Wapenaar C., Hermann P., Verschuur D., Berkhout A., 1990. Decomposition of multicomponent seismic data into primary P- and S-wave responses, Geophys. Prospect. , 38, 633– 662. https://doi.org/10.1111/j.1365-2478.1990.tb01867.x Google Scholar CrossRef Search ADS   APPENDIX A: P/S SEPARATION OF LAND SEISMIC DATA Robertsson & Curtis (2002) derived the following expressions to extract the upgoing P- and S-wave potentials by taking the divergence and curl of up/down separated land seismic data:   \begin{equation} (\nabla \cdot v^U) = \frac{k_\alpha ^2}{k_\beta ^2} \left(\frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y}\right) + F_{v_z}^{(\nabla \cdot v) } \ast v_z, \end{equation} (A1)  \begin{equation} (\nabla \times v^U)_x = -\frac{\partial v_z}{\partial y} + F_{v_x}^{(\nabla \times v)_x } \ast v_x + F_{v_y}^{(\nabla \times v)_x} \ast v_y, \end{equation} (A2)  \begin{equation} (\nabla \times v^U)_y = - \frac{\partial v_z}{\partial x} + F_{v_x}^{(\nabla \times v)_y } \ast v_x + F_{v_y}^{(\nabla \times v)_y} \ast v_y, \end{equation} (A3)where in the wavenumber–frequency domain the filters can be expressed as (Robertsson & Curtis 2002):   \begin{equation} \widetilde{F}_{v_z}^{(\nabla \cdot v) } = ik_\alpha ^2 \frac{\left(1-2k_\beta ^{-2} \left(k_x^2 + k_y^2\right)\right)}{2k_z^{(\alpha )}}, \end{equation} (A4)  \begin{equation} \widetilde{F}_{v_x}^{(\nabla \times v)_x } = i \frac{k_x k_y}{2k_z^{(\beta )}}, \end{equation} (A5)  \begin{equation} \widetilde{F}_{v_y}^{(\nabla \times v)_x } = i \frac{k_\beta ^2 - 2k_y^2 - k_x^2}{2k_z^{(\beta )}}, \end{equation} (A6)  \begin{equation} \widetilde{F}_{v_x}^{(\nabla \times v)_y } = i \frac{k_\beta ^2 - 2k_x^2 - k_y^2}{2k_z^{(\beta )}}, \end{equation} (A7)  \begin{equation} \widetilde{F}_{v_y}^{(\nabla \times v)_y } = \widetilde{F}_{v_x}^{(\nabla \times v)_x }. \end{equation} (A8)Similarly as with the up/down separation technique we are seeking expressions that require local measurements of the wavefield only. The Taylor expansion of eqs (A4)–(A8) up to first-order spatial gradients yields (Robertsson & Curtis 2002):   \begin{equation} \widetilde{F}_{v_z}^{(\nabla \cdot v) } \approx i\frac{k_\alpha }{2}, \end{equation} (A9)  \begin{equation} \widetilde{F}_{v_x}^{(\nabla \times v)_x } \approx 0, \end{equation} (A10)  \begin{equation} \widetilde{F}_{v_y}^{(\nabla \times v)_x } \approx i\frac{k_\beta }{2}, \end{equation} (A11)  \begin{equation} \widetilde{F}_{v_x}^{(\nabla \times v)_y } \approx i\frac{k_\beta }{2}, \end{equation} (A12)  \begin{equation} \widetilde{F}_{v_y}^{(\nabla \times v)_y } \approx 0. \end{equation} (A13)Inserting eqs (A9)–(A13) into eqs (A1)–(A3) and transforming back to the space–frequency domain results in the following expressions for the upgoing P- and S-wave potentials (Robertsson & Curtis 2002):   \begin{equation} (\nabla \cdot v^U) \approx \frac{\beta ^2}{\alpha ^2} \frac{\partial v_x}{\partial x} + \frac{i\omega }{2 \alpha } v_z, \end{equation} (A14)  \begin{equation} (\nabla \times v^U)_x \approx - \frac{\partial v_z}{\partial y} + \frac{i\omega }{2 \beta } v_y, \end{equation} (A15)  \begin{equation} (\nabla \times v^U)_y \approx - \frac{\partial v_z}{\partial x} + \frac{i\omega }{2 \beta } v_x. \end{equation} (A16) APPENDIX B: P/S SEPARATION OF OBS DATA The upgoing P- and S-wave potential fields at the seabed are given by (Amundsen et al. 2000):   \begin{equation} (\nabla \cdot v^U) = -\frac{1}{2} p + \frac{i \omega \rho }{k_\beta ^{2}} \left(\frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} \right) + F_{v_z}^{(\nabla \cdot v) } \ast v_z, \end{equation} (B1)  \begin{eqnarray} (\nabla \times v^U)_x &=& F_p^{(\nabla \times v)_x} \ast p + F_{v_x}^{(\nabla \times v)_x } \ast v_x + F_{v_y}^{(\nabla \times v)_x } \ast v_y \nonumber\\ && -\, F_{v_z}^{(\nabla \times v)_x } \ast v_z, \end{eqnarray} (B2)  \begin{eqnarray} (\nabla \times v^U)_y &=& F_p^{(\nabla \times v)_y} \ast p + F_{v_x}^{(\nabla \times v)_y } \ast v_x + F_{v_y}^{(\nabla \times v)_y } \ast v_y \nonumber\\ && -\, F_{v_z}^{(\nabla \times v)_y } \ast v_z, \end{eqnarray} (B3)with filters in the wavenumber–frequency domain (Amundsen et al. 2000):   \begin{equation} \widetilde{F}_{v_z}^{(\nabla \cdot v)} = \frac{\omega \rho }{2k_z^{(\alpha )}} \left(1 - \frac{k_x^2 + k_y^2}{k_\beta ^2}\right), \end{equation} (B4)  \begin{equation} \widetilde{F}_p^{(\nabla \times v)_x} = \frac{k_y}{2k_z^{(\beta )}}, \end{equation} (B5)  \begin{equation} \widetilde{F}_{v_x}^{(\nabla \times v)_x} = \frac{\omega \rho }{2k_z^{(\beta )}} \left(\frac{k_x k_y}{k_\beta ^2}\right), \end{equation} (B6)  \begin{equation} \widetilde{F}_{v_y}^{(\nabla \times v)_x} = \frac{\omega \rho }{2k_z^{(\beta )}} \left(1 - \frac{k_x^2 + 2k_y^2}{k_\beta ^2}\right), \end{equation} (B7)  \begin{equation} \widetilde{F}_{v_z}^{(\nabla \times v)_x} = \omega \rho \frac{k_y}{k_\beta ^2}, \end{equation} (B8)  \begin{equation} \widetilde{F}_p^{(\nabla \times v)_y} = \frac{k_x}{2k_z^{(\beta )}}, \end{equation} (B9)  \begin{equation} \widetilde{F}_{v_x}^{(\nabla \times v)_y} = \frac{\omega \rho }{2k_z^{(\beta )}} \left(1 - \frac{2k_x^2 + k_y^2}{k_\beta ^2}\right), \end{equation} (B10)  \begin{equation} \widetilde{F}_{v_y}^{(\nabla \times v)_y} = \frac{\omega \rho }{2k_z^{(\beta )}} \left(\frac{k_x k_y}{k_\beta ^2}\right), \end{equation} (B11)  \begin{equation} \widetilde{F}_{v_z}^{(\nabla \times v)_y} = \omega \rho \frac{k_x}{k_\beta ^2}. \end{equation} (B12)The Taylor approximations to eqs (B4)–(B12) up to first-order spatial gradients yields:   \begin{equation} \widetilde{F}_{v_z}^{(\nabla \cdot v)} \approx \frac{\omega \rho }{2 k_\alpha }, \end{equation} (B13)  \begin{equation} \widetilde{F}_p^{(\nabla \times v)_x} \approx \frac{k_y}{2 k_\beta }, \end{equation} (B14)  \begin{equation} \widetilde{F}_{v_x}^{(\nabla \times v)_x} \approx 0, \end{equation} (B15)  \begin{equation} \widetilde{F}_{v_y}^{(\nabla \times v)_x} \approx \frac{\omega \rho }{2 k_\beta }, \end{equation} (B16)  \begin{equation} \widetilde{F}_{v_z}^{(\nabla \times v)_x} \approx \omega \rho \frac{k_y}{k_\beta ^2}, \end{equation} (B17)  \begin{equation} \widetilde{F}_p^{(\nabla \times v)_y} \approx \frac{k_x}{2 k_\beta }, \end{equation} (B18)  \begin{equation} \widetilde{F}_{v_x}^{(\nabla \times v)_y} \approx \frac{\omega \rho }{2 k_\beta }, \end{equation} (B19)  \begin{equation} \widetilde{F}_{v_y}^{(\nabla \times v)_y} \approx 0, \end{equation} (B20)  \begin{equation} \widetilde{F}_{v_z}^{(\nabla \times v)_y} \approx \omega \rho \frac{k_x}{k_\beta ^2}. \end{equation} (B21)The spatial wavefield gradient-based upgoing P- and S-wave potentials are obtained by inserting eqs (B13)–(B21) into eqs (B1)–(B3) and transforming back to the space–frequency domain:   \begin{equation} (\nabla \cdot v^U) \approx -\frac{1}{2}p + \frac{1}{i\omega } \rho \beta ^2 \left(\frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y}\right) + \frac{\alpha \rho }{2} v_z, \end{equation} (B22)  \begin{equation} (\nabla \times v^U)_x \approx \frac{\beta }{2 i \omega } \frac{\partial p}{\partial y} + \frac{\rho \beta }{2} v_y - \frac{\rho \beta ^2}{i \omega } \frac{\partial v_z}{\partial y}, \end{equation} (B23)  \begin{equation} (\nabla \times v^U)_y \approx \frac{\beta }{2 i \omega } \frac{\partial p}{\partial x} + \frac{\rho \beta }{2} v_x - \frac{\rho \beta ^2}{i \omega } \frac{\partial v_z}{\partial x}. \end{equation} (B24) © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

Spatial wavefield gradient-based seismic wavefield separation

Loading next page...
 
/lp/ou_press/spatial-wavefield-gradient-based-seismic-wavefield-separation-SJkqjiuWlP
Publisher
The Royal Astronomical Society
Copyright
© The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx499
Publisher site
See Article on Publisher Site

Abstract

Summary Measurements of the horizontal and vertical components of particle motion combined with estimates of the spatial gradients of the seismic wavefield enable seismic data to be acquired and processed using single dedicated multicomponent stations (e.g. rotational sensors) and/or small receiver groups instead of large receiver arrays. Here, we present seismic wavefield decomposition techniques that use spatial wavefield gradient data to separate land and ocean bottom data into their upgoing/downgoing and P/S constituents. Our method is based on the elastodynamic representation theorem with the derived filters requiring local measurements of the wavefield and its spatial gradients only. We demonstrate with synthetic data and a land seismic field data example that combining translational measurements with spatial wavefield gradient estimates allows separating seismic data recorded either at the Earth’s free-surface or at the sea bottom into upgoing/downgoing and P/S wavefield constituents for typical incidence angle ranges of body waves. A key finding is that the filter application only requires knowledge of the elastic properties exactly at the recording locations and is valid for a wide elastic property range. Numerical modelling, Body waves, Rotational seismology, Wave propagation 1 INTRODUCTION The decomposition of multicomponent land and ocean bottom seismic (OBS) recordings into their upgoing/downgoing and P/S wavefield constituents is a key processing step for multicomponent data analysis. Important applications of wavefield separation techniques are water layer multiple removal and true amplitude imaging. However, wavefield separation of seismic data remains one of the main challenges in seismic data processing to date. The simplest conventional approach to multicomponent data wavefield separation is based on the assumption that the vertical component measurements correspond to P-wave recordings, whereas S-waves are recorded solely on the horizontal components. However, this assumption breaks down at non-normal incidence angles for which P- and S-waves are recorded both on the vertical and horizontal components. More sophisticated approaches to the separation of wave modes on multicomponent recordings have been proposed using array-based filtering techniques (e.g. Dankbaar 1985; Wapenaar et al. 1990; Wapenaar & Haimé 1990; Amundsen et al. 2000). Dankbaar (1985) derived P/S separation filters as a function of the geophone receiving characteristics. Wapenaar et al. (1990) formulated a wave equation-based approach to decompose multicomponent surface seismic data into their primary P- and S-wave components. Wapenaar & Haimé (1990) used the elastodynamic representation theorem to derive upgoing/downgoing and P/S wavefield decomposition filters. Amundsen et al. (2000) extended this approach to multicomponent OBS recordings for receiver-side water layer multiple removal and P/S separation. Array-based wavefield separation filters typically require long arrays of densely spaced multicomponent stations, which significantly increases the field effort and cost and is not feasible in areas with limited access (e.g. the ocean bottom). Additionally, the application of the filters requires knowledge of the near-surface elastic properties. Lateral variations in elastic properties and statics are likely to take place along such extended receiver lines and can result in incorrect separation results. These limitations of array-based decomposition techniques generally render them less practical for field applications. The use of spatial wavefield gradient data as a new observable in seismic exploration has gained increasing attention in recent years and could open up new ways of acquiring and processing multicomponent seismic data (Igel et al. 2012; Schmelzbach et al. 2016b). If multicomponent receivers are closely spaced, not only the vector wavefield but also the spatial gradients of the seismic wavefield can be estimated from the array data (e.g. Spudich et al. 1995; Langston 2007a,b; Liu & Holt 2015; Porter et al. 2016; Sollberger et al. 2016a). Wavefield separation techniques based on spatial wavefield gradient data have been proposed by Robertsson & Muyzert (1999) and Robertsson & Curtis (2002), using spatial wavefield gradient estimates from spatially compact groups of three-component (3C) receivers. Robertsson & Muyzert (1999) proposed a compact volume acquisition design consisting of only four stations to separate the curl- and divergence-free components of the wavefield, but their method requires the presence of a receiver at depth to compute the vertical spatial derivatives. A major step forward was the realization by Robertsson & Curtis (2002) that the free-surface boundary condition allows expressing the vertical derivatives of the particle velocity components in terms of horizontal derivatives. Hence, recording stations at depth to estimate vertical gradients are not required to isolate the upgoing P- and S-waves of land surface seismic data. In multicomponent OBS exploration, spatial wavefield gradient data have mainly been used for deghosting of pressure data (e.g. Osen et al. 2002; Amundsen et al. 2013). Barak et al. (2013, 2015) recently discussed the potential benefit of using rotational data acquired on the sea bottom for enhanced wavefield characterization. In this paper, we present spatial wavefield gradient-based techniques to separate land and OBS recordings into their (1) upgoing/downgoing and (2) P/S wavefield constituents. Up/down separation of OBS data is important for water layer multiple suppression, whereas in the land seismic case we discuss up/down decomposition as an important step towards P/S separation for imaging. In contrast to Robertsson & Curtis (2002), we present filters to isolate the horizontal and vertical particle velocity components of the upgoing S- and P-waves, respectively, which represent more straightforward quantities for interpretation instead of the originally proposed P- and S-wave potentials. For OBS data, we extend the seabed wavefield decomposition filters by Amundsen et al. (2000) to become applicable to local measurements of the wavefield and its spatial gradients only. Furthermore, we investigate the effect of strong near-surface velocity gradients on the wavefield separation results and make the important observation that the P- and S-wave velocities required for our technique are only those exactly at the receiver locations. Finally, we apply our method to field data acquired on land, which, to the best of our knowledge, is the first application of spatial wavefield gradient-based filters to land seismic data. 2 METHODOLOGY Our wavefield separation filters are derived from the elastodynamic representation theorem, which allows to calculate the seismic wavefield inside a volume from recordings along the surface enclosing the volume (Wapenaar & Haimé 1990; Amundsen et al. 2000). Adjusting the surface geometry with corresponding boundary conditions to either the free-surface (zero vertical traction) or the fluid/solid boundary at the ocean bottom (zero tangential stress) and introducing analytical expressions for the Green’s tensor in a homogeneous medium allow for the derivation of formulations for the upgoing wavefield at an infinitesimal distance below the recording interface. 2.1 Wavefield separation of land seismic data using spatial wavefield gradients 2.1.1 Up/down separation Land seismic data acquired at the Earth’s free-surface comprise of superimposed recordings of upgoing waves (P and S) and their downgoing reflections and mode-conversions (Aki & Richards 2002). As a consequence, the recorded composite ground motion differs in polarization from the incident waves dependent on the incidence angle and the medium properties. Up/down separation of land seismic data allows to extract the true amplitude upgoing wavefield and can therefore significantly improve, for example, amplitude-versus-offset analysis and full-waveform inversion. Expressions for the upgoing horizontal ($$v_x^{U}$$ and $$v_y^{U}$$) and vertical ($$v_z^{U}$$) particle velocity components that are exact for homogeneous and isotropic media at an infinitesimal distance below the Earth’s surface are given by (Wapenaar & Haimé 1990):   \begin{equation} v_x^{U} = \frac{1}{2}\left(v_x + F^{v_x}_{v_z} \ast v_z\right), \end{equation} (1)  \begin{equation} v_y^{U} = \frac{1}{2}\left(v_y + F^{v_y}_{v_z} \ast v_z\right) , \end{equation} (2)  \begin{equation} v_z^{U} = \frac{1}{2}\left(v_z - F^{v_z}_{v_x} \ast v_x - F^{v_z}_{v_y} \ast v_y\right) , \end{equation} (3)where vx and vy, and vz are the total horizontal and vertical particle velocities, respectively, recorded at the free-surface. The superscripts of the filter F denote the wavefield component to be extracted, whereas the subscripts denote the wavefield quantity on which the filter acts. Spatiotemporal convolutions are denoted by an asterisk. In the wavenumber–frequency domain, the filters can be expressed as (Wapenaar & Haimé 1990):   \begin{equation} \widetilde{F}^{v_x}_{v_z} = \frac{k_x}{k_z^{(\alpha )}} \left(1-\frac{2}{k_\beta ^2}\left(k_x^2 + k_y^2 + k_z^{(\alpha )}k_z^{(\beta )}\right)\right), \end{equation} (4)  \begin{equation} \widetilde{F}^{v_y}_{v_z} = \frac{k_y}{k_z^{(\alpha )}} \left(1-\frac{2}{k_\beta ^2}\left(k_x^2 + k_y^2 + k_z^{(\alpha )}k_z^{(\beta )}\right)\right) , \end{equation} (5)  \begin{equation} \widetilde{F}^{v_z}_{v_x} = \frac{k_x}{k_z^{(\beta )}} \left(1-\frac{2}{k_\beta ^2}\left(k_x^2 + k_y^2 + k_z^{(\alpha )}k_z^{(\beta )}\right)\right), \end{equation} (6)  \begin{equation} \widetilde{F}^{v_z}_{v_y} = \frac{k_y}{k_z^{(\beta )}} \left(1-\frac{2}{k_\beta ^2}\left(k_x^2 + k_y^2 + k_z^{(\alpha )}k_z^{(\beta )}\right)\right), \end{equation} (7)where $$k_z^{(\alpha )} = k_\alpha \sqrt{1-\frac{k_x^2 + k_y^2}{k_\alpha ^2}}$$ and $$k_z^{(\beta )} = k_\beta \sqrt{1-\frac{k_x^2 + k_y^2}{k_\beta ^2}}$$ are the unsigned vertical wavenumbers for P- and S-waves, respectively, and kx and ky are the horizontal wavenumbers. The P- and S-wavenumbers are denoted by kα = ω/α and kβ = ω/β, respectively, where α is the P-wave velocity, β the S-wave velocity and ω the angular frequency. The implementation of eqs (4)–(7) would require horizontal wavenumber estimates conventionally derived from long receiver arrays. Instead, we seek formulations that only contain first-order horizontal spatial derivatives of the wavefield (corresponding to kx and ky) and that can be estimated using spatially compact receiver arrays. This is achieved by keeping the first term in the Taylor approximations around kx = 0 and ky = 0 to eqs (4)–(7):   \begin{equation} \widetilde{F}^{v_x}_{v_z} \approx \frac{k_x}{k_\alpha }\left(1 - \frac{2\beta }{\alpha }\right) , \end{equation} (8)  \begin{equation} \widetilde{F}^{v_y}_{v_z} \approx \frac{k_y}{k_\alpha }\left(1 - \frac{2\beta }{\alpha }\right) , \end{equation} (9)  \begin{equation} \widetilde{F}^{v_z}_{v_x} \approx \frac{k_x}{k_\beta }\left(1 - \frac{2\beta }{\alpha }\right), \end{equation} (10)  \begin{equation} \widetilde{F}^{v_z}_{v_y} \approx \frac{k_y}{k_\beta }\left(1 - \frac{2\beta }{\alpha }\right). \end{equation} (11) Inserting eqs (8)–(11) into eqs (1)–(3) and transforming the expressions back to the space–frequency domain yields the following expressions for the upgoing particle velocity components:   \begin{equation} v_x^U \approx \frac{1}{2}\left(v_x + \frac{1}{i\omega }(\alpha -2\beta )\frac{\partial v_z}{\partial x}\right), \end{equation} (12)  \begin{equation} v_y^U \approx \frac{1}{2}\left(v_y + \frac{1}{i\omega }(\alpha -2\beta )\frac{\partial v_z}{\partial y}\right), \end{equation} (13)  \begin{equation} v_z^U \approx \frac{1}{2}\left(v_z - \frac{1}{i\omega }\left(\beta -\frac{2\beta ^2}{\alpha }\right)\left(\frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y}\right)\right). \end{equation} (14) Eqs (12)–(14) combine the total recorded wavefield with spatial gradient estimates scaled by the seismic velocities to suppress the downgoing reflections and mode conversions generated at the Earth’s surface. It is also worthwhile to note that the first terms in eqs (12)–(14), which are equal to half of the total recorded wavefield, correspond to the conventional interpretation of multicomponent data assuming that vertical component measurements correspond to P-waves, whereas the horizontal components correspond to S-waves. This only holds for vertically incident waves. In contrast, eqs (12)–(14) better approximate the upgoing wavefield for non-normal incidence angles. We test the performance of our up/down separation filters on synthetic data generated using a 2-D finite-difference (FD) code (Virieux 1986). A force source emitting a 50 Hz Ricker wavelet was placed at 100 m depth in a homogeneous medium with P- and S-wave velocities of 1800 and 600 m s−1, respectively. The incident P- and S-wave arrivals were recorded at the free-surface. The free-surface condition was implemented by imaging the stresses and particle velocities on several grid points above the free-surface (Levander 1988; Crase 1990), which we found to be in good agreement with an analytical solution (Sherwood 1958). A reference solution of the upgoing wavefield at the free-surface was derived using the FD-injection method presented by Robertsson et al. (2015) that allows decomposing FD-computed seismograms into their upgoing and downgoing constituents to within machine precision. Fig. 1 displays the recorded horizontal and vertical particle velocity components and the up/down separated results using eqs (12) and (14), respectively. To assess the performance of the spatial wavefield gradient filters we calculated the relative root-mean-square error (RMSE) against incidence angle between the estimated upgoing wavefield and the FD-injection reference solution for the isolated P- or S-wave arrivals at the free-surface (Fig. 2). The dashed lines correspond to taking half of the total recorded wavefield (denoted here as the conventional method), whereas the solid lines correspond to the error from using our spatial gradient-based eqs (12) and (14). Fig. 2(a) shows that our spatial gradient filters significantly improve the approximation to the horizontal component of the upgoing P-wave compared to the conventional approach, with errors below 10 per cent up to incidence angles of about 30°. In contrast, simply taking half of the total wavefield to approximate the vertical component of the upgoing P-wavefield already results in errors below 10 per cent up to 50° incidence (Fig. 2b) and our spatial gradient filters do not significantly reduce the error. In case of an incident P-wave recorded on the vertical component, the amplitudes of the upgoing P-wave and its downgoing reflections have similar magnitudes, while the contributions of P-to-S conversions to the vertical component recordings are very small and predominantly project onto the horizontal component recording (Aki & Richards 2002). This observation explains why taking half of the total vertical component recording is already a good approximation to the vertical component of the upgoing P-wave, whereas our spatial wavefield gradient-based filters are required to better estimate the horizontal component of the upgoing P-wave. Figure 1. View largeDownload slide Up/down separation of a synthetic land seismic data example. (a) Total horizontal component recording. (b) Total vertical component recording. (c) Estimated upgoing horizontal component. (d) Estimated upgoing vertical component. Figure 1. View largeDownload slide Up/down separation of a synthetic land seismic data example. (a) Total horizontal component recording. (b) Total vertical component recording. (c) Estimated upgoing horizontal component. (d) Estimated upgoing vertical component. Figure 2. View largeDownload slide Relative RMSE between the FD-injection reference solution and the estimated particle velocity components of isolated upgoing P- and S-waves plotted against incidence angle for the synthetic land seismic data example in Fig. 1. The dashed line corresponds to the error from using the conventional method (taking half of the total wavefield), while the solid line represents the error from using spatial wavefield gradients (eqs 12 and 14). (a) Error for the horizontal component of the upgoing P-wave. (b) Error for the vertical component of the upgoing P-wave. (c) Error for the horizontal component of the upgoing S-wave. (d) Error for the vertical component of the upgoing S-wave. Figure 2. View largeDownload slide Relative RMSE between the FD-injection reference solution and the estimated particle velocity components of isolated upgoing P- and S-waves plotted against incidence angle for the synthetic land seismic data example in Fig. 1. The dashed line corresponds to the error from using the conventional method (taking half of the total wavefield), while the solid line represents the error from using spatial wavefield gradients (eqs 12 and 14). (a) Error for the horizontal component of the upgoing P-wave. (b) Error for the vertical component of the upgoing P-wave. (c) Error for the horizontal component of the upgoing S-wave. (d) Error for the vertical component of the upgoing S-wave. In case of an incident S-wave, the spatial wavefield gradient filters do not provide a better approximation to the upgoing horizontal component compared to taking half of the total wavefield (Fig. 2c), whereas our spatial gradient filters enable to extract more accurate amplitude information on the vertical component of the upgoing S-wave in comparison to the conventional solution (Fig. 2d) with errors below 10 per cent up to about 20° incidence angle. The reason why the incidence angle range over which the upgoing S-wave can be reliably estimated is much smaller than for a P-wave is due to the free-surface effect. In particular, S-wave arrivals at the free-surface generate horizontally propagating S-to-P mode conversions for incidence angles beyond the critical angle and these waveforms cannot be removed by our approximated filter expressions. The critical angle is a function of the near-surface P- and S-wave velocities and is about 20° in our example. However, we anticipate that in case of field data, due to strong vertical velocity gradients in the shallow subsurface, body waves will mostly arrive at steep incidence angles close to vertical and well below the critical angle for S-waves. 2.1.2 P/S separation The up/down separation of land seismic data removes the effect of the free-surface, but for imaging purposes an additional P/S separation step is desirable. Robertsson & Curtis (2002) derived expressions to estimate the upgoing P- and S-wave potentials using spatially compact receiver groups (see Appendix A). In contrast to Robertsson & Curtis (2002), we prefer to work with particle velocities instead of wavefield potentials to obtain the horizontal component of the upgoing S-wave constituents and the vertical component of the upgoing P-wave constituents. Similar to the up/down separation technique, we seek expressions that remove undesired wavefield constituents from the total recorded particle velocity components. This can be achieved by normalizing eqs (A14)–(A16) as follows:   \begin{equation} v^{S^U}_x \approx \frac{\beta }{i\omega } (\nabla \times v^U)_y \approx \frac{1}{2} \left(v_x - \frac{1}{i \omega } 2 \beta \frac{\partial v_z}{\partial x}\right), \end{equation} (15)  \begin{equation} v^{S^U}_y \approx \frac{\beta }{i\omega } (\nabla \times v^U)_x \approx \frac{1}{2} \left(v_y - \frac{1}{i \omega } 2 \beta \frac{\partial v_z}{\partial y}\right), \end{equation} (16)  \begin{equation} v^{P^U}_z \approx \frac{\alpha }{i\omega } (\nabla \cdot v^U) \approx \frac{1}{2} \left(v_z + \frac{1}{i\omega } \frac{2\beta ^2}{\alpha } \left(\frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} \right)\right). \end{equation} (17) Eqs (15)–(17) separate the recorded particle velocity components into their P- and S-wave constituents and additionally remove the free-surface effect. For example, the spatial wavefield gradient term ∂vx/∂x + ∂vy/∂y in eq. (17) is proportional to the divergence at the free-surface (Robertsson & Curtis 2002). The divergence measured at the free-surface is sensitive to S-to-P conversions and properly scaled by the medium parameters it effectively removes the total vertical component S-wave constituent. In case of an incident P-wave, this divergence term does not only contain the upgoing P-waves but also the downgoing reflected P-waves and is close to zero (destructive interference) when the amplitudes of both constituents have similar magnitudes. Because the vertical components of the upgoing and downgoing P-waves are about equal in amplitude for a considerable range of incidence angles (Aki & Richards 2002), the vertical component of the upgoing P-wave is retained by scaling the total vertical component P-wave arrival by a factor 1/2 (see Fig. 2b). To demonstrate the P/S separation technique on synthetic land seismic data we use the same 2-D model as for the up/down separation. Fig. 3 displays the total horizontal and vertical particle velocity recordings and the P/S separated results using eqs (15) and (17). Fig. 4 shows the relative RMSE versus incidence angle between the P/S separated results using spatial wavefield gradients (eqs 15 and 17) and the FD-injection reference solution for an upgoing P-wave arrival recorded on the vertical component and an upgoing S-wave arrival recorded on the horizontal component (solid lines). For comparison, we also plotted the error from using the conventional method (dashed lines), which assumes P-waves to be recorded on the vertical component and S-waves to be recorded on the horizontal component only. Our spatial wavefield gradient method significantly improves the P/S separated results for the upgoing horizontal and vertical components, reducing the errors below 10 per cent up to about 25° and 20° incidence, respectively. Although, the P-wave arrival recorded on the horizontal component is fully removed up to 50° incidence, the main error source arises from extracting the horizontal component of the upgoing S-wave beyond the critical angle. Similarly, the vertical component of the incident S-wave is only suppressed up to the critical angle. Nevertheless, the improved P/S separated results demonstrate that our spatial wavefield gradient filters result in good quality recordings for P/S imaging and are particularly attractive to enhance S-wave analyses on horizontal component recordings. Figure 3. View largeDownload slide P/S separation of a synthetic land seismic data example. (a) Total horizontal component recording. (b) Total vertical component recording. (c) Estimated upgoing horizontal component with the P-wave arrival suppressed. (d) Estimated upgoing vertical component with the S-wave arrival suppressed. Figure 3. View largeDownload slide P/S separation of a synthetic land seismic data example. (a) Total horizontal component recording. (b) Total vertical component recording. (c) Estimated upgoing horizontal component with the P-wave arrival suppressed. (d) Estimated upgoing vertical component with the S-wave arrival suppressed. Figure 4. View largeDownload slide Relative RMSE between the FD-injection reference solution and the estimated P/S separation results for the synthetic land seismic data example in Fig. 3. The dashed line corresponds to the error from using the conventional method (P- and S-waves only recorded on the vertical and horizontal components, respectively), while the solid line represents the error from using spatial wavefield gradients (eqs 15 and 17). (a) Error for the P-wave suppression on the upgoing horizontal component. (b) Error for the S-wave suppression on the upgoing vertical component. Figure 4. View largeDownload slide Relative RMSE between the FD-injection reference solution and the estimated P/S separation results for the synthetic land seismic data example in Fig. 3. The dashed line corresponds to the error from using the conventional method (P- and S-waves only recorded on the vertical and horizontal components, respectively), while the solid line represents the error from using spatial wavefield gradients (eqs 15 and 17). (a) Error for the P-wave suppression on the upgoing horizontal component. (b) Error for the S-wave suppression on the upgoing vertical component. 2.2 Wavefield separation of OBS data using spatial wavefield gradients 2.2.1 Up/down separation Up/down separation based on the elastodynamic representation theorem and using spatial wavefield gradients can be extended to OBS data for receiver-side water layer multiple removal. Amundsen et al. (2000) derived analytical expressions for the upgoing horizontal and vertical wavefields that are exact for homogeneous and isotropic media at an infinitesimal distance below the seabed:   \begin{equation} v_x^{U} = \frac{1}{2}\left(v_x + F^{v_x}_{v_z} \ast v_z\right), \end{equation} (18)  \begin{equation} v_y^{U} = \frac{1}{2}\left(v_y + F^{v_y}_{v_z} \ast v_z\right), \end{equation} (19)  \begin{equation} v_z^{U} = \frac{1}{2}\left(v_z - F^{v_z}_{v_x} \ast v_x - F^{v_z}_{v_y} \ast v_y - F^{v_z}_{p} \ast p\right). \end{equation} (20) The filters operating on the particle velocity components are identical to those derived for the land seismic case (eqs 4–7). Eq. (20) contains an additional spatial filter term $$F^{v_z}_{p}$$ that operates on the pressure recordings p. In the wavenumber–frequency domain this filter is given by (Amundsen et al. 2000):   \begin{equation} \widetilde{F}^{v_z}_{p} = \frac{1}{\rho \omega k^{(\beta )}_z}\left(k^2_x + k^2_y + k^{(\alpha )}_z k^{(\beta )}_z \right), \end{equation} (21) where ρ is the density of the seafloor. In order to adapt the up/down separation technique of OBS data to local measurements of the wavefield only, we follow the same approach as for the land seismic case by expanding the filters up to first-order spatial gradients. For $$\widetilde{F}^{v_z}_{p}$$ this results in   \begin{equation} \widetilde{F}^{v_z}_{p} \approx \frac{k_\alpha }{\rho \omega }. \end{equation} (22) Inserting the truncated filter expressions (8)–(11) and (22) into eqs (18)–(20) and transforming back to the space–frequency domain results in the following expressions for the upgoing particle velocity components at an infinitesimal distance below the seabed:   \begin{equation} v_x^U \approx \frac{1}{2}\left(v_x + \frac{1}{i\omega }(\alpha -2\beta )\frac{\partial v_z}{\partial x}\right), \end{equation} (23)  \begin{equation} v_y^U \approx \frac{1}{2}\left(v_y + \frac{1}{i\omega }(\alpha -2\beta )\frac{\partial v_z}{\partial y}\right), \end{equation} (24)  \begin{equation} v_z^U \approx \frac{1}{2}\left(v_z -\frac{1}{i\omega } \left(\beta -\frac{2\beta ^2}{\alpha }\right)\left(\frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y}\right) - \frac{1}{\alpha \rho } p\right). \end{equation} (25) Eqs (23)–(25) suppress all downgoing events, including the direct wave, receiver ghosts and accompanying water layer reverberations (receiver-side multiples) in addition to remaining downgoing reflections and mode conversions generated at the seabed. We demonstrate the spatial gradient-based multiple attenuation technique on synthetic FD modelled vertical particle velocity recordings at the seabed. Our 2-D model consists of a 100 m thick water layer bounded above by a free-surface. Below the water layer is a homogeneous medium with a P-wave velocity of 1800 m s−1, an S-wave velocity of 600 m s−1 and a density of 1600 m kg−3. An explosive point source emitting a 50 Hz Ricker wavelet was placed at 5 m depth below the sea surface. In order to generate overlapping downgoing and upgoing arrivals recorded by a multicomponent OBS sensor, a second explosive point source was located at 150 m below the seafloor. Figs 5(a) and (b) show the total as well as the estimated upgoing vertical particle velocity recordings, obtained by evaluating eq. (25). The downgoing direct wave (black arrow) is effectively suppressed, while the upgoing arrivals (green arrows) are preserved, even in case of full overlap of the upgoing and downgoing events as observed for example at about 100 m offset. The later arriving water layer reverberations (grey arrows) are also effectively suppressed. Figure 5. View largeDownload slide Receiver-side water layer multiple suppression of a synthetic OBS data example. (a) Total vertical component recording. (b) Estimated upgoing vertical component. The downgoing direct wave is denoted by the black arrow, whereas the upgoing waves are indicated by the green arrows. Downgoing water layer reverberations are highlighted by the grey arrows. Figure 5. View largeDownload slide Receiver-side water layer multiple suppression of a synthetic OBS data example. (a) Total vertical component recording. (b) Estimated upgoing vertical component. The downgoing direct wave is denoted by the black arrow, whereas the upgoing waves are indicated by the green arrows. Downgoing water layer reverberations are highlighted by the grey arrows. The performance of eq. (25) to suppress the downgoing direct wave is shown in Fig. 6, which displays the relative RMSE between the filtered data and the FD-injection reference solution for the upgoing wavefield (equal to zero) (solid line). The dashed line represents the error from using the conventional multiple suppression approach, which corresponds to combining the particle velocity and pressure recordings only (PZ summation, Amundsen 1993) and ignoring the spatial gradient term in eq. (25). Both approaches effectively suppress vertical component water layer multiples with errors below 10 per cent up to incidence angles of about 40°, although our spatial wavefield gradient-based technique does not significantly improve the results compared to the conventional PZ summation. Figure 6. View largeDownload slide Relative RMSE between the FD-injection reference solution and the estimated upgoing vertical component plotted against incidence angle for the synthetic OBS data example in Fig. 5. The dashed line corresponds to the error from using the conventional method (PZ summation), while the solid line represents the error from using spatial wavefield gradients (eq. 25). Figure 6. View largeDownload slide Relative RMSE between the FD-injection reference solution and the estimated upgoing vertical component plotted against incidence angle for the synthetic OBS data example in Fig. 5. The dashed line corresponds to the error from using the conventional method (PZ summation), while the solid line represents the error from using spatial wavefield gradients (eq. 25). 2.2.2 P/S separation P/S separation based on the elastodynamic representation theorem using spatial wavefield gradients can be extended to OBS data and the corresponding derivations are given in Appendix B. Similar to the land seismic case, we prefer to interpret the horizontal particle velocity component of the upgoing S-wave constituents and the vertical particle velocity component of the upgoing P-wave constituents, which can be achieved by normalizing the P- and S-wave potential equations (B22)–(B24) as follows:   \begin{equation} v_x^{S^U} \approx \frac{1}{\rho \beta } (\nabla \times v^U)_y \approx \frac{1}{2} \left(v_x - \frac{1}{i\omega }2 \beta \frac{\partial v_z}{\partial x} + \frac{1}{i \omega } \frac{1}{\rho } \frac{\partial p}{\partial x}\right), \end{equation} (26)  \begin{equation} v_y^{S^U} \approx \frac{1}{\rho \beta } (\nabla \times v^U)_x \approx \frac{1}{2} \left(v_y - \frac{1}{i \omega } 2 \beta \frac{\partial v_z}{\partial y} + \frac{1}{i \omega } \frac{1}{\rho } \frac{\partial p}{\partial y}\right), \end{equation} (27)  \begin{eqnarray} v_z^{P^U} \approx \frac{1}{\alpha \rho } (\nabla \cdot v) \approx \frac{1}{2} \left(v_z + \frac{1}{i\omega } \frac{2 \beta ^2}{\alpha } \left(\frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y}\right) - \frac{1}{\alpha \rho }p \right).\nonumber\\ \end{eqnarray} (28) In addition to extracting the upgoing P- or S-waves, eqs (26)–(28) also suppress all downgoing events such as the receiver-side water layer multiples. To illustrate the P/S separation method, OBS data were modelled using the same 2-D subsurface model as for the up/down separation but now using a force source emitting a 50 Hz Ricker wavelet at 100 m below the seafloor to generate both P- and S-waves. Fig. 7 shows the total horizontal and vertical particle velocity recordings at the seafloor and the P/S separated results using eqs (26) and (28). Our technique effectively suppresses the horizontal component of the P-wave and the vertical component of the S-wave up to moderate incidence angles. The sea surface multiple (receiver ghost), which is observed on the vertical component at about 0.21 s, has also been removed by the combined up/down and P/S separation. Figure 7. View largeDownload slide P/S separation of a synthetic OBS data example. (a) Total horizontal component recording. (b) Total vertical component recording. (c) Estimated upgoing horizontal component with the P-wave arrival suppressed. (d) Estimated upgoing vertical component with the S-wave arrival suppressed. Figure 7. View largeDownload slide P/S separation of a synthetic OBS data example. (a) Total horizontal component recording. (b) Total vertical component recording. (c) Estimated upgoing horizontal component with the P-wave arrival suppressed. (d) Estimated upgoing vertical component with the S-wave arrival suppressed. The performance of the seabed P/S separation filters (solid line, eqs 26 and 28) as well as the conventional method (dashed line) against FD-injection as a function of incidence angle is displayed in Fig. 8 and shows very similar behaviour as for the land seismic case (Fig. 4). Figure 8. View largeDownload slide Relative RMSE between the FD-injection reference solution and the estimated P/S separation results for the synthetic OBS data example in Fig. 7. The dashed line corresponds to the error from using the conventional method (P- and S-waves only recorded on the vertical and horizontal components, respectively), while the solid line represents the error from using spatial wavefield gradients (eqs 26 and 28). (a) Error for the P-wave suppression on the upgoing horizontal component. (b) Error for the S-wave suppression on the upgoing vertical component. Figure 8. View largeDownload slide Relative RMSE between the FD-injection reference solution and the estimated P/S separation results for the synthetic OBS data example in Fig. 7. The dashed line corresponds to the error from using the conventional method (P- and S-waves only recorded on the vertical and horizontal components, respectively), while the solid line represents the error from using spatial wavefield gradients (eqs 26 and 28). (a) Error for the P-wave suppression on the upgoing horizontal component. (b) Error for the S-wave suppression on the upgoing vertical component. 3 VELOCITIES FOR WAVEFIELD SEPARATION The application of wavefield separation techniques using spatial filters requires the knowledge of the ‘near-surface’ P- and S-wave velocities (e.g. Dankbaar 1985; Robertsson & Curtis 2002). It has been unclear whether these ‘near-surface’ velocities are the velocities at or below the free-surface, possibly averaged over some volume, and whether they are frequency dependent or not. Robertsson & Curtis (2002) speculated that the required P- and S-wave velocities likely are those experienced by the wavefield, whereby longer wavelengths are influenced by the elastic properties of a larger volume. In order to experimentally investigate the velocity dependence of our wavefield separation filters, we apply our techniques to a synthetic data example with a strong vertical velocity gradient in the shallow subsurface. The P- and S-wave velocities in our synthetic model linearly increase from 1800 to 4800 m s−1 and from 600 to 1800 m s−1, respectively, between the receivers at the free-surface and a depth of 100 m. An explosive source emitting a 50 Hz Ricker wavelet was excited 100 m below the surface. We used an iterative frequency-domain Gauss–Newton inversion scheme (Pratt et al. 1998) to determine the velocities that minimize the RMSE between the estimated upgoing horizontal and vertical wavefield obtained using eqs (12)–(14) and the reference upgoing wavefield obtained by means of FD-injection. The inversion approach consistently yields the P- and S-wave velocities exactly at the free-surface for frequencies above 30 Hz (Fig. 9). Furthermore, the accuracy of the wavefield separation results in case of a strong vertical velocity gradient is the same as in case of a homogeneous subsurface (see Figs 2a and b). Below 30 Hz, the optimum velocity seems to be frequency-dependent with the inverted P-wave velocities being lower than the minimum P-wave velocity in our subsurface model, whereas the inverted S-wave velocities are higher. The frequency dependence can be attributed to near-field effects: for the lowest frequencies in our example, the plane wave assumption inherent to the up/down separation algorithm is not met. In the far field, the velocities required to separate seismic data into their desired wavefield constituents using spatial wavefield gradients are exactly those at the receiver locations. Our synthetic data example demonstrates that highly accurate results can be achieved in the presence of a vertical velocity gradient. Bearing in mind that our wavefield decomposition filters were derived under a homogeneous space Green’s function assumption, this is perhaps not an entirely intuitive result. Figure 9. View largeDownload slide Inverted P- and S-wave velocities as a function of frequency in order to minimize the error of the up/down wavefield separation technique compared to the FD-injection reference solution in case of a strong vertical velocity gradient in the velocity model. Figure 9. View largeDownload slide Inverted P- and S-wave velocities as a function of frequency in order to minimize the error of the up/down wavefield separation technique compared to the FD-injection reference solution in case of a strong vertical velocity gradient in the velocity model. Our results are consistent with observations recently made by Robertsson et al. (2015). Robertsson et al. (2015) found that in order to isolate the upgoing/downgoing constituents of seismic recordings using the FD-injection method, only knowledge of the material properties exactly at the recording location is needed. This observation has so far not been fully appreciated, but has significant implications for field data processing. The shallow subsurface is typically characterized by strong velocity variations and the fact that only the knowledge of the P- and S-wave velocities at the recording locations are required significantly simplifies wavefield separation in practice. 4 FIELD DATA EXAMPLE We applied our spatial wavefield gradient-based technique to a multicomponent land seismic data set acquired in northern Switzerland to characterize the geometry and infill of a Quaternary valley (Schmelzbach et al. 2014; Sollberger et al. 2014). Near-surface studies benefit from the enhanced resolution and improved lithological characterization achieved by S-wave and converted-wave exploration (Helbig & Mesdag 1982). We therefore test our P/S separation filters to suppress P-waves on the horizontal particle velocity component recordings for S-wave imaging. The survey consisted of a 250 m long 2-D line of 3C geophones spaced at 1.5 m intervals. Seismic data were generated by striking a prismatic wedge (3C vector source, Schmelzbach et al. 2016a) with a sledgehammer at 3 m intervals along the recording line. Fig. 10(a) displays a horizontal component recording of a vertically directed source. Based on hodogram analyses and arrival identifications previously carried out by Sollberger et al. (2016b), the P-wave first arrival is marked by the red arrow, whereas the S-wave first arrival is highlighted by the blue arrow. At later times, the shot gather is mainly dominated by Rayleigh waves (green arrow). Figure 10. View largeDownload slide P/S separation of field data. The P- and S-wave first arrivals are indicated by the red and blue arrows, respectively, whereas a Rayleigh wave is indicated by the green arrow. (a) Total horizontal component recording. (b) Estimated upgoing horizontal component with the P-wave arrival suppressed. Figure 10. View largeDownload slide P/S separation of field data. The P- and S-wave first arrivals are indicated by the red and blue arrows, respectively, whereas a Rayleigh wave is indicated by the green arrow. (a) Total horizontal component recording. (b) Estimated upgoing horizontal component with the P-wave arrival suppressed. The P/S separation of the recorded horizontal particle velocity was carried out using eq. (15). The horizontal spatial wavefield gradients of the vertical particle velocities were calculated by finite-differencing the recordings from three neighbouring stations and centring the estimate at the middle receiver. The S-wave velocity was obtained from the traveltime tomography results in Schmelzbach et al. (2014) and is about 250 m s−1 at the surface. Fig. 10(b) shows that our spatial wavefield gradient-based P/S separation filters effectively suppress the undesired P-wave arrival to obtain a horizontal component recording that only contains S-waves. To verify the wavefield separation results we carried out the P/S separation on synthetic data that were generated from the traveltime tomography results in Schmelzbach et al. (2014). Fig. 11 displays the total and P/S separated horizontal component synthetic data and shows an overall very good agreement to both the raw and separated field data (Fig. 10). Hodogram analyses and snapshots of the FD-simulated wavefield confirm that the isolated body wave arrival in Fig. 11(b) is indeed an S-wave arrival. Figure 11. View largeDownload slide P/S separation of synthetic FD data generated from the traveltime tomography results in Schmelzbach et al. (2014). The P- and S-wave first arrivals are indicated by the red and blue arrows, respectively, whereas a Rayleigh wave is indicated by the green arrow. (a) Total horizontal component recording. (b) Estimated upgoing horizontal component with the P-wave arrival suppressed. Figure 11. View largeDownload slide P/S separation of synthetic FD data generated from the traveltime tomography results in Schmelzbach et al. (2014). The P- and S-wave first arrivals are indicated by the red and blue arrows, respectively, whereas a Rayleigh wave is indicated by the green arrow. (a) Total horizontal component recording. (b) Estimated upgoing horizontal component with the P-wave arrival suppressed. 5 DISCUSSION The techniques presented in this paper were tested on synthetic data examples with a typical near-surface P-to-S-wave velocity ratio equal to 3. It is important to stress that additional synthetic data tests showed that our observations generally hold in case of realistic near-surface Poisson’s ratios (between 0.2 and 0.45), as also demonstrated by the field data example. The spatial gradient data in our synthetic and field data examples were estimated by finite-differencing the particle velocity and pressure recordings from neighbouring stations, which in practice require spatially compact receiver groups. Recently, reliable and portable rotational sensors have been made commercially available for seismic purposes (Pierson et al. 2016; Bernauer et al. 2017). Such sensors allow to directly measure the spatial wavefield gradient quantities ∂vz/∂x and ∂vz/∂y, which correspond to the rotational motion around the two horizontal axes at the free-surface. These spatial gradients of the vertical particle velocity component are required for the up/down and P/S separation of land and seabed horizontal particle velocity component data. Collocated measurements of particle velocities and spatial wavefield gradients (e.g. 6C sensors; Sollberger et al. 2017) allow for single station wavefield separation using the expressions derived in this study, which is very attractive for field applications. 6 CONCLUSIONS We present a wavefield decomposition technique combining local measurements of particle velocity and spatial wavefield gradient data to separate multicomponent land and seabed seismic data into their upgoing/downgoing and P/S constituents. Generally, our wavefield separation filters involving first-order derivatives are more effective in suppressing P-waves than S-waves and, hence, are of particular interest for S-wave and P-to-S-converted wave analysis for steep to moderate incidence angles. The filter performance as a function of incidence angle could be improved by including higher-order derivatives, which would require larger arrays and the gradient estimation may be more prone to errors compared to estimating first-order derivatives. Spatial gradients may be derived from small sensor arrays or, given that the horizontal gradient of the vertical particle velocity corresponds at the free-surface to a rotational motion, some spatial gradients can directly be measured with rotational sensors. Therefore, our filters are particularly interesting for the application to single station multicomponent measurements with dedicated sensors at sites where one or only a few stations can be installed such as at the sea bottom or in extraterrestrial seismology. Synthetic data tests revealed that the filter performance only shows a weak dependence on the subsurface elastic parameter variations. Hence, the derived wavefield separation approach has a wide potential application range in terms of environment. An interesting observation is that the filters require as input parameters the elastic parameters found exactly at the receivers. This observation may have implications for other wave equation-based techniques that either depend on or are used to infer on material properties. As only the receiver elastic properties need to be known, the estimation of these material parameters using standard geophysical techniques (e.g. traveltime tomography, surface wave tomography) is straightforward and facilitates the application of our wavefield separation approaches to field data. In principle, the presented wavefield separation technique can easily be adopted to estimate the material properties from identified isolated arrivals. ACKNOWLEDGEMENTS This work was supported by Swiss National Science Foundation grant 200021_156996 and the sponsors of the CARNEVAL consortium (Nagra, OMV and Schlumberger Gould Research). We wish to thank editor Herve Chauris and two anonymous reviewers for constructive reviews that improved the paper. REFERENCES Aki K., Richards P.G., 2002. Quantitative Seimology , University Science Books. Amundsen L., 1993. Wavenumber-based filtering of marine point-source data, Geophysics , 58( 9), 1335– 1348. https://doi.org/10.1190/1.1443516 Google Scholar CrossRef Search ADS   Amundsen L., Ikelle L., Martin J., 2000. Multiple attenuation and P/S splitting of multicomponent OBC data at a heterogeneous sea floor, Wave Motion , 32, 67– 78. https://doi.org/10.1016/S0165-2125(99)00047-5 Google Scholar CrossRef Search ADS   Amundsen L., Ikelle L., Robertsson J.O.A., Westerdahl H., 2013. Wavefield decomposition of seabed node marine recordings, in 83rd Annual International Meeting, SEG, Expanded Abstracts , pp. 146– 150. Barak O., Brune R., Herkenhoff F., Dash R., Rector J., Ronen S., 2013. Seven-component seismic data, in 83rd Annual International Meeting, SEG, Expanded Abstracts , pp. 5151– 5155. Barak O., Key K., Constable S., Milligan P., Ronen S., 2015. Acquiring rotation data on the ocean bottom without rotation sensors, in 85th Annual International Meeting, SEG, Expanded Abstracts , pp. 2148– 2152. Bernauer F. et al.  , 2017. BlueSeis3A - full characterization of a 3C broadband rotational ground motion sensor for seismology, in 19th EGU General Assembly Conference Abstracts , p. 15512. Crase E., 1990. High-order (space and time) finite-difference modeling of the elastic wave equation, in 60th Annual International Meeting, SEG, Expanded Abstracts , pp. 987– 991. Dankbaar J., 1985. Separation of P- and S-waves, Geophys. Prospect. , 33, 970– 986. https://doi.org/10.1111/j.1365-2478.1985.tb00792.x Google Scholar CrossRef Search ADS   Helbig K., Mesdag C.S., 1982. The potential of shear wave observations, Geophys. Prospect. , 30, 413– 431. https://doi.org/10.1111/j.1365-2478.1982.tb01314.x Google Scholar CrossRef Search ADS   Igel H., Brokesova J., Evans J.R., Zembaty Z., 2012. Preface to the special issue on “Advances in rotational seismology: instrumentation, theory, observations, and engineering”, J. Seismol. , 16( 4), 571– 572. https://doi.org/10.1007/s10950-012-9307-6 Google Scholar CrossRef Search ADS   Langston C., 2007a. Spatial Gradient Analysis for Linear Seismic Arrays, Bull. seism. Soc. Am. , 97( 1), 265– 280. https://doi.org/10.1785/0120060100 Google Scholar CrossRef Search ADS   Langston C., 2007b. Wave gradiometry in two dimensions, Bull. seism. Soc. Am. , 97( 2), 401– 416. https://doi.org/10.1785/0120060138 Google Scholar CrossRef Search ADS   Levander A., 1988. Fourth-order finite-difference P–SV seismograms, Geophysics , 53, 1425– 1436. https://doi.org/10.1190/1.1442422 Google Scholar CrossRef Search ADS   Liu Y., Holt W.E., 2015. Wave gradiometry and its link with Helmholtz equation solutions applied to USArray in the eastern U.S., J. geophys. Res. , 120( 8), 5717– 5746. https://doi.org/10.1002/2015JB011982 Google Scholar CrossRef Search ADS   Osen A., Amundsen L., Reitan A., 2002. Toward optimal spatial filters for demultiple and wavefield splitting of ocean-bottom seismic data, Geophysics , 67, 1983– 1990. https://doi.org/10.1190/1.1527098 Google Scholar CrossRef Search ADS   Pierson R., Laughlin D., Brune R., 2016. Advances in rotational seismic measurements, in 86th Annual International Meeting, SEG, Expanded Abstracts , pp. 2263– 2267. Porter R., Liu Y., Holt W.E., 2016. Lithospheric records of orogeny within the continental U.S., Geophys. Res. Lett. , 43( 1), 144– 153. https://doi.org/10.1002/2015GL066950 Google Scholar CrossRef Search ADS   Pratt R., Shin C., Hicks G., 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion, Geophys. J. Int. , 133, 341– 362. https://doi.org/10.1046/j.1365-246X.1998.00498.x Google Scholar CrossRef Search ADS   Robertsson J., Curtis A., 2002. Wavefield separation using densely deployed three-component single-sensor groups in land surface-seismic recordings, Geophysics , 67( 5), 1624– 1633. https://doi.org/10.1190/1.1512809 Google Scholar CrossRef Search ADS   Robertsson J., Muyzert E., 1999. Wavefield separation using a volume distribution of three component recordings, Geophys. Res. Lett. , 26( 18), 2821– 2824. https://doi.org/10.1029/1999GL010472 Google Scholar CrossRef Search ADS   Robertsson J., van Manen D.-J., Schmelzbach C., Van Renterghem C., Amundsen L., 2015. Finite-difference modelling of wavefield constituents, Geophys. J. Int. , 203( 2), 1334– 1342. https://doi.org/10.1093/gji/ggv379 Google Scholar CrossRef Search ADS   Schmelzbach C. et al.  , 2014. Multi-method geophysical imaging of a Quaternary valley in northern Switzerland, in 84th Annual International Meeting, SEG, Expanded Abstracts , pp. 2083– 2087. Schmelzbach C., Sollberger D., Greenhalgh S.A., Horstmeyer H., Maurer H., Robertsson J.O.A., 2016a. 9C seismic data acquisition for near-surface applications: recording, waveform reciprocity and 4C rotation, in 78th Conference and Exhibition , EAGE, Workshop Programme, WS04 B03. Schmelzbach C., Sollberger D., Van Renterghem C., Häusler M., Robertsson J., Greenhalgh S., 2016b. Seismic spatial wavefield gradient and rotational rate measurements as new observables in land seismic exploration, in 18th EGU General Assembly Conference Abstracts , p. 3787. Sherwood J., 1958. Elastic wave propagation in a semi-infinite solid medium, Proc. Phys. Soc. , 71, 207– 219. https://doi.org/10.1088/0370-1328/71/2/308 Google Scholar CrossRef Search ADS   Sollberger D., Schmelzbach C., Horstmeyer H., Reiser F., Rabenstein L., Maurer H., Robertsson J.O.A., Greenhalgh S., 2014. Mapping near surface-elastic parameters of Quaternary sediments using multicomponent seismic techniques, in 16th EGU General Assembly Conference Abstracts , p. 4014. Sollberger D., Schmelzbach C., Robertsson J., Greenhalgh S., Nakamura Y., Khan A., 2016a. The shallow elastic structure of the lunar crusts: new insights from seismic wavefield gradient analysis, Geophys. Res. Lett. , 43, 1– 10. https://doi.org/10.1002/2016GL070883 Google Scholar CrossRef Search ADS   Sollberger D., Schmelzbach C., Van Renterghem C., Robertsson J., Greenhalgh S., 2016b. Single-component elastic wavefield separation at the free surface using source- and receiver-side gradients, in 86th Annual International Meeting, SEG, Expanded Abstracts , pp. 2268– 2273. Sollberger D., Schmelzbach C., Van Renterghem C., Robertsson J., Greenhalgh S., 2017. Automated, six-component, single-station ground-roll identification and suppression by combined processing of translational and rotational ground motion, in 87th Annual International Meeting, SEG, Expanded Abstracts , pp. 5064– 5068. Spudich P., Steck L.K., Hellweg M., Fletcher J.B., Baker L.M., 1995. Transient stresses at Parkfield, California, produced by the M 7.4 Landers earthquake of June 28, 1992: observations from the UPSAR dense seismograph array, J. geophys. Res. , 100, 675– 690. https://doi.org/10.1029/94JB02477 Google Scholar CrossRef Search ADS   Virieux J., 1986. P–SV wave propagation in heterogeneous media: velocity-stress finite-difference method, Geophysics , 51( 4), 889– 901. https://doi.org/10.1190/1.1442147 Google Scholar CrossRef Search ADS   Wapenaar C., Haimé G., 1990. Elastic extrapolation of primary seismic P- and S-waves, Geophys. Prospect. , 38, 23– 60. https://doi.org/10.1111/j.1365-2478.1990.tb01833.x Google Scholar CrossRef Search ADS   Wapenaar C., Hermann P., Verschuur D., Berkhout A., 1990. Decomposition of multicomponent seismic data into primary P- and S-wave responses, Geophys. Prospect. , 38, 633– 662. https://doi.org/10.1111/j.1365-2478.1990.tb01867.x Google Scholar CrossRef Search ADS   APPENDIX A: P/S SEPARATION OF LAND SEISMIC DATA Robertsson & Curtis (2002) derived the following expressions to extract the upgoing P- and S-wave potentials by taking the divergence and curl of up/down separated land seismic data:   \begin{equation} (\nabla \cdot v^U) = \frac{k_\alpha ^2}{k_\beta ^2} \left(\frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y}\right) + F_{v_z}^{(\nabla \cdot v) } \ast v_z, \end{equation} (A1)  \begin{equation} (\nabla \times v^U)_x = -\frac{\partial v_z}{\partial y} + F_{v_x}^{(\nabla \times v)_x } \ast v_x + F_{v_y}^{(\nabla \times v)_x} \ast v_y, \end{equation} (A2)  \begin{equation} (\nabla \times v^U)_y = - \frac{\partial v_z}{\partial x} + F_{v_x}^{(\nabla \times v)_y } \ast v_x + F_{v_y}^{(\nabla \times v)_y} \ast v_y, \end{equation} (A3)where in the wavenumber–frequency domain the filters can be expressed as (Robertsson & Curtis 2002):   \begin{equation} \widetilde{F}_{v_z}^{(\nabla \cdot v) } = ik_\alpha ^2 \frac{\left(1-2k_\beta ^{-2} \left(k_x^2 + k_y^2\right)\right)}{2k_z^{(\alpha )}}, \end{equation} (A4)  \begin{equation} \widetilde{F}_{v_x}^{(\nabla \times v)_x } = i \frac{k_x k_y}{2k_z^{(\beta )}}, \end{equation} (A5)  \begin{equation} \widetilde{F}_{v_y}^{(\nabla \times v)_x } = i \frac{k_\beta ^2 - 2k_y^2 - k_x^2}{2k_z^{(\beta )}}, \end{equation} (A6)  \begin{equation} \widetilde{F}_{v_x}^{(\nabla \times v)_y } = i \frac{k_\beta ^2 - 2k_x^2 - k_y^2}{2k_z^{(\beta )}}, \end{equation} (A7)  \begin{equation} \widetilde{F}_{v_y}^{(\nabla \times v)_y } = \widetilde{F}_{v_x}^{(\nabla \times v)_x }. \end{equation} (A8)Similarly as with the up/down separation technique we are seeking expressions that require local measurements of the wavefield only. The Taylor expansion of eqs (A4)–(A8) up to first-order spatial gradients yields (Robertsson & Curtis 2002):   \begin{equation} \widetilde{F}_{v_z}^{(\nabla \cdot v) } \approx i\frac{k_\alpha }{2}, \end{equation} (A9)  \begin{equation} \widetilde{F}_{v_x}^{(\nabla \times v)_x } \approx 0, \end{equation} (A10)  \begin{equation} \widetilde{F}_{v_y}^{(\nabla \times v)_x } \approx i\frac{k_\beta }{2}, \end{equation} (A11)  \begin{equation} \widetilde{F}_{v_x}^{(\nabla \times v)_y } \approx i\frac{k_\beta }{2}, \end{equation} (A12)  \begin{equation} \widetilde{F}_{v_y}^{(\nabla \times v)_y } \approx 0. \end{equation} (A13)Inserting eqs (A9)–(A13) into eqs (A1)–(A3) and transforming back to the space–frequency domain results in the following expressions for the upgoing P- and S-wave potentials (Robertsson & Curtis 2002):   \begin{equation} (\nabla \cdot v^U) \approx \frac{\beta ^2}{\alpha ^2} \frac{\partial v_x}{\partial x} + \frac{i\omega }{2 \alpha } v_z, \end{equation} (A14)  \begin{equation} (\nabla \times v^U)_x \approx - \frac{\partial v_z}{\partial y} + \frac{i\omega }{2 \beta } v_y, \end{equation} (A15)  \begin{equation} (\nabla \times v^U)_y \approx - \frac{\partial v_z}{\partial x} + \frac{i\omega }{2 \beta } v_x. \end{equation} (A16) APPENDIX B: P/S SEPARATION OF OBS DATA The upgoing P- and S-wave potential fields at the seabed are given by (Amundsen et al. 2000):   \begin{equation} (\nabla \cdot v^U) = -\frac{1}{2} p + \frac{i \omega \rho }{k_\beta ^{2}} \left(\frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} \right) + F_{v_z}^{(\nabla \cdot v) } \ast v_z, \end{equation} (B1)  \begin{eqnarray} (\nabla \times v^U)_x &=& F_p^{(\nabla \times v)_x} \ast p + F_{v_x}^{(\nabla \times v)_x } \ast v_x + F_{v_y}^{(\nabla \times v)_x } \ast v_y \nonumber\\ && -\, F_{v_z}^{(\nabla \times v)_x } \ast v_z, \end{eqnarray} (B2)  \begin{eqnarray} (\nabla \times v^U)_y &=& F_p^{(\nabla \times v)_y} \ast p + F_{v_x}^{(\nabla \times v)_y } \ast v_x + F_{v_y}^{(\nabla \times v)_y } \ast v_y \nonumber\\ && -\, F_{v_z}^{(\nabla \times v)_y } \ast v_z, \end{eqnarray} (B3)with filters in the wavenumber–frequency domain (Amundsen et al. 2000):   \begin{equation} \widetilde{F}_{v_z}^{(\nabla \cdot v)} = \frac{\omega \rho }{2k_z^{(\alpha )}} \left(1 - \frac{k_x^2 + k_y^2}{k_\beta ^2}\right), \end{equation} (B4)  \begin{equation} \widetilde{F}_p^{(\nabla \times v)_x} = \frac{k_y}{2k_z^{(\beta )}}, \end{equation} (B5)  \begin{equation} \widetilde{F}_{v_x}^{(\nabla \times v)_x} = \frac{\omega \rho }{2k_z^{(\beta )}} \left(\frac{k_x k_y}{k_\beta ^2}\right), \end{equation} (B6)  \begin{equation} \widetilde{F}_{v_y}^{(\nabla \times v)_x} = \frac{\omega \rho }{2k_z^{(\beta )}} \left(1 - \frac{k_x^2 + 2k_y^2}{k_\beta ^2}\right), \end{equation} (B7)  \begin{equation} \widetilde{F}_{v_z}^{(\nabla \times v)_x} = \omega \rho \frac{k_y}{k_\beta ^2}, \end{equation} (B8)  \begin{equation} \widetilde{F}_p^{(\nabla \times v)_y} = \frac{k_x}{2k_z^{(\beta )}}, \end{equation} (B9)  \begin{equation} \widetilde{F}_{v_x}^{(\nabla \times v)_y} = \frac{\omega \rho }{2k_z^{(\beta )}} \left(1 - \frac{2k_x^2 + k_y^2}{k_\beta ^2}\right), \end{equation} (B10)  \begin{equation} \widetilde{F}_{v_y}^{(\nabla \times v)_y} = \frac{\omega \rho }{2k_z^{(\beta )}} \left(\frac{k_x k_y}{k_\beta ^2}\right), \end{equation} (B11)  \begin{equation} \widetilde{F}_{v_z}^{(\nabla \times v)_y} = \omega \rho \frac{k_x}{k_\beta ^2}. \end{equation} (B12)The Taylor approximations to eqs (B4)–(B12) up to first-order spatial gradients yields:   \begin{equation} \widetilde{F}_{v_z}^{(\nabla \cdot v)} \approx \frac{\omega \rho }{2 k_\alpha }, \end{equation} (B13)  \begin{equation} \widetilde{F}_p^{(\nabla \times v)_x} \approx \frac{k_y}{2 k_\beta }, \end{equation} (B14)  \begin{equation} \widetilde{F}_{v_x}^{(\nabla \times v)_x} \approx 0, \end{equation} (B15)  \begin{equation} \widetilde{F}_{v_y}^{(\nabla \times v)_x} \approx \frac{\omega \rho }{2 k_\beta }, \end{equation} (B16)  \begin{equation} \widetilde{F}_{v_z}^{(\nabla \times v)_x} \approx \omega \rho \frac{k_y}{k_\beta ^2}, \end{equation} (B17)  \begin{equation} \widetilde{F}_p^{(\nabla \times v)_y} \approx \frac{k_x}{2 k_\beta }, \end{equation} (B18)  \begin{equation} \widetilde{F}_{v_x}^{(\nabla \times v)_y} \approx \frac{\omega \rho }{2 k_\beta }, \end{equation} (B19)  \begin{equation} \widetilde{F}_{v_y}^{(\nabla \times v)_y} \approx 0, \end{equation} (B20)  \begin{equation} \widetilde{F}_{v_z}^{(\nabla \times v)_y} \approx \omega \rho \frac{k_x}{k_\beta ^2}. \end{equation} (B21)The spatial wavefield gradient-based upgoing P- and S-wave potentials are obtained by inserting eqs (B13)–(B21) into eqs (B1)–(B3) and transforming back to the space–frequency domain:   \begin{equation} (\nabla \cdot v^U) \approx -\frac{1}{2}p + \frac{1}{i\omega } \rho \beta ^2 \left(\frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y}\right) + \frac{\alpha \rho }{2} v_z, \end{equation} (B22)  \begin{equation} (\nabla \times v^U)_x \approx \frac{\beta }{2 i \omega } \frac{\partial p}{\partial y} + \frac{\rho \beta }{2} v_y - \frac{\rho \beta ^2}{i \omega } \frac{\partial v_z}{\partial y}, \end{equation} (B23)  \begin{equation} (\nabla \times v^U)_y \approx \frac{\beta }{2 i \omega } \frac{\partial p}{\partial x} + \frac{\rho \beta }{2} v_x - \frac{\rho \beta ^2}{i \omega } \frac{\partial v_z}{\partial x}. \end{equation} (B24) © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

Journal

Geophysical Journal InternationalOxford University Press

Published: Mar 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Unlimited reading

Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere.

Stay up to date

Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates.

Organize your research

It’s easy to organize your research with our built-in tools.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

Monthly Plan

  • Read unlimited articles
  • Personalized recommendations
  • No expiration
  • Print 20 pages per month
  • 20% off on PDF purchases
  • Organize your research
  • Get updates on your journals and topic searches

$49/month

Start Free Trial

14-day Free Trial

Best Deal — 39% off

Annual Plan

  • All the features of the Professional Plan, but for 39% off!
  • Billed annually
  • No expiration
  • For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588

$360/year

billed annually
Start Free Trial

14-day Free Trial