# Parsimonious surface wave interferometry

Parsimonious surface wave interferometry Summary To decrease the recording time of a 2-D seismic survey from a few days to one hour or less, we present a parsimonious surface wave interferometry method. Interferometry allows for the creation of a large number of virtual shot gathers from just two reciprocal shot gathers by crosscoherence of trace pairs. Then, the virtual surface waves can be inverted for the S-wave velocity model by wave-equation dispersion inversion (WD). Synthetic and field data tests suggest that parsimonious WD (PWD) gives S-velocity tomograms that are comparable to those obtained from a conventional survey with a shot at each receiver. The limitation of PWD is that the virtual data lose some information so that the resolution of the S-velocity tomogram can be modestly lower than that of the S-velocity tomogram inverted from a conventional survey. Waveform inversion, Computational seismology, Seismic interferometry, Surface waves and free oscillations INTRODUCTION Surface wave inversion is an established method for estimating S-wave velocity models in near-surface, earthquake, and engineering seismology (Xia et al. 1999; Shapiro et al. 2005; Lin et al. 2008; Socco et al. 2010; Yilmaz 2015; Li et al. 2017a). For conventional surface wave surveys along a line of geophones, receivers are deployed along the survey line and common shot gathers (CSGs) are recorded at selected inline source positions. To significantly reduce recording time, we recall the concept of a parsimonious refraction survey where (Fu et al. 2017; Hanafy et al. 2017) reciprocal sources are located only at the beginning and the end locations of the receiver line as shown in Fig. 1. The two reciprocal shot gathers are used to interferometrically generate many virtual shot gathers with virtual refractions (Wapenaar et al. 2010; Schuster 2016). For refraction arrivals, we define this procedure as parsimonious refraction interferometry because it uses a stationary-phase principle to decompose the long-offset refraction data into virtual refraction data with short source-receiver offsets. That is, only two reciprocal shot gathers with 2N traces can be used to generate O(N2) virtual traces where there are N shot gathers with a virtual source at each of the N geophones along the survey line. These virtual shot gathers can be inverted for the P-wave velocity model by traveltime tomography. Figure 1. View largeDownload slide Velocity model where the reciprocal sources are at s and $$\bar{{\bf s}}$$, and the solid white ray paths represent the surface waves excited by the actual source at s. For D(ω)sg and $${D(\omega )_{s\bar{g}}}$$, the phase in D(ω)sg cancels that in $$D(\omega )_{s\bar{g}}$$ along the common ray paths marked by the shortest white line. The resulting phase is associated with the dashed yellow ray path associated with the virtual surface wave trace $$R({\bf s}, \omega )_{g\bar{g}}$$ excited by the virtual source (blue star) at g and recorded by the geophone at $$\bar{\bf g}$$. Figure 1. View largeDownload slide Velocity model where the reciprocal sources are at s and $$\bar{{\bf s}}$$, and the solid white ray paths represent the surface waves excited by the actual source at s. For D(ω)sg and $${D(\omega )_{s\bar{g}}}$$, the phase in D(ω)sg cancels that in $$D(\omega )_{s\bar{g}}$$ along the common ray paths marked by the shortest white line. The resulting phase is associated with the dashed yellow ray path associated with the virtual surface wave trace $$R({\bf s}, \omega )_{g\bar{g}}$$ excited by the virtual source (blue star) at g and recorded by the geophone at $$\bar{\bf g}$$. In a typical exploration survey, the surface waves are considered to be noise so the goal is to remove the surface waves (Dong et al. 2006; Halliday & Curtis 2008; Xue et al. 2009). In contrast, for many engineering surveys, the surface waves are considered to be useful signal so the goal is to invert the surface wave dispersion curves (Park et al. 1998, 1999; Hayashi & Suzuki 2004; Socco et al. 2010; Gouédard et al. 2012; Boiero et al. 2013; Yilmaz 2015). We now extend parsimonious refraction interferometry to surface waves, where trace pairs from the two reciprocal shot gathers are correlated to give approximately N virtual shot gathers that contain virtual surface waves, where N is the number of geophones between the two shots. These virtual surface waves are inverted by wave-equation dispersion inversion (WD; Li & Schuster 2016; Li et al. 2017a) to give the S-velocity tomogram. We denote this combination of interferometry and inversion methods for surface waves as parsimonious WD (PWD). The next section describes the theory of PWD. Then, we present the results of applying PWD to both synthetic and seismic field data. The final section presents the summary and conclusions. THEORY OF PARSIMONIOUS WAVE-EQUATION INVERSION Assume a point source at each end of the survey line that overlies the irregularly layered medium in Fig. 1, where surface waves propagate along the free surface. Here, we assume an isotropic medium with no attenuation. There are velocity variations in the upper medium and N evenly spaced receivers (black triangles) on the recording surface between the two actual sources; these two actual sources at the ends of the survey line are denoted as reciprocal sources (red stars). The actual trace D(ω)sg in the frequency domain is excited by the source at s and recorded at g by a vertical particle-velocity geophone; the point source is a vertical-displacement at s oscillating at frequency ω. The key idea is to crosscorrelate a pair of traces at different receivers $$\bar{{{\bf g}}}$$ and g to obtain a virtual trace at g excited by a virtual source (blue star) at $$\bar{{{\bf g}}}$$. If this is done for all geophone pairs then O(N) virtual shot gathers will be generated, with a virtual shot at each geophone and N virtual traces for each shot gather. The leftmost (rightmost) reciprocal shot generates the right-going (left-going) surface waves. In our examples we assume that the fundamental Rayleigh mode is isolated and its dispersion curve is inverted for the subsurface shear-velocity model. Parsimonious surface-wave interferometry We now present the theory of surface wave interferometry in the context of a parsimonious experiment where there is a line of vertical-component geophones and there are just two vertical-displacement shots, one at each end of the line. We assume an isotropic layered model with a shear wave velocity smoothly varying along the lateral direction and there is no attenuation in the model. Haney & Nakahara (2014, 2016) derived expressions for surface wave Green’s tensors that include near-field behaviour. Here, in order to describe the theory of surface wave interferometry in a simple way, the Rayleigh-wave approximation Green’s function (Snieder et al. 2002) can be described by   \begin{eqnarray} G({\bf g}|{\bf s})&=&\sum _\eta A_\eta \frac{e^{i\theta (k_\eta , {\bf s},{\bf g})+i\pi / 4}}{\sqrt{0.5 \pi k_\eta q_\eta }}, \end{eqnarray} (1)where the receivers at g are in the far-field region of the harmonic point source at s. The summation is over different propagation modes, the horizontal free surface is at z = 0, and kη is the wavenumber of the integer-valued mode denoted by η. For simplicity of exposition we will assume the phase term θ(kη, s, g) = kηX and the geometrical spreading factor to be qη = X, where X = |s − g| is the horizontal distance between the vertical-displacement point source at s = (xs, ys, 0) and the vertical-component geophone at g = (xg, yg, 0). The term Aη accounts for the source amplitude and radiation pattern for a trace at g excited by a point source at s. In general, the phase θ(kη, s, g) and geometrical spreading qη terms are computed by ray tracing for a laterally varying velocity (Tanimoto 1990). We now assume that the surface waves are filtered so that there is only a single mode η = 1 propagating according to equation 1. For a linear recording line with an endline shot at s, the spectral product $$G(\bar{{{\bf g}}}|{\bf s})G({\bf g}|{\bf s})^*$$ of two traces, one recorded at the near-shot position g and the other at the far-shot position $$\bar{{{\bf g}}}$$ along the recording line, becomes   \begin{eqnarray} \Phi (\bar{{{\bf g}}}, {\bf g})_{{\bf s}}=G(\bar{{{\bf g}}}|{\bf s})G({\bf g}|{\bf s})^*&=&A_1 \skew7\bar{A}_1 \frac{e^{ik_1 |{\bf g}-\bar{{{\bf g}}}|}}{\sqrt{0.5 \pi k_1\overline{X}}\sqrt{0.5 \pi k_1 X}}, \end{eqnarray} (2)where X = |s − g| ($$\overline{X}=|{\bf s}-\bar{{{\bf g}}}|$$), and $$\skew7\bar{A}_1$$ (A1) describe the source and radiation characteristics for the trace recorded at $$\bar{{{\bf g}}}$$ (g). The arrival in the virtual trace $$\Phi (\bar{{{\bf g}}}, {\bf g})_{\bf s}$$ has the phase characteristics of a Rayleigh wave (with mode η = 1) excited by a vertical ground displacement at g and recorded by a vertical-component geophone at $$\bar{{{\bf g}}}$$. Therefore, if the recorded traces G(g|s) and $$G(\bar{{{\bf g}}}|{\bf s})$$ are filtered so that they only contain a single mode of the one-way propagating Rayleigh wave, then we conclude that the phase-velocity curve C(ω)η computed from the virtual shot gather $$\Phi (\bar{{{\bf g}}}, {\bf g})_{\bf s}$$ will be identical to that of an actual shot gather $$G(\bar{{{\bf g}}}|{\bf g})$$. The magnitudes of the virtual and actual spectra will differ in the C(ω)η − ω domain but the trajectory of the two curves will be the same. Only an accurate estimate of the phase is needed to give an accurate trajectory of C(ω)η, so that WD can provide an accurate shear-velocity tomogram. If the filtering is incomplete so that there are two modes in the filtered data   \begin{eqnarray} G({\bf g}|{\bf s})&=& A_1 \frac{e^{ik_1X+i\pi /4}}{\sqrt{0.5 \pi k_1 X}}+ A_2 \frac{e^{ik_2X+i\pi /4}}{\sqrt{0.5 \pi k_2 X}} , \end{eqnarray} (3)then the correlated traces yield the virtual trace   \begin{eqnarray} &&{G(\bar{{{\bf g}}}|{{\bf s}})G({\bf g}|{\bf s})^*= \overbrace{A_1 \skew4\bar{A}_1 \frac{e^{ik_1 |{\bf g}\!-\!\bar{{{\bf g}}}|}}{\sqrt{0.5 \pi k_1\overline{X}}\sqrt{0.5 \pi k_1 X}}}^{\text{mode} 1}} \nonumber\\ &&+ \overbrace{A_2 \skew4\bar{A}_2 \frac{e^{ik_2 |{\bf g}-\bar{{{\bf g}}}|}}{\sqrt{0.5 \pi k_2\overline{X}}\sqrt{0.5 \pi k_2 X}}}^{\text{mode} 2}\nonumber \\ &&+ \overbrace{A_1 \skew4\bar{A}_2 \frac{e^{ik_1 |{\bf s}\!-\!\bar{{{\bf g}}}|-ik_2|{\bf s}\!-\! {\bf g}|}}{\sqrt{0.5 \pi k_1 \overline{X}}\sqrt{0.5 \pi k_2 X}} \!+\! \skew4\bar{A}_1 A_2 \frac{e^{ik_2 |{\bf s}-\bar{{{\bf g}}}|-ik_1|{\bf s}- {\bf g}|}}{\sqrt{0.5 \pi k_1 X}\sqrt{0.5 \pi k_2 \overline{X}}}}^{\text{crosstalk}}.\nonumber\\ \end{eqnarray} (4)In this case, the dispersion curves associated with the virtual shot gather $$G(\bar{{{\bf g}}}|{\bf s})G({{\bf g}}|{\bf s})^*$$ will contain the correct C(ω)η for modes η = 1 and η = 2, but they will also contain the noisy crosstalk modes. Such noise will hinder an accurate estimation of C(ω)η. We conclude that, prior to crosscorrelation, it is important to filter the shot gather so that it only contains one mode propagating in one direction. In practice, this mode is the fundamental mode of the Rayleigh wave because it is the strongest arrival in the records. Workflow of parsimonious wave-equation inversion The workflow for PWD can be summarized by the following three steps. For one shot gather, crosscoherence of the traces at receivers A = g and B = $$\bar{{{\bf g}}}$$ generates a virtual seismogram where the receiver at the A = g position is now the new virtual source position shown in Fig. 2. Virtual surface waves can be reconstructed at each receiver location by summation over sources and by using the crosscoherence of the normalized traces recorded by the receivers at g and $$\bar{{{\bf g}}}$$. In the frequency domain this operation is represented as (Xue et al. 2009; Nakata et al. 2013):   \begin{eqnarray} R({\bf s}, \omega )_{g\bar{g}}=\sum _{s} \frac{D(\omega )_{\bar{g}s}D(\omega )^*_{gs}}{|D(\omega )_{gs}||D(\omega )_{\bar{g}s}|+\lambda }, \end{eqnarray} (5)where $$R({\bf s}, \omega )_{g\bar{g}}$$ represents the time-symmetric virtual trace in the frequency domain with right-going waves recorded at $$\bar{{{\bf g}}}$$ for the virtual source at g, and λ is the damping parameter that is from 1 per cent to 10 per cent of the signal amplitude (Nakata et al. 2013). Eq. (5) sums over each source deployed at the left end of the line, and for our examples we find that it is sufficient to use just one source at the end of the line. The recorded data are FK filtered so only right-going surface waves are contained in D(ω)gs and $$D(\omega )_{\bar{g}s}$$ if s is at the left end of the line, and if s is at the right end of the line then the FK filter is used to eliminate left-going waves in D(ω)gs and $$D(\omega )_{\bar{g}s}$$. Here, g is several wavelengths from $$\bar{{{\bf g}}}$$ in order for the surface waves to develop. Halliday & Curtis (2008) suggest that the fundamental mode of the Rayleigh waves can be reliably recovered by interferometry for inline sources on the surface if there is no crosstalk from other modes. Therefore, the PWD method will concentrate on recovering and inverting the fundamental modes of Rayleigh waves for the N virtual shot gathers created from just two reciprocal shot gathers with their shot locations at the two ends of the survey line. To mitigate the crosstalk between different modes, the reciprocal shot gather can be transformed into the phase-velocity and frequency domains by a frequency-domain fast Fourier transform (FFT) and a linear Radon transform (LRT). Then, the fundamental Rayleigh wave can be windowed to eliminate higher-order modes, and the result will be a dispersion curve that is devoid of higher-mode energy. This windowed result can be inverse transformed back into the space–time domain to give a shot gather that only contains the fundamental mode of the Rayleigh waves (Luo et al. 2009; Hu et al. 2016). The filtered reciprocal shot gather can then be used to generate the virtual shot gathers by crosscoherence of trace pairs and summation over the two reciprocal sources. Use an LRT method to transform the filtered virtual shot gathers $$R({\bf s}, \omega )_{g\bar{g}}$$ into the C(g, ω) − ω domain, where C(g, ω) is the phase velocity for the virtual shot gather with the virtual shot at g. The fundamental Rayleigh dispersion curve C(g, ω)obs is picked for each shot gather with the virtual source at g. These curves are used to generate the associated wavenumber dispersion curve by the formula κ(g, ω)obs = ω/C(g, ω)obs. The dispersion curves are used to give the misfit function   \begin{eqnarray} \epsilon =\frac{1}{2} \sum _{\omega } \sum _g (\kappa ({\bf g}, \omega )^{\text{pred}} - \kappa ({\bf g}, \omega )^{\text{obs}})^2, \end{eqnarray} (6)where κ(g, ω)pred is the predicted surface wave dispersion curve computed from a 2-D finite-difference solution to the elastic wave equation with a vertical point source at g. The S-wave velocity is iteratively updated by the iterative steepest-descent formula of Li & Schuster (2016) and Li et al. (2017a):   \begin{eqnarray} s({\bf x})^{(k+1)} &=& s({\bf x})^{(k)} - \alpha \sum _{\omega } \sum _g \Delta \kappa ({\bf g},\omega ) \frac{{\partial } \kappa ({\bf g},\omega )^{\text{pred}} }{{\partial } s({\bf x})}, \end{eqnarray} (7)where Δκ(g, ω) = κ(g, ω)pred − κ(g, ω)obs is the residual at the kth iteration, and s(x)(k) is the shear slowness at the kth iteration. Here, α is the step length by any backtracking line-search method (Nocedal & Wright 1999) and the k superscript denotes the kth iteration. For pedagogical simplicity, we assume a single shot gather $$\hat{D}({{\bf g}},\omega )$$ for a virtual or real source at s and geophones at g, and the notation for the source location is silent. The misfit function will include an additional summation over different shot gathers if more than one virtual shot gather is used. Mathematical details for deriving the Fréchet derivative $$\frac{{\partial } \kappa (\omega )}{{\partial } s({\bf x})}$$ are given in Li et al. (2017a), where   \begin{eqnarray} &&{\frac{{\partial } \kappa }{{\partial } s(\mathbf {x})} = \frac{4s(\mathbf {x})\rho ({\bf x})}{A} \text{Real} \Big\lbrace \overbrace{\frac{{\partial } D({\bf x},\omega )_{x}}{{\partial } x}}^{\text{source}=f({{\bf x}},\mathbf {s},\omega )_1}}\nonumber \\ && \int \overbrace{\hat{D}({{\bf g}},\omega )_{\text{obs}}^*\left(\frac{{\partial } G({\bf g}|{\bf x})_{zz}}{{\partial } z}\right)\text{d}x_{g}}^{\text{backprojected data}=B(\mathbf {x},\mathbf {s},\omega )_{1}^*} +\overbrace{\frac{{\partial } D({\bf x},\omega )_{z}}{{\partial } z}}^{\text{source}=f({{\bf x}},\mathbf {s},\omega )_2}\nonumber \\ && \int \overbrace{\hat{D}({{\bf g}},\omega )_{\text{obs}}^*\frac{{\partial } G({\bf g}|{\bf x})_{xz}}{{\partial } x}\text{d}x_{g}}^{\text{backprojected data}=B(\mathbf {x},\mathbf {s},\omega )_{2}^*}\nonumber \\ && -\frac{1}{2} \overbrace{\left(\frac{{\partial } D({\bf x},\omega )_{z}}{{\partial } x}+\frac{{\partial } D({\bf x},\omega )_{x}}{{\partial } z}\right)}^{\text{source}={f(\mathbf {x},\mathbf {s},\omega )_3}}\nonumber \\ && \int \overbrace{\hat{D}({{\bf g}},\omega )_{\text{obs}}^*\left(\frac{{\partial } G({\bf g}|{\bf x})_{xz}}{{\partial } z} +\frac{{\partial } G({\bf g}|{\bf x})_{zz}}{{\partial } x}\right)}^{\text{backprojected data}={ B(\mathbf {x},\mathbf {s},\omega )_{3}^*}}\text{d}x_{g} \Big\rbrace. \end{eqnarray} (8)Here, $$\hat{D}({{\bf g}},\omega )_{\text{obs}}^*$$ is the weighted conjugated data for the vertical component displacement recorded at g, G(g|x)ij is the Green’s tensor for the ith component of the displacement field at (x, z) for a jth-component source at (xg, 0), s(x) is the shear-slowness model for the background medium, and A is the normalization term. Here, f(x, s, ω)i for i ∈ (1, 2, 3) is the source field at $$\bf x$$ that originates at $$\bf s$$, and B(x, s, ω)i represents the backprojected field at $$\bf x$$ for a source at s; it is assumed that this field is recorded at the surface and backpropagated by the adjoint of the forward modelling operator. The interpretation of equations (7) and (8) is that the traces and source wavelet are weighted by terms proportional to the wavenumber residual Δκ(ω) and backpropagated into the medium to update the slowness of the S-velocity model. Fig. 3 summarizes the workflow for the PWD method. Figure 2. View largeDownload slide Crosscoherence of two traces at different source–receiver offsets yields a trace (on the right) with a virtual shot at the A receiver. Figure 2. View largeDownload slide Crosscoherence of two traces at different source–receiver offsets yields a trace (on the right) with a virtual shot at the A receiver. Figure 3. View largeDownload slide The workflow for the parsimonious surface wave WD method, where illumination compensation is applied to the gradient as a pre-conditioner. Here, RT denotes the Radon transform and the comparison between the conventional data and virtual data is only used for assessing the accuracy of PWD with synthetic data. In practice, the conventional data for shots at all geophones is usually not available. Figure 3. View largeDownload slide The workflow for the parsimonious surface wave WD method, where illumination compensation is applied to the gradient as a pre-conditioner. Here, RT denotes the Radon transform and the comparison between the conventional data and virtual data is only used for assessing the accuracy of PWD with synthetic data. In practice, the conventional data for shots at all geophones is usually not available. PARSIMONIOUS WAVE-EQUATION DISPERSION INVERSION TESTS In this section, the PWD method is used to invert for the S-velocity model from three different sets of synthetic data and one field data set recorded near the Gulf of Aqaba in Saudi Arabia. We also compared the PWD tomogram (inverted from N virtual shot gathers) with the conventional WD tomogram (inverted from N actual shot gathers). In each example the observed surface wave data are generated using time-domain finite-difference solutions to the 2-D elastic wave equation (Levander 1988). Stair-step velocity model The geometry of a conventional seismic survey and a parsimonious survey with a stair-step S-wave velocity model are shown in Fig. 4. The 2-D model size is 20 × 120 m2, the S-wave velocities in the top and bottom layers are 500 ms−1 and 800 ms−1, respectively, the P-wave velocity Vp is constrained to be $$V_p=\sqrt{3} V_s$$ and the source wavelet is a Ricker wavelet with a peak frequency of 30 Hz. The conventional survey consists of 60 shots (red stars) shooting into 60 receivers (white triangles) evenly deployed at 2 m intervals with a 2 m shot interval. To test the parsimonious method, we used two reciprocal sources located at the beginning and the end of the survey line (black stars) while the data are recorded at all receiver locations. Figure 4. View largeDownload slide The acquisition geometry of the conventional and parsimonious surveys. Figure 4. View largeDownload slide The acquisition geometry of the conventional and parsimonious surveys. We computed 60 virtual shot gathers from two reciprocal shot gathers by the cross-coherence method. Fig. 5(a) displays one virtual shot gather (red traces), which almost matches the corresponding shot gather in the conventional survey data (black traces), except their amplitudes have some differences. These amplitude differences will not significantly affect the result of S-wave inversion because the WD method inverts the residuals associated with the phase velocity (Li & Hanafy 2016; Li & Schuster 2016; Li et al. 2017a). Figure 5. View largeDownload slide (a) Actual shot gather (red) and virtual shot gather (black) computed by crosscoherence, and (b) the spectra of single traces. Figure 5. View largeDownload slide (a) Actual shot gather (red) and virtual shot gather (black) computed by crosscoherence, and (b) the spectra of single traces. Fig. 5(b) shows the spectra computed from the actual traces and the virtual traces computed by the crosscoherence method, where the useful frequencies in the original traces (red spectrum) are from 10 to 80 Hz while those in the virtual traces computed by crosscorrelation (black spectrum) are only from 15 to 75 Hz. However, in a noisy environment the signals above a critical frequency are too noisy so a high-cut filter should be applied to the virtual traces computed by crosscoherence. The high-cut frequency is the critical frequency identified in the data. A phase-velocity model (see Fig. 6) is computed from both the actual shot gathers and the virtual shot gathers using the LRT. The phase-velocity model computed from the 60 actual shot gathers (Fig. 6a) closely resembles the model inverted from the 60 virtual shot gathers (i.e. the parsimonious data) (Fig. 6b). The PWD and conventional WD tomograms are shown in Figs 7(c) and (d). Both tomograms are mostly consistent with the true S-wave velocity model in Fig. 7(a). For comparison, Fig. 7(b) shows the inaccurate S-velocity tomogram inverted from just two reciprocal shot gathers (without using the crosscoherence operation to generate the virtual traces). The normalized data residual plotted against iteration number in Fig. 8 shows that the PWD method has a convergence rate similar to that of the conventional WD method, but the conventional WD tomogram has a smaller final residual. Figure 6. View largeDownload slide Phase-velocity models computed from the (a) conventional data and (b) virtual data. Figure 6. View largeDownload slide Phase-velocity models computed from the (a) conventional data and (b) virtual data. Figure 7. View largeDownload slide (a) True S-wave velocity model; (b) WD tomogram inverted from two reciprocal shot gathers; tomograms computed by the (c) conventional WD and (d) parsimonious WD methods. Figure 7. View largeDownload slide (a) True S-wave velocity model; (b) WD tomogram inverted from two reciprocal shot gathers; tomograms computed by the (c) conventional WD and (d) parsimonious WD methods. Figure 8. View largeDownload slide RMS residual plotted against iteration number for the two WD tomograms in Figs 7(c) and (d), where the red and black lines are for the conventional and parsimonious residuals, respectively. Figure 8. View largeDownload slide RMS residual plotted against iteration number for the two WD tomograms in Figs 7(c) and (d), where the red and black lines are for the conventional and parsimonious residuals, respectively. Local-anomalies model The next test is for a velocity model with local-velocity anomalies shown in Fig. 9(a). This model is based on geological observations in a shallow open-pit mine with mineral deposits. The 2-D model size is 20 × 170 m2 and the source wavelet is a Ricker wavelet with a peak frequency of 30 Hz. We follow the same workflow shown in Fig. 3, where the input data for conventional WD consists of 40 shot gathers with 85 traces and receivers located every 2 m, with the shots evenly distributed every 4 m on the ground surface. The parsimonious data consist of two CSGs generated by the shots at the two ends of the survey line. These two reciprocal CSGs are used to generate 40 virtual shots gathers by crosscoherence interferometry. Figure 9. View largeDownload slide (a) True S-wave velocity model, (b) initial S-wave velocity model, (c) conventional WD tomogram, and (d) parsimonious WD tomogram. Figure 9. View largeDownload slide (a) True S-wave velocity model, (b) initial S-wave velocity model, (c) conventional WD tomogram, and (d) parsimonious WD tomogram. The S-velocity tomograms inverted by the conventional WD and PWD methods are shown in Figs 9(c) and (d). Both tomograms are in reasonable agreement with the true model, especially with the high-velocity anomaly delineated by the white dashed lines. After 31 iterations (Fig. 10), the normalized data misfit residuals of PWD and conventional WD decreased to the approximate values of 0.2 and 0.1, respectively. Figure 10. View largeDownload slide RMS residual as a function of iteration number for the conventional (red) and parsimonious (black) data computed for the Fig. 10a model. Figure 10. View largeDownload slide RMS residual as a function of iteration number for the conventional (red) and parsimonious (black) data computed for the Fig. 10a model. Aqaba field data test PWD is now tested on a field data set recorded next to the Gulf of Aqaba (see Fig. 11) in Saudi Arabia (Hanafy et al. 2014). For the conventional data survey, a total of 120 common shot gathers (CSGs) are recorded, and each shot gather has 120 traces at equal shot and receiver intervals of 2.5 m. The total length of the profile is 297.5 m and a 90 kg weight drop is used as the seismic source, with 10 to 15 stacks at each shot location. We select three shot gathers (shot locations in the middle and at the two ends of the survey line) to generate 120 virtual shot gathers by crosscoherence. Fig. 12 shows that the surface waves in the virtual shot gather compare well with those in the actual shot gather for the same source location. Figure 11. View largeDownload slide (a) Google map shows the survey location of the Aqaba experiment, and (b) photo shows an earthquake rupture (indicated by the red dashed line) on the surface. The black arrow indicates the location of the seismic survey line. Figure 11. View largeDownload slide (a) Google map shows the survey location of the Aqaba experiment, and (b) photo shows an earthquake rupture (indicated by the red dashed line) on the surface. The black arrow indicates the location of the seismic survey line. Figure 12. View largeDownload slide Conventional (black) and virtual (red) trace comparisons at different shot locations. Figure 12. View largeDownload slide Conventional (black) and virtual (red) trace comparisons at different shot locations. Following the workflow shown in Fig. 3, the conventional and PWD tomograms are computed and shown in Figs 13(b) and (c). The S-velocity tomograms generally agree with one another for the long-wavelength features. According to the P-velocity tomogram in Fig. 13(a), the actual fault position is at x = 150 m, which is at the same position indicated in the PWD tomogram. It is evident that PWD provides an image similar to that computed from the data recorded by a conventional seismic survey with 120 shots. However, the conventional survey required more than an order-of-magnitude increase in survey time compared to the parsimonious survey. Figure 13. View largeDownload slide (a) P-velocity inverted from first-arrival traveltimes tomogram, S-velocity tomograms computed from the (b) conventional and (c) parsimonious surface wave data. Figure 13. View largeDownload slide (a) P-velocity inverted from first-arrival traveltimes tomogram, S-velocity tomograms computed from the (b) conventional and (c) parsimonious surface wave data. Time-lapse monitoring test Seismic time-lapse techniques are often applied to large areas, such as hydrocarbon reservoirs (Landro 2001) or geologic storage reservoirs (Carcione et al. 2006; Picotti et al. 2012). For near-surface geophysics, the potential application of seismic methods for tracking time-variant processes is summarized by Jefferson et al. (1998). As an example, seismic surface wave techniques have been proposed for the detection of varying levels of the water table of shallow aquifers, heavy-oil in oil fields (Pasquet et al. 2015; Dubos-Sallee et al. 2015; Ikeda et al. 2016; Pasquet et al. 2016), or monitoring climate effects on earthworks (Bergamo et al. 2016a,b). A significant limitation of a small-scale time-lapse seismic survey, for example, 240 shot gathers, is that it takes at least 16 working hours for completion if the source excitation at each of the 240 source points requires more than 5 min of effort. This is a major problem if the geological properties of the subsurface target, such as in a fluid-injection experiment, rapidly vary in time. In comparison, acquisition of parsimonious data is fast and efficient because it only requires several shot gathers where the source can be far from the fluid-injection site. In this section, we will use synthetic data to evaluate the feasibility of using time-lapse parsimonious WD to monitor the subsurface fluid distribution after dumping a viscous fluid on to the surface of the recording survey. In this case we assume the viscous fluid increases the S-velocity of the porous rock. Synthetic model test The time-lapse fluid model shown in Fig. 14(a) is a two-layered model where the fluid is leaking from the surface into the model. A 2-D time-domain finite-difference method is used to compute the vertical-component elastic seismograms. The model size is 25 × 120 m2 where the grid spacing and time sampling are 1m and 0.25 ms, respectively, and the source time-history is a Ricker wavelet peaked at 30 Hz. Here, we assume the fluid increases the S-velocity of the subsurface while in practice we find that the S-velocity can be lowered after injecting water into porous sand. Figure 14. View largeDownload slide S -wave velocity models corresponding to the (a) time-lapse fluid-injection tests, (b) conventional WD tomogram, (c) parsimonious WD tomogram and (d) velocity changes of the parsimonious S-velocity tomograms before and after fluid injection at different times. Figure 14. View largeDownload slide S -wave velocity models corresponding to the (a) time-lapse fluid-injection tests, (b) conventional WD tomogram, (c) parsimonious WD tomogram and (d) velocity changes of the parsimonious S-velocity tomograms before and after fluid injection at different times. First, we generate the conventional-acquisition shot gathers which consist of 30 CSGs with sources located every 4 m along the recording line. Each shot is recorded by 60 receivers spaced at 2 m intervals. The conventional WD S-velocity tomogram, shown in Fig. 14(b), is in good agreement with the true fluid-injection velocity model. Now, we test the PWD approach, where we only use two reciprocal shot gathers with shots located at x = 0 m and x = 120 m. In both shot gathers the data are recorded at 60 receivers spaced at 2 m intervals along the ground surface. These reciprocal shot gathers are used to generate 30 virtual shot gathers by crosscoherence of traces. The corresponding S-velocity tomograms are shown in Fig. 14(c), where the tomograms show a clear time-lapse velocity anomaly that is almost the same as seen in the conventional WD tomogram. However, the conventional WD tomogram has a modestly higher resolution than the parsimonious data tomogram. This lowered resolution is possibly caused by imperfect filtering so that mode coupling artefacts appear in the correlation. The difference between the parsimonious S-velocity tomograms computed from data recorded before and after fluid injection is shown in Fig. 14(d). This clearly reveals the velocity changes as a function of time. DISCUSSION A limitation of PWD is that it typically provides tomograms that have modestly lower resolution than the conventional tomograms. The partial cause is that the virtual data are created by correlating two traces with one another, so that the spectrum of the virtual source wavelet is the square of the wavelet spectrum. This tends to emphasize the central frequencies at the expense of downweighting the higher and lower frequencies, as illustrated in Fig. 5. Thus, crosscoherence or spectral deconvolution is required, but its effectiveness can be degraded by noise in the data and the choice of the damping parameter. In addition, mode mixing can be activated by crosscoherence, so ideally higher-order modes should be eliminated prior to crosscoherence of traces (Hu et al. 2016). Another possibility is to design matching filters that match the traces from the parsimonious data to the conventional data. These matching filters can then be applied to the time lapse virtual data to deconvolve the source wavelet. This is a subject of future research. There are two other limitations. The first one is that PWD theory assumes no attenuation in the data. In this case attenuation should be estimated from the data and, if possible, the data should be corrected for attenuation effects. The second limitation is that the accuracy of the estimated dispersion curve depends on the signal-to-ratio of the data and the window selection to isolate the surface waves of interest. Thus, sensitivity tests should be carried out to determine the window design that gives the most robust estimate of the dispersion curves. Similar to a multi-frequency FWI strategy (Masoni et al. 2016), we can utilize multi-offset dispersion curves to invert for the S-velocity at different depths. Shot-offset traces have high horizontal resolution at the expense of decreased accuracy at the lower frequencies, whereas wide-offset data have the opposite characteristics (Strobbia & Foti 2006). So, the combined short-to-wide offset data with gradually increasing frequency ranges can be used to optimize the accuracy of the reconstructed S-velocity model. Details are provided in Li et al. (2017b). CONCLUSIONS We presented the theory of PWD where a dense set of virtual surface wave data is interferometrically generated from a pair of reciprocal shot gathers, and the shear-velocity model is computed by WD. This method is denoted as PWD and assumes an isotropic elastic medium with no attenuation. We compare the S-velocity tomograms generated by PWD and conventional WD using synthetic and seismic field data sets. Our results show that the PWD can provide nearly the same S-velocity tomograms as conventional WD, but with slightly less resolution. Our synthetic-data results also suggest that time-lapse parsimonious surface wave monitoring can be a robust and reliable technique for monitoring small-scale and rapidly varying subsurface conditions in the subsurface. Possible applications of PWD are for environmental monitoring of water-injection experiments for hazard remediation, archaeological sites where sources are not permitted over sensitive areas, characterizing the processes in erupting geysers or volcanoes, and engineering surveys. The PWD methodology can be adapted to inversion of dispersion curves computed from Love waves and trapped P-waves in wave guides. ACKNOWLEDGEMENTS We thank the sponsors for supporting the Consortium of Subsurface Imaging and Fluid Modelling (CSIM). We also thank KAUST for providing funding by the CRG grant OCRF-2014-CRG3-2300. For computer time, this research used the resources of the IT Research Computing Group and the Supercomputing Laboratory at KAUST. REFERENCES Bergamo P.B., Dashwood S., Uhlemann R., Swift J.E., Chambers D., Gunn A., Donohue S., 2016a. Time-lapse monitoring of climate effects on earthworks using surface waves, Geophysics , 81, EN1– EN15. Google Scholar CrossRef Search ADS   Bergamo P.B., Dashwood S., Uhlemann R., Swift J.E., Chambers D., Gunn A., Donohue S., 2016b. Time-lapse monitoring of fluid-induced geophysical property variations within an unstable earthwork using P-wave refraction, Geophysics , 81, EN17– EN27. Google Scholar CrossRef Search ADS   Boiero D., Wiarda E., Vermeer P., 2013. Surface-and guided-wave inversion for near-surface modelling in land and shallow marine seismic data, Leading Edge , 32, 638– 646. Google Scholar CrossRef Search ADS   Carcione J.M., Picotti S., Gei D., Rossi G., 2006. Physics and seismic modelling for monitoring CO2 storage, Pure appl. Geophys. , 163, 175– 207. Google Scholar CrossRef Search ADS   Dong S., He R., Schuster G.T., 2006. Interferometric prediction and least squares subtraction of surface waves, SEG Technical Program Expanded Abstracts 2006 , Society of Exploration Geophysicists, pp. 2783– 2786. Dubos-Sallee N., Rasolofosaon P., Barbouteau S., Charlety J., Duret F., 2015. Time-lapse surface wave analysis using noise correlation on an onshore heavy-oil field, SEG Technical Program Expanded Abstracts 2015 , Society of Exploration Geophysicists, pp. 5479– 5484. Fu L., Hanafy S., Schuster G.T., 2017. Parsimonious wave equation refraction interferometry and inversion, Geophys. Prospect. , doi:10.1111/1365–2478.12488. Gouédard P., Yao H., Ernst F., van der Hilst R.D., 2012. Surface wave Eikonal tomography in heterogeneous media using exploration data, Geophys. J. Int. , 191, 781– 788. Google Scholar CrossRef Search ADS   Halliday D., Curtis A., 2008. Seismic interferometry, surface waves and source distribution, Geophys. J. Int. , 175, 1067– 1087. Google Scholar CrossRef Search ADS   Hanafy S.M., Jonsson S., Kger Y., 2014. Imaging normal faults in alluvial fans using geophysical techniques: field example from the coast of Gulf of Aqaba, Saudi Arabia, SEG Technical Program Expanded Abstracts 2014 , Society of Exploration Geophysicists, pp. 4670– 4674. Hanafy S.M., Li J., Schuster G., 2017. Time-lapse monitoring of subsurface fluid flow using parsimonious seismic interferometry, SAGEEP Technical Program Expanded Abstracts 2017 , SAGEEP, pp. 1– 5. Haney M.M., Nakahara H., 2014. Surface-wave Greens tensors in the near field, Bull. seism. Soc. Am. , 104, 1578– 1586. Google Scholar CrossRef Search ADS   Haney M.M., Nakahara H., 2016. Surface-wave Greens tensors in the near field, Bull. seism. Soc. Am. , 106, 816– 818. Google Scholar CrossRef Search ADS   Hayashi K., Suzuki H., 2004. CMP cross-correlation analysis of multi-channel surface-wave data, Explor. Geophys. , 35, 7– 13. Google Scholar CrossRef Search ADS   Hu Y., Cheng F., Bin M., Shen C., 2016. Waveform inversion of Rayleigh-wave in time domain for shallow shear-wave velocity based on the waveform separation, in AGU Annual Fall Meeting Abstracts , San Francisco, AGU, pp. 1– 5. Ikeda T., Tsuji T., Takanashi M., Kurosawa I., Nakatsukasa M., White D., Worth K., Roberts B., 2016. Time-lapse monitoring of shallow subsurface in the aquistore CO2 storage site from surface-wave analysis using a continuous and controlled seismic source, SEG Technical Program Expanded Abstracts 2016 , Society of Exploration Geophysicists, pp. 5479– 5484. Jefferson R.D., Steeples D.W., Black R.A., Carr T., 1998. Effects of soil-moisture content on shallow-seismic data, Geophysics , 63, 1357– 1362. Google Scholar CrossRef Search ADS   Landro M., 2001. Discrimination between pressure and fluid saturation changes from time-lapse seismic data, Geophysics , 66, 836– 844. Google Scholar CrossRef Search ADS   Levander A.R., 1988. Fourth-order finite-difference P–SV seismograms, Geophysics , 53, 1425– 1436. Google Scholar CrossRef Search ADS   Li J., Hanafy S., 2016. Skeletonized inversion of surface waves: active source versus controlled noise comparison, Interpretation , 3, 11– 19. Li J., Schuster G., 2016. Skeletonized wave equation of surface wave dispersion inversion, SEG Technical Program Expanded Abstracts 2016 , Society of Exploration Geophysicists, pp. 3630– 3635. Li J., Feng Z., Schuster G.T., 2017a. Wave-equation dispersion inversion, Geophys. J. Int. , 208, 1567– 1578. Google Scholar CrossRef Search ADS   Li J., Lin F., Schuster G.T., 2017b. Wave equation dispersion inversion with irregular topography, Geophys. J. Int. , in press. Lin F.-C., Moschetti M.P., Ritzwoller M.H., 2008. Surface wave tomography of the Western United States from ambient seismic noise: Rayleigh and Love wave phase velocity maps, Geophys. J. Int. , 173, 281– 298. Google Scholar CrossRef Search ADS   Luo Y., Xia J., Miller R.D., Xu Y., Liu J., Liu Q., 2009. Rayleigh-wave mode separation by high-resolution linear Radon transform, Geophys. J. Int. , 179, 254– 264. Google Scholar CrossRef Search ADS   Masoni I., Boelle J.-L., Brossier R., Virieux J., 2016. Layer stripping FWI for surface waves, SEG Technical Program Expanded Abstracts 2016 , Society of Exploration Geophysicists, pp. 1369– 1373. Nakata N., Snieder R., Kuroda S., Ito S., Aizawa T., Kunimi T., 2013. Monitoring a building using deconvolution interferometry. I: Earthquake-data analysis, Bull. seism. Soc. Am. , 103, 1662– 1678. Google Scholar CrossRef Search ADS   Nocedal J., Wright S., 1999, Numerical Optimization , Springer Series in Operations Research, Springer. Google Scholar CrossRef Search ADS   Park C.B., Miller R.D., Xia J., 1998. Imaging dispersion curves of surface waves on multi-channel record, SEG Expanded Abstracts , pp. 1377– 1380. Park C.B., Miller R.D., Xia J., 1999. Multichannel analysis of surface waves, Geophysics , 64, 800– 808. Google Scholar CrossRef Search ADS   Pasquet S., Bodet L., Bergamo P., Guérin R., Martin R., Mourgues R., Tournat V., 2016. Small-scale seismic monitoring of varying water levels in granular media, Vadose Zone J. , 15, 7– 15. Google Scholar CrossRef Search ADS   Pasquet S., Bodet L., Dhemaied A., Mouhri A., Vitale Q., Rejiba F., Flipo N., Guérin R., 2015. Detecting different water table levels in a shallow aquifer with combined P-, surface and SH-wave surveys: Insights from Vp/Vs or Poisson’s ratios, J. Appl. Geophys. , 113, 38– 50. Google Scholar CrossRef Search ADS   Picotti S., Carcione J.M., Gei D., Rossi G., Santos J.E., 2012. Seismic modelling to monitor CO2 geological storage: The Atzbach-Schwanenstadt gas field, J. geophys. Res. , 117, 103– 121. Google Scholar CrossRef Search ADS   Schuster G.T., 2016. Seismic interferometry, in Encyclopedia of Exploration Geophysics , pp. Q1– Q41, eds Grechka V., Wapenaar K., Society of Exploration Geophysicists. Shapiro N.M., Campillo M., Stehly L., Ritzwoller M.H., 2005. High-resolution surface-wave tomography from ambient seismic noise, Science , 307, 1615– 1618. Google Scholar CrossRef Search ADS PubMed  Snieder R., Grêt A., Douma H., Scales J., 2002. Coda wave interferometry for estimating nonlinear behaviour in seismic velocity, Science , 295, 2253– 2255. Google Scholar CrossRef Search ADS PubMed  Socco L.V., Foti S., Boiero D., 2010. Surface-wave analysis for building near-surface velocity models established approaches and new perspectives, Geophysics , 75, A83– A102. Google Scholar CrossRef Search ADS   Strobbia C., Foti S., 2006. Multi-offset phase analysis of surface wave data (MOPA), J. Appl. Geophys. , 59, 300– 313. Google Scholar CrossRef Search ADS   Tanimoto T., 1990. Modelling curved surface wave paths: membrane surface wave synthetics, Geophys. J. Int. , 102, 89– 100. Google Scholar CrossRef Search ADS   Wapenaar K., Slob E., Snieder R., Curtis A., 2010. Tutorial on seismic interferometry: Part 2–Underlying theory and new advances, Geophysics , 75, A211– A227. Xia J.R., Miller D., Park C.B., 1999. Estimation of near-surface shear-wave velocity by inversion of Rayleigh waves, Geophysics , 64, 691– 700. Google Scholar CrossRef Search ADS   Xue Y., Dong S., Schuster G.T., 2009. Interferometric prediction and subtraction of surface waves with a nonlinear local filter, Geophysics , 74, SI1– SI8. Google Scholar CrossRef Search ADS   Yilmaz Ö., 2015. Engineering Seismology with Applications to Geotechnical Engineering , Society of Exploration Geophysicists. Google Scholar CrossRef Search ADS   © 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

# Parsimonious surface wave interferometry

, Volume 212 (3) – Mar 1, 2018
10 pages

/lp/ou_press/parsimonious-surface-wave-interferometry-fBECE1aLTl
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx467
Publisher site
See Article on Publisher Site

### Abstract

Summary To decrease the recording time of a 2-D seismic survey from a few days to one hour or less, we present a parsimonious surface wave interferometry method. Interferometry allows for the creation of a large number of virtual shot gathers from just two reciprocal shot gathers by crosscoherence of trace pairs. Then, the virtual surface waves can be inverted for the S-wave velocity model by wave-equation dispersion inversion (WD). Synthetic and field data tests suggest that parsimonious WD (PWD) gives S-velocity tomograms that are comparable to those obtained from a conventional survey with a shot at each receiver. The limitation of PWD is that the virtual data lose some information so that the resolution of the S-velocity tomogram can be modestly lower than that of the S-velocity tomogram inverted from a conventional survey. Waveform inversion, Computational seismology, Seismic interferometry, Surface waves and free oscillations INTRODUCTION Surface wave inversion is an established method for estimating S-wave velocity models in near-surface, earthquake, and engineering seismology (Xia et al. 1999; Shapiro et al. 2005; Lin et al. 2008; Socco et al. 2010; Yilmaz 2015; Li et al. 2017a). For conventional surface wave surveys along a line of geophones, receivers are deployed along the survey line and common shot gathers (CSGs) are recorded at selected inline source positions. To significantly reduce recording time, we recall the concept of a parsimonious refraction survey where (Fu et al. 2017; Hanafy et al. 2017) reciprocal sources are located only at the beginning and the end locations of the receiver line as shown in Fig. 1. The two reciprocal shot gathers are used to interferometrically generate many virtual shot gathers with virtual refractions (Wapenaar et al. 2010; Schuster 2016). For refraction arrivals, we define this procedure as parsimonious refraction interferometry because it uses a stationary-phase principle to decompose the long-offset refraction data into virtual refraction data with short source-receiver offsets. That is, only two reciprocal shot gathers with 2N traces can be used to generate O(N2) virtual traces where there are N shot gathers with a virtual source at each of the N geophones along the survey line. These virtual shot gathers can be inverted for the P-wave velocity model by traveltime tomography. Figure 1. View largeDownload slide Velocity model where the reciprocal sources are at s and $$\bar{{\bf s}}$$, and the solid white ray paths represent the surface waves excited by the actual source at s. For D(ω)sg and $${D(\omega )_{s\bar{g}}}$$, the phase in D(ω)sg cancels that in $$D(\omega )_{s\bar{g}}$$ along the common ray paths marked by the shortest white line. The resulting phase is associated with the dashed yellow ray path associated with the virtual surface wave trace $$R({\bf s}, \omega )_{g\bar{g}}$$ excited by the virtual source (blue star) at g and recorded by the geophone at $$\bar{\bf g}$$. Figure 1. View largeDownload slide Velocity model where the reciprocal sources are at s and $$\bar{{\bf s}}$$, and the solid white ray paths represent the surface waves excited by the actual source at s. For D(ω)sg and $${D(\omega )_{s\bar{g}}}$$, the phase in D(ω)sg cancels that in $$D(\omega )_{s\bar{g}}$$ along the common ray paths marked by the shortest white line. The resulting phase is associated with the dashed yellow ray path associated with the virtual surface wave trace $$R({\bf s}, \omega )_{g\bar{g}}$$ excited by the virtual source (blue star) at g and recorded by the geophone at $$\bar{\bf g}$$. In a typical exploration survey, the surface waves are considered to be noise so the goal is to remove the surface waves (Dong et al. 2006; Halliday & Curtis 2008; Xue et al. 2009). In contrast, for many engineering surveys, the surface waves are considered to be useful signal so the goal is to invert the surface wave dispersion curves (Park et al. 1998, 1999; Hayashi & Suzuki 2004; Socco et al. 2010; Gouédard et al. 2012; Boiero et al. 2013; Yilmaz 2015). We now extend parsimonious refraction interferometry to surface waves, where trace pairs from the two reciprocal shot gathers are correlated to give approximately N virtual shot gathers that contain virtual surface waves, where N is the number of geophones between the two shots. These virtual surface waves are inverted by wave-equation dispersion inversion (WD; Li & Schuster 2016; Li et al. 2017a) to give the S-velocity tomogram. We denote this combination of interferometry and inversion methods for surface waves as parsimonious WD (PWD). The next section describes the theory of PWD. Then, we present the results of applying PWD to both synthetic and seismic field data. The final section presents the summary and conclusions. THEORY OF PARSIMONIOUS WAVE-EQUATION INVERSION Assume a point source at each end of the survey line that overlies the irregularly layered medium in Fig. 1, where surface waves propagate along the free surface. Here, we assume an isotropic medium with no attenuation. There are velocity variations in the upper medium and N evenly spaced receivers (black triangles) on the recording surface between the two actual sources; these two actual sources at the ends of the survey line are denoted as reciprocal sources (red stars). The actual trace D(ω)sg in the frequency domain is excited by the source at s and recorded at g by a vertical particle-velocity geophone; the point source is a vertical-displacement at s oscillating at frequency ω. The key idea is to crosscorrelate a pair of traces at different receivers $$\bar{{{\bf g}}}$$ and g to obtain a virtual trace at g excited by a virtual source (blue star) at $$\bar{{{\bf g}}}$$. If this is done for all geophone pairs then O(N) virtual shot gathers will be generated, with a virtual shot at each geophone and N virtual traces for each shot gather. The leftmost (rightmost) reciprocal shot generates the right-going (left-going) surface waves. In our examples we assume that the fundamental Rayleigh mode is isolated and its dispersion curve is inverted for the subsurface shear-velocity model. Parsimonious surface-wave interferometry We now present the theory of surface wave interferometry in the context of a parsimonious experiment where there is a line of vertical-component geophones and there are just two vertical-displacement shots, one at each end of the line. We assume an isotropic layered model with a shear wave velocity smoothly varying along the lateral direction and there is no attenuation in the model. Haney & Nakahara (2014, 2016) derived expressions for surface wave Green’s tensors that include near-field behaviour. Here, in order to describe the theory of surface wave interferometry in a simple way, the Rayleigh-wave approximation Green’s function (Snieder et al. 2002) can be described by   \begin{eqnarray} G({\bf g}|{\bf s})&=&\sum _\eta A_\eta \frac{e^{i\theta (k_\eta , {\bf s},{\bf g})+i\pi / 4}}{\sqrt{0.5 \pi k_\eta q_\eta }}, \end{eqnarray} (1)where the receivers at g are in the far-field region of the harmonic point source at s. The summation is over different propagation modes, the horizontal free surface is at z = 0, and kη is the wavenumber of the integer-valued mode denoted by η. For simplicity of exposition we will assume the phase term θ(kη, s, g) = kηX and the geometrical spreading factor to be qη = X, where X = |s − g| is the horizontal distance between the vertical-displacement point source at s = (xs, ys, 0) and the vertical-component geophone at g = (xg, yg, 0). The term Aη accounts for the source amplitude and radiation pattern for a trace at g excited by a point source at s. In general, the phase θ(kη, s, g) and geometrical spreading qη terms are computed by ray tracing for a laterally varying velocity (Tanimoto 1990). We now assume that the surface waves are filtered so that there is only a single mode η = 1 propagating according to equation 1. For a linear recording line with an endline shot at s, the spectral product $$G(\bar{{{\bf g}}}|{\bf s})G({\bf g}|{\bf s})^*$$ of two traces, one recorded at the near-shot position g and the other at the far-shot position $$\bar{{{\bf g}}}$$ along the recording line, becomes   \begin{eqnarray} \Phi (\bar{{{\bf g}}}, {\bf g})_{{\bf s}}=G(\bar{{{\bf g}}}|{\bf s})G({\bf g}|{\bf s})^*&=&A_1 \skew7\bar{A}_1 \frac{e^{ik_1 |{\bf g}-\bar{{{\bf g}}}|}}{\sqrt{0.5 \pi k_1\overline{X}}\sqrt{0.5 \pi k_1 X}}, \end{eqnarray} (2)where X = |s − g| ($$\overline{X}=|{\bf s}-\bar{{{\bf g}}}|$$), and $$\skew7\bar{A}_1$$ (A1) describe the source and radiation characteristics for the trace recorded at $$\bar{{{\bf g}}}$$ (g). The arrival in the virtual trace $$\Phi (\bar{{{\bf g}}}, {\bf g})_{\bf s}$$ has the phase characteristics of a Rayleigh wave (with mode η = 1) excited by a vertical ground displacement at g and recorded by a vertical-component geophone at $$\bar{{{\bf g}}}$$. Therefore, if the recorded traces G(g|s) and $$G(\bar{{{\bf g}}}|{\bf s})$$ are filtered so that they only contain a single mode of the one-way propagating Rayleigh wave, then we conclude that the phase-velocity curve C(ω)η computed from the virtual shot gather $$\Phi (\bar{{{\bf g}}}, {\bf g})_{\bf s}$$ will be identical to that of an actual shot gather $$G(\bar{{{\bf g}}}|{\bf g})$$. The magnitudes of the virtual and actual spectra will differ in the C(ω)η − ω domain but the trajectory of the two curves will be the same. Only an accurate estimate of the phase is needed to give an accurate trajectory of C(ω)η, so that WD can provide an accurate shear-velocity tomogram. If the filtering is incomplete so that there are two modes in the filtered data   \begin{eqnarray} G({\bf g}|{\bf s})&=& A_1 \frac{e^{ik_1X+i\pi /4}}{\sqrt{0.5 \pi k_1 X}}+ A_2 \frac{e^{ik_2X+i\pi /4}}{\sqrt{0.5 \pi k_2 X}} , \end{eqnarray} (3)then the correlated traces yield the virtual trace   \begin{eqnarray} &&{G(\bar{{{\bf g}}}|{{\bf s}})G({\bf g}|{\bf s})^*= \overbrace{A_1 \skew4\bar{A}_1 \frac{e^{ik_1 |{\bf g}\!-\!\bar{{{\bf g}}}|}}{\sqrt{0.5 \pi k_1\overline{X}}\sqrt{0.5 \pi k_1 X}}}^{\text{mode} 1}} \nonumber\\ &&+ \overbrace{A_2 \skew4\bar{A}_2 \frac{e^{ik_2 |{\bf g}-\bar{{{\bf g}}}|}}{\sqrt{0.5 \pi k_2\overline{X}}\sqrt{0.5 \pi k_2 X}}}^{\text{mode} 2}\nonumber \\ &&+ \overbrace{A_1 \skew4\bar{A}_2 \frac{e^{ik_1 |{\bf s}\!-\!\bar{{{\bf g}}}|-ik_2|{\bf s}\!-\! {\bf g}|}}{\sqrt{0.5 \pi k_1 \overline{X}}\sqrt{0.5 \pi k_2 X}} \!+\! \skew4\bar{A}_1 A_2 \frac{e^{ik_2 |{\bf s}-\bar{{{\bf g}}}|-ik_1|{\bf s}- {\bf g}|}}{\sqrt{0.5 \pi k_1 X}\sqrt{0.5 \pi k_2 \overline{X}}}}^{\text{crosstalk}}.\nonumber\\ \end{eqnarray} (4)In this case, the dispersion curves associated with the virtual shot gather $$G(\bar{{{\bf g}}}|{\bf s})G({{\bf g}}|{\bf s})^*$$ will contain the correct C(ω)η for modes η = 1 and η = 2, but they will also contain the noisy crosstalk modes. Such noise will hinder an accurate estimation of C(ω)η. We conclude that, prior to crosscorrelation, it is important to filter the shot gather so that it only contains one mode propagating in one direction. In practice, this mode is the fundamental mode of the Rayleigh wave because it is the strongest arrival in the records. Workflow of parsimonious wave-equation inversion The workflow for PWD can be summarized by the following three steps. For one shot gather, crosscoherence of the traces at receivers A = g and B = $$\bar{{{\bf g}}}$$ generates a virtual seismogram where the receiver at the A = g position is now the new virtual source position shown in Fig. 2. Virtual surface waves can be reconstructed at each receiver location by summation over sources and by using the crosscoherence of the normalized traces recorded by the receivers at g and $$\bar{{{\bf g}}}$$. In the frequency domain this operation is represented as (Xue et al. 2009; Nakata et al. 2013):   \begin{eqnarray} R({\bf s}, \omega )_{g\bar{g}}=\sum _{s} \frac{D(\omega )_{\bar{g}s}D(\omega )^*_{gs}}{|D(\omega )_{gs}||D(\omega )_{\bar{g}s}|+\lambda }, \end{eqnarray} (5)where $$R({\bf s}, \omega )_{g\bar{g}}$$ represents the time-symmetric virtual trace in the frequency domain with right-going waves recorded at $$\bar{{{\bf g}}}$$ for the virtual source at g, and λ is the damping parameter that is from 1 per cent to 10 per cent of the signal amplitude (Nakata et al. 2013). Eq. (5) sums over each source deployed at the left end of the line, and for our examples we find that it is sufficient to use just one source at the end of the line. The recorded data are FK filtered so only right-going surface waves are contained in D(ω)gs and $$D(\omega )_{\bar{g}s}$$ if s is at the left end of the line, and if s is at the right end of the line then the FK filter is used to eliminate left-going waves in D(ω)gs and $$D(\omega )_{\bar{g}s}$$. Here, g is several wavelengths from $$\bar{{{\bf g}}}$$ in order for the surface waves to develop. Halliday & Curtis (2008) suggest that the fundamental mode of the Rayleigh waves can be reliably recovered by interferometry for inline sources on the surface if there is no crosstalk from other modes. Therefore, the PWD method will concentrate on recovering and inverting the fundamental modes of Rayleigh waves for the N virtual shot gathers created from just two reciprocal shot gathers with their shot locations at the two ends of the survey line. To mitigate the crosstalk between different modes, the reciprocal shot gather can be transformed into the phase-velocity and frequency domains by a frequency-domain fast Fourier transform (FFT) and a linear Radon transform (LRT). Then, the fundamental Rayleigh wave can be windowed to eliminate higher-order modes, and the result will be a dispersion curve that is devoid of higher-mode energy. This windowed result can be inverse transformed back into the space–time domain to give a shot gather that only contains the fundamental mode of the Rayleigh waves (Luo et al. 2009; Hu et al. 2016). The filtered reciprocal shot gather can then be used to generate the virtual shot gathers by crosscoherence of trace pairs and summation over the two reciprocal sources. Use an LRT method to transform the filtered virtual shot gathers $$R({\bf s}, \omega )_{g\bar{g}}$$ into the C(g, ω) − ω domain, where C(g, ω) is the phase velocity for the virtual shot gather with the virtual shot at g. The fundamental Rayleigh dispersion curve C(g, ω)obs is picked for each shot gather with the virtual source at g. These curves are used to generate the associated wavenumber dispersion curve by the formula κ(g, ω)obs = ω/C(g, ω)obs. The dispersion curves are used to give the misfit function   \begin{eqnarray} \epsilon =\frac{1}{2} \sum _{\omega } \sum _g (\kappa ({\bf g}, \omega )^{\text{pred}} - \kappa ({\bf g}, \omega )^{\text{obs}})^2, \end{eqnarray} (6)where κ(g, ω)pred is the predicted surface wave dispersion curve computed from a 2-D finite-difference solution to the elastic wave equation with a vertical point source at g. The S-wave velocity is iteratively updated by the iterative steepest-descent formula of Li & Schuster (2016) and Li et al. (2017a):   \begin{eqnarray} s({\bf x})^{(k+1)} &=& s({\bf x})^{(k)} - \alpha \sum _{\omega } \sum _g \Delta \kappa ({\bf g},\omega ) \frac{{\partial } \kappa ({\bf g},\omega )^{\text{pred}} }{{\partial } s({\bf x})}, \end{eqnarray} (7)where Δκ(g, ω) = κ(g, ω)pred − κ(g, ω)obs is the residual at the kth iteration, and s(x)(k) is the shear slowness at the kth iteration. Here, α is the step length by any backtracking line-search method (Nocedal & Wright 1999) and the k superscript denotes the kth iteration. For pedagogical simplicity, we assume a single shot gather $$\hat{D}({{\bf g}},\omega )$$ for a virtual or real source at s and geophones at g, and the notation for the source location is silent. The misfit function will include an additional summation over different shot gathers if more than one virtual shot gather is used. Mathematical details for deriving the Fréchet derivative $$\frac{{\partial } \kappa (\omega )}{{\partial } s({\bf x})}$$ are given in Li et al. (2017a), where   \begin{eqnarray} &&{\frac{{\partial } \kappa }{{\partial } s(\mathbf {x})} = \frac{4s(\mathbf {x})\rho ({\bf x})}{A} \text{Real} \Big\lbrace \overbrace{\frac{{\partial } D({\bf x},\omega )_{x}}{{\partial } x}}^{\text{source}=f({{\bf x}},\mathbf {s},\omega )_1}}\nonumber \\ && \int \overbrace{\hat{D}({{\bf g}},\omega )_{\text{obs}}^*\left(\frac{{\partial } G({\bf g}|{\bf x})_{zz}}{{\partial } z}\right)\text{d}x_{g}}^{\text{backprojected data}=B(\mathbf {x},\mathbf {s},\omega )_{1}^*} +\overbrace{\frac{{\partial } D({\bf x},\omega )_{z}}{{\partial } z}}^{\text{source}=f({{\bf x}},\mathbf {s},\omega )_2}\nonumber \\ && \int \overbrace{\hat{D}({{\bf g}},\omega )_{\text{obs}}^*\frac{{\partial } G({\bf g}|{\bf x})_{xz}}{{\partial } x}\text{d}x_{g}}^{\text{backprojected data}=B(\mathbf {x},\mathbf {s},\omega )_{2}^*}\nonumber \\ && -\frac{1}{2} \overbrace{\left(\frac{{\partial } D({\bf x},\omega )_{z}}{{\partial } x}+\frac{{\partial } D({\bf x},\omega )_{x}}{{\partial } z}\right)}^{\text{source}={f(\mathbf {x},\mathbf {s},\omega )_3}}\nonumber \\ && \int \overbrace{\hat{D}({{\bf g}},\omega )_{\text{obs}}^*\left(\frac{{\partial } G({\bf g}|{\bf x})_{xz}}{{\partial } z} +\frac{{\partial } G({\bf g}|{\bf x})_{zz}}{{\partial } x}\right)}^{\text{backprojected data}={ B(\mathbf {x},\mathbf {s},\omega )_{3}^*}}\text{d}x_{g} \Big\rbrace. \end{eqnarray} (8)Here, $$\hat{D}({{\bf g}},\omega )_{\text{obs}}^*$$ is the weighted conjugated data for the vertical component displacement recorded at g, G(g|x)ij is the Green’s tensor for the ith component of the displacement field at (x, z) for a jth-component source at (xg, 0), s(x) is the shear-slowness model for the background medium, and A is the normalization term. Here, f(x, s, ω)i for i ∈ (1, 2, 3) is the source field at $$\bf x$$ that originates at $$\bf s$$, and B(x, s, ω)i represents the backprojected field at $$\bf x$$ for a source at s; it is assumed that this field is recorded at the surface and backpropagated by the adjoint of the forward modelling operator. The interpretation of equations (7) and (8) is that the traces and source wavelet are weighted by terms proportional to the wavenumber residual Δκ(ω) and backpropagated into the medium to update the slowness of the S-velocity model. Fig. 3 summarizes the workflow for the PWD method. Figure 2. View largeDownload slide Crosscoherence of two traces at different source–receiver offsets yields a trace (on the right) with a virtual shot at the A receiver. Figure 2. View largeDownload slide Crosscoherence of two traces at different source–receiver offsets yields a trace (on the right) with a virtual shot at the A receiver. Figure 3. View largeDownload slide The workflow for the parsimonious surface wave WD method, where illumination compensation is applied to the gradient as a pre-conditioner. Here, RT denotes the Radon transform and the comparison between the conventional data and virtual data is only used for assessing the accuracy of PWD with synthetic data. In practice, the conventional data for shots at all geophones is usually not available. Figure 3. View largeDownload slide The workflow for the parsimonious surface wave WD method, where illumination compensation is applied to the gradient as a pre-conditioner. Here, RT denotes the Radon transform and the comparison between the conventional data and virtual data is only used for assessing the accuracy of PWD with synthetic data. In practice, the conventional data for shots at all geophones is usually not available. PARSIMONIOUS WAVE-EQUATION DISPERSION INVERSION TESTS In this section, the PWD method is used to invert for the S-velocity model from three different sets of synthetic data and one field data set recorded near the Gulf of Aqaba in Saudi Arabia. We also compared the PWD tomogram (inverted from N virtual shot gathers) with the conventional WD tomogram (inverted from N actual shot gathers). In each example the observed surface wave data are generated using time-domain finite-difference solutions to the 2-D elastic wave equation (Levander 1988). Stair-step velocity model The geometry of a conventional seismic survey and a parsimonious survey with a stair-step S-wave velocity model are shown in Fig. 4. The 2-D model size is 20 × 120 m2, the S-wave velocities in the top and bottom layers are 500 ms−1 and 800 ms−1, respectively, the P-wave velocity Vp is constrained to be $$V_p=\sqrt{3} V_s$$ and the source wavelet is a Ricker wavelet with a peak frequency of 30 Hz. The conventional survey consists of 60 shots (red stars) shooting into 60 receivers (white triangles) evenly deployed at 2 m intervals with a 2 m shot interval. To test the parsimonious method, we used two reciprocal sources located at the beginning and the end of the survey line (black stars) while the data are recorded at all receiver locations. Figure 4. View largeDownload slide The acquisition geometry of the conventional and parsimonious surveys. Figure 4. View largeDownload slide The acquisition geometry of the conventional and parsimonious surveys. We computed 60 virtual shot gathers from two reciprocal shot gathers by the cross-coherence method. Fig. 5(a) displays one virtual shot gather (red traces), which almost matches the corresponding shot gather in the conventional survey data (black traces), except their amplitudes have some differences. These amplitude differences will not significantly affect the result of S-wave inversion because the WD method inverts the residuals associated with the phase velocity (Li & Hanafy 2016; Li & Schuster 2016; Li et al. 2017a). Figure 5. View largeDownload slide (a) Actual shot gather (red) and virtual shot gather (black) computed by crosscoherence, and (b) the spectra of single traces. Figure 5. View largeDownload slide (a) Actual shot gather (red) and virtual shot gather (black) computed by crosscoherence, and (b) the spectra of single traces. Fig. 5(b) shows the spectra computed from the actual traces and the virtual traces computed by the crosscoherence method, where the useful frequencies in the original traces (red spectrum) are from 10 to 80 Hz while those in the virtual traces computed by crosscorrelation (black spectrum) are only from 15 to 75 Hz. However, in a noisy environment the signals above a critical frequency are too noisy so a high-cut filter should be applied to the virtual traces computed by crosscoherence. The high-cut frequency is the critical frequency identified in the data. A phase-velocity model (see Fig. 6) is computed from both the actual shot gathers and the virtual shot gathers using the LRT. The phase-velocity model computed from the 60 actual shot gathers (Fig. 6a) closely resembles the model inverted from the 60 virtual shot gathers (i.e. the parsimonious data) (Fig. 6b). The PWD and conventional WD tomograms are shown in Figs 7(c) and (d). Both tomograms are mostly consistent with the true S-wave velocity model in Fig. 7(a). For comparison, Fig. 7(b) shows the inaccurate S-velocity tomogram inverted from just two reciprocal shot gathers (without using the crosscoherence operation to generate the virtual traces). The normalized data residual plotted against iteration number in Fig. 8 shows that the PWD method has a convergence rate similar to that of the conventional WD method, but the conventional WD tomogram has a smaller final residual. Figure 6. View largeDownload slide Phase-velocity models computed from the (a) conventional data and (b) virtual data. Figure 6. View largeDownload slide Phase-velocity models computed from the (a) conventional data and (b) virtual data. Figure 7. View largeDownload slide (a) True S-wave velocity model; (b) WD tomogram inverted from two reciprocal shot gathers; tomograms computed by the (c) conventional WD and (d) parsimonious WD methods. Figure 7. View largeDownload slide (a) True S-wave velocity model; (b) WD tomogram inverted from two reciprocal shot gathers; tomograms computed by the (c) conventional WD and (d) parsimonious WD methods. Figure 8. View largeDownload slide RMS residual plotted against iteration number for the two WD tomograms in Figs 7(c) and (d), where the red and black lines are for the conventional and parsimonious residuals, respectively. Figure 8. View largeDownload slide RMS residual plotted against iteration number for the two WD tomograms in Figs 7(c) and (d), where the red and black lines are for the conventional and parsimonious residuals, respectively. Local-anomalies model The next test is for a velocity model with local-velocity anomalies shown in Fig. 9(a). This model is based on geological observations in a shallow open-pit mine with mineral deposits. The 2-D model size is 20 × 170 m2 and the source wavelet is a Ricker wavelet with a peak frequency of 30 Hz. We follow the same workflow shown in Fig. 3, where the input data for conventional WD consists of 40 shot gathers with 85 traces and receivers located every 2 m, with the shots evenly distributed every 4 m on the ground surface. The parsimonious data consist of two CSGs generated by the shots at the two ends of the survey line. These two reciprocal CSGs are used to generate 40 virtual shots gathers by crosscoherence interferometry. Figure 9. View largeDownload slide (a) True S-wave velocity model, (b) initial S-wave velocity model, (c) conventional WD tomogram, and (d) parsimonious WD tomogram. Figure 9. View largeDownload slide (a) True S-wave velocity model, (b) initial S-wave velocity model, (c) conventional WD tomogram, and (d) parsimonious WD tomogram. The S-velocity tomograms inverted by the conventional WD and PWD methods are shown in Figs 9(c) and (d). Both tomograms are in reasonable agreement with the true model, especially with the high-velocity anomaly delineated by the white dashed lines. After 31 iterations (Fig. 10), the normalized data misfit residuals of PWD and conventional WD decreased to the approximate values of 0.2 and 0.1, respectively. Figure 10. View largeDownload slide RMS residual as a function of iteration number for the conventional (red) and parsimonious (black) data computed for the Fig. 10a model. Figure 10. View largeDownload slide RMS residual as a function of iteration number for the conventional (red) and parsimonious (black) data computed for the Fig. 10a model. Aqaba field data test PWD is now tested on a field data set recorded next to the Gulf of Aqaba (see Fig. 11) in Saudi Arabia (Hanafy et al. 2014). For the conventional data survey, a total of 120 common shot gathers (CSGs) are recorded, and each shot gather has 120 traces at equal shot and receiver intervals of 2.5 m. The total length of the profile is 297.5 m and a 90 kg weight drop is used as the seismic source, with 10 to 15 stacks at each shot location. We select three shot gathers (shot locations in the middle and at the two ends of the survey line) to generate 120 virtual shot gathers by crosscoherence. Fig. 12 shows that the surface waves in the virtual shot gather compare well with those in the actual shot gather for the same source location. Figure 11. View largeDownload slide (a) Google map shows the survey location of the Aqaba experiment, and (b) photo shows an earthquake rupture (indicated by the red dashed line) on the surface. The black arrow indicates the location of the seismic survey line. Figure 11. View largeDownload slide (a) Google map shows the survey location of the Aqaba experiment, and (b) photo shows an earthquake rupture (indicated by the red dashed line) on the surface. The black arrow indicates the location of the seismic survey line. Figure 12. View largeDownload slide Conventional (black) and virtual (red) trace comparisons at different shot locations. Figure 12. View largeDownload slide Conventional (black) and virtual (red) trace comparisons at different shot locations. Following the workflow shown in Fig. 3, the conventional and PWD tomograms are computed and shown in Figs 13(b) and (c). The S-velocity tomograms generally agree with one another for the long-wavelength features. According to the P-velocity tomogram in Fig. 13(a), the actual fault position is at x = 150 m, which is at the same position indicated in the PWD tomogram. It is evident that PWD provides an image similar to that computed from the data recorded by a conventional seismic survey with 120 shots. However, the conventional survey required more than an order-of-magnitude increase in survey time compared to the parsimonious survey. Figure 13. View largeDownload slide (a) P-velocity inverted from first-arrival traveltimes tomogram, S-velocity tomograms computed from the (b) conventional and (c) parsimonious surface wave data. Figure 13. View largeDownload slide (a) P-velocity inverted from first-arrival traveltimes tomogram, S-velocity tomograms computed from the (b) conventional and (c) parsimonious surface wave data. Time-lapse monitoring test Seismic time-lapse techniques are often applied to large areas, such as hydrocarbon reservoirs (Landro 2001) or geologic storage reservoirs (Carcione et al. 2006; Picotti et al. 2012). For near-surface geophysics, the potential application of seismic methods for tracking time-variant processes is summarized by Jefferson et al. (1998). As an example, seismic surface wave techniques have been proposed for the detection of varying levels of the water table of shallow aquifers, heavy-oil in oil fields (Pasquet et al. 2015; Dubos-Sallee et al. 2015; Ikeda et al. 2016; Pasquet et al. 2016), or monitoring climate effects on earthworks (Bergamo et al. 2016a,b). A significant limitation of a small-scale time-lapse seismic survey, for example, 240 shot gathers, is that it takes at least 16 working hours for completion if the source excitation at each of the 240 source points requires more than 5 min of effort. This is a major problem if the geological properties of the subsurface target, such as in a fluid-injection experiment, rapidly vary in time. In comparison, acquisition of parsimonious data is fast and efficient because it only requires several shot gathers where the source can be far from the fluid-injection site. In this section, we will use synthetic data to evaluate the feasibility of using time-lapse parsimonious WD to monitor the subsurface fluid distribution after dumping a viscous fluid on to the surface of the recording survey. In this case we assume the viscous fluid increases the S-velocity of the porous rock. Synthetic model test The time-lapse fluid model shown in Fig. 14(a) is a two-layered model where the fluid is leaking from the surface into the model. A 2-D time-domain finite-difference method is used to compute the vertical-component elastic seismograms. The model size is 25 × 120 m2 where the grid spacing and time sampling are 1m and 0.25 ms, respectively, and the source time-history is a Ricker wavelet peaked at 30 Hz. Here, we assume the fluid increases the S-velocity of the subsurface while in practice we find that the S-velocity can be lowered after injecting water into porous sand. Figure 14. View largeDownload slide S -wave velocity models corresponding to the (a) time-lapse fluid-injection tests, (b) conventional WD tomogram, (c) parsimonious WD tomogram and (d) velocity changes of the parsimonious S-velocity tomograms before and after fluid injection at different times. Figure 14. View largeDownload slide S -wave velocity models corresponding to the (a) time-lapse fluid-injection tests, (b) conventional WD tomogram, (c) parsimonious WD tomogram and (d) velocity changes of the parsimonious S-velocity tomograms before and after fluid injection at different times. First, we generate the conventional-acquisition shot gathers which consist of 30 CSGs with sources located every 4 m along the recording line. Each shot is recorded by 60 receivers spaced at 2 m intervals. The conventional WD S-velocity tomogram, shown in Fig. 14(b), is in good agreement with the true fluid-injection velocity model. Now, we test the PWD approach, where we only use two reciprocal shot gathers with shots located at x = 0 m and x = 120 m. In both shot gathers the data are recorded at 60 receivers spaced at 2 m intervals along the ground surface. These reciprocal shot gathers are used to generate 30 virtual shot gathers by crosscoherence of traces. The corresponding S-velocity tomograms are shown in Fig. 14(c), where the tomograms show a clear time-lapse velocity anomaly that is almost the same as seen in the conventional WD tomogram. However, the conventional WD tomogram has a modestly higher resolution than the parsimonious data tomogram. This lowered resolution is possibly caused by imperfect filtering so that mode coupling artefacts appear in the correlation. The difference between the parsimonious S-velocity tomograms computed from data recorded before and after fluid injection is shown in Fig. 14(d). This clearly reveals the velocity changes as a function of time. DISCUSSION A limitation of PWD is that it typically provides tomograms that have modestly lower resolution than the conventional tomograms. The partial cause is that the virtual data are created by correlating two traces with one another, so that the spectrum of the virtual source wavelet is the square of the wavelet spectrum. This tends to emphasize the central frequencies at the expense of downweighting the higher and lower frequencies, as illustrated in Fig. 5. Thus, crosscoherence or spectral deconvolution is required, but its effectiveness can be degraded by noise in the data and the choice of the damping parameter. In addition, mode mixing can be activated by crosscoherence, so ideally higher-order modes should be eliminated prior to crosscoherence of traces (Hu et al. 2016). Another possibility is to design matching filters that match the traces from the parsimonious data to the conventional data. These matching filters can then be applied to the time lapse virtual data to deconvolve the source wavelet. This is a subject of future research. There are two other limitations. The first one is that PWD theory assumes no attenuation in the data. In this case attenuation should be estimated from the data and, if possible, the data should be corrected for attenuation effects. The second limitation is that the accuracy of the estimated dispersion curve depends on the signal-to-ratio of the data and the window selection to isolate the surface waves of interest. Thus, sensitivity tests should be carried out to determine the window design that gives the most robust estimate of the dispersion curves. Similar to a multi-frequency FWI strategy (Masoni et al. 2016), we can utilize multi-offset dispersion curves to invert for the S-velocity at different depths. Shot-offset traces have high horizontal resolution at the expense of decreased accuracy at the lower frequencies, whereas wide-offset data have the opposite characteristics (Strobbia & Foti 2006). So, the combined short-to-wide offset data with gradually increasing frequency ranges can be used to optimize the accuracy of the reconstructed S-velocity model. Details are provided in Li et al. (2017b). CONCLUSIONS We presented the theory of PWD where a dense set of virtual surface wave data is interferometrically generated from a pair of reciprocal shot gathers, and the shear-velocity model is computed by WD. This method is denoted as PWD and assumes an isotropic elastic medium with no attenuation. We compare the S-velocity tomograms generated by PWD and conventional WD using synthetic and seismic field data sets. Our results show that the PWD can provide nearly the same S-velocity tomograms as conventional WD, but with slightly less resolution. Our synthetic-data results also suggest that time-lapse parsimonious surface wave monitoring can be a robust and reliable technique for monitoring small-scale and rapidly varying subsurface conditions in the subsurface. Possible applications of PWD are for environmental monitoring of water-injection experiments for hazard remediation, archaeological sites where sources are not permitted over sensitive areas, characterizing the processes in erupting geysers or volcanoes, and engineering surveys. The PWD methodology can be adapted to inversion of dispersion curves computed from Love waves and trapped P-waves in wave guides. ACKNOWLEDGEMENTS We thank the sponsors for supporting the Consortium of Subsurface Imaging and Fluid Modelling (CSIM). We also thank KAUST for providing funding by the CRG grant OCRF-2014-CRG3-2300. For computer time, this research used the resources of the IT Research Computing Group and the Supercomputing Laboratory at KAUST. REFERENCES Bergamo P.B., Dashwood S., Uhlemann R., Swift J.E., Chambers D., Gunn A., Donohue S., 2016a. Time-lapse monitoring of climate effects on earthworks using surface waves, Geophysics , 81, EN1– EN15. Google Scholar CrossRef Search ADS   Bergamo P.B., Dashwood S., Uhlemann R., Swift J.E., Chambers D., Gunn A., Donohue S., 2016b. Time-lapse monitoring of fluid-induced geophysical property variations within an unstable earthwork using P-wave refraction, Geophysics , 81, EN17– EN27. Google Scholar CrossRef Search ADS   Boiero D., Wiarda E., Vermeer P., 2013. Surface-and guided-wave inversion for near-surface modelling in land and shallow marine seismic data, Leading Edge , 32, 638– 646. Google Scholar CrossRef Search ADS   Carcione J.M., Picotti S., Gei D., Rossi G., 2006. Physics and seismic modelling for monitoring CO2 storage, Pure appl. Geophys. , 163, 175– 207. Google Scholar CrossRef Search ADS   Dong S., He R., Schuster G.T., 2006. Interferometric prediction and least squares subtraction of surface waves, SEG Technical Program Expanded Abstracts 2006 , Society of Exploration Geophysicists, pp. 2783– 2786. Dubos-Sallee N., Rasolofosaon P., Barbouteau S., Charlety J., Duret F., 2015. Time-lapse surface wave analysis using noise correlation on an onshore heavy-oil field, SEG Technical Program Expanded Abstracts 2015 , Society of Exploration Geophysicists, pp. 5479– 5484. Fu L., Hanafy S., Schuster G.T., 2017. Parsimonious wave equation refraction interferometry and inversion, Geophys. Prospect. , doi:10.1111/1365–2478.12488. Gouédard P., Yao H., Ernst F., van der Hilst R.D., 2012. Surface wave Eikonal tomography in heterogeneous media using exploration data, Geophys. J. Int. , 191, 781– 788. Google Scholar CrossRef Search ADS   Halliday D., Curtis A., 2008. Seismic interferometry, surface waves and source distribution, Geophys. J. Int. , 175, 1067– 1087. Google Scholar CrossRef Search ADS   Hanafy S.M., Jonsson S., Kger Y., 2014. Imaging normal faults in alluvial fans using geophysical techniques: field example from the coast of Gulf of Aqaba, Saudi Arabia, SEG Technical Program Expanded Abstracts 2014 , Society of Exploration Geophysicists, pp. 4670– 4674. Hanafy S.M., Li J., Schuster G., 2017. Time-lapse monitoring of subsurface fluid flow using parsimonious seismic interferometry, SAGEEP Technical Program Expanded Abstracts 2017 , SAGEEP, pp. 1– 5. Haney M.M., Nakahara H., 2014. Surface-wave Greens tensors in the near field, Bull. seism. Soc. Am. , 104, 1578– 1586. Google Scholar CrossRef Search ADS   Haney M.M., Nakahara H., 2016. Surface-wave Greens tensors in the near field, Bull. seism. Soc. Am. , 106, 816– 818. Google Scholar CrossRef Search ADS   Hayashi K., Suzuki H., 2004. CMP cross-correlation analysis of multi-channel surface-wave data, Explor. Geophys. , 35, 7– 13. Google Scholar CrossRef Search ADS   Hu Y., Cheng F., Bin M., Shen C., 2016. Waveform inversion of Rayleigh-wave in time domain for shallow shear-wave velocity based on the waveform separation, in AGU Annual Fall Meeting Abstracts , San Francisco, AGU, pp. 1– 5. Ikeda T., Tsuji T., Takanashi M., Kurosawa I., Nakatsukasa M., White D., Worth K., Roberts B., 2016. Time-lapse monitoring of shallow subsurface in the aquistore CO2 storage site from surface-wave analysis using a continuous and controlled seismic source, SEG Technical Program Expanded Abstracts 2016 , Society of Exploration Geophysicists, pp. 5479– 5484. Jefferson R.D., Steeples D.W., Black R.A., Carr T., 1998. Effects of soil-moisture content on shallow-seismic data, Geophysics , 63, 1357– 1362. Google Scholar CrossRef Search ADS   Landro M., 2001. Discrimination between pressure and fluid saturation changes from time-lapse seismic data, Geophysics , 66, 836– 844. Google Scholar CrossRef Search ADS   Levander A.R., 1988. Fourth-order finite-difference P–SV seismograms, Geophysics , 53, 1425– 1436. Google Scholar CrossRef Search ADS   Li J., Hanafy S., 2016. Skeletonized inversion of surface waves: active source versus controlled noise comparison, Interpretation , 3, 11– 19. Li J., Schuster G., 2016. Skeletonized wave equation of surface wave dispersion inversion, SEG Technical Program Expanded Abstracts 2016 , Society of Exploration Geophysicists, pp. 3630– 3635. Li J., Feng Z., Schuster G.T., 2017a. Wave-equation dispersion inversion, Geophys. J. Int. , 208, 1567– 1578. Google Scholar CrossRef Search ADS   Li J., Lin F., Schuster G.T., 2017b. Wave equation dispersion inversion with irregular topography, Geophys. J. Int. , in press. Lin F.-C., Moschetti M.P., Ritzwoller M.H., 2008. Surface wave tomography of the Western United States from ambient seismic noise: Rayleigh and Love wave phase velocity maps, Geophys. J. Int. , 173, 281– 298. Google Scholar CrossRef Search ADS   Luo Y., Xia J., Miller R.D., Xu Y., Liu J., Liu Q., 2009. Rayleigh-wave mode separation by high-resolution linear Radon transform, Geophys. J. Int. , 179, 254– 264. Google Scholar CrossRef Search ADS   Masoni I., Boelle J.-L., Brossier R., Virieux J., 2016. Layer stripping FWI for surface waves, SEG Technical Program Expanded Abstracts 2016 , Society of Exploration Geophysicists, pp. 1369– 1373. Nakata N., Snieder R., Kuroda S., Ito S., Aizawa T., Kunimi T., 2013. Monitoring a building using deconvolution interferometry. I: Earthquake-data analysis, Bull. seism. Soc. Am. , 103, 1662– 1678. Google Scholar CrossRef Search ADS   Nocedal J., Wright S., 1999, Numerical Optimization , Springer Series in Operations Research, Springer. Google Scholar CrossRef Search ADS   Park C.B., Miller R.D., Xia J., 1998. Imaging dispersion curves of surface waves on multi-channel record, SEG Expanded Abstracts , pp. 1377– 1380. Park C.B., Miller R.D., Xia J., 1999. Multichannel analysis of surface waves, Geophysics , 64, 800– 808. Google Scholar CrossRef Search ADS   Pasquet S., Bodet L., Bergamo P., Guérin R., Martin R., Mourgues R., Tournat V., 2016. Small-scale seismic monitoring of varying water levels in granular media, Vadose Zone J. , 15, 7– 15. Google Scholar CrossRef Search ADS   Pasquet S., Bodet L., Dhemaied A., Mouhri A., Vitale Q., Rejiba F., Flipo N., Guérin R., 2015. Detecting different water table levels in a shallow aquifer with combined P-, surface and SH-wave surveys: Insights from Vp/Vs or Poisson’s ratios, J. Appl. Geophys. , 113, 38– 50. Google Scholar CrossRef Search ADS   Picotti S., Carcione J.M., Gei D., Rossi G., Santos J.E., 2012. Seismic modelling to monitor CO2 geological storage: The Atzbach-Schwanenstadt gas field, J. geophys. Res. , 117, 103– 121. Google Scholar CrossRef Search ADS   Schuster G.T., 2016. Seismic interferometry, in Encyclopedia of Exploration Geophysics , pp. Q1– Q41, eds Grechka V., Wapenaar K., Society of Exploration Geophysicists. Shapiro N.M., Campillo M., Stehly L., Ritzwoller M.H., 2005. High-resolution surface-wave tomography from ambient seismic noise, Science , 307, 1615– 1618. Google Scholar CrossRef Search ADS PubMed  Snieder R., Grêt A., Douma H., Scales J., 2002. Coda wave interferometry for estimating nonlinear behaviour in seismic velocity, Science , 295, 2253– 2255. Google Scholar CrossRef Search ADS PubMed  Socco L.V., Foti S., Boiero D., 2010. Surface-wave analysis for building near-surface velocity models established approaches and new perspectives, Geophysics , 75, A83– A102. Google Scholar CrossRef Search ADS   Strobbia C., Foti S., 2006. Multi-offset phase analysis of surface wave data (MOPA), J. Appl. Geophys. , 59, 300– 313. Google Scholar CrossRef Search ADS   Tanimoto T., 1990. Modelling curved surface wave paths: membrane surface wave synthetics, Geophys. J. Int. , 102, 89– 100. Google Scholar CrossRef Search ADS   Wapenaar K., Slob E., Snieder R., Curtis A., 2010. Tutorial on seismic interferometry: Part 2–Underlying theory and new advances, Geophysics , 75, A211– A227. Xia J.R., Miller D., Park C.B., 1999. Estimation of near-surface shear-wave velocity by inversion of Rayleigh waves, Geophysics , 64, 691– 700. Google Scholar CrossRef Search ADS   Xue Y., Dong S., Schuster G.T., 2009. Interferometric prediction and subtraction of surface waves with a nonlinear local filter, Geophysics , 74, SI1– SI8. Google Scholar CrossRef Search ADS   Yilmaz Ö., 2015. Engineering Seismology with Applications to Geotechnical Engineering , Society of Exploration Geophysicists. Google Scholar CrossRef Search ADS   © 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

## 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
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### 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. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations