# Superresolution near-field imaging with surface waves

Superresolution near-field imaging with surface waves Summary We present the theory for near-field superresolution imaging with surface waves and time reverse mirrors (TRMs). Theoretical formulae and numerical results show that applying the TRM operation to surface waves in an elastic half-space can achieve superresolution imaging of subwavelength scatterers if they are located less than about 1/2 of the shear wavelength from the source line. We also show that the TRM operation for a single frequency is equivalent to natural migration, which uses the recorded data to approximate the Green’s functions for migration, and only costs O(N4) algebraic operations for post-stack migration compared to O(N6) operations for natural pre-stack migration. Here, we assume the sources and receivers are on an N × N grid and there are N2 trial image points on the free surface. Our theoretical predictions of superresolution are validated with tests on synthetic data. The field-data tests suggest that hidden faults at the near surface can be detected with subwavelength imaging of surface waves by using the TRM operation if they are no deeper than about 1/2 the dominant shear wavelength. Image processing, Surface waves and free oscillations 1 INTRODUCTION Superresolution imaging is defined as estimating the location of an object to a resolution much better than λ/2, where λ is the effective wavelength of the wavefield. For optical wavefields, De Fornel (2001) and Jia et al. (2010) utilized evanescent energy in the near-field for their high-resolution optical superlenses. Similarly, Lerosey et al. (2007) demonstrated that superresolution imaging is also possible with microwaves recorded in the far-field region of subwavelength scatterers, when the scatterers are near the source. Their microwave experiment with subwavelength scatterers near the source showed that imaging resolution could be better than λ/30 by the use of a time reversal mirror (TRM) operation (Fink 1993, 2008) applied to the recorded microwave data. Schuster et al. (2012) analysed the TRM operation in the context of seismic scattering and theoretically showed that near-field seismic superresolution is enabled by the evanescent field associated with the geometrical-spreading factor. Their theoretical treatment of the evanescent field was restricted to body-wave arrivals in the acoustic approximation. Both simulations and field experiments (Dutta et al. 2016) validated the theoretical predictions. Their imaging method, denoted as the seismic scanning tunnelling macroscope (SSTM) operation, was extended to elastic body waves by Tarhini et al. (2017), and the superresolution properties were validated by tests on both synthetic data and field data. The TRM operation is equivalent to natural migration (Schuster 2002; Schuster et al. 2012; AlTheyab et al. 2016) because it uses the recorded Green’s functions for migration and so makes no assumptions about the velocity model. These ‘natural’ Green’s functions eliminate the need for knowing the velocity model. AlTheyab et al. (2016) and Liu et al. (2017) used natural pre-stack migration of scattered surface waves to estimate the locations of hidden faults that were within λ/3 of the free surface. Migrating surface waves is a robust imaging strategy because surface waves are more than an order-of-magnitude stronger than most body-wave reflections. There are two important questions posed by the AlTheyab et al. (2016) study with natural migration. Can surface waves be used to achieve superresolution imaging? The answer is yes as we will show both theoretically and numerically. How is post-stack migration related to the natural pre-stack migration of AlTheyab et al. (2016)? As will be shown, TRM imaging of data generated by simultaneously exploding sources is equivalent to natural post-stack migration while Altheyab’s imaging is natural pre-stack migration. This paper now provides the details to the above answers and compares the benefits and liabilities of natural pre-stack migration and natural post-stack migration. In our paper, migration of data to the trial source locations using a natural Green’s function will also be denoted as TRM migration. Both synthetic data and field data examples are provided to validate the effectiveness of superresolution imaging with TRM migration of surface waves. This paper is organized into four sections. Following the introduction we present the theory for TRM migration of surface waves in a homogeneous elastic medium. Formulae predict superresolution imaging with surface waves recorded in the near-field region of the subwavelength scatterer. Instead of the 1/r singularity near the scatterer, the superresolution properties of surface waves result from a −ln (r) singularity; where r is the distance between the near-field scatterer and trial image point. This weaker singularity for surface waves is not surprising because they spread out along a 2-D membrane compared to an expanding sphere for body waves. The theory section is followed by numerical results for migrating surface waves with TRM migration applied to both synthetic and field data. The last section presents the conclusions. 2 THEORY We will first present the geometrical spreading term for the Rayleigh-wave Green’s function for an elastic half-space, which is asymptotically proportional to $$1/ \sqrt{X}$$ and −ln (X) in the far-field and near-field regions of the source, respectively. Here, X is the distance between the source and receiver on the horizontal free surface. We will then use these Green’s functions in the equations for TRMs. 2.1 Rayleigh Green’s Function Zabolotskaya et al. (2012) shows that the surface-wave Green’s function G(g|s)R for a volume point source s ∈ B on the free surface of an elastic half-space and a z-component record at g ∈ B consists of a weighted zeroth-order Hankel function of the second kind:   \begin{eqnarray} G({\bf g}|{\bf s})^R =\alpha H_0^{(2)}(|k_R X|), \end{eqnarray} (1)where   \begin{eqnarray} \alpha =-\frac{\pi ^2 \bar{s} \eta u_0^2}{\rho c_l^2} ik_R^2(1+\eta \xi _l) \quad \mbox{and}\quad X=||{\bf s}-{\bf x}||. \end{eqnarray} (2)Here, R stands for the Rayleigh wave component of the Greens function, ρ is the density, $$\bar{s}$$ is the source strength, $$c_l=\sqrt{\frac{\lambda +2\mu }{\rho }}$$, $$c_t=\sqrt{\frac{\mu }{\rho }}$$, ξ = cR/ct, cR is the Rayleigh velocity, kR = ω/cR, $$\xi _t=\sqrt{1-\xi ^2}$$, $$\eta =-2\frac{\sqrt{1-\xi ^2}}{2-\xi ^2}$$, $$\xi _l=\sqrt{1-\frac{\xi ^2 c_t^2}{c_l^2}}$$, and $$u_0=\frac{1}{2^{1/2} \pi }[\xi _t+1/\xi _t +\eta ^2(\xi _l+1/\xi _l)+4 \eta ]^{-1/2}$$. Scalar wave equation and homogeneous half-space model. Weighting the Rayleigh Green’s function by −iπ/α gives the scaled Green’s function   \begin{eqnarray} G({\bf g}|{\bf s})_{\text{dir}}&=& -\frac{i \pi }{\alpha } G({\bf g}|{\bf s})^R=-i\pi H_0^{(2)}(k_RX), \end{eqnarray} (3)for a homogeneous 3-D half-space model, where ‘dir’ denotes the direct Green’s function. An alternative interpretation is that eq. (3) solves the Helmholtz equation for a line source in a 2-D homogeneous medium (Feshbach & Morse 1953):   \begin{eqnarray} \left[\nabla _1^2 + {\omega ^2}{s_o^2}\right]\overbrace{ \frac{-i \pi }{\alpha }G({\bf x}|{\bf s})^R}^{G({\bf x}|{\bf s})_{\text{dir}}} &=&-4 \pi \delta ({\bf s}-{\bf x}), \end{eqnarray} (4)where $$\nabla _1^2=(\frac{\partial ^2}{\partial x^2},\frac{\partial ^2}{\partial y^2})$$ is the 2-D Laplacian, so = 1/cR is the slowness of the Rayleigh-wave phase velocity at the frequency ω, and s is the location of the harmonic line source. Eq. (4) says that the scaled Rayleigh-wave propagation in a homogeneous half-space at z = 0 is the same as that for an acoustic wave propagating in a 2-D membrane (Tanimoto 1990). The asymptotic far-field (kX ≫ 1) and near-field (kX ≪ 1) representations (Feshbach & Morse 1953)1 of the scaled Rayleigh Green’s function in eq. (3) are   \begin{eqnarray} \lim _{kX \gg 1} G(kX)_{\text{dir}} \sim -\frac{e^{-i(kX+\pi /4)}}{\sqrt{0.5 k X/\pi^{\vphantom{1}} }}, \end{eqnarray} (5)  \begin{eqnarray} \lim _{kX \ll 1}G(kX)_{\text{dir}} \sim -{2} \ln (k X), \end{eqnarray} (6)where kR is replaced by the generic wavenumber symbol k. Eq. (5) is related to Rayleigh-wave propagation in a layered medium. For example, Yomogida & Aki (1985) shows in the far-field approximation that the Rayleigh waves in a smoothly varying background medium2 have a far-field decay similar to that of eq. (5). However, in a smoothly varying medium the lateral distance X in $$1/\sqrt{X}$$ is replaced by a geometrical spreading factor determined by ray tracing on the free surface with the appropriate phase-velocity distribution (Tanimoto 1990). The near-field surface-waves in a layered medium should have a -ln (kX)-like behaviour in the near-field region of the point source (see eq. 6) as long as the shallowest interface is more than a wavelength from the source at the free surface. Scalar wave equation and point-scatterer+half-space model. Assume a point scatterer at xo on the free surface such that the slowness model $$s({\bf x})=s_o +\frac{ 2 \pi \delta s }{s_o}\delta ({\bf x}-{\bf x}_o)$$, where |δs| ≪ 1 is the slowness perturbation. In this case the Green’s function satisfies   \begin{eqnarray} \left[\nabla _1^2 + {\omega ^2}\left({s_o+\frac{ 2 \pi \delta s }{s_o} \delta ({\bf x}-{\bf x}_o)}\right)^2\right]G({\bf x}|{\bf s}) &=&-4 \pi \delta ({\bf s}-{\bf x}), \end{eqnarray} (7)or after rearranging   \begin{eqnarray} \left[\nabla _1^2 + {\omega ^2}{s_o}^2\right]G({\bf x}|{\bf s}) &\approx &-4 \pi \delta ({\bf s}-{\bf x}) -4 \pi \omega ^2 \delta s \delta ({\bf x}-{\bf x}_o)G({\bf x}|{\bf s}),\nonumber\\ \end{eqnarray} (8)where the second-order perturbation term in |δs|2ω2 ≪ 1 is harmlessly ignored. The solution to this equation is found by multiplying by G(x΄|x)dir and integrating x over the free surface to get   \begin{eqnarray} G({\bf x}|{\bf s}) &=&G({\bf x}|{\bf s})_{\text{dir}} + \delta s\omega ^2 G({\bf x}|{\bf x}_o)_{\text{dir}} G({\bf x}_o|{\bf s}),\nonumber \\ &\approx &\overbrace{G({\bf x}|{\bf s})_{\text{dir}}}^{G({\bf x}|{\bf s})_{\text{dir}}=\text{direct Rayleigh wave}} +\overbrace{\omega ^2 \delta s G({\bf x}|{\bf x}_o)_{\text{dir}}G({\bf x}_o|{\bf s})_{\text{dir}}}^{G({\bf x}|{\bf s})_{\text{scatt}}=\text{scattered Rayleigh wave}}, \end{eqnarray} (9)where the scattering is assumed to be small in magnitude so G(xo|s) ≈ G(xo|s)dir under the Born approximation. Rayleigh waves and an arbitrary source distribution. Now assume an arbitrary volume source distribution f (s) on the free surface of the half-space model with a point scatterer near the free surface. In this case the vertical component of the scaled Rayleigh wavefield ϕ(x) satisfies   \begin{eqnarray} \left[\nabla _1^2 + {\omega ^2}{s({\bf x})^2}\right] \phi ({\bf x}) &=&-4 \pi f({\bf s}), \end{eqnarray} (10)where $$s({\bf x})=s_o +\frac{ 2 \pi \delta s }{s_o}\delta ({\bf x}-{\bf x}_o)$$, x = (x, y, 0), and f (s) represents the source distribution at s = (xs, ys, 0) harmonically oscillating with frequency ω and magnitude |f (s)|. Here, the source distribution is confined to the area on the free surface denoted by Bs. The 2-D Green’s function in eq. (9) can be used to solve for ϕ(x) in eq. (10) to give   \begin{eqnarray} {\bf g}\in B_g;\,\, \,\,\,\phi ({\bf g}) &=& \int _{B_s} G({\bf g}|{\bf s}_o) f({\bf s}_o) \text{d}{\bf s}_o^2, \end{eqnarray} (11)where $$\text{d}{\bf s}_o^2$$ indicates $$\text{d}x_{s_o}\text{d}y_{s_o}$$. This equation describes the vertical-component of Rayleigh waves excited by an arbitrary volume source distribution on the surface of a homogeneous elastic half-space with a point scatterer near the free surface at xo. 2.2 Back-scattered surface-wave migration by TRMs The starting point for migrating the back-scattered surface waves is to define the background medium to be the model in Fig. 1 where there is an embedded scatterer at xo near the free surface. The usual goal of the TRM operation (Fink 1993, 2006, 2008) is to image the localized source distribution described by f (s), not the distribution of scatterers as in exploration geophysics (Yilmaz 2001). However, we will see that a subwavelength scatterer (such as a near-surface fault discontinuity) near the source line will manifest its presence by a sharp needle-like image in the migration image. Detecting these scatterers, such as small faults, by migrating the back-scattered surface waves with superresolution imaging is the main goal of our work. Figure 1. View largeDownload slide Plan view of the half-space model where the source, receiver and scatterer are on the free surface. The scatterer at x0 is in the near-field region of the source at s ∈ Bs, and the receivers g ∈ Bg are in the far-field region. The distance between the source at s and the scatterer at x0 is ε. The trial image point at s΄ is along the line denoted by Bs. Figure 1. View largeDownload slide Plan view of the half-space model where the source, receiver and scatterer are on the free surface. The scatterer at x0 is in the near-field region of the source at s ∈ Bs, and the receivers g ∈ Bg are in the far-field region. The distance between the source at s and the scatterer at x0 is ε. The trial image point at s΄ is along the line denoted by Bs. Multiplying eq. (11) by the conjugated Green’s function G(g|s΄)* and integrating over the geophone coordinates g ∈ Bg and the angular-frequency range denoted by Ωo gives the estimated source distribution   \begin{eqnarray} {\bf s}^{\prime } \in B_s;\,\,\,f({\bf s}^{\prime })^{est} &=&\int _{\Omega _o} \int _{B_g} G({\bf g}|{\bf s}^{\prime })^*\phi ({\bf g})\text{d}{\bf g}^2 \text{d} \omega , \end{eqnarray} (12)  \begin{eqnarray} =\int _{\Omega _o} \int _{B_g}\,G({\bf g}|{\bf s}^{\prime })^* \overbrace{{\left [ \int _{B_s} G({\bf g}|{\bf s}_{o})f({\bf s}_{o})d{\bf s}_{o}^{2}\right ] }} ^{m({\bf s},{\bf s}^{\prime })_{\text{dir}-\text{dir}} } \, {\rm d}{\bf g^{2}} {\rm d}\omega , \end{eqnarray} (13)where dg2 indicates dxgdyg and f (s΄)est is the estimated source distribution at the trial image point s΄. Here, the total Green’s function G(g|so) = G(g|so)dir + G(g|so)scatt is the sum of the direct surface-wave and scattered surface-wave Green’s functions in eq. (9). The recorded shot gather is used for the conjugated Green’s function G(g|s΄)*, which naturally backprojects the data from g to s΄ and takes into account all scattering effects in the medium. Unlike standard migration used in the oil industry, no velocity model is needed for the backpropagation operation. The term ϕ(g) in eq. (12) is denoted as post-stack data because the entire source distribution f (so) explodes at time zero and the resulting surface waves are recorded at g to give $$\phi ({\bf g})=\int _{B_s}G({\bf g}|{\bf s}_o)f({\bf s}_o) {}\text{d}{\bf s}_o^2$$. This is equivalent to the exploding reflector concept in exploration geophysics that is used to generate post-stack data (Yilmaz 2001). The key difference is that instead of imaging the source distribution at depth, eq. (13) images the source distribution on the surface. In practice, ϕ(g) is computed by taking all of the shot gathers recorded on the surface, filtering or muting all but the back-scattered surface wave arrivals, and inexpensively stacking the traces from all the sources recorded at g to get ϕ(g). This concept is illustrated in Fig. 2(a) for post-stack forward modelling and post-stack migration in Fig. 2(c). Here, large-amplitude anomalies with high resolution in the scattered migration image at Bs indicate the nearby presence of scatterers. Figure 2. View largeDownload slide Diagrams depict rays associated with (a) post-stack modelling, (b) pre-stack modelling, (c) natural post-stack migration and (d) natural pre-stack migration for a subwavelength scatterer at x in the near-field region of s. The post-stack trace ϕ(g) recorded at g is generated by simultaneously firing all of the shots (red stars); ϕ(g) is then migrated [depicted by dashed arrows in (c)] to the trial image points at s. In contrast, pre-stack traces recorded at g in (b) are generated by firing a point source at s, and then migrating the scattered surface waves in the traces to the trial image points at the source locations s on the source line. This process must be repeated for all the source positions, and then the pre-stack migration images are stacked together. The wiggly rays correspond to the evanescent surface waves that are in the near-field region of the sources, the solid straight rays correspond to surface waves in the far-field region of the geophones, and the dashed lines in (c) and (d) denote the ray paths associated with backpropagated wavefields. Figure 2. View largeDownload slide Diagrams depict rays associated with (a) post-stack modelling, (b) pre-stack modelling, (c) natural post-stack migration and (d) natural pre-stack migration for a subwavelength scatterer at x in the near-field region of s. The post-stack trace ϕ(g) recorded at g is generated by simultaneously firing all of the shots (red stars); ϕ(g) is then migrated [depicted by dashed arrows in (c)] to the trial image points at s. In contrast, pre-stack traces recorded at g in (b) are generated by firing a point source at s, and then migrating the scattered surface waves in the traces to the trial image points at the source locations s on the source line. This process must be repeated for all the source positions, and then the pre-stack migration images are stacked together. The wiggly rays correspond to the evanescent surface waves that are in the near-field region of the sources, the solid straight rays correspond to surface waves in the far-field region of the geophones, and the dashed lines in (c) and (d) denote the ray paths associated with backpropagated wavefields. If the source distribution is restricted to be a point source such that f (so) = δ(s − so), where the point source is in the near-field of the scatterer at xo in Fig. 1, then eq. (13) yields the surface-wave TRM image m(s, s΄) that is similar to the one for body waves (Schuster et al. 2012):   \begin{eqnarray} &&{{\bf s},{\bf s}^{\prime } \in B_s;\,\overbrace{ m({\bf s},{\bf s}^{\prime })}^{f({\bf s}^{\prime })^{\text{est}}}= \int _{\Omega _o} \int _{B_g} G({\bf s}^{\prime }|{\bf g})^*G({\bf g}|{\bf s}) \text{d}{\bf g}^2 \text{d}\omega ,}\nonumber \\ &=&\int _{\Omega _o} \int _{B_g}\big[\, \overbrace{G({\bf s}^{\prime }|{\bf g})_{\text{dir}}^*G({\bf g}|{\bf s})_{\text{dir}}}^{m({\bf s},{\bf s}^{\prime })_{\text{dir}-\text{dir}}} +\overbrace{G({\bf s}^{\prime }|{\bf g})_{\text{scatt}}^* G({\bf g}|{\bf s})_{\text{scatt}}}^{m({\bf s},{\bf s}^{\prime })_{\text{scatt}-\text{scatt}}} \nonumber \\ \nonumber \\ &&+\overbrace{G({\bf s}^{\prime }|{\bf g})_{\text{scatt}}^*G({\bf g}|{\bf s})_{\text{dir}}+ G({\bf s}^{\prime }|{\bf g})_{\text{dir}}^*G({\bf g}|{\bf s})_{\text{scatt}}}^{m({\bf s},{\bf s}^{\prime })_{\text{dir}-\text{scatt}}}\,\big]{}\text{d}{\bf g}^2 \text{d}\omega , \end{eqnarray} (14)where m(s, s΄) defines, for the actual point source at s, the migration image at the trial image location s΄ ∈ Bs. The next sections examine the contributions of the integrands in eq. (14) to superresolution imaging. Direct–direct images. For the Fig. 1 model, the direct-direct imaging term in eq. (14) can be analysed by inserting the far-field direct-wave Green’s function in eq. (5) into the m(s, s΄)dir–dir term in eq. (14) to give   \begin{eqnarray} \lim _{^{\scriptscriptstyle k|{\bf s}-{\bf g}| \gg 1}_{ \scriptscriptstyle k|{\bf s}^{\prime }-{\bf g}| \gg 1}} m({\bf s},{\bf s}^{\prime })_{\text{dir-dir}}\sim \frac{\pi }{0.5} \int _{\Omega _0} \left[ \int _{B_g} \frac{ e^{-ik(|{\bf s}-{\bf g}|-|{\bf s}^{\prime }-{\bf g}|) }}{k\sqrt{|{\bf s}-{\bf g}|}\sqrt{|{\bf s}^{\prime }-{\bf g}|}} {}\text{d}{\bf g}^2 \right]\text{d}\omega . \nonumber \\ \end{eqnarray} (15)Here, the far-field geophones are assumed to be many wavelengths from the sources so that k|s΄ − g|, k|s − g| ≫ 1. The far-field trial image points s΄ ∈ Bs produce kernels in eq. (15) that are too smooth to create sharp superresolution features in the migration image. Scattered–scattered images. For the single-scatterer model in Fig. 1, Appendix A shows that, in the near-field approximation k|s − xo|, k|s΄ − xo| ≪ 1, the scattered-scattered kernel m(s, s΄)scatt–scatt in eq. (14) becomes   \begin{eqnarray} \lim _{^{\scriptscriptstyle k|{\bf x}_o-{\bf g}| \gg 1}_{ \scriptscriptstyle k|{\bf s}^{\prime }-{\bf x}_o| \ll 1}}m({\bf s},{\bf s}^{\prime })_{\scriptscriptstyle \text{scatt-scatt}} &\sim& \alpha + \eta [\ln (|{\bf s}^{\prime }-{\bf x}_o|)+\ln (|{\bf s}-{\bf x}_o|)]\nonumber \\ &&+\,\gamma \ln (|{\bf s}-{\bf x}_o|)\ln (|{\bf s}^{\prime }-{\bf x}_o|)\,, \end{eqnarray} (16)where α, η and γ are constants. Unlike the smooth kernel $$\frac{ 1 }{\sqrt{|{\bf s}-{\bf g}|}\sqrt{|{\bf s}^{\prime }-{\bf g}|}}$$ in eq. (15), the logarithmic singularity ln (|s΄ − xo|) in eq. (16) provides the superresolution characteristics in the migration image for k|s΄ − xo| ≪ 1. Even though the original goal is to image the source distribution f (s) in eq. (13), subwavelength scatterers at xo in the near-field region of the source line Bs are indirectly detected by a strong sharp peak in the migration image m(s, s΄) for s΄ near xo. This is illustrated by the black curve in Fig. 3 which plots $$\text{Im}[H_0^{(2)}(kr)]$$ against −0.6 < kx < 0.6 for different values of the minimum distance ε > 0 between the source line and scatterer at xo. The sharply peaked black curve is what gives rise to superresolution imaging for ε < .001λ, where the black double-sided arrow indicates horizontal resolution better than 1/12λ. On the other hand, if k|x − xo| > 1 then the migration image resolution will be wavelength limited as illustrated by the blue curve in Fig. 3. Figure 3. View largeDownload slide Plot of the imaginary part of the Hankel function $$H_0^{(2)}(kr)$$ against kx, where $$r=\sqrt{(x-x_o)^2+\epsilon ^2}$$, x is the horizontal distance of the trial image point s΄ in Fig. 1 and ε is the minimum distance between the source line and the point scatterer at xo. For the black curve, the lateral resolution limit is indicated by the double-sided black arrow. Figure 3. View largeDownload slide Plot of the imaginary part of the Hankel function $$H_0^{(2)}(kr)$$ against kx, where $$r=\sqrt{(x-x_o)^2+\epsilon ^2}$$, x is the horizontal distance of the trial image point s΄ in Fig. 1 and ε is the minimum distance between the source line and the point scatterer at xo. For the black curve, the lateral resolution limit is indicated by the double-sided black arrow. Direct–scattered images. The direct-scattered kernels in eq. (14) consist of a smooth kernel G(g|s΄)dir (for |s΄ − g| ≫ λ) and a weakly singular kernel G(xo|s΄)scatt (for |s΄ − xo| ≪ λ) for geophones in the far-field region of s and s΄. The weakly singular kernel will provide superresolution peaks in the migration image for the trial image point s΄ near xo. Halo filtering. The goal with TRM is to image near-surface faults with superresolution, and not be distracted by the strong peaky images at each of the source locations. Therefore, traces at g within the near-field region of s΄ and s, i.e. |g − s|, |g − s΄| < 0.5λ, should be muted in order to avoid the direct Rayleigh wave singularities ln |s − g| or ln |s΄ − g| in m(s|s΄)dir–scatt and m(s|s΄)dir–dir. Otherwise there will be sharp strong features in the migration image that will masquerade as false signatures of scatterers. This muting of near-field traces will be accomplished by the halo filter h(s, s΄, g),   $$h({\bf s},{\bf s}^{\prime },{\bf g})=\left\lbrace \begin{array}{cc}1&\quad {\rm if} \,| {\bf s}-{\bf g}|, |{\bf s}^{\prime }-{\bf g}|>0.5 \lambda \\ 0&\quad {\rm if} \,| {\bf s}-{\bf g}|,|{\bf s}^{\prime }-{\bf g}|\le 0.5 \lambda \end{array} ,\right.$$ (17)which can be inserted into the integrand of eq. (13). 2.3 Natural pre-stack migration Eq. (13) is the surface wave TRM equation for imaging the source distribution f (s). It can be considered as a natural post-stack migration operation because ϕ(g) can be interpreted as post-stack data3 associated with the exploding source distribution f (s), and the integration weighted by the extrapolation kernel by G(g|s΄)* is only over the geophone positions on the surface (see Figs 2a and 2c). In contrast, the ray diagram associated with natural pre-stack migration used by AlTheyab et al. (2016) is illustrated in Fig. 2(d) and is mathematically described by   \begin{eqnarray} m ({\bf x}) &=& \int _{{\bf x}_s \in B} \int _{{\bf x}_r\in B}\int 2\omega ^2 W(\omega )^{*}G({\bf x}|{\bf x}_{s})^{*}G({\bf x}|{\bf x}_{r})^{*}\nonumber\\ &&\times \, \overbrace{u({\bf x}_s,{\bf x}_r)}^{\text{pre-stack data}}\text{d}\omega \text{d}x_r \text{d}y_r \text{d}x_s \text{d}y_s, \end{eqnarray} (18)where m(x) is the image at point x, ω is the angular frequency, W(ω) represents the source-wavelet spectrum and * denotes complex conjugation. x, xs and xr are, respectively, the trial-image, source and receiver positions at the surface locations in the set B; G(x|xs) is the Green’s function for the vertical-component source at xs and receiver at x, and G(x|xr) is the Green’s function for a vertical-component particle-velocity recording that only contains the transmitted wavefield without back-scattering. Here, the pre-stack modelling operation is illustrated by the ray diagram in Fig. 2(b), where a single source excites scattered waves that are recorded by the receivers. In this case, the natural pre-stack migration operation for a single frequency is a 4-D integral, two of the integrals are over the source coordinates and two are over the receiver coordinates. Therefore, if there is a grid of N × N sources shooting into a grid of N × N receivers the computational cost will be O(N4) algebraic operations per image point for natural pre-stack migration (Fig. 2d). This means that the total cost will be O(N6) operations for all N2 image points on the surface. In comparison, natural post-stack migration (Fig. 2c) will only cost O(N2) operations for each image point per frequency. In this case the total cost of natural post-stack migration is O(N4) operations for all N2 image points. 3 NUMERICAL RESULTS The surface-wave TRM operations will be tested on both synthetic elastic data and field data collected in Saudi Arabia. The synthetic data are generated by computing solutions to the elastic wave equation, and the input data for TRM are the muted back-scattered surface waves recorded at each vertical-displacement geophone. These tests demonstrate both the effectiveness and limitations of using formula (13) for TRM migration of surface waves. They also validate that post-stack migration of scattered surface waves can lead to superresolution imaging of subwavelength scatterers in the near-field region of the source location. The workflow for implementing post-stack migration of surface waves is illustrated in Fig. 4, which is summarized in the following six steps: Isolate the back-scattered surface waves from other arrivals. In our examples, a muting window is designed according to the estimated phase velocity of the recorded surface waves to mute everything except the back-scattered surface waves. Halo filtering is applied to the recorded data to obtain the Green’s function G(g|s΄)*, where the near-field traces within one dominant wavelength of the source are muted. The halo filter in eq. (17) is applied to the back-scattered surface waves, and stacking over the source index forms the post-stack back-scattered data ϕ(g) for receiver g. Migrate the post-stack back-scattered data ϕ(g) to compute the migration image f (s΄)0 by using eq. (13). A low-pass filter is applied to the post-stack back-scattered data, and then the migration image f (s΄)1 is computed with eq. (13). The enhanced image is obtained by multiplication of the images f (s΄)0 and f (s΄)1. Figure 4. View largeDownload slide Workflow for implementing natural post-stack migration of surface waves. Figure 4. View largeDownload slide Workflow for implementing natural post-stack migration of surface waves. 3.1 2-D buried scatterer A finite-difference method is used to generate 2-D elastic data for the scatterer model in Fig. 5. The model is discretized into 81 × 241 grid-points with the grid size of 0.5 m, where only one subwavelength scatterer is buried at the depth of 3.5 m, which is less than 1/3 of the dominant wavelength. That is, it is within the near-field region of the source and receiver locations directly above the scatterer. There are 120 receivers (z component) located on the surface with a 1 m geophone interval, and there are 120 sources shooting, with one source located at each geophone. The source wavelet is a Ricker wavelet with the peak frequency of 25 Hz, and the sources consist of point displacements in the z-direction. Figure 5. View largeDownload slide The shear velocity model for the 2-D simulations. The background shear velocity is 800 m s−1, the dashed line is the imaginary fault location, and only one subwavelength scatterer (indicated by the grey circle) with shear velocity of 400 m s−1 is embedded around the depth of 3 m. Figure 5. View largeDownload slide The shear velocity model for the 2-D simulations. The background shear velocity is 800 m s−1, the dashed line is the imaginary fault location, and only one subwavelength scatterer (indicated by the grey circle) with shear velocity of 400 m s−1 is embedded around the depth of 3 m. The CSGs (with source location (x, z) = (0,0) m) for the total, direct-transmitted, scattered and muted back-scattered surface waves are shown in Fig. 6. The total (Fig. 6a) and direct-transmitted (Fig. 6b) waves are generated for the scatterer and scatterer-free models in Fig. 5, respectively; the scattered waves (Fig. 6c) are obtained by the subtraction of the total and direct-transmitted waves; and the back-scattered surface waves (Fig. 6d) are muted from the total waves by using an appropriate muting window. These shot gathers are input into migration eq. (13) to give the single-source migration images in Figs 7(a1)–(a4) for different locations of the actual source location s along the free surface. As expected, the scatterer near the source location of sx = 37 m should be imaged in the total (Fig. 7a1), scattered-wave (Fig. 7a3) and muted back-scattered (Fig. 7a4) wave images, and undetectable in the direct–direct image (Fig. 7a2). Note the skinny yellowish band along the diagonals in Figs 7(a1) and (a4). This band of strong migration amplitudes is generated by migration of the direct Rayleigh wave from the source to a near-field geophone. Figure 6. View largeDownload slide Common shot gathers (with source location (x, z) = (0,0) m) for the (a) total, (b) direct-transmitted, (c) scattered and (d) muted back-scattered waves. Traces are for the vertical component of particle velocity. Figure 6. View largeDownload slide Common shot gathers (with source location (x, z) = (0,0) m) for the (a) total, (b) direct-transmitted, (c) scattered and (d) muted back-scattered waves. Traces are for the vertical component of particle velocity. Figure 7. View largeDownload slide Magnified single-source migration images without low-pass filtering for the (a1) total, (a2) direct, (a3) scattered and (a4) muted back-scattered waves. The second row contains the same images except the input data are low-pass filtered between 4 Hz and 18 Hz. The third row contains the enhanced migration images computed from eq. (19). Figure 7. View largeDownload slide Magnified single-source migration images without low-pass filtering for the (a1) total, (a2) direct, (a3) scattered and (a4) muted back-scattered waves. The second row contains the same images except the input data are low-pass filtered between 4 Hz and 18 Hz. The third row contains the enhanced migration images computed from eq. (19). Superresolution behaviour due to scatterers in the near-field of the source is frequency independent (Schuster et al. 2012), which means the same resolution should be achieved when the data are filtered by a low-pass filter. Here, a low-pass filter (4–18 Hz) is applied to the shot gathers in Fig. 6, and then input into the TRM eq. (13) to give the single-source migration images mf (s, s΄) in Figs 7(b1)–(b4). The scatterers near the source location of sx = 37 m are still imaged in the total (Fig. 7b1) and scattered-wave (Fig. 7b3) images. In contrast, the scatterer is slightly smeared in the muted back-scattered image (Fig. 7b4) and undetectable in the direct–direct image (Fig. 7b2). The enhanced TRM image mE(s, s΄) can be obtained by taking the product of the raw TRM image m0(s, s΄) with the low-pass filtered TRM image mf (s, s΄):   $$m_{E}({\bf s},{\bf s}^{\prime }) = m_{0}({\bf s},{\bf s}^{\prime }) \times m_{f}({\bf s},{\bf s}^{\prime }).$$ (19) These enhanced TRM images are shown in Figs 7(c1)–(c4) and more clearly reveal the location of the scatterer without the distraction of artefacts along the diagonal. The migration images m(s, s΄) at the shot location sx = 37 m for different migration kernels are shown in Fig. 8, where the trial shot location is along the horizontal axis and the normalized amplitude is along the vertical axis. Here, the black double arrow represents the dominant shear wavelength, and the blue, red and black curves represent, respectively, the TRM profiles without low-pass filtering, with low-pass filtering data, and the corresponding enhanced image profile. It is obvious that the needle-like profiles of the enhanced images mE(s, s΄)scatt–scatt (black curves in Figs 8c and d) achieve superresolution imaging better than λ/10. However, the resolution of the image from the back-scattered waves (Fig. 8d) is slightly reduced compared with that of the image from the scattered waves (Fig. 8c). In practice, we will mute all events except the back-scattered surface waves to compute the TRM image and use halo filtering to avoid geophone-related singularities in the migration kernels. An alternative to muting is FK filtering. Figure 8. View largeDownload slide The TRM profiles taken from the single-source migration images with the shot offset of 37 m for (a) total, (b) direct, (c) scattered and (d) back-scattered waves. The blue, red and black curves represent, respectively, the TRM profiles without low-pass filtering, with low-pass filtering, and the enhanced image in each subplot. The black double arrow indicates the dominant shear wavelength. Figure 8. View largeDownload slide The TRM profiles taken from the single-source migration images with the shot offset of 37 m for (a) total, (b) direct, (c) scattered and (d) back-scattered waves. The blue, red and black curves represent, respectively, the TRM profiles without low-pass filtering, with low-pass filtering, and the enhanced image in each subplot. The black double arrow indicates the dominant shear wavelength. 3.2 3-D fault model The TRM migration is now tested on 3-D synthetic data associated with the 3-D fault model, where Fig. 9(a) shows the depth slice at z = 2 m. There is an imaginary fault along the diagonal position, and some discretized subwavelength scatterers (the black dots) are embedded along the fault, the shear velocity of the scatterers is 400 m s−1. The model is discretized into 81 × 241 × 241 grid-points in the z, y and x directions, with the grid size of 0.5 m, and the background shear velocity is 600 m s−1. Figure 9. View largeDownload slide (a) Depth slice of the 3-D fault model (z = 2 m), where the background shear velocity is 600 m s−1, there are subwavelength scatterers (the black dots) located along the fault position, and the shear velocity of the scatterers is 400 m s−1 and b) the corresponding enhanced post-stack TRM migration image. The dominant background shear velocity wavelength is 24 m, which is indicated by the double-sided red arrow. Figure 9. View largeDownload slide (a) Depth slice of the 3-D fault model (z = 2 m), where the background shear velocity is 600 m s−1, there are subwavelength scatterers (the black dots) located along the fault position, and the shear velocity of the scatterers is 400 m s−1 and b) the corresponding enhanced post-stack TRM migration image. The dominant background shear velocity wavelength is 24 m, which is indicated by the double-sided red arrow. Figure 10. View largeDownload slide (a) Google map showing the location of the Qademah-fault seismic experiment. (b) Geophone layout for the Qadema-fault seismic experiment, where the green triangles represent the geophones. There are 12 parallel lines, each line has 24 geophones, the inline geophone interval is 5 m and the crossline offset is 10 m. Shots are located at each geophone, and a total of 288 shot gathers are acquired in total. The red dashed line represents the interpreted fault location. Figure 10. View largeDownload slide (a) Google map showing the location of the Qademah-fault seismic experiment. (b) Geophone layout for the Qadema-fault seismic experiment, where the green triangles represent the geophones. There are 12 parallel lines, each line has 24 geophones, the inline geophone interval is 5 m and the crossline offset is 10 m. Shots are located at each geophone, and a total of 288 shot gathers are acquired in total. The red dashed line represents the interpreted fault location. There are 3600 receivers (z-component) evenly distributed on the surface, with a 2 m spacing interval in both the x- and y-directions. The data are generated by computing solutions to the 3-D elastic wave equation with a Ricker source wavelet having the peak frequency of 25 Hz. Point sources are placed at each geophone location with a z-component displacement to generate 3600 shot gathers, each with 3600 traces. An appropriate muting window is designed, based on the phase velocity of the surface wave for each shot, to extract the back-scattered surface waves. In addition, a Halo filter is used to mute the near-field traces within one dominant wavelength of the source. The muted back-scattered surface waves are used to compute the post-stack TRM migration image in Fig. 9(b). Here, the workflow in Fig. 4 is used to compute the post-stack image and the strong migration amplitudes clearly reveal the locations of the scatterers at subwavelength resolution. 3.3 3-D Qademah fault data A 3-D land survey was carried out along the Red Sea coast over the Qademah fault system, about 30 km north of the KAUST campus (Hanafy 2015). There are 288 vertical-component geophones arranged in 12 parallel lines, and each line has 24 geophones. The inline geophone interval is 5 m and the crossline interval is 10 m. The geophone geometry layout is shown in Fig. 10(a), where one shot is fired at each geophone location for a total of 288 shot gathers. The source is generated by a 200 lb hammer striking a metal plate on the ground. The red dashed line represents the interpreted fault location from geological evidence (Roobol & Kadi 2008). The post-stack TRM operation can be used to image near-field subwavelength scatterers,where the scatterers at shallow depth are sensitive to higher frequency data, while deeper scatterers are mostly sensitive to lower frequency data. First, we apply different narrowband filters with different peak frequencies (from 13 to 41 Hz, each with a 4 Hz increase in the frequency band) to the raw Qademah data. A CSG (no. 120) filtered by a narrowband filter with the peak frequency of 25 Hz is shown in Fig. 11(a), where the red dashed curves are computed from the estimated phase velocity, and they are used to mute the back-scattered surface waves from the total wavefield. The muted back-scattered surface waves shown in Fig. 11(b) are the input into the TRM operation in eq. (13) for computation of the post-stack TRM image. For comparison, the data in Fig. 11(b) are also used to calculate the natural pre-stack migration images (Liu et al. 2017). Figure 11. View largeDownload slide Common shot gather no. 120 from the Qademah fault data filtered by a bandpass filter with 25 Hz peak frequency, (a) the total wavefield, where the red dashed curve is used to design a muting window for muting the back-scattered waves, and (b) the muted back-scattered waves. Figure 11. View largeDownload slide Common shot gather no. 120 from the Qademah fault data filtered by a bandpass filter with 25 Hz peak frequency, (a) the total wavefield, where the red dashed curve is used to design a muting window for muting the back-scattered waves, and (b) the muted back-scattered waves. According to Hanafy (2015)’s interpretation, the fault ‘F1’ (at an offset ≈90 m in the inline direction) in Fig. 12(c) is the north-south Qademah fault, which is believed to be a normal rotational fault resulting from the opening of the Red Sea, the fault ‘F2’ in the westside of ‘F1’ is interpreted as an antithetic fault. The area between faults ‘F1’ and ‘F2’ are composed of loose sediments and colluvial rocks. The surface waves scattered by these colluvial rocks with size less than the dominant shear wavelength are to be imaged by the TRM operation. Figure 12. View largeDownload slide (a) The natural post-stack migration (TRM) image, (b) natural pre-stack migration (Liu et al. 2017) (NPM) image for the Qademah data with different dominant frequencies and (c) Rayleigh phase-velocity tomogram (Hanafy 2015) at different frequencies. Here, the red double-sided arrows indicate the approximate dominant shear wavelength at different frequencies. The dashed pink curves in (a) labelled with ‘S1’, ‘S2’ and ‘S3’ indicate three segments of subwavelength scatterers between faults ‘F1’ and ‘F2’. The dashed pink curves in (b) indicate the location of Qademah fault ‘F1’. These two white dashed curves in (c) labelled with ‘F1’ and ‘F2’ represent the Qademah fault and another small antithetic fault, respectively. Figure 12. View largeDownload slide (a) The natural post-stack migration (TRM) image, (b) natural pre-stack migration (Liu et al. 2017) (NPM) image for the Qademah data with different dominant frequencies and (c) Rayleigh phase-velocity tomogram (Hanafy 2015) at different frequencies. Here, the red double-sided arrows indicate the approximate dominant shear wavelength at different frequencies. The dashed pink curves in (a) labelled with ‘S1’, ‘S2’ and ‘S3’ indicate three segments of subwavelength scatterers between faults ‘F1’ and ‘F2’. The dashed pink curves in (b) indicate the location of Qademah fault ‘F1’. These two white dashed curves in (c) labelled with ‘F1’ and ‘F2’ represent the Qademah fault and another small antithetic fault, respectively. Fig. 12(a) only shows the post-stack TRM images (from the top-to-bottom panels) calculated from the filtered back-scattered data with dominant frequencies of 29, 25 and 21 Hz, respectively. The dashed pink curve (on the top panel) labelled with ‘S1’ indicates the location of the first segment of the fault, and the estimated depth of these subwavelength scatterers is around 5 m, which is roughly estimated as one third of the corresponding shear wavelength (Stokoe & Nazarian 1985). The post-stack TRM image in the middle panel of Fig. 12(a) is calculated from 25 Hz back-scattered data, where the dashed pink curve labelled with ‘S2’ indicates the location of the second segment of the fault. We found that there is an overlap between the ‘S1’ and ‘S2’ images, which indicates that the scatterers which belong to ‘S1’ are also imaged by 25 Hz data. The strong amplitudes in ‘S2’ reveal that the approximate depth of these subwavelength scatterers is around 6.7 m. The bottom panel in Fig. 12(a) shows the post-stack TRM image calculated from 21 Hz back-scattered data, the dashed pink curve labelled as ‘S3’ indicates the location of the third segment of the fault, and the strong amplitudes suggest that the approximate depth of these subwavelength scatterers is around 10 m. The red double-sided arrows represent the dominant shear wavelength in each panel, where the horizontal resolution of the fault images is less than half the dominant shear wavelength. For comparison, the natural pre-stack migration (NPM) images calculated by eq. (18) are computed by Liu et al. (2017) and shown in Fig. 12(b). The NPM images are computed from the filtered back-scattered data with dominant frequencies of 29, 25 and 21 Hz, respectively, for the top-to-bottom panels. The dashed pink curve in each panel indicates the location of the Qademah Fault (‘F1’). Figs 13(a)–(c) show the TRM and NPM image profiles along ‘A1 − B1’, ‘A2 − B2’ and ‘A3 − B3’, respectively. Here, the red-dashed and blue-solid curves represent the TRM and NPM profiles, respectively. And the black double-sided arrows indicate the dominant shear wavelength in each subplot. We can see the TRM image profile has achieved a resolution better than one half of the dominant shear wavelength, while the resolution of NPM is around one dominant shear wavelength. This difference in resolution limits for pre-stack and post-stack migration are roughly in agreement with the resolution limits predicted by (Chen & Schuster 1999). Figure 13. View largeDownload slide The profile comparisons between the post-stack TRM (Fig. 12a) and NPM (Fig. 12b) images along (a) ‘A1 − B1’, (b) ‘A2 − B2’ and (c) ‘A3 − B3’. The red dashed and blue solid curves represent image results from TRM and NPM, respectively. The black double-side arrows indicate the dominant shear wavelength in each subplot, and amplitude has been normalized. Figure 13. View largeDownload slide The profile comparisons between the post-stack TRM (Fig. 12a) and NPM (Fig. 12b) images along (a) ‘A1 − B1’, (b) ‘A2 − B2’ and (c) ‘A3 − B3’. The red dashed and blue solid curves represent image results from TRM and NPM, respectively. The black double-side arrows indicate the dominant shear wavelength in each subplot, and amplitude has been normalized. 4 SUMMARY AND CONCLUSIONS We theoretically and numerically analysed the efficacy of surface-wave TRM migration. Theoretical formulae and synthetic data results show that surface waves can be used to achieve superresolution imaging of subwavelength scatterers if they are located less than about 1/3 of a shear wavelength from the source line. We also show that the TRM operation can be used for natural post-stack migration, and only costs O(N4) algebraic operation per frequency at the N2 trial image points on the surface. This compares to O(N6) operations for natural pre-stack migration at a single frequency. Here, we assume the sources and receivers are on a N × N grid on the free surface and there are N2 trial image points on the surface. These theoretical formulae show, for the first time, the connections between pre-stack natural migration, post-stack natural migration and the TRM operations. It also shows the potential for achieving superresolution imaging of scattered surface waves in the near-field region of the source. The field-data tests suggest that hidden faults can be detected with the post-stack TRM operation if they are no deeper than about 1/2 the dominant wavelength. These tests also show images with resolution that is about 1/3 that of the dominant wavelength. The most significant challenge in natural post-stack migration is to process the data so that only the back-scattered surface waves are used. This is a subject of ongoing research. ACKNOWLEDGEMENTS The research reported in this publication was supported by the King Abdullah University of Science and Technology (KAUST) in Thuwal, Saudi Arabia and the CRG grant. We are grateful to the sponsors of the Center for Subsurface Imaging and Modeling (CSIM) Consortium for their financial support. For computer time, this research used the resources of the Supercomputing Laboratory at KAUST and the IT Research Computing Group. We thank them for providing the computational resources required for carrying out this work. We greatly appreciate the constructive comments and suggestions from editor Herve Chauris, and two anonymous reviewers, which helped improve this paper. Footnotes 1 See eq. 7.2.18 in Feshbach & Morse (1953), where $$H_0^{(1)}(kX)=H_0^{(2)}(kX)^*$$. 2 That is, the horizontal variations are small on the scale of the largest wavelength. 3 The post-stack traces ϕ(g) result from exploding the many different sources at s ∈ Bs with strength f (s). REFERENCES AlTheyab A., Lin F.-C., Schuster G.T., 2016. Imaging near-surface heterogeneities by natural migration of backscattered surface waves, Geophys. J. Int. , 204( 2), 1332– 1341. Google Scholar CrossRef Search ADS   Chen J., Schuster G.T., 1999. Resolution limits of migrated images, Geophysics , 64( 4), 1046– 1053. Google Scholar CrossRef Search ADS   De Fornel F., 2001. Evanescent Waves: from Newtonian Optics to Atomic Optics , Vol. 73, Springer Science & Business Media. Dutta G., AlTheyab A., Tarhini A., Hanafy S., Schuster G.T., 2016. Extracting 220 Hz information from 55 Hz field data by near-field superresolution imaging, Geophys. J. Int. , 206, 1194– 1203. Google Scholar CrossRef Search ADS   Feshbach H., Morse P.M., 1953. Methods of Theoretical Physics , McGraw-Hill Interamericana. Fink M., 1993. Time-reversal mirrors, J. Phys. D: Appl. Phys. , 26( 9), 1333– 1350. Google Scholar CrossRef Search ADS   Fink M., 2006. Time-reversal acoustics in complex environments, Geophysics , 71( 4), SI151– SI164. Google Scholar CrossRef Search ADS   Fink M., 2008. Time-reversal waves and super resolution, J. Phys.: Conf. Ser. , 124, 012004, pp. 1–29. Google Scholar CrossRef Search ADS   Hanafy S.M., 2015. Mapping the Qademah fault with traveltime, surface-wave, and resistivity tomograms, SEG Technical Program Expanded Abstracts , pp. 3347– 3351. Jia H., Ke M., Hao R., Ye Y., Liu F., Liu Z., 2010. Subwavelength imaging by a simple planar acoustic superlens, Appl. Phys. Lett. , 97( 17), 173507, doi:10.1063/1.3507893. Google Scholar CrossRef Search ADS   Lerosey G., de Rosny J., Tourin A., Fink M., 2007. Focusing beyond the diffraction limit with far-field time reversal, Science , 315( 5815), 1120– 1122. Google Scholar CrossRef Search ADS PubMed  Liu Z., Altheyab A., Hanafy S., Schuster G., 2017. Imaging near-surface heterogeneities by natural migration of surface waves: field data test, Geophysics , 82( 3), 1– 9. Google Scholar CrossRef Search ADS   Roobol M., Kadi K., 2008. Cenozoic faulting in the Rabigh area, central west Saudi Arabia (including the sites of King Abdullah Economic City and King Abdullah University for science and technology), Saudi Geological Survey Technical Report SGS-TR-2008-6, 1(250,000). Schuster G.T., 2002. Reverse-time migration= generalized diffraction stack migration, in SEG Technical Program Expanded Abstracts 2002 , pp. 1280– 1283, Society of Exploration Geophysicists. Schuster G.T., Hanafy S., Huang Y., 2012. Theory and feasibility tests for a seismic scanning tunnelling macroscope, Geophys. J. Int. , 190( 3), 1593– 1606. Google Scholar CrossRef Search ADS   Stokoe K.H., Nazarian S., 1985. Use of Rayleigh waves in liquefaction studies, in Measurement and Use of Shear Wave Velocity for Evaluating Dynamic Soil Properties , p. 17, ASCE. Tanimoto T., 1990. Modelling curved surface wave paths: membrane surface wave synthetics, Geophys. J. Int. , 102( 1), 89– 100. Google Scholar CrossRef Search ADS   Tarhini A., Guo B., Dutta G., Schuster G.T., 2017. Extension of seismic scanning tunnelling macroscope to elastic waves, Pure Appl. Geophys. , doi:10.1007/s00024-017-1692-x. Yilmaz Ö., 2001. Seismic Data Analysis , Society of Exploration Geophysicists. Google Scholar CrossRef Search ADS   Yomogida K., Aki K., 1985. Waveform synthesis of surface waves in a laterally heterogeneous earth by the Gaussian beam method, J. geophys. Res. , 90( B9), 7665– 7688. Google Scholar CrossRef Search ADS   Zabolotskaya E.A., Ilinskii Y.A., Hay T.A., Hamilton M.F., 2012. Green’s functions for a volume source in an elastic half-space, J. acoust. Soc. Am. , 131( 3), 1831– 1842. Google Scholar CrossRef Search ADS PubMed  APPENDIX A: SCATTERED-SCATTERED IMAGE For the Fig. 1 model, assume that the scatterer at xo is in the near-field region of both s and s΄, and it is also in the far-field region of the geophone at g. Using eqs (5) and (6) and the definitions in eq. (9) gives the asymptotic expression for the scattered Green’s functions in eq. (14):   \begin{eqnarray} \lim _{^{\scriptscriptstyle k|{\bf x}_o-{\bf g}| \gg 1}_{ \scriptscriptstyle k|{\bf s}-{\bf x}_o| \ll 1}}G({\bf g}|{\bf s})_{\text{scatt}}&=&-\pi ^2 R\omega ^2 \,\,\lim _{^{\scriptscriptstyle k|{\bf x}_o-{\bf g}| \gg 1}_{ \scriptscriptstyle k|{\bf s}-{\bf x}_o| \ll 1}}H_o^{(2)}(k|{\bf s}-{\bf x}_o|)H_o^{(2)}(k|{\bf x}_o-{\bf g}|), \nonumber\\ &\sim & -2R \omega ^2 \ln ( k|{\bf s}-{\bf x}_o|) \frac{e^{-i(k|{\bf g}-{\bf x}_o|+\pi /4)}}{\sqrt{0.5 k |{\bf x}_o-{\bf g}|/\pi }}, \end{eqnarray} (A1)  \begin{eqnarray} \lim _{^{\scriptscriptstyle k|{\bf x}_o-{\bf g}| \gg 1}_{ \scriptscriptstyle k|{\bf s}^{\prime }-{\bf x}_o| \ll 1}}G({\bf s}^{\prime }|{\bf g})_{\text{scatt}}^* &\sim & -2R \omega ^2\ln ( k|{\bf s}^{\prime }-{\bf x}_o|) \frac{e^{i(k|{\bf g}-{\bf x}_o|+\pi /4)}}{\sqrt{0.5 k |{\bf x}_o-{\bf g}|/\pi }}, \end{eqnarray} (A2)where R = δs is the scattering coefficient, k = ω/c, and c is the Rayleigh-wave phase velocity. Conveniently setting $$R=\frac{1}{\sqrt{8 \pi }}$$, the scattered-scattered image m(s, s΄)scatt–scatt in eq. (14) becomes   \begin{eqnarray} m({\bf s},{\bf s}^{\prime })_{\text{scatt--scatt}}&=& \int _{\Omega _0} \left[ \int _{B_g} G({\bf s}^{\prime }|{\bf g})_{\text{scatt}}^*G({\bf g}|{\bf s})_{\text{scatt}} \text{d}{\bf g}^2 \right]\text{d}\omega ,\nonumber \\ &=& \beta \int _{\Omega _k} k^3 { \ln (k |{\bf s}^{\prime }-{\bf x}_o|)\ln (k |{\bf s}-{\bf x}_o|) } \text{d}k, \end{eqnarray} (A3)where the integration over frequency has been transformed to an integration over the wavenumber region denoted by Ωk and $$\beta =\int _{B_g} \frac{c^5}{ |{\bf x}_o-{\bf g}|}\text{d}{\bf g}^2$$. Here, s, s΄ ∈ Bs and the integration over geophone positions is along the set of points denoted by Bg. Eq. (A3) can be expressed as   \begin{eqnarray*} m({\bf s},{\bf s}^{\prime })_{\scriptscriptstyle \text{scatt--scatt}}&=&\beta \int _k {\text{d}k}\,k^3\lbrace {{(\ln k)^2 +\ln k[\ln (|{\bf s}^{\prime }-{\bf x}_o|)+\ln (|{\bf s}-{\bf x}_o|)]+\ln (|{\bf s}^{\prime }-{\bf x}_o|)\ln ( |{\bf s}-{\bf x}_o|)}{}}\rbrace ,\nonumber \\ &=&\alpha + \eta [\ln (|{\bf s}^{\prime }-{\bf x}_o|)+\ln (|{\bf s}-{\bf x}_o|)]+\gamma \ln (|{\bf s}-{\bf x}_o|)\ln (|{\bf s}^{\prime }-{\bf x}_o|)\,, \end{eqnarray*} where α = β∫kk3(ln k)2dk, η = β∫kk3ln kdk, and γ = β∫k3dk are finite-valued numbers that depend on the range of the integration over wavenumbers. The weakly singular ln |s΄ − xo| term is responsible for near-field superresolution imaging of surface waves for trial image points s΄ near the scatterer at xo. © 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

# Superresolution near-field imaging with surface waves

, Volume 212 (2) – Feb 1, 2018
12 pages

/lp/ou_press/superresolution-near-field-imaging-with-surface-waves-oSIHilNcBj
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx466
Publisher site
See Article on Publisher Site

### Abstract

Summary We present the theory for near-field superresolution imaging with surface waves and time reverse mirrors (TRMs). Theoretical formulae and numerical results show that applying the TRM operation to surface waves in an elastic half-space can achieve superresolution imaging of subwavelength scatterers if they are located less than about 1/2 of the shear wavelength from the source line. We also show that the TRM operation for a single frequency is equivalent to natural migration, which uses the recorded data to approximate the Green’s functions for migration, and only costs O(N4) algebraic operations for post-stack migration compared to O(N6) operations for natural pre-stack migration. Here, we assume the sources and receivers are on an N × N grid and there are N2 trial image points on the free surface. Our theoretical predictions of superresolution are validated with tests on synthetic data. The field-data tests suggest that hidden faults at the near surface can be detected with subwavelength imaging of surface waves by using the TRM operation if they are no deeper than about 1/2 the dominant shear wavelength. Image processing, Surface waves and free oscillations 1 INTRODUCTION Superresolution imaging is defined as estimating the location of an object to a resolution much better than λ/2, where λ is the effective wavelength of the wavefield. For optical wavefields, De Fornel (2001) and Jia et al. (2010) utilized evanescent energy in the near-field for their high-resolution optical superlenses. Similarly, Lerosey et al. (2007) demonstrated that superresolution imaging is also possible with microwaves recorded in the far-field region of subwavelength scatterers, when the scatterers are near the source. Their microwave experiment with subwavelength scatterers near the source showed that imaging resolution could be better than λ/30 by the use of a time reversal mirror (TRM) operation (Fink 1993, 2008) applied to the recorded microwave data. Schuster et al. (2012) analysed the TRM operation in the context of seismic scattering and theoretically showed that near-field seismic superresolution is enabled by the evanescent field associated with the geometrical-spreading factor. Their theoretical treatment of the evanescent field was restricted to body-wave arrivals in the acoustic approximation. Both simulations and field experiments (Dutta et al. 2016) validated the theoretical predictions. Their imaging method, denoted as the seismic scanning tunnelling macroscope (SSTM) operation, was extended to elastic body waves by Tarhini et al. (2017), and the superresolution properties were validated by tests on both synthetic data and field data. The TRM operation is equivalent to natural migration (Schuster 2002; Schuster et al. 2012; AlTheyab et al. 2016) because it uses the recorded Green’s functions for migration and so makes no assumptions about the velocity model. These ‘natural’ Green’s functions eliminate the need for knowing the velocity model. AlTheyab et al. (2016) and Liu et al. (2017) used natural pre-stack migration of scattered surface waves to estimate the locations of hidden faults that were within λ/3 of the free surface. Migrating surface waves is a robust imaging strategy because surface waves are more than an order-of-magnitude stronger than most body-wave reflections. There are two important questions posed by the AlTheyab et al. (2016) study with natural migration. Can surface waves be used to achieve superresolution imaging? The answer is yes as we will show both theoretically and numerically. How is post-stack migration related to the natural pre-stack migration of AlTheyab et al. (2016)? As will be shown, TRM imaging of data generated by simultaneously exploding sources is equivalent to natural post-stack migration while Altheyab’s imaging is natural pre-stack migration. This paper now provides the details to the above answers and compares the benefits and liabilities of natural pre-stack migration and natural post-stack migration. In our paper, migration of data to the trial source locations using a natural Green’s function will also be denoted as TRM migration. Both synthetic data and field data examples are provided to validate the effectiveness of superresolution imaging with TRM migration of surface waves. This paper is organized into four sections. Following the introduction we present the theory for TRM migration of surface waves in a homogeneous elastic medium. Formulae predict superresolution imaging with surface waves recorded in the near-field region of the subwavelength scatterer. Instead of the 1/r singularity near the scatterer, the superresolution properties of surface waves result from a −ln (r) singularity; where r is the distance between the near-field scatterer and trial image point. This weaker singularity for surface waves is not surprising because they spread out along a 2-D membrane compared to an expanding sphere for body waves. The theory section is followed by numerical results for migrating surface waves with TRM migration applied to both synthetic and field data. The last section presents the conclusions. 2 THEORY We will first present the geometrical spreading term for the Rayleigh-wave Green’s function for an elastic half-space, which is asymptotically proportional to $$1/ \sqrt{X}$$ and −ln (X) in the far-field and near-field regions of the source, respectively. Here, X is the distance between the source and receiver on the horizontal free surface. We will then use these Green’s functions in the equations for TRMs. 2.1 Rayleigh Green’s Function Zabolotskaya et al. (2012) shows that the surface-wave Green’s function G(g|s)R for a volume point source s ∈ B on the free surface of an elastic half-space and a z-component record at g ∈ B consists of a weighted zeroth-order Hankel function of the second kind:   \begin{eqnarray} G({\bf g}|{\bf s})^R =\alpha H_0^{(2)}(|k_R X|), \end{eqnarray} (1)where   \begin{eqnarray} \alpha =-\frac{\pi ^2 \bar{s} \eta u_0^2}{\rho c_l^2} ik_R^2(1+\eta \xi _l) \quad \mbox{and}\quad X=||{\bf s}-{\bf x}||. \end{eqnarray} (2)Here, R stands for the Rayleigh wave component of the Greens function, ρ is the density, $$\bar{s}$$ is the source strength, $$c_l=\sqrt{\frac{\lambda +2\mu }{\rho }}$$, $$c_t=\sqrt{\frac{\mu }{\rho }}$$, ξ = cR/ct, cR is the Rayleigh velocity, kR = ω/cR, $$\xi _t=\sqrt{1-\xi ^2}$$, $$\eta =-2\frac{\sqrt{1-\xi ^2}}{2-\xi ^2}$$, $$\xi _l=\sqrt{1-\frac{\xi ^2 c_t^2}{c_l^2}}$$, and $$u_0=\frac{1}{2^{1/2} \pi }[\xi _t+1/\xi _t +\eta ^2(\xi _l+1/\xi _l)+4 \eta ]^{-1/2}$$. Scalar wave equation and homogeneous half-space model. Weighting the Rayleigh Green’s function by −iπ/α gives the scaled Green’s function   \begin{eqnarray} G({\bf g}|{\bf s})_{\text{dir}}&=& -\frac{i \pi }{\alpha } G({\bf g}|{\bf s})^R=-i\pi H_0^{(2)}(k_RX), \end{eqnarray} (3)for a homogeneous 3-D half-space model, where ‘dir’ denotes the direct Green’s function. An alternative interpretation is that eq. (3) solves the Helmholtz equation for a line source in a 2-D homogeneous medium (Feshbach & Morse 1953):   \begin{eqnarray} \left[\nabla _1^2 + {\omega ^2}{s_o^2}\right]\overbrace{ \frac{-i \pi }{\alpha }G({\bf x}|{\bf s})^R}^{G({\bf x}|{\bf s})_{\text{dir}}} &=&-4 \pi \delta ({\bf s}-{\bf x}), \end{eqnarray} (4)where $$\nabla _1^2=(\frac{\partial ^2}{\partial x^2},\frac{\partial ^2}{\partial y^2})$$ is the 2-D Laplacian, so = 1/cR is the slowness of the Rayleigh-wave phase velocity at the frequency ω, and s is the location of the harmonic line source. Eq. (4) says that the scaled Rayleigh-wave propagation in a homogeneous half-space at z = 0 is the same as that for an acoustic wave propagating in a 2-D membrane (Tanimoto 1990). The asymptotic far-field (kX ≫ 1) and near-field (kX ≪ 1) representations (Feshbach & Morse 1953)1 of the scaled Rayleigh Green’s function in eq. (3) are   \begin{eqnarray} \lim _{kX \gg 1} G(kX)_{\text{dir}} \sim -\frac{e^{-i(kX+\pi /4)}}{\sqrt{0.5 k X/\pi^{\vphantom{1}} }}, \end{eqnarray} (5)  \begin{eqnarray} \lim _{kX \ll 1}G(kX)_{\text{dir}} \sim -{2} \ln (k X), \end{eqnarray} (6)where kR is replaced by the generic wavenumber symbol k. Eq. (5) is related to Rayleigh-wave propagation in a layered medium. For example, Yomogida & Aki (1985) shows in the far-field approximation that the Rayleigh waves in a smoothly varying background medium2 have a far-field decay similar to that of eq. (5). However, in a smoothly varying medium the lateral distance X in $$1/\sqrt{X}$$ is replaced by a geometrical spreading factor determined by ray tracing on the free surface with the appropriate phase-velocity distribution (Tanimoto 1990). The near-field surface-waves in a layered medium should have a -ln (kX)-like behaviour in the near-field region of the point source (see eq. 6) as long as the shallowest interface is more than a wavelength from the source at the free surface. Scalar wave equation and point-scatterer+half-space model. Assume a point scatterer at xo on the free surface such that the slowness model $$s({\bf x})=s_o +\frac{ 2 \pi \delta s }{s_o}\delta ({\bf x}-{\bf x}_o)$$, where |δs| ≪ 1 is the slowness perturbation. In this case the Green’s function satisfies   \begin{eqnarray} \left[\nabla _1^2 + {\omega ^2}\left({s_o+\frac{ 2 \pi \delta s }{s_o} \delta ({\bf x}-{\bf x}_o)}\right)^2\right]G({\bf x}|{\bf s}) &=&-4 \pi \delta ({\bf s}-{\bf x}), \end{eqnarray} (7)or after rearranging   \begin{eqnarray} \left[\nabla _1^2 + {\omega ^2}{s_o}^2\right]G({\bf x}|{\bf s}) &\approx &-4 \pi \delta ({\bf s}-{\bf x}) -4 \pi \omega ^2 \delta s \delta ({\bf x}-{\bf x}_o)G({\bf x}|{\bf s}),\nonumber\\ \end{eqnarray} (8)where the second-order perturbation term in |δs|2ω2 ≪ 1 is harmlessly ignored. The solution to this equation is found by multiplying by G(x΄|x)dir and integrating x over the free surface to get   \begin{eqnarray} G({\bf x}|{\bf s}) &=&G({\bf x}|{\bf s})_{\text{dir}} + \delta s\omega ^2 G({\bf x}|{\bf x}_o)_{\text{dir}} G({\bf x}_o|{\bf s}),\nonumber \\ &\approx &\overbrace{G({\bf x}|{\bf s})_{\text{dir}}}^{G({\bf x}|{\bf s})_{\text{dir}}=\text{direct Rayleigh wave}} +\overbrace{\omega ^2 \delta s G({\bf x}|{\bf x}_o)_{\text{dir}}G({\bf x}_o|{\bf s})_{\text{dir}}}^{G({\bf x}|{\bf s})_{\text{scatt}}=\text{scattered Rayleigh wave}}, \end{eqnarray} (9)where the scattering is assumed to be small in magnitude so G(xo|s) ≈ G(xo|s)dir under the Born approximation. Rayleigh waves and an arbitrary source distribution. Now assume an arbitrary volume source distribution f (s) on the free surface of the half-space model with a point scatterer near the free surface. In this case the vertical component of the scaled Rayleigh wavefield ϕ(x) satisfies   \begin{eqnarray} \left[\nabla _1^2 + {\omega ^2}{s({\bf x})^2}\right] \phi ({\bf x}) &=&-4 \pi f({\bf s}), \end{eqnarray} (10)where $$s({\bf x})=s_o +\frac{ 2 \pi \delta s }{s_o}\delta ({\bf x}-{\bf x}_o)$$, x = (x, y, 0), and f (s) represents the source distribution at s = (xs, ys, 0) harmonically oscillating with frequency ω and magnitude |f (s)|. Here, the source distribution is confined to the area on the free surface denoted by Bs. The 2-D Green’s function in eq. (9) can be used to solve for ϕ(x) in eq. (10) to give   \begin{eqnarray} {\bf g}\in B_g;\,\, \,\,\,\phi ({\bf g}) &=& \int _{B_s} G({\bf g}|{\bf s}_o) f({\bf s}_o) \text{d}{\bf s}_o^2, \end{eqnarray} (11)where $$\text{d}{\bf s}_o^2$$ indicates $$\text{d}x_{s_o}\text{d}y_{s_o}$$. This equation describes the vertical-component of Rayleigh waves excited by an arbitrary volume source distribution on the surface of a homogeneous elastic half-space with a point scatterer near the free surface at xo. 2.2 Back-scattered surface-wave migration by TRMs The starting point for migrating the back-scattered surface waves is to define the background medium to be the model in Fig. 1 where there is an embedded scatterer at xo near the free surface. The usual goal of the TRM operation (Fink 1993, 2006, 2008) is to image the localized source distribution described by f (s), not the distribution of scatterers as in exploration geophysics (Yilmaz 2001). However, we will see that a subwavelength scatterer (such as a near-surface fault discontinuity) near the source line will manifest its presence by a sharp needle-like image in the migration image. Detecting these scatterers, such as small faults, by migrating the back-scattered surface waves with superresolution imaging is the main goal of our work. Figure 1. View largeDownload slide Plan view of the half-space model where the source, receiver and scatterer are on the free surface. The scatterer at x0 is in the near-field region of the source at s ∈ Bs, and the receivers g ∈ Bg are in the far-field region. The distance between the source at s and the scatterer at x0 is ε. The trial image point at s΄ is along the line denoted by Bs. Figure 1. View largeDownload slide Plan view of the half-space model where the source, receiver and scatterer are on the free surface. The scatterer at x0 is in the near-field region of the source at s ∈ Bs, and the receivers g ∈ Bg are in the far-field region. The distance between the source at s and the scatterer at x0 is ε. The trial image point at s΄ is along the line denoted by Bs. Multiplying eq. (11) by the conjugated Green’s function G(g|s΄)* and integrating over the geophone coordinates g ∈ Bg and the angular-frequency range denoted by Ωo gives the estimated source distribution   \begin{eqnarray} {\bf s}^{\prime } \in B_s;\,\,\,f({\bf s}^{\prime })^{est} &=&\int _{\Omega _o} \int _{B_g} G({\bf g}|{\bf s}^{\prime })^*\phi ({\bf g})\text{d}{\bf g}^2 \text{d} \omega , \end{eqnarray} (12)  \begin{eqnarray} =\int _{\Omega _o} \int _{B_g}\,G({\bf g}|{\bf s}^{\prime })^* \overbrace{{\left [ \int _{B_s} G({\bf g}|{\bf s}_{o})f({\bf s}_{o})d{\bf s}_{o}^{2}\right ] }} ^{m({\bf s},{\bf s}^{\prime })_{\text{dir}-\text{dir}} } \, {\rm d}{\bf g^{2}} {\rm d}\omega , \end{eqnarray} (13)where dg2 indicates dxgdyg and f (s΄)est is the estimated source distribution at the trial image point s΄. Here, the total Green’s function G(g|so) = G(g|so)dir + G(g|so)scatt is the sum of the direct surface-wave and scattered surface-wave Green’s functions in eq. (9). The recorded shot gather is used for the conjugated Green’s function G(g|s΄)*, which naturally backprojects the data from g to s΄ and takes into account all scattering effects in the medium. Unlike standard migration used in the oil industry, no velocity model is needed for the backpropagation operation. The term ϕ(g) in eq. (12) is denoted as post-stack data because the entire source distribution f (so) explodes at time zero and the resulting surface waves are recorded at g to give $$\phi ({\bf g})=\int _{B_s}G({\bf g}|{\bf s}_o)f({\bf s}_o) {}\text{d}{\bf s}_o^2$$. This is equivalent to the exploding reflector concept in exploration geophysics that is used to generate post-stack data (Yilmaz 2001). The key difference is that instead of imaging the source distribution at depth, eq. (13) images the source distribution on the surface. In practice, ϕ(g) is computed by taking all of the shot gathers recorded on the surface, filtering or muting all but the back-scattered surface wave arrivals, and inexpensively stacking the traces from all the sources recorded at g to get ϕ(g). This concept is illustrated in Fig. 2(a) for post-stack forward modelling and post-stack migration in Fig. 2(c). Here, large-amplitude anomalies with high resolution in the scattered migration image at Bs indicate the nearby presence of scatterers. Figure 2. View largeDownload slide Diagrams depict rays associated with (a) post-stack modelling, (b) pre-stack modelling, (c) natural post-stack migration and (d) natural pre-stack migration for a subwavelength scatterer at x in the near-field region of s. The post-stack trace ϕ(g) recorded at g is generated by simultaneously firing all of the shots (red stars); ϕ(g) is then migrated [depicted by dashed arrows in (c)] to the trial image points at s. In contrast, pre-stack traces recorded at g in (b) are generated by firing a point source at s, and then migrating the scattered surface waves in the traces to the trial image points at the source locations s on the source line. This process must be repeated for all the source positions, and then the pre-stack migration images are stacked together. The wiggly rays correspond to the evanescent surface waves that are in the near-field region of the sources, the solid straight rays correspond to surface waves in the far-field region of the geophones, and the dashed lines in (c) and (d) denote the ray paths associated with backpropagated wavefields. Figure 2. View largeDownload slide Diagrams depict rays associated with (a) post-stack modelling, (b) pre-stack modelling, (c) natural post-stack migration and (d) natural pre-stack migration for a subwavelength scatterer at x in the near-field region of s. The post-stack trace ϕ(g) recorded at g is generated by simultaneously firing all of the shots (red stars); ϕ(g) is then migrated [depicted by dashed arrows in (c)] to the trial image points at s. In contrast, pre-stack traces recorded at g in (b) are generated by firing a point source at s, and then migrating the scattered surface waves in the traces to the trial image points at the source locations s on the source line. This process must be repeated for all the source positions, and then the pre-stack migration images are stacked together. The wiggly rays correspond to the evanescent surface waves that are in the near-field region of the sources, the solid straight rays correspond to surface waves in the far-field region of the geophones, and the dashed lines in (c) and (d) denote the ray paths associated with backpropagated wavefields. If the source distribution is restricted to be a point source such that f (so) = δ(s − so), where the point source is in the near-field of the scatterer at xo in Fig. 1, then eq. (13) yields the surface-wave TRM image m(s, s΄) that is similar to the one for body waves (Schuster et al. 2012):   \begin{eqnarray} &&{{\bf s},{\bf s}^{\prime } \in B_s;\,\overbrace{ m({\bf s},{\bf s}^{\prime })}^{f({\bf s}^{\prime })^{\text{est}}}= \int _{\Omega _o} \int _{B_g} G({\bf s}^{\prime }|{\bf g})^*G({\bf g}|{\bf s}) \text{d}{\bf g}^2 \text{d}\omega ,}\nonumber \\ &=&\int _{\Omega _o} \int _{B_g}\big[\, \overbrace{G({\bf s}^{\prime }|{\bf g})_{\text{dir}}^*G({\bf g}|{\bf s})_{\text{dir}}}^{m({\bf s},{\bf s}^{\prime })_{\text{dir}-\text{dir}}} +\overbrace{G({\bf s}^{\prime }|{\bf g})_{\text{scatt}}^* G({\bf g}|{\bf s})_{\text{scatt}}}^{m({\bf s},{\bf s}^{\prime })_{\text{scatt}-\text{scatt}}} \nonumber \\ \nonumber \\ &&+\overbrace{G({\bf s}^{\prime }|{\bf g})_{\text{scatt}}^*G({\bf g}|{\bf s})_{\text{dir}}+ G({\bf s}^{\prime }|{\bf g})_{\text{dir}}^*G({\bf g}|{\bf s})_{\text{scatt}}}^{m({\bf s},{\bf s}^{\prime })_{\text{dir}-\text{scatt}}}\,\big]{}\text{d}{\bf g}^2 \text{d}\omega , \end{eqnarray} (14)where m(s, s΄) defines, for the actual point source at s, the migration image at the trial image location s΄ ∈ Bs. The next sections examine the contributions of the integrands in eq. (14) to superresolution imaging. Direct–direct images. For the Fig. 1 model, the direct-direct imaging term in eq. (14) can be analysed by inserting the far-field direct-wave Green’s function in eq. (5) into the m(s, s΄)dir–dir term in eq. (14) to give   \begin{eqnarray} \lim _{^{\scriptscriptstyle k|{\bf s}-{\bf g}| \gg 1}_{ \scriptscriptstyle k|{\bf s}^{\prime }-{\bf g}| \gg 1}} m({\bf s},{\bf s}^{\prime })_{\text{dir-dir}}\sim \frac{\pi }{0.5} \int _{\Omega _0} \left[ \int _{B_g} \frac{ e^{-ik(|{\bf s}-{\bf g}|-|{\bf s}^{\prime }-{\bf g}|) }}{k\sqrt{|{\bf s}-{\bf g}|}\sqrt{|{\bf s}^{\prime }-{\bf g}|}} {}\text{d}{\bf g}^2 \right]\text{d}\omega . \nonumber \\ \end{eqnarray} (15)Here, the far-field geophones are assumed to be many wavelengths from the sources so that k|s΄ − g|, k|s − g| ≫ 1. The far-field trial image points s΄ ∈ Bs produce kernels in eq. (15) that are too smooth to create sharp superresolution features in the migration image. Scattered–scattered images. For the single-scatterer model in Fig. 1, Appendix A shows that, in the near-field approximation k|s − xo|, k|s΄ − xo| ≪ 1, the scattered-scattered kernel m(s, s΄)scatt–scatt in eq. (14) becomes   \begin{eqnarray} \lim _{^{\scriptscriptstyle k|{\bf x}_o-{\bf g}| \gg 1}_{ \scriptscriptstyle k|{\bf s}^{\prime }-{\bf x}_o| \ll 1}}m({\bf s},{\bf s}^{\prime })_{\scriptscriptstyle \text{scatt-scatt}} &\sim& \alpha + \eta [\ln (|{\bf s}^{\prime }-{\bf x}_o|)+\ln (|{\bf s}-{\bf x}_o|)]\nonumber \\ &&+\,\gamma \ln (|{\bf s}-{\bf x}_o|)\ln (|{\bf s}^{\prime }-{\bf x}_o|)\,, \end{eqnarray} (16)where α, η and γ are constants. Unlike the smooth kernel $$\frac{ 1 }{\sqrt{|{\bf s}-{\bf g}|}\sqrt{|{\bf s}^{\prime }-{\bf g}|}}$$ in eq. (15), the logarithmic singularity ln (|s΄ − xo|) in eq. (16) provides the superresolution characteristics in the migration image for k|s΄ − xo| ≪ 1. Even though the original goal is to image the source distribution f (s) in eq. (13), subwavelength scatterers at xo in the near-field region of the source line Bs are indirectly detected by a strong sharp peak in the migration image m(s, s΄) for s΄ near xo. This is illustrated by the black curve in Fig. 3 which plots $$\text{Im}[H_0^{(2)}(kr)]$$ against −0.6 < kx < 0.6 for different values of the minimum distance ε > 0 between the source line and scatterer at xo. The sharply peaked black curve is what gives rise to superresolution imaging for ε < .001λ, where the black double-sided arrow indicates horizontal resolution better than 1/12λ. On the other hand, if k|x − xo| > 1 then the migration image resolution will be wavelength limited as illustrated by the blue curve in Fig. 3. Figure 3. View largeDownload slide Plot of the imaginary part of the Hankel function $$H_0^{(2)}(kr)$$ against kx, where $$r=\sqrt{(x-x_o)^2+\epsilon ^2}$$, x is the horizontal distance of the trial image point s΄ in Fig. 1 and ε is the minimum distance between the source line and the point scatterer at xo. For the black curve, the lateral resolution limit is indicated by the double-sided black arrow. Figure 3. View largeDownload slide Plot of the imaginary part of the Hankel function $$H_0^{(2)}(kr)$$ against kx, where $$r=\sqrt{(x-x_o)^2+\epsilon ^2}$$, x is the horizontal distance of the trial image point s΄ in Fig. 1 and ε is the minimum distance between the source line and the point scatterer at xo. For the black curve, the lateral resolution limit is indicated by the double-sided black arrow. Direct–scattered images. The direct-scattered kernels in eq. (14) consist of a smooth kernel G(g|s΄)dir (for |s΄ − g| ≫ λ) and a weakly singular kernel G(xo|s΄)scatt (for |s΄ − xo| ≪ λ) for geophones in the far-field region of s and s΄. The weakly singular kernel will provide superresolution peaks in the migration image for the trial image point s΄ near xo. Halo filtering. The goal with TRM is to image near-surface faults with superresolution, and not be distracted by the strong peaky images at each of the source locations. Therefore, traces at g within the near-field region of s΄ and s, i.e. |g − s|, |g − s΄| < 0.5λ, should be muted in order to avoid the direct Rayleigh wave singularities ln |s − g| or ln |s΄ − g| in m(s|s΄)dir–scatt and m(s|s΄)dir–dir. Otherwise there will be sharp strong features in the migration image that will masquerade as false signatures of scatterers. This muting of near-field traces will be accomplished by the halo filter h(s, s΄, g),   $$h({\bf s},{\bf s}^{\prime },{\bf g})=\left\lbrace \begin{array}{cc}1&\quad {\rm if} \,| {\bf s}-{\bf g}|, |{\bf s}^{\prime }-{\bf g}|>0.5 \lambda \\ 0&\quad {\rm if} \,| {\bf s}-{\bf g}|,|{\bf s}^{\prime }-{\bf g}|\le 0.5 \lambda \end{array} ,\right.$$ (17)which can be inserted into the integrand of eq. (13). 2.3 Natural pre-stack migration Eq. (13) is the surface wave TRM equation for imaging the source distribution f (s). It can be considered as a natural post-stack migration operation because ϕ(g) can be interpreted as post-stack data3 associated with the exploding source distribution f (s), and the integration weighted by the extrapolation kernel by G(g|s΄)* is only over the geophone positions on the surface (see Figs 2a and 2c). In contrast, the ray diagram associated with natural pre-stack migration used by AlTheyab et al. (2016) is illustrated in Fig. 2(d) and is mathematically described by   \begin{eqnarray} m ({\bf x}) &=& \int _{{\bf x}_s \in B} \int _{{\bf x}_r\in B}\int 2\omega ^2 W(\omega )^{*}G({\bf x}|{\bf x}_{s})^{*}G({\bf x}|{\bf x}_{r})^{*}\nonumber\\ &&\times \, \overbrace{u({\bf x}_s,{\bf x}_r)}^{\text{pre-stack data}}\text{d}\omega \text{d}x_r \text{d}y_r \text{d}x_s \text{d}y_s, \end{eqnarray} (18)where m(x) is the image at point x, ω is the angular frequency, W(ω) represents the source-wavelet spectrum and * denotes complex conjugation. x, xs and xr are, respectively, the trial-image, source and receiver positions at the surface locations in the set B; G(x|xs) is the Green’s function for the vertical-component source at xs and receiver at x, and G(x|xr) is the Green’s function for a vertical-component particle-velocity recording that only contains the transmitted wavefield without back-scattering. Here, the pre-stack modelling operation is illustrated by the ray diagram in Fig. 2(b), where a single source excites scattered waves that are recorded by the receivers. In this case, the natural pre-stack migration operation for a single frequency is a 4-D integral, two of the integrals are over the source coordinates and two are over the receiver coordinates. Therefore, if there is a grid of N × N sources shooting into a grid of N × N receivers the computational cost will be O(N4) algebraic operations per image point for natural pre-stack migration (Fig. 2d). This means that the total cost will be O(N6) operations for all N2 image points on the surface. In comparison, natural post-stack migration (Fig. 2c) will only cost O(N2) operations for each image point per frequency. In this case the total cost of natural post-stack migration is O(N4) operations for all N2 image points. 3 NUMERICAL RESULTS The surface-wave TRM operations will be tested on both synthetic elastic data and field data collected in Saudi Arabia. The synthetic data are generated by computing solutions to the elastic wave equation, and the input data for TRM are the muted back-scattered surface waves recorded at each vertical-displacement geophone. These tests demonstrate both the effectiveness and limitations of using formula (13) for TRM migration of surface waves. They also validate that post-stack migration of scattered surface waves can lead to superresolution imaging of subwavelength scatterers in the near-field region of the source location. The workflow for implementing post-stack migration of surface waves is illustrated in Fig. 4, which is summarized in the following six steps: Isolate the back-scattered surface waves from other arrivals. In our examples, a muting window is designed according to the estimated phase velocity of the recorded surface waves to mute everything except the back-scattered surface waves. Halo filtering is applied to the recorded data to obtain the Green’s function G(g|s΄)*, where the near-field traces within one dominant wavelength of the source are muted. The halo filter in eq. (17) is applied to the back-scattered surface waves, and stacking over the source index forms the post-stack back-scattered data ϕ(g) for receiver g. Migrate the post-stack back-scattered data ϕ(g) to compute the migration image f (s΄)0 by using eq. (13). A low-pass filter is applied to the post-stack back-scattered data, and then the migration image f (s΄)1 is computed with eq. (13). The enhanced image is obtained by multiplication of the images f (s΄)0 and f (s΄)1. Figure 4. View largeDownload slide Workflow for implementing natural post-stack migration of surface waves. Figure 4. View largeDownload slide Workflow for implementing natural post-stack migration of surface waves. 3.1 2-D buried scatterer A finite-difference method is used to generate 2-D elastic data for the scatterer model in Fig. 5. The model is discretized into 81 × 241 grid-points with the grid size of 0.5 m, where only one subwavelength scatterer is buried at the depth of 3.5 m, which is less than 1/3 of the dominant wavelength. That is, it is within the near-field region of the source and receiver locations directly above the scatterer. There are 120 receivers (z component) located on the surface with a 1 m geophone interval, and there are 120 sources shooting, with one source located at each geophone. The source wavelet is a Ricker wavelet with the peak frequency of 25 Hz, and the sources consist of point displacements in the z-direction. Figure 5. View largeDownload slide The shear velocity model for the 2-D simulations. The background shear velocity is 800 m s−1, the dashed line is the imaginary fault location, and only one subwavelength scatterer (indicated by the grey circle) with shear velocity of 400 m s−1 is embedded around the depth of 3 m. Figure 5. View largeDownload slide The shear velocity model for the 2-D simulations. The background shear velocity is 800 m s−1, the dashed line is the imaginary fault location, and only one subwavelength scatterer (indicated by the grey circle) with shear velocity of 400 m s−1 is embedded around the depth of 3 m. The CSGs (with source location (x, z) = (0,0) m) for the total, direct-transmitted, scattered and muted back-scattered surface waves are shown in Fig. 6. The total (Fig. 6a) and direct-transmitted (Fig. 6b) waves are generated for the scatterer and scatterer-free models in Fig. 5, respectively; the scattered waves (Fig. 6c) are obtained by the subtraction of the total and direct-transmitted waves; and the back-scattered surface waves (Fig. 6d) are muted from the total waves by using an appropriate muting window. These shot gathers are input into migration eq. (13) to give the single-source migration images in Figs 7(a1)–(a4) for different locations of the actual source location s along the free surface. As expected, the scatterer near the source location of sx = 37 m should be imaged in the total (Fig. 7a1), scattered-wave (Fig. 7a3) and muted back-scattered (Fig. 7a4) wave images, and undetectable in the direct–direct image (Fig. 7a2). Note the skinny yellowish band along the diagonals in Figs 7(a1) and (a4). This band of strong migration amplitudes is generated by migration of the direct Rayleigh wave from the source to a near-field geophone. Figure 6. View largeDownload slide Common shot gathers (with source location (x, z) = (0,0) m) for the (a) total, (b) direct-transmitted, (c) scattered and (d) muted back-scattered waves. Traces are for the vertical component of particle velocity. Figure 6. View largeDownload slide Common shot gathers (with source location (x, z) = (0,0) m) for the (a) total, (b) direct-transmitted, (c) scattered and (d) muted back-scattered waves. Traces are for the vertical component of particle velocity. Figure 7. View largeDownload slide Magnified single-source migration images without low-pass filtering for the (a1) total, (a2) direct, (a3) scattered and (a4) muted back-scattered waves. The second row contains the same images except the input data are low-pass filtered between 4 Hz and 18 Hz. The third row contains the enhanced migration images computed from eq. (19). Figure 7. View largeDownload slide Magnified single-source migration images without low-pass filtering for the (a1) total, (a2) direct, (a3) scattered and (a4) muted back-scattered waves. The second row contains the same images except the input data are low-pass filtered between 4 Hz and 18 Hz. The third row contains the enhanced migration images computed from eq. (19). Superresolution behaviour due to scatterers in the near-field of the source is frequency independent (Schuster et al. 2012), which means the same resolution should be achieved when the data are filtered by a low-pass filter. Here, a low-pass filter (4–18 Hz) is applied to the shot gathers in Fig. 6, and then input into the TRM eq. (13) to give the single-source migration images mf (s, s΄) in Figs 7(b1)–(b4). The scatterers near the source location of sx = 37 m are still imaged in the total (Fig. 7b1) and scattered-wave (Fig. 7b3) images. In contrast, the scatterer is slightly smeared in the muted back-scattered image (Fig. 7b4) and undetectable in the direct–direct image (Fig. 7b2). The enhanced TRM image mE(s, s΄) can be obtained by taking the product of the raw TRM image m0(s, s΄) with the low-pass filtered TRM image mf (s, s΄):   $$m_{E}({\bf s},{\bf s}^{\prime }) = m_{0}({\bf s},{\bf s}^{\prime }) \times m_{f}({\bf s},{\bf s}^{\prime }).$$ (19) These enhanced TRM images are shown in Figs 7(c1)–(c4) and more clearly reveal the location of the scatterer without the distraction of artefacts along the diagonal. The migration images m(s, s΄) at the shot location sx = 37 m for different migration kernels are shown in Fig. 8, where the trial shot location is along the horizontal axis and the normalized amplitude is along the vertical axis. Here, the black double arrow represents the dominant shear wavelength, and the blue, red and black curves represent, respectively, the TRM profiles without low-pass filtering, with low-pass filtering data, and the corresponding enhanced image profile. It is obvious that the needle-like profiles of the enhanced images mE(s, s΄)scatt–scatt (black curves in Figs 8c and d) achieve superresolution imaging better than λ/10. However, the resolution of the image from the back-scattered waves (Fig. 8d) is slightly reduced compared with that of the image from the scattered waves (Fig. 8c). In practice, we will mute all events except the back-scattered surface waves to compute the TRM image and use halo filtering to avoid geophone-related singularities in the migration kernels. An alternative to muting is FK filtering. Figure 8. View largeDownload slide The TRM profiles taken from the single-source migration images with the shot offset of 37 m for (a) total, (b) direct, (c) scattered and (d) back-scattered waves. The blue, red and black curves represent, respectively, the TRM profiles without low-pass filtering, with low-pass filtering, and the enhanced image in each subplot. The black double arrow indicates the dominant shear wavelength. Figure 8. View largeDownload slide The TRM profiles taken from the single-source migration images with the shot offset of 37 m for (a) total, (b) direct, (c) scattered and (d) back-scattered waves. The blue, red and black curves represent, respectively, the TRM profiles without low-pass filtering, with low-pass filtering, and the enhanced image in each subplot. The black double arrow indicates the dominant shear wavelength. 3.2 3-D fault model The TRM migration is now tested on 3-D synthetic data associated with the 3-D fault model, where Fig. 9(a) shows the depth slice at z = 2 m. There is an imaginary fault along the diagonal position, and some discretized subwavelength scatterers (the black dots) are embedded along the fault, the shear velocity of the scatterers is 400 m s−1. The model is discretized into 81 × 241 × 241 grid-points in the z, y and x directions, with the grid size of 0.5 m, and the background shear velocity is 600 m s−1. Figure 9. View largeDownload slide (a) Depth slice of the 3-D fault model (z = 2 m), where the background shear velocity is 600 m s−1, there are subwavelength scatterers (the black dots) located along the fault position, and the shear velocity of the scatterers is 400 m s−1 and b) the corresponding enhanced post-stack TRM migration image. The dominant background shear velocity wavelength is 24 m, which is indicated by the double-sided red arrow. Figure 9. View largeDownload slide (a) Depth slice of the 3-D fault model (z = 2 m), where the background shear velocity is 600 m s−1, there are subwavelength scatterers (the black dots) located along the fault position, and the shear velocity of the scatterers is 400 m s−1 and b) the corresponding enhanced post-stack TRM migration image. The dominant background shear velocity wavelength is 24 m, which is indicated by the double-sided red arrow. Figure 10. View largeDownload slide (a) Google map showing the location of the Qademah-fault seismic experiment. (b) Geophone layout for the Qadema-fault seismic experiment, where the green triangles represent the geophones. There are 12 parallel lines, each line has 24 geophones, the inline geophone interval is 5 m and the crossline offset is 10 m. Shots are located at each geophone, and a total of 288 shot gathers are acquired in total. The red dashed line represents the interpreted fault location. Figure 10. View largeDownload slide (a) Google map showing the location of the Qademah-fault seismic experiment. (b) Geophone layout for the Qadema-fault seismic experiment, where the green triangles represent the geophones. There are 12 parallel lines, each line has 24 geophones, the inline geophone interval is 5 m and the crossline offset is 10 m. Shots are located at each geophone, and a total of 288 shot gathers are acquired in total. The red dashed line represents the interpreted fault location. There are 3600 receivers (z-component) evenly distributed on the surface, with a 2 m spacing interval in both the x- and y-directions. The data are generated by computing solutions to the 3-D elastic wave equation with a Ricker source wavelet having the peak frequency of 25 Hz. Point sources are placed at each geophone location with a z-component displacement to generate 3600 shot gathers, each with 3600 traces. An appropriate muting window is designed, based on the phase velocity of the surface wave for each shot, to extract the back-scattered surface waves. In addition, a Halo filter is used to mute the near-field traces within one dominant wavelength of the source. The muted back-scattered surface waves are used to compute the post-stack TRM migration image in Fig. 9(b). Here, the workflow in Fig. 4 is used to compute the post-stack image and the strong migration amplitudes clearly reveal the locations of the scatterers at subwavelength resolution. 3.3 3-D Qademah fault data A 3-D land survey was carried out along the Red Sea coast over the Qademah fault system, about 30 km north of the KAUST campus (Hanafy 2015). There are 288 vertical-component geophones arranged in 12 parallel lines, and each line has 24 geophones. The inline geophone interval is 5 m and the crossline interval is 10 m. The geophone geometry layout is shown in Fig. 10(a), where one shot is fired at each geophone location for a total of 288 shot gathers. The source is generated by a 200 lb hammer striking a metal plate on the ground. The red dashed line represents the interpreted fault location from geological evidence (Roobol & Kadi 2008). The post-stack TRM operation can be used to image near-field subwavelength scatterers,where the scatterers at shallow depth are sensitive to higher frequency data, while deeper scatterers are mostly sensitive to lower frequency data. First, we apply different narrowband filters with different peak frequencies (from 13 to 41 Hz, each with a 4 Hz increase in the frequency band) to the raw Qademah data. A CSG (no. 120) filtered by a narrowband filter with the peak frequency of 25 Hz is shown in Fig. 11(a), where the red dashed curves are computed from the estimated phase velocity, and they are used to mute the back-scattered surface waves from the total wavefield. The muted back-scattered surface waves shown in Fig. 11(b) are the input into the TRM operation in eq. (13) for computation of the post-stack TRM image. For comparison, the data in Fig. 11(b) are also used to calculate the natural pre-stack migration images (Liu et al. 2017). Figure 11. View largeDownload slide Common shot gather no. 120 from the Qademah fault data filtered by a bandpass filter with 25 Hz peak frequency, (a) the total wavefield, where the red dashed curve is used to design a muting window for muting the back-scattered waves, and (b) the muted back-scattered waves. Figure 11. View largeDownload slide Common shot gather no. 120 from the Qademah fault data filtered by a bandpass filter with 25 Hz peak frequency, (a) the total wavefield, where the red dashed curve is used to design a muting window for muting the back-scattered waves, and (b) the muted back-scattered waves. According to Hanafy (2015)’s interpretation, the fault ‘F1’ (at an offset ≈90 m in the inline direction) in Fig. 12(c) is the north-south Qademah fault, which is believed to be a normal rotational fault resulting from the opening of the Red Sea, the fault ‘F2’ in the westside of ‘F1’ is interpreted as an antithetic fault. The area between faults ‘F1’ and ‘F2’ are composed of loose sediments and colluvial rocks. The surface waves scattered by these colluvial rocks with size less than the dominant shear wavelength are to be imaged by the TRM operation. Figure 12. View largeDownload slide (a) The natural post-stack migration (TRM) image, (b) natural pre-stack migration (Liu et al. 2017) (NPM) image for the Qademah data with different dominant frequencies and (c) Rayleigh phase-velocity tomogram (Hanafy 2015) at different frequencies. Here, the red double-sided arrows indicate the approximate dominant shear wavelength at different frequencies. The dashed pink curves in (a) labelled with ‘S1’, ‘S2’ and ‘S3’ indicate three segments of subwavelength scatterers between faults ‘F1’ and ‘F2’. The dashed pink curves in (b) indicate the location of Qademah fault ‘F1’. These two white dashed curves in (c) labelled with ‘F1’ and ‘F2’ represent the Qademah fault and another small antithetic fault, respectively. Figure 12. View largeDownload slide (a) The natural post-stack migration (TRM) image, (b) natural pre-stack migration (Liu et al. 2017) (NPM) image for the Qademah data with different dominant frequencies and (c) Rayleigh phase-velocity tomogram (Hanafy 2015) at different frequencies. Here, the red double-sided arrows indicate the approximate dominant shear wavelength at different frequencies. The dashed pink curves in (a) labelled with ‘S1’, ‘S2’ and ‘S3’ indicate three segments of subwavelength scatterers between faults ‘F1’ and ‘F2’. The dashed pink curves in (b) indicate the location of Qademah fault ‘F1’. These two white dashed curves in (c) labelled with ‘F1’ and ‘F2’ represent the Qademah fault and another small antithetic fault, respectively. Fig. 12(a) only shows the post-stack TRM images (from the top-to-bottom panels) calculated from the filtered back-scattered data with dominant frequencies of 29, 25 and 21 Hz, respectively. The dashed pink curve (on the top panel) labelled with ‘S1’ indicates the location of the first segment of the fault, and the estimated depth of these subwavelength scatterers is around 5 m, which is roughly estimated as one third of the corresponding shear wavelength (Stokoe & Nazarian 1985). The post-stack TRM image in the middle panel of Fig. 12(a) is calculated from 25 Hz back-scattered data, where the dashed pink curve labelled with ‘S2’ indicates the location of the second segment of the fault. We found that there is an overlap between the ‘S1’ and ‘S2’ images, which indicates that the scatterers which belong to ‘S1’ are also imaged by 25 Hz data. The strong amplitudes in ‘S2’ reveal that the approximate depth of these subwavelength scatterers is around 6.7 m. The bottom panel in Fig. 12(a) shows the post-stack TRM image calculated from 21 Hz back-scattered data, the dashed pink curve labelled as ‘S3’ indicates the location of the third segment of the fault, and the strong amplitudes suggest that the approximate depth of these subwavelength scatterers is around 10 m. The red double-sided arrows represent the dominant shear wavelength in each panel, where the horizontal resolution of the fault images is less than half the dominant shear wavelength. For comparison, the natural pre-stack migration (NPM) images calculated by eq. (18) are computed by Liu et al. (2017) and shown in Fig. 12(b). The NPM images are computed from the filtered back-scattered data with dominant frequencies of 29, 25 and 21 Hz, respectively, for the top-to-bottom panels. The dashed pink curve in each panel indicates the location of the Qademah Fault (‘F1’). Figs 13(a)–(c) show the TRM and NPM image profiles along ‘A1 − B1’, ‘A2 − B2’ and ‘A3 − B3’, respectively. Here, the red-dashed and blue-solid curves represent the TRM and NPM profiles, respectively. And the black double-sided arrows indicate the dominant shear wavelength in each subplot. We can see the TRM image profile has achieved a resolution better than one half of the dominant shear wavelength, while the resolution of NPM is around one dominant shear wavelength. This difference in resolution limits for pre-stack and post-stack migration are roughly in agreement with the resolution limits predicted by (Chen & Schuster 1999). Figure 13. View largeDownload slide The profile comparisons between the post-stack TRM (Fig. 12a) and NPM (Fig. 12b) images along (a) ‘A1 − B1’, (b) ‘A2 − B2’ and (c) ‘A3 − B3’. The red dashed and blue solid curves represent image results from TRM and NPM, respectively. The black double-side arrows indicate the dominant shear wavelength in each subplot, and amplitude has been normalized. Figure 13. View largeDownload slide The profile comparisons between the post-stack TRM (Fig. 12a) and NPM (Fig. 12b) images along (a) ‘A1 − B1’, (b) ‘A2 − B2’ and (c) ‘A3 − B3’. The red dashed and blue solid curves represent image results from TRM and NPM, respectively. The black double-side arrows indicate the dominant shear wavelength in each subplot, and amplitude has been normalized. 4 SUMMARY AND CONCLUSIONS We theoretically and numerically analysed the efficacy of surface-wave TRM migration. Theoretical formulae and synthetic data results show that surface waves can be used to achieve superresolution imaging of subwavelength scatterers if they are located less than about 1/3 of a shear wavelength from the source line. We also show that the TRM operation can be used for natural post-stack migration, and only costs O(N4) algebraic operation per frequency at the N2 trial image points on the surface. This compares to O(N6) operations for natural pre-stack migration at a single frequency. Here, we assume the sources and receivers are on a N × N grid on the free surface and there are N2 trial image points on the surface. These theoretical formulae show, for the first time, the connections between pre-stack natural migration, post-stack natural migration and the TRM operations. It also shows the potential for achieving superresolution imaging of scattered surface waves in the near-field region of the source. The field-data tests suggest that hidden faults can be detected with the post-stack TRM operation if they are no deeper than about 1/2 the dominant wavelength. These tests also show images with resolution that is about 1/3 that of the dominant wavelength. The most significant challenge in natural post-stack migration is to process the data so that only the back-scattered surface waves are used. This is a subject of ongoing research. ACKNOWLEDGEMENTS The research reported in this publication was supported by the King Abdullah University of Science and Technology (KAUST) in Thuwal, Saudi Arabia and the CRG grant. We are grateful to the sponsors of the Center for Subsurface Imaging and Modeling (CSIM) Consortium for their financial support. For computer time, this research used the resources of the Supercomputing Laboratory at KAUST and the IT Research Computing Group. We thank them for providing the computational resources required for carrying out this work. We greatly appreciate the constructive comments and suggestions from editor Herve Chauris, and two anonymous reviewers, which helped improve this paper. Footnotes 1 See eq. 7.2.18 in Feshbach & Morse (1953), where $$H_0^{(1)}(kX)=H_0^{(2)}(kX)^*$$. 2 That is, the horizontal variations are small on the scale of the largest wavelength. 3 The post-stack traces ϕ(g) result from exploding the many different sources at s ∈ Bs with strength f (s). REFERENCES AlTheyab A., Lin F.-C., Schuster G.T., 2016. Imaging near-surface heterogeneities by natural migration of backscattered surface waves, Geophys. J. Int. , 204( 2), 1332– 1341. Google Scholar CrossRef Search ADS   Chen J., Schuster G.T., 1999. Resolution limits of migrated images, Geophysics , 64( 4), 1046– 1053. Google Scholar CrossRef Search ADS   De Fornel F., 2001. Evanescent Waves: from Newtonian Optics to Atomic Optics , Vol. 73, Springer Science & Business Media. Dutta G., AlTheyab A., Tarhini A., Hanafy S., Schuster G.T., 2016. Extracting 220 Hz information from 55 Hz field data by near-field superresolution imaging, Geophys. J. Int. , 206, 1194– 1203. Google Scholar CrossRef Search ADS   Feshbach H., Morse P.M., 1953. Methods of Theoretical Physics , McGraw-Hill Interamericana. Fink M., 1993. Time-reversal mirrors, J. Phys. D: Appl. Phys. , 26( 9), 1333– 1350. Google Scholar CrossRef Search ADS   Fink M., 2006. Time-reversal acoustics in complex environments, Geophysics , 71( 4), SI151– SI164. Google Scholar CrossRef Search ADS   Fink M., 2008. Time-reversal waves and super resolution, J. Phys.: Conf. Ser. , 124, 012004, pp. 1–29. Google Scholar CrossRef Search ADS   Hanafy S.M., 2015. Mapping the Qademah fault with traveltime, surface-wave, and resistivity tomograms, SEG Technical Program Expanded Abstracts , pp. 3347– 3351. Jia H., Ke M., Hao R., Ye Y., Liu F., Liu Z., 2010. Subwavelength imaging by a simple planar acoustic superlens, Appl. Phys. Lett. , 97( 17), 173507, doi:10.1063/1.3507893. Google Scholar CrossRef Search ADS   Lerosey G., de Rosny J., Tourin A., Fink M., 2007. Focusing beyond the diffraction limit with far-field time reversal, Science , 315( 5815), 1120– 1122. Google Scholar CrossRef Search ADS PubMed  Liu Z., Altheyab A., Hanafy S., Schuster G., 2017. Imaging near-surface heterogeneities by natural migration of surface waves: field data test, Geophysics , 82( 3), 1– 9. Google Scholar CrossRef Search ADS   Roobol M., Kadi K., 2008. Cenozoic faulting in the Rabigh area, central west Saudi Arabia (including the sites of King Abdullah Economic City and King Abdullah University for science and technology), Saudi Geological Survey Technical Report SGS-TR-2008-6, 1(250,000). Schuster G.T., 2002. Reverse-time migration= generalized diffraction stack migration, in SEG Technical Program Expanded Abstracts 2002 , pp. 1280– 1283, Society of Exploration Geophysicists. Schuster G.T., Hanafy S., Huang Y., 2012. Theory and feasibility tests for a seismic scanning tunnelling macroscope, Geophys. J. Int. , 190( 3), 1593– 1606. Google Scholar CrossRef Search ADS   Stokoe K.H., Nazarian S., 1985. Use of Rayleigh waves in liquefaction studies, in Measurement and Use of Shear Wave Velocity for Evaluating Dynamic Soil Properties , p. 17, ASCE. Tanimoto T., 1990. Modelling curved surface wave paths: membrane surface wave synthetics, Geophys. J. Int. , 102( 1), 89– 100. Google Scholar CrossRef Search ADS   Tarhini A., Guo B., Dutta G., Schuster G.T., 2017. Extension of seismic scanning tunnelling macroscope to elastic waves, Pure Appl. Geophys. , doi:10.1007/s00024-017-1692-x. Yilmaz Ö., 2001. Seismic Data Analysis , Society of Exploration Geophysicists. Google Scholar CrossRef Search ADS   Yomogida K., Aki K., 1985. Waveform synthesis of surface waves in a laterally heterogeneous earth by the Gaussian beam method, J. geophys. Res. , 90( B9), 7665– 7688. Google Scholar CrossRef Search ADS   Zabolotskaya E.A., Ilinskii Y.A., Hay T.A., Hamilton M.F., 2012. Green’s functions for a volume source in an elastic half-space, J. acoust. Soc. Am. , 131( 3), 1831– 1842. Google Scholar CrossRef Search ADS PubMed  APPENDIX A: SCATTERED-SCATTERED IMAGE For the Fig. 1 model, assume that the scatterer at xo is in the near-field region of both s and s΄, and it is also in the far-field region of the geophone at g. Using eqs (5) and (6) and the definitions in eq. (9) gives the asymptotic expression for the scattered Green’s functions in eq. (14):   \begin{eqnarray} \lim _{^{\scriptscriptstyle k|{\bf x}_o-{\bf g}| \gg 1}_{ \scriptscriptstyle k|{\bf s}-{\bf x}_o| \ll 1}}G({\bf g}|{\bf s})_{\text{scatt}}&=&-\pi ^2 R\omega ^2 \,\,\lim _{^{\scriptscriptstyle k|{\bf x}_o-{\bf g}| \gg 1}_{ \scriptscriptstyle k|{\bf s}-{\bf x}_o| \ll 1}}H_o^{(2)}(k|{\bf s}-{\bf x}_o|)H_o^{(2)}(k|{\bf x}_o-{\bf g}|), \nonumber\\ &\sim & -2R \omega ^2 \ln ( k|{\bf s}-{\bf x}_o|) \frac{e^{-i(k|{\bf g}-{\bf x}_o|+\pi /4)}}{\sqrt{0.5 k |{\bf x}_o-{\bf g}|/\pi }}, \end{eqnarray} (A1)  \begin{eqnarray} \lim _{^{\scriptscriptstyle k|{\bf x}_o-{\bf g}| \gg 1}_{ \scriptscriptstyle k|{\bf s}^{\prime }-{\bf x}_o| \ll 1}}G({\bf s}^{\prime }|{\bf g})_{\text{scatt}}^* &\sim & -2R \omega ^2\ln ( k|{\bf s}^{\prime }-{\bf x}_o|) \frac{e^{i(k|{\bf g}-{\bf x}_o|+\pi /4)}}{\sqrt{0.5 k |{\bf x}_o-{\bf g}|/\pi }}, \end{eqnarray} (A2)where R = δs is the scattering coefficient, k = ω/c, and c is the Rayleigh-wave phase velocity. Conveniently setting $$R=\frac{1}{\sqrt{8 \pi }}$$, the scattered-scattered image m(s, s΄)scatt–scatt in eq. (14) becomes   \begin{eqnarray} m({\bf s},{\bf s}^{\prime })_{\text{scatt--scatt}}&=& \int _{\Omega _0} \left[ \int _{B_g} G({\bf s}^{\prime }|{\bf g})_{\text{scatt}}^*G({\bf g}|{\bf s})_{\text{scatt}} \text{d}{\bf g}^2 \right]\text{d}\omega ,\nonumber \\ &=& \beta \int _{\Omega _k} k^3 { \ln (k |{\bf s}^{\prime }-{\bf x}_o|)\ln (k |{\bf s}-{\bf x}_o|) } \text{d}k, \end{eqnarray} (A3)where the integration over frequency has been transformed to an integration over the wavenumber region denoted by Ωk and $$\beta =\int _{B_g} \frac{c^5}{ |{\bf x}_o-{\bf g}|}\text{d}{\bf g}^2$$. Here, s, s΄ ∈ Bs and the integration over geophone positions is along the set of points denoted by Bg. Eq. (A3) can be expressed as   \begin{eqnarray*} m({\bf s},{\bf s}^{\prime })_{\scriptscriptstyle \text{scatt--scatt}}&=&\beta \int _k {\text{d}k}\,k^3\lbrace {{(\ln k)^2 +\ln k[\ln (|{\bf s}^{\prime }-{\bf x}_o|)+\ln (|{\bf s}-{\bf x}_o|)]+\ln (|{\bf s}^{\prime }-{\bf x}_o|)\ln ( |{\bf s}-{\bf x}_o|)}{}}\rbrace ,\nonumber \\ &=&\alpha + \eta [\ln (|{\bf s}^{\prime }-{\bf x}_o|)+\ln (|{\bf s}-{\bf x}_o|)]+\gamma \ln (|{\bf s}-{\bf x}_o|)\ln (|{\bf s}^{\prime }-{\bf x}_o|)\,, \end{eqnarray*} where α = β∫kk3(ln k)2dk, η = β∫kk3ln kdk, and γ = β∫k3dk are finite-valued numbers that depend on the range of the integration over wavenumbers. The weakly singular ln |s΄ − xo| term is responsible for near-field superresolution imaging of surface waves for trial image points s΄ near the scatterer at xo. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Feb 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 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. ### 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

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

14-day Free Trial