Passive body-wave interferometric imaging with directionally constrained migration

10.1093/gji/ggy306/5060754

Publication date

2018

Document Version

Final published version

Published in

Geophysical Journal International

Citation (APA)

Almagro Vidal, C., van der Neut, J., Verdel, A., Hartstra, I., & Wapenaar, K. (2018). Passive body-wave

interferometric imaging with directionally constrained migration. Geophysical Journal International, 215(2),

1022-1036. https://doi.org/10.1093/gji/ggy306/5060754

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

Geophys. J. Int. (2018)215, 1022–1036 doi: 10.1093/gji/ggy306 Advance Access publication 2018 July 27

GJI Seismology

Passive body-wave interferometric imaging with directionally

constrained migration

Carlos Almagro Vidal,

Joost van der Neut,

Arie Verdel,

Iris Eline Hartstra

and

Kees Wapenaar

1_{Department of Geoscience and Engineering, Faculty of Civil Engineering and Geosciences, Delft University of Technology, 2628 CN Delft, The Netherlands.} E-mail:carlos.almagrovidal@gmail.com

2_{TNO Sustainable Geo Energy, 3584 CB, Utrecht, The Netherlands}

Accepted 2018 July 25. Received 2018 July 18; in original form 2018 January 23

S U M M A R Y

Passive seismic interferometry enables the estimation of the reflection response of the subsur-face using passive receiver recordings at the sursubsur-face from sources located deep in the Earth. Interferometric imaging makes use of this retrieved reflection response in order to study the subsurface. Successful interferometric imaging relies on the availability of passive recordings from sufficient sources in the subsurface. Ideally, these sources should be homogeneously dis-tributed, which is unlikely to happen in practical applications. Incomplete source distributions result in the retrieval of inaccurate reflection responses, containing artefacts which can disturb the interferometric imaging process. We propose an alternative imaging method for passive data based on illumination diagnosis and directionally constrained migration. In this method, passive responses from single transient sources are cross-correlated individually, and the dom-inant radiation direction from each virtual source is estimated. The correlated responses are imaged individually, thereby limiting the source wavefield to the dominant radiation direction of the virtual source. This constraint enables the construction of accurate images from indi-vidual sources with a significantly reduced amount of migrated interferometric artefacts. We also show that the summation of all individual imaging results improves the subsurface image by constructive interference, while migrated crosstalk and artefacts experience cancellation. This process, called Image Interferometry, shows that in case of limited subsurface illumi-nation the interferometric integration can be applied in the image domain rather than in the virtual reflection-response domain, thus eliminating the need for the retrieval of the reflection response as an intermediate step.

Key words: Seismic Interferometry; Body waves; Crustal imaging.

1 I N T R O D U C T I O N

Seismic interferometry (SI) aims to reconstruct the impulsive response between receivers, as if one of the receivers were a source (Schuster

2001; Weaver & Lobkis2001; Campillo & Paul2003; Wapenaar2003; Snieder2004; Schuster2009; Galetti & Curtis2012). The retrieval of the impulsive response between the receiver pairs results from the superposition of the correlations of the receiver recordings over individual contributions from sources surrounding the receivers. SI can be applied to surface waves as well as to body waves. For body waves, SI has been successfully applied in passive seismics with naturally occurring passive sources (Roux et al.2005; Nishida2013). The recovery of the impulse response from the correlation responses stems from the constructive interference of events in stationary phase regions, and the cancellation of the remaining correlated events corresponding to non-physical arrivals (Snieder et al.2006). Examples of SI applications with body waves to passive seismics are traveltime tomography (Nakata et al.2015; Olivier et al.2015) and retrieval of reflection events (Rickett & Claerbout1999; Abe et al.2007; Draganov et al.2007; Poli et al.2012; Bou´e et al.2013; Lin et al.2013). The retrieval of reflections with deep passive sources enables the study and imaging of the subsurface by treating the correlation responses as virtual-source records and consecutively employ them for depth migration (Tonegawa et al.2009; Draganov et al.2010; Ruigrok et al.2010). To accomplish this, it is required to have passive sources illuminating the receivers uniformly over all angles. We refer to this procedure as conventional passive SI imaging.

For the case of complete illumination, an alternative to the conventional passive SI imaging entails migrating the correlated responses due to individual passive sources, followed by summation of the migration results, thus obviating the requirement to construct a complete 1022 CThe Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.

perform lithospheric imaging of the Moho.

2 C O R R E L AT I O N F U N C T I O N

For transient sources, Wapenaar & Fokkema (2006) introduce a relation in acoustic media to retrieve the Green’s function ˆG(x_A, x0, ω)

between a receiver at x_Aand a virtual source at x0 from recordings at these two positions of a continuous distribution of passive sources

(individually located at xB). In a passive seismic configuration with the receiver locations at the free surface, the retrieved Green’s function ˆ

G corresponds to the reflection response of the medium, ˆR3  :

{ ˆR3  (xA, x0, ω)}| ˆS(ω)|2≈  xB ρ(xB) cP(xB)vˆ3obs(xA, xB, ω)  ˆ vobs 3 (x0, xB, ω)∗dxB, (1)

where R stands for the real part,{}∗denotes complex conjugation, ω is the angular frequency and ˆ indicates the wavefield is in the frequency domain. In eq. (1), the observed wavefield ˆvobs

3 (x0, xB, ω) = ˆG3,3(x0, xB, ω) ˆS(ω), is the vertical particle-velocity Green’s function

ˆ

G_3,3(x0, xB, ω) due to a vertical point-force at xB, multiplied by the Fourier transform of the source function ˆS(ω). The product on the right

hand side of eq. (1) corresponds to a cross-correlation in the time domain. The integration over x_Bis defined by the distribution of passive point-force sources in the medium, andρ and cPstand for the mass density and acoustic velocity of the medium at the locations of these

passive sources. The result of the integration consists of the real part of the desired impulsive reflection response ˆR3  (xA, x0, ω) (this is the

representation for a vertical point-force source at x0and vertical particle-velocity wavefield at xA), multiplied by the power spectrum of the

source function| ˆS(ω)|2_{= ˆS(ω){ ˆS(ω)}}∗_.

The correct estimation of the reflection response with eq. (1) requires the correlation responses of records from uniformly distributed passive sources with the same spectrum ˆS(ω), which illuminate the receivers from all possible angles. In many cases, passive sources are sparsely distributed and clustered. In that case, we carry out the approximation by discretizing eq. (1):

{ ˆR3  (xA, x0, ω)}| ˆS(ω)|2∝



x_B

ˆ

Cx_B(xA, x0, ω) , (2)

where ˆCx_B(xA, x0, ω) stands for the correlation function of a single passive source at xB:

ˆ Cx_B(xA, x0, ω) = ˆvobs3 (xA, xB, ω)  ˆ vobs 3 (x0, xB, ω) _∗ . (3)

Eqs (2) and (3) ignore scaling factorsρ and cPat the source locations xBbecause they are not known in practice.

Fig.1shows an example of applying eq. (2) in an acoustic model (Fig.1a) with three different source distributions. Figs1(c), (d) and (e) show the interferometric results for the three different source distributions displayed on top of their corresponding results. The case in Fig.1(c) produces a complete retrieval of the reflection response due to the homogeneous source distribution in the subsurface. Within the cone limited by the first arrivals (indicated by the red dashed line), this result resembles the reference response in Fig.1(b). In this scenario, the estimation of the reflected events by constructive interference is optimal while the source distribution in the subsurface is sufficiently dense to achieve an acceptable cancellation of correlation artefacts from different stationary points. The case in Fig.1(d) presents the retrieved reflection result obtained with a limited amount of sources, clustered in one region of the subsurface. The retrieved result shows the constructive interference is restricted to the reflection events that can be obtained for that limited illumination range, but the destructive interference still manages to eliminate the correlation events not in stationary phase within the array. The case in Fig.1(e) displays the cross-correlation result obtained from a single source, where no constructive nor destructive interference can be achieved. The correlated events seem to match with the reflections from scenario1(b), yet they show incorrect arrival times since they are not in stationary phase with the source–receivers geometry. In this study, we assume the integration of individual passive recordings in eq. (2) is not attainable, due to a lack of passive sources. Hence, the focus lies on the migration of the data in the correlation function ˆCx_Bthat is in stationary phase, and in minimizing correlated artefacts

1024 C. Almagro Vidal et al.

Figure 1. Three different interferometric results using SI by cross-correlation. (a) Acoustic velocity model with receivers at the free surface (x_A, indicated by yellow triangles) and passive sources (x_B, indicated by red dots). (b) Reference reflection response R₃ due to an actual source at x0= 6000 m (red star). (c) Retrieved virtual-source reflection response at the same location (hollow red star) from a SI scenario with equipartitioned illumination from the subsurface. (d) Same result from a limited clustered source distribution. (e) Result obtained from a single source in the subsurface. The red dashed line delimits the direct/first arrival. The respective source distributions are displayed on top.

3 D I R E C T I O N A L C O N S T R A I N T S

The difficulty of working with incomplete distributions of passive sources in the subsurface is that part of the necessary information to retrieve a proper reflection response is missing. The migration of the correlation functions turns into an incomplete process since it fails to suppress the correlated artefacts in the image result. In this section we aim to obtain additional information that could serve to constrain the migration process.

The scattered field originated from free-surface reflections due to a passive source in the subsurface brings information about the reflectors in the medium. In Figs2(a)–(c) we illustrate the process of retrieving a reflection response between two receivers. The specular ray from the source (the direct arrival to the first receiver, to become virtual source at x0) defines the direction in which the correct retrieved

reflection ray can be found. For each passive-/virtual-source pair, in laterally invariant media, there exists a unique ray parameter that defines this specular ray.

Almagro Vidal et al. (2014) introduced a method to determine the dominant ray parameter of a correlation function at a specific virtual-source location. The original aim of the method was a qualitative analysis of ambient-noise recordings: to separate noise recordings which are dominated by surface waves from those suitable for the retrieval of body-wave reflections. This ray parameter analysis also provides a quantitative analysis of the illumination characteristics of the passive source. The correlation function Cx_B features a source function

around zero time lag, which quantifies the illumination characteristics of the passive source. We name this section the virtual-source function. Since direct arrivals are generally the most energetic events, they dominate the virtual-source function. The illumination diagnosis over the virtual-source function determines the specular-ray path of the direct wavefield from the passive source with respect to the virtual-source location (Fig.2b).

Figure 2. Illustration of the reflection-response retrieval by passive SI with one reflector, and its relation to directionally constrained migration. The receivers

are shown with yellow triangles and the passive source with a red star. (a) A receiver at xArecords a field originating from a subsurface source (xB) after being scattered by a reflector. A receiver at x0records the direct field from the source. The specular ray from the passive source passes along these receivers. (b) The cross-correlation of the response at xAwith the one at x0retrieves the reflection response at xAas if a source were located at receiver x0(red triangle). The locations of the passive source and the virtual source define a unique ray parameter (px0

xB). (c) The value of this ray parameter defines the “specular reflection”

direction from the free surface by this virtual source. In order to find the reflector location in stationary phase, only this ray parameter, and not the location of the passive source x_B, is needed.

In Almagro Vidal et al. (2014) the analysis of the ray parameter distribution of the virtual-source function is described with a linear-slant stack on the time-domain correlation function Cx_B(xA, x0, t) at t = τ + p · (xH,A− xH,0), forτ = 0:

 Cx_B(x0, p) =  x_A Cx_B  xA, x0, p · (xH,A− xH,0)  dxA, (4)

where p is the ray parameter vector, xH,A correspond to the horizontal coordinates of xA and Cx_B is the illumination distribution of the

virtual-source function for each virtual source x0. However, when the distance of the passive source to the acquisition array is of the same

order of magnitude as the array aperture, a linear slant-stack does not suffice and a parabolic approximation is required for better precision on the ray parameter analysis of the virtual-source radiation (van der Neut et al.2011). The dominant ray parameter that delimits the illumination direction of the wavefield at the virtual-source location is defined as:

px0

xB= arg max p

Cx_B(x0, p). (5)

A display of the illumination distribution of a virtual-source function Cx_B is shown in Figs3(g), (h) and (i) (with their respective dominant

ray parameter px0

x_B), corresponding to the parabolic slant-stack applied on the virtual-source functions in Figs3(d), (e) and (f), respectively.

All results correspond to the model scenario described in Figs3(a), (b) and (c), for the same virtual-source location.

4 M I G R AT I O N S C H E M E

Since from the cross-correlation we aim to obtain correct reflections for a specific ray parameter only, we require a directionally constrained migration scheme. The method we propose here is an adaptation from the work of Popov et al. (2010), where, similar to other migration methods, the imaging condition is defined by the correlation of a forward wavefield with the backprojection of the receiver wavefield; in our configuration, the fields emitted from the virtual-source and receiver locations, respectively. This method uses high-frequency asymptotics of Gaussian beams to reconstruct the Green’s functions of the medium. The summation of the beams of different directions approximates the wavefield in the medium. Every individual Gaussian beam is defined by its ray centred coordinates s(x) and n(x) of any location x of the medium in the proximity of the beam ( ˇCerven´y et al.1982). In Popov et al. (2010), the Green’s function between location x0and any point

x in the 3-D medium is represented as the integration of individual Gaussian beams ( ˆuG B_{) over different directions (described by azimuthal}

and polar anglesθ and φ). This is expressed as: ˆ G(x, x0, ω) =  _π 0  2π 0 ˆ (θ, φ, ω) ˆuG B (s, n, x0, θ, φ, ω) dθ dφ, (6)

where the ray centred coordinates s and n define the observation location x associated with the beams passing its proximity. ˆ defines the initial amplitudes of the Gaussian beams (Popov1982). The behaviour of the Gaussian beam can be controlled by their width and curvature. These parameters are defined at the receiver locations, following the construction of Nowack et al. (2006). An adequate estimation of the beam width can be found in Hill (1990).

1026 C. Almagro Vidal et al.

Figure 3. Results obtained from illumination diagnosis of the virtual source at x0= 2500 m from three different sources located at 3000 m depth, at xB= 1250, 3000 and 4750 m, respectively. (a), (b) and (c) Passive source location (red star) and receiver array (yellow triangles) displayed inside the acoustic velocity model. Red triangle represents the virtual source. (d), (e) and (f) Virtual-source functions from CxB(xA, x0, t) at x0= 2500 m from each individual source of (a), (b) and (c). Illumination diagnosis analyzes the energy distribution in ray parameters of the virtual-source radiation (green lines), in order to identify the dominant direction of the most energetic arrivals from the passive source (blue line) to the acquisition array. (g), (h) and (i) ray parameter distribution 

Cx_B(x0, p) with the dominant direction represented by pxx0_B(blue square).

For the passive-seismic case with isotropic illumination, the forward wavefield should radiate in all angles in the migration process. However, for the migration of the correlation function of a single source CxB(xA, x0, t), the forward wavefield is to be limited to the

dominant illumination direction only. Using the results from the illumination diagnosis previously described, we aim to constrain the illumination of the forward wavefield by imposing the radiation pattern of the virtual-source function. Making use of the medium velocity cP(x0) at the virtual-source location, we convert the coordinates of the virtual-source function from ray parameters into angular directions:



Cx_B(x0, θ, φ) = Cx_B(x0, p), using the horizontal-slowness coordinates of the ray parameter p(θ, φ) = (cos(θ) sin(φ)_c

P ,

sin(θ) sin(φ)

cP ). The approximated Green’s function due to a directionally constrained virtual-source located at x0evaluated at x is weighted according to the radiation pattern

Figure 4. Wavefields and imaging result at instant t= 0.78 s for the model in Fig.3(a). The passive source is the red star and the receivers are the yellow triangles. The red triangle represents the virtual-source located at x0= 2500 m. (a) Forward wavefield guided in a single direction from the virtual-source location (red line). (b) Backprojection of the receiver wavefield with correlation function. (c) Partial image result using the correlation imaging condition between Dx_Band Ux_Bat instant t= 0.78 s.

described by the normalized ray parameter distribution of the virtual-source function Cx_B(as the ones shown in Figs3g, h and i):

ˆ GG B_xB(x, x0, ω) =  _π 0  2π 0  Cx_B(x0, θ, φ) ˆ(θ, φ, ω) ˆuG B(s, n, x0, θ, φ, ω) dθ dφ. (7)

This equation can be simplified by constraining it to the direction in which the ray parameter distribution attains a maximum, px0

xB: ˆ GG B x_B(x, x0, ω) ∝ ˆ(pxx0_B, ω) ˆu G B_{(s, n, x} 0, pxx0_B, ω). (8) Therefore, ˆGG B

x_B(x, x0, ω) is now constructed by a single Gaussian beam in the direction of the direct arrival from the passive source, pxx0_B.

The forward or downgoing wavefield (D_x

B, Fig.4a) at the instant t is generated for the virtual-source position x0by using the Green’s function approximation of eq. (8):

D_xB(x, x0, t) ≈ 1 π   _∞ 0 ˆ GG B_xB (x, x0, ω) ˆSx_B(ω) e−iωtdω , (9)

where ˆSx_B(ω) stands for the source function of the corresponding passive source. The source function can be estimated from the direct arrival

of the wavefield, depending on the transient behaviour of the passive source. If this is not the case, an approximation can be obtained by isolating the virtual-source function from Cx_Band using it as source function.

For the backprojection of the receiver wavefield, we build the asymptotic form of the correlation function using the Gaussian beam approximation from eq. (6) and adapt the Kirchhoff integral to the boundary defined by the receivers at x_A:

ˆ CG B x_B (s, n, x0, θ, φ, ω) = _ˆ (θ, φ, ω)∗  xA  ˆuG B_{(s, n, x} A, θ, φ, ω)∗CˆxB(xA, x0, ω) dxA. (10)

Therefore, the receiver or upgoing wavefield (U_x

B, Fig.4b) at an instant t is calculated at the locations x by summing the Gaussian beam form of the correlation function CG B

x_B in all directions: U_x B(x, x0, t) = −2 π  _π 0  2π 0   ∞ 0 ˆ CG B x_B (s, n, x0, θ, φ, ω) e−iωtdω dθ dφ. (11)

The upgoing wavefield U_xBcontains the autocorrelation of the source signal provided by the correlation function Cx_B(xA, x0, t).

The estimation of the backprojection of the correlation function has been described here following the Gaussian beam summation method (Popov et al.2010). However, unlike the forward wavefield, the construction of the receiver wavefield is not necessarily constrained to this specific method of wavefield reconstruction.

The zero-time-lag correlation of the two wavefields D_x

Band Ux_Bset the imaging condition to the image result (IxB, Fig.4c):

IxB(x, x0)=

 T 0

D_x

B(x, x0, t) Ux_B(x, x0, t) dt, (12)

where IxB(x, x0) is the partial image produced by the passive source xBand illuminated by the virtual source at x0. The contribution of every

virtual source completes the final image: Ix_B(x)=



x0

Ix_B(x, x0) dx0. (13)

The result obtained in Ix_Bidentifies that part of the medium that can be reliably imaged for the limited ray parameter provided by the single

1028 C. Almagro Vidal et al.

Figure 5. (a), (b) and (c) Partial image results Ix_B(x) obtained from sources xB1, xB2and xB3, following scenarios in Figs3a,3b and3c. (d), (e) and (f) Same

as in (a), (b) and (c) using pre-stack depth migration (WEM) on conventional passive SI imaging.

5 S Y N T H E T I C R E S U L T S

We use the 2-D acoustic scenarios depicted in Figs3(a), (b) and (c); The three scenarios share the same acoustic model and an array with 41 receivers at the free surface (yellow triangles, both x_Aand x0between 2000 and 4000 m, with 50 m spacing), and a different single passive

source in each of the cases (red stars). In these results no taper is applied to the array edges. To obtain the migration results from each passive source we use the correlation imaging condition described in eq. (12) and integrate over all the virtual sources as in eq. (13).

Fig.5displays the synthetic results for the three different source locations. Depending on the location of the passive source, the results show that only the parts of the reflectors that contain specular-reflection points are imaged.

Figs5(a), (b) and (c) show the three respective partial images obtained by using the directionally constrained migration scheme. These results can be compared to the respective results in Figs5(d), (e) and (f), which were produced when applying a conventional pre-stack depth migration (wave-equation migration, WEM) to the virtual-source records retrieved for all 41 receiver positions. Imaging results from using this conventional scheme show migrated artefacts from the correlation function that are not in stationary phase within the array and thus are not specular-reflection points. This is observable when comparing the synclinal structure at 1400 m depth in Figs5(b) and (e), or the area above the shallowest reflector in Figs5(a) and (d), and Figs5(c) and (f).

6 I M A G E I N T E R F E R O M E T RY

Conventional passive SI imaging retrieves the reflection response prior to applying an imaging scheme. This method integrates the correlation results of all passive sources x_B, from a well-sampled distribution of passive sources in the subsurface:

R₃  (x_A, x0, t) ∗ Sac(t)∝ H(t)  xB vobs 3 (xA, xB, t) ⊗ v obs 3 (x0, xB, t) dxB, (14)

where⊗ symbolises the cross-correlation product, H(t) is the Heaviside function and Sac(t) is an average of the auto-correlation of the

respective passive sourcesSac(t)= Sx_B(t)⊗ Sx_B(t)



. Once the reflection response has been retrieved, a standard active imaging technique is applied to the virtual-source reflection responses (assuming well sampled receivers xAand virtual sources x0):

I (x)∝  x0  t G (x, x0, t) ·   x_A G(x, x_A, t) ⊗  R₃  (x_A, x0, t) ∗ Sac(t) dx_A  dt dx0, (15)

Figure 6. Depiction of different orders of surface-related multiples and their dominant ray parameters with respect to the virtual source x0, and the implication for directionally constrained migration. (a) Surface multiple scattering from the passive source at x_B(red star) arrives at the surface receiver at x0(yellow triangle) with different dominant ray parameters: px0

x_B, pxM01, p

x0

M2, etc. The corresponding specular-reflection points to these events P, M1and M2, are located

along their corresponding specular rays and dominant ray parameters. (b) The direction constraint limits the migration to a single ray parameter (px0

x_B) that

images the primary energy into its correct location in the medium. For ray angles other than normal incidence, this constraint causes the multiple energy to be imaged at different locations (M’1and M’2), thus reducing the total imprint of the free-surface multiples in the partial image result.

where G stands for the Green’s function of the medium with respect to the receiver locations. Following Schuster (2001), we change the order in which the integrals are put into effect: In order to obtain an image result due to an individual passive source at x_Bwe rewrite eq. (15):

IxB(x)∝  x0  t G (x, x0, t) ×   xA G(x, x_A, t) ⊗  H (t)vobs 3 (xA, xB, t) ⊗ v obs 3 (x0, xB, t) dxA  dt dx0, (16)

where we image first and subsequently integrate over the passive sources: I (x)=



x_B

Ix_B(x) dxB. (17)

This procedure of interchanging the integral order has previously been applied in Artman (2006), where he combined the observed wavefields vobs

3 with wavefield extrapolation operators to build the upgoing and downgoing imaging wavefields separately first, and subsequently

correlated them in the image domain.

We employ the term Image Interferometry (II, Schuster et al.2004) for the integration of the partial image results due to individual sources in the subsurface. The partial images Ix_Bcomplement one another and the result of this integration produces an output equivalent to

the one in eq. (15). In the case of utilizing directionally constrained migration, the downgoing wavefield G(x, x0, t) in eq. (16) is replaced by

the directionally customized downgoing wavefield D_x

B(x, x0, t) from eq. (9). With this substitution, II profits from partial image results of

the subsurface with fewer artefacts, therefore producing a more reliable final image result of the subsurface. By applying II with directionally constrained migration, the retrieval of the reflection response as an intermediate result is obviated and, although it still relies on a densely sampled array of receivers, it allows to have a sparsely sampled passive-source distribution.

In II with directionally constrained migration, only the correlated events from primaries prevail in the final image result due to constructive interference. The imaged events identified as correlation artefacts are suppressed in the final image result since they do not coincide at the same time/depth in the different partial image results. In the case of an oblique incidence in the directional constraint, coherent events other than primaries, such as free-surface multiples, also are reduced. In Fig.6(a) we depict how from the virtual-source location every free-surface

1030 C. Almagro Vidal et al.

Figure 7. Comparison of migrated results. (a) Image Interferometry due to the directionally constrained migrated results of the three individual passive sources.

(b) Prestack depth migration of conventional passive SI imaging using pre-stack depth migration (WEM) with the same three passive sources.

multiple defines a specific ray parameter. The directionally constrained migration only considers the dominant direction of the virtual-source function, and obviates the imaging directions that would correspond to the surface-related multiples. Hence, free-surface multiples are imaged along the migration path of primaries (see Fig.6b). Unless the migration path is close to normal incidence, free-surface multiples are wrongly imaged at different locations for every virtual source x0and passive source xB, thus reducing their imprint in the final image result.

In II we integrate eq. (17) by summing together the individual partial images Ix_Bwith weight factors:

I (x)=

x_BWxB(p x0

x_B, Sx_B) Ix_B(x), (18)

where the weights Wx_Bserve to balance the strength and contribution of events from different angles. This weighting process may additionally

include frequency balancing to correct for the different frequency spectra the sources may have. Also, in case source functions have varying frequency content, strength and transient behaviour, it may help to use other definitions of the correlation function based on alternative SI methods (such as deconvolution (Mehta et al.2007; Vasconcelos & Snieder2008) and multidimensional deconvolution (Wapenaar et al.

2011; Nakata et al.2014; Hartstra et al.2017).

Fig.7(a) is the result of stacking the partial image results in Figs5(a), (b) and (c). The artefacts from imaging the individual sources are now largely suppressed by destructive interference. Fig.7(b) is the conventional pre-stack migration result using WEM with the virtual-source records retrieved from cross-correlating sequentially first and adding consecutively the three individual passive sources (following eqs14and

15). The latter result is the same as stacking the partial image results of Figs5(d), (e) and (f). In Fig.7(b), the migration of events of the correlation function that are not in stationary phase leaves incorrect events clearly visible in the centre of the synclinal structure at 1400 m depth. In Fig.7(a), the imprint of migration artefacts is reduced in the area above the shallow reflector and between the latter and the deep reflector.

7 F I E L D - D AT A E X A M P L E

Nowack et al. (2010) imaged the lithosphere with Gaussian beams by using teleseismic body-wave events. Their approach stacked multiple teleseismic events from limited azimuths and applied an adapted version of the Gaussian beam migration scheme presented by Hill (2001). In our implementation, we employ regional earthquakes independently in order to image the lithosphere using the directionally constrained migration, which is based on the Gaussian beam summation method by Popov et al. (2010).

7.1 Veracruz-Oaxaca line (VEOX)

We tested the migration scheme on the recordings of local earthquakes in Southern Mexico, obtained by the passive acquisition array “Veracruz-Oaxaca line, VEOX” (Caltech 2010). This array contained 45 three-component broad-band receivers deployed in line across the isthmus of Tehuantepec (ca. 280 km long). The VEOX line has been employed to study the subduction of the Cocos Plate under the North American Plate (Kim et al.2011; Chen & Clayton2012). Tomography along the VEOX line showed the lithospheric structure of the North American Plate, the subduction of the Cocos Plate beneath it, and its interaction with a subduction plate from the North (Melgar & P´erez-Campos2011).

The tectonic complexity of the region explains the presence of seismic tremors throughout the area with source locations at a broad range of depths (Kim et al.2012). We identify the active seismic tremors during the continuous recording period of the array (March 2007 to January 2009) with the aid of the earthquake browser of the Incorporated Research Institutions for Seismology (IRIS). Since the array covers

Figure 8. (a), (b) and (c) Geographical setting of the field data, with epicentre and depth of the passive source (red star) and receivers employed for migration.

(d), (e) and (f) Earthquake recordings. Earthquake names based on IRIS ID. Red shaded area represents the S-wave arrivals that are removed before processing.

a large aperture, we focus on seismic tremors in the region with magnitudes larger than 4 Mw. To produce the field data results we employed

earthquakes that featured favourable radiation conditions in order to avoid dealing with polarity reversals. Also, because the receiver array has an irregular spacing and suffered from strong noise variations, elastic wavefield separation is not performed. Therefore, we concentrated on earthquakes with source locations at depths larger than 120 km and separated the P-wave scattering wavefield from the S-wave arrivals with a smooth time-gating. We work with the vertical particle velocity in order to neglect the P to S conversions occurring before the S-wave arrivals. However, we are aware that this is a crude approximation.

We employ a laterally invariant velocity model with values for the crustal and upper mantle based on the standard velocity model iasp91 (Kennett & Engdahl1991). The use of this model to conduct the migration compromises the image interferometry result. The simplicity of the velocity model neglects the complexity of the crustal structure in depth, and we ignore whether its velocity/depth values correspond to the crystalline basement in this region of the Earth. Also, the heterogeneities of the crustal lithosphere along the receiver line are ignored. The region has positive topography features at both Southern and Northern limits that could produce strong variations in the actual velocity of the lithosphere. Since we do not know the nature of these possible velocity changes, we cannot account for them in the model.

Thus far, we employed SI by cross-correlation to construct the correlation function that we employ for migration. In this section though, in order to improve resolution, we construct an alternative function to the correlation function by applying an array based source-signal deconvolution to the passive recording. We use this approach as an alternative to the trace-by-trace deconvolution since it balances the SNR over the array response and the removal of the source function is applied to the whole array simultaneously. The alternative function Mx_Bis

1032 C. Almagro Vidal et al.

Figure 9. Same as in Fig.8, but for different earthquake locations and their respective recordings.

the result of implementing the following multidimensional deconvolution formula (Hartstra et al.2017):

ˆ Mx_B =  ˆ vobs 3 − ˆv obs 3,dir  ˆ vobs 3 † ˆ vobs 3  ˆ vobs 3 †₊ 2 I −1 . (19)

In this equation, the recorded responses of the receivers in the array ˆvobs

3 are organized as a column vector ˆv obs

3 . Moreover, ˆv

obs

3,dir is a

time-windowed estimate of the direct arrival of the observed wavefield,{}†represents transposition and complex conjugation, I is the identity matrix and is a stabilization factor. The employment of the array deconvolution approach describes an ill-conditioned problem. The inversion in eq. (19) removes the source function of the seismic tremor and improves the time resolution of the correlation result. We applied a fourth order bandpass Butterworth filter between 0.001 and 2.0−4.5 Hz, depending on the passive recording, before estimating the correlation function. After deconvolution we utilized the same filter for all the correlation functions with a band-pass from 0.001 to 0.9 Hz. All virtual sources were employed without the use of array tapers.

The earthquake settings mapped in Figs8(b) and (c) correspond to tremors with epicentres oriented along the array, which is optimal for our 2-D imaging approach and assuming a laterally invariant medium. The earthquake settings in Fig.8(a) and Figs9(a), (b) and (c) show the epicentres detached from the receiver lines. The dominant frequency from the corresponding correlation functions was not higher than 0.3 Hz. This allowed the Moho reflection to be near the edge but within the Fresnel zone of the array, and thus we can apply our 2-D imaging approach to this arrival. In Figs8(d), (e), (f), and9(d), (e) and (f), the seismograms are displayed with the smooth time-windows employed to remove the S-wave arrivals depicted in red. For every seismic tremor, we work with different amounts of receivers due to the fact that some

Figure 10. (a), (b) and (c) Image results obtained using directionally constrained migration from the earthquakes displayed in Figs8d, e and f, respectively. In orange, the interpretation of the basin of Veracruz. The red line depicts the interpretation of the Moho. The blue line in (b) shows the interpretation of converted-wave migration artefact. For comparison, (d) shows the receiver function (RF) result in the region, from Kim et al. (2011).

Figure 11. Same as in Fig.10, but for the earthquakes displayed in Figs9(d), (e) and (f).

stations were inoperative during the recording time or were discarded because of extreme variations in the SNR between receivers. While conventional reflection-response retrieval relies on having the same receivers operative for all passive source recordings, the directionally constrained migration scheme we propose does not.

The migration results are shown in Figs10and11. At 15 km depth, results Figs10(a) and (c) depict the same reflector from lateral position 150–220 km along the line (orange line), which could determine the bounds of the sedimentary basin of Veracruz. Between 30 and 50 km depth a strong dipping reflector is imaged (red line) that we could interpret as the Moho. However, the imprecision of the results, possibly caused by the inaccuracy of the velocity model, and the low SNR of the tremor recordings, prevents confirming whether the bottom of the crust has been mapped accurately in depth.

Displays in Fig.11, show the image results of earthquakes that featured a higher SNR and permitted a clearer migration result. In all three resultsFigs 11(a), (b) and (c) the same strong dipping reflector at 30–50 km depth is identifiable, with small variations in depth among the results, probably due to the inaccurate velocity model employed during the migration. Nevertheless, the inaccuracy of the velocity model and in this case also the mislocation of the hypocentre of the tremors with respect to the array leave doubts whether this reflector is properly positioned in the image results. This circumstance dissuaded us from producing the result of stacking the migration results: The image-interferometry result would not yield a proper constructive interference and the final result would have been difficult to interpret.

1034 C. Almagro Vidal et al.

Figure 12. Numerical validation experiment: (a) Synthetic-earthquake recording. Red shaded area depicts the S-wave arrivals that are removed before

processing. (b) Compressional-velocity model. Receiver locations at the surface and double-couple source represented by the “beach-ball”. (c) Image result obtained using directionally constrained migration from the synthetic earthquake displayed in (a), using the same laterally invariant velocity model as employed for the field-data results. In red, the interpretation of the synthetic Moho. The blue line depicts the interpretation of the Spconverted-wave artefact.

7.2 Numerical validation of field data results

In order to validate the results of imaging the Moho we perform a numerical modelling based on the recording of earthquake 2482604 (Figs9a, d and11a). For the elastic 2-D forward modelling we use the velocity model displayed in Fig.12(b). This velocity model is inspired by the geometry of the lithosphere on the results of Melgar & P´erez-Campos (2011). We use exactly the same velocity values of the lithosphere as from the standard velocity model iasp91 (Kennett & Engdahl1991). Although very simple, the only purpose of this synthetic model is to migrate the synthetic reflections which correspond to the 60 km of the model. We isolate the signal from the direct arrival of the recording and employ it as source function. The location of the source is set at the estimated depth along the subduction slab and we use a double-couple model as source mechanism. The receiver array imitates the geometry with irregular spacing in the recording of the reference earthquake.

During the migration of the synthetic-earthquake recording we applied the same processing steps as described for the field-data results. Fig.12(a) shows the seismograms of the synthetic result and the smooth time window to remove the S-wave arrivals. The result of applying the directionally constrained migration to the synthetic earthquake recording is shown in Fig.12(c). The main imaged reflector corresponds to the synthetic Moho at depths between 30 and 50 km. The source mechanism was oriented such that the radiation characteristics of the source delimit the section of the reflector that can be imaged, which could explain the lack of continuity of the corresponding reflector in the field-data results. The imaged artefacts with negative amplitudes (white features, see blue arrow and line in Fig.12c) correspond to the correlation of the direct P and the S-P converted wave (Sp). This artefact could also be affecting the field-data imaging results. The strength

of this artefact varies with the source location, the source mechanism and the difference between the P- and S-wave velocity in the mantle. For instance, the imaged reflector in Fig.10(b) could be strongly affected by such artefact since it is obtained from the shallowest of the analyzed earthquakes (see blue arrow and line). The observed P-wave direct arrivals of the earthquakes employed for the results in Figs11(a), (b) and (c) feature strong amplitudes at short offsets and no polarity changes along the recordings. Such source radiation characteristics limit the presence of the converted waves to large offsets in the vertical particle-velocity recordings. Should the source mechanism be oriented featuring maximum shear-wave amplitudes at short offset, we would observe stronger converted-wave arrivals. This synthetic result confirms the potential of imaging lithospheric structures such as the Moho by applying directionally constrained migration to regional earthquake recordings, although converted waves in the P-coda can affect the interpretation.

8 D I S C U S S I O N

Since in passive seismics the occurrence of sources in the subsurface is beyond our control and may happen over long periods of time, through the introduced migration scheme we propose to apply the virtual-source and receiver integrations of the migration process first, and leave the interferometric integration (i.e. the summation over the earthquake sources) for the last. This makes it possible to obtain an increasingly

either velocity model, depending on the P- or S field used to obtain the correlation function. This implementation is expected to produce independent results (PP, SP, etc.) with coinciding imaged specular-reflector points (per wave type) and stronger destructive interference for correlation artefacts.

Finally, we emphasize the importance of correctly using the limited information provided by the correlation function. When analyzing a single correlated event in the correlation function, generally two sections can be distinguished: the section which corresponds to specular reflection points and thus is in stationary phase within the array, and the section which does not. When employing the correct velocity model, the backprojection of the correlated event with directionally constrained migration concentrates on imaging only that section of the correlated event that contains specular reflector points. However, the identification of the receiver pair (virtual source and receiver) of the section of the correlated event that corresponds to a specular reflection remains uncertain. In order to resolve the receiver pair identification, a future development could exploit the information brought by midpoint interferometry. This technique replaces the integration over multiple passive sources by the integration over receiver pairs for a single source recording. This analysis has shown positive results for laterally invariant media (Ruigrok & Almagro Vidal2013). The extraction of this information would allow the application of one-way traveltime tomography between three terms (virtual-source position, receiver location and specular-reflection point) in order to update the initial velocity model of the subsurface.

9 C O N C L U S I O N S

We presented a passive migration method for generating partial reflection images from a limited number of subsurface sources. Our scheme takes the illumination characteristics of the passive sources into account. It uses this information to image only energy in stationary phase for the corresponding virtual source, thus limiting the migration of correlated energy that would only contribute to migration artefacts. In case of limited and irregular passive-source distributions, the scheme produces better results than conventional SI imaging.

Under specific circumstances, the explicit reconstruction of the Green’s function as an intermediate step is not necessary to image the subsurface. The contribution from an individual passive source can resolve reflector geometries in the subsurface. This could be further improved with the eventual addition of images from other passive sources. This process of adding images produced by individual passive sources enhances and complements the imaging of reflectors, and produces cancellation of the already limited migration artefacts and non-physical correlated events. By this process, we are postponing the use of the interferometric integration (i.e. the summation over the passive sources) to the image domain, thus obviating the explicit reconstruction of the complete reflection response.

The application of this method on field data has allowed us to obtain imaging results of the lithosphere and interpret the Moho. Although consistent with the interpretation of previous studies, the small variations of the depth of the imaged interface between different partial images emphasize the importance of a model of the subsurface with accurate velocity values before performing successful interferometric imaging.

A C K N O W L E D G E M E N T S

The facilities of IRIS Data Services, and specifically the IRIS Data Management Center, were used for access to the waveforms and related metadata used in this study. IRIS Data Services are funded through the Seismological Facilities for the Advancement of Geoscience and EarthScope (SAGE) Proposal of the National Science Foundation under Cooperative Agreement EAR-1261681. Arie Verdel’s paper contribution has been partly funded by the European Union’s Horizon 2020 research and innovation programme under grant agreement number 727550 (GEMex). We would like to thank the editor Martin Schimmel and reviewers Carlos da Costa Filho and Robert Nowack for providing very constructive comments that helped improve this work. We thank Jan Thorbecke for providing the numerical codes to produce the synthetic results (https://janth.home.xs4all.nl/). We are grateful to Kasper Van Wijk for his assistance with the IRIS Data Management Center. We also thank Elmer Ruigrok for thoughtful discussions and comments on this work.

R E F E R E N C E S

Abe, S., Kurashimo, E., Sato, H., Hirata, N., Iwasaki, T. & Kawanaka, T., 2007. Interferometric seismic imaging of crustal struc-ture using scattered teleseismic waves, Geophys. Res. Lett., 34(19), doi:10.1029/2007GL030633.

Almagro Vidal, C., Draganov, D., van der Neut, J., Drijkoningen, G. & Wapenaar, K., 2014. Retrieval of reflections from ambient noise using illumination diagnosis, J. Geophys. Int., 198(3), 1572–1584.

Artman, B., 2006. Imaging passive seismic data, Geophysics, 71(4), SI177– SI187.

1036 C. Almagro Vidal et al.

Bou´e, P., Poli, P., Campillo, M., Pedersen, H., Briand, X. & Roux, P., 2013. Teleseismic correlations of ambient seismic noise for deep global imaging of the Earth, J. Geophys. Int., 194(2), 844–848.

Caltech, 2010. Veracruz-Oaxaca Subduction Experiment Dataset. Campillo, M. & Paul, A., 2003. Long-range correlations in the diffuse

seis-mic coda, Science, 299(5606), 547–549. ˇ

Cerven´y, V., Popov, M.M. & Pˇsenˇc´ık, I., 1982. Computation of wave fields in inhomogeneous media-Gaussian beam approach, J. Geophys. Int., 70(1), 109–128.

Chen, T. & Clayton, R.W., 2012. Structure of central and southern Mexico from velocity and attenuation tomography, J. Geophys. Res., 117(B9), B09302, doi:10.1029/2012JB009233.

Draganov, D., Wapenaar, K., Mulder, W., Singer, J. & Verdel, A., 2007. Retrieval of reflections from seismic background-noise measurements,

Geophys. Res. Lett., 34(4), L04305–1-4, doi:10.1029/2006/GL028735.

Draganov, D., Campman, X., Thorbecke, J., Verdel, A. & Wapenaar, K., 2010. Event-driven seismic interferometry with ambient seismic noise, in Proceedings of the 72nd EAGE Conference & Exhibition, European Association of Geoscientists & Engineers, Barcelona, Spain.

Galetti, E. & Curtis, A., 2012. Generalised receiver functions and seismic interferometry, Tectonophysics, 532, 1–26.

Hartstra, I., Almagro Vidal, C. & Wapenaar, K., 2017. Full-field multidi-mensional deconvolution to retrieve body-wave reflections from sparse passive sources, J. Geophys. Int., 210(2), 609–620.

Hill, N.R., 1990. Gaussian beam migration, Geophysics, 55(11), 1416–1428. Hill, N.R., 2001. Prestack gaussian-beam depth migration, Geophysics,

66(4), 1240–1250.

Kennett, B. & Engdahl, E., 1991. Traveltimes for global earthquake location and phase identification, J. Geophys. Int., 105(2), 429–465.

Kim, Y., Clayton, R.W. & Keppie, F., 2011. Evidence of a collision between the Yucat´an Block and Mexico in the Miocene, J. Geophys. Int., 187(2), 989–1000.

Kim, Y., Clayton, R.W. & Jackson, J.M., 2012. Distribution of hydrous min-erals in the subduction system beneath Mexico, Earth Planet. Sci. Lett.,

341, 58–67.

Lin, F.-C., Li, D., Clayton, R.W. & Hollis, D., 2013. High-resolution 3d shallow crustal structure in Long Beach, California: Application of am-bient noise tomography on a dense seismic array. noise tomography with a dense array, Geophysics, 78(4), Q45–Q56.

Mehta, K., Snieder, R. & Graizer, V., 2007. Extraction of near-surface prop-erties for a lossy layered medium using the propagator matrix, J. Geophys.

Int., 169(1), 271–280.

Melgar, D. & P´erez-Campos, X., 2011. Imaging the Moho and subducted oceanic crust at the Isthmus of Tehuantepec, Mexico, from receiver func-tions, Pure Appl. Geophys., 168(8-9), 1449–1460.

Nakata, N., Snieder, R. & Behm, M., 2014. Body-wave interferometry using regional earthquakes with multidimensional deconvolution after wave-field decomposition at free surface, J. Geophys. Int., 199(2), 1125–1137. Nakata, N., Chang, J.P., Lawrence, J.F. & Bou´e, P., 2015. Body wave ex-traction and tomography at Long Beach, California, with ambient-noise interferometry, J. Geophys. Res., 120(2), 1159–1173.

Nishida, K., 2013. Global propagation of body waves revealed by cross-correlation analysis of seismic hum, Geophys. Res. Lett., 40(9), 1691– 1696.

Nowack, R.L., Dasgupta, S., Schuster, G.T. & Sheng, J.-M., 2006. Correla-tion migraCorrela-tion using Gaussian beams of scattered teleseismic body waves,

Bull. Seismol. Soc. Am., 96(1), 1–10.

Nowack, R.L., Chen, W.P. & Tseng, T.L., 2010. Application of Gaussian-beam migration to multiscale imaging of the lithosphere beneath the Hi-CLIMB array in Tibet, Bull. Seismol. Soc. Am., 100(4), 1743–1754.

Olivier, G., Brenguier, F., Campillo, M., Lynch, R. & Roux, P., 2015. Body-wave reconstruction from ambient seismic noise correlations in an under-ground mine, Geophysics, 80(3), KS11–KS25.

Poli, P., Pedersen, H. & Campillo, M., 2012. Emergence of body waves from cross-correlation of short period seismic noise, J. Geophys. Int., 188(2), 549–558.

Popov, M.M., 1982. A new method of computation of wave fields using Gaussian beams, Wave Motion, 4(1), 85–97.

Popov, M.M., Semtchenok, N.M., Popov, P.M. & Verdel, A.R., 2010. Depth migration by the Gaussian beam summation method, Geophysics, 75(2), S81–S93.

Rickett, J. & Claerbout, J., 1999. Acoustic daylight imaging via spectral fac-torization: Helioseismology and reservoir monitoring, The Leading Edge,

18(8), 957–960.

Roux, P., Sabra, K.G., Gerstoft, P., Kuperman, W. & Fehler, M.C., 2005. P-waves from cross-correlation of seismic noise, Geophys. Res. Lett.,

32(19), L19303.

Ruigrok, E. & Almagro Vidal, C., 2013. Body-wave receiver-pair seismic interferometry, in Proceedings of the 75th EAGE Conference & Exhibition

incorporating SPE EUROPEC, European Association of Geoscientists &

Engineers, London, United Kingdom.

Ruigrok, E., Campman, X., Draganov, D. & Wapenaar, K., 2010. High-resolution lithospheric imaging with seismic interferometry, J. Geophys.

Int., 183(1), 339–357.

Schuster, G.T., 2001. Theory of daylight/interferometric imaging-tutorial, in

Proceedings of the 63rd EAGE Conference & Exhibition, European

As-sociation of Geoscientists & Engineers , Amsterdam, The Netherlands. Schuster, G.T., 2009. Seismic Interferometry, Cambridge University Press,

ISBN-13:9780521871242.

Schuster, G.T., Yu, J., Sheng, J. & Rickett, J., 2004. Interferometric/daylight seismic imaging, J. Geophys. Int., 157(2), 838–852.

Snieder, R., 2004. Extracting the Green’s function from the correlation of coda waves: a derivation based on stationary phase, Phys. Rev. E, 69(4), 046610.

Snieder, R., Wapenaar, K. & Larner, K., 2006. Spurious multi-ples in seismic interferometry of primaries, Geophysics, 71(4), SI111–SI124.

Tonegawa, T., Nishida, K., Watanabe, T. & Shiomi, K., 2009. Seismic in-terferometry of teleseismic S-wave coda for retrieval of body waves: an application to the Philippine Sea slab underneath the Japanese Islands, J.

Geophys. Int., 178(3), 1574–1586.

van der Neut, J.R., Ruigrok, E.G., Draganov, D. & Wapenaar, K., 2011. Re-trieval of reflections from ambient noise using the incident fields (point-spread function) as a diagnostic tool, in Proceedings of the 73rd EAGE

Conference & Exhibition, European Association of Geoscientists &

En-gineers , Vienna, Austria.

Vasconcelos, I. & Snieder, R., 2008. Interferometry by deconvolution: Part 1 - theory for acoustic waves and numerical examples, Geophysics, 73(3), S115–S128.

Wapenaar, K., 2003. Synthesis of an inhomogeneous medium from its acous-tic transmission response, Geophysics, 68(5), 1756–1759.

Wapenaar, K. & Fokkema, J., 2006. Green’s function representations for seismic interferometry, Geophysics, 71(4), SI33–SI46.

Wapenaar, K., van der Neut, J., Ruigrok, E., Draganov, D., Hunziker, J., Slob, E., Thorbecke, J. & Snieder, R., 2011. Seismic interferometry by crosscorrelation and by multidimensional deconvolution: a systematic comparison, J. Geophys. Int., 185(3), 1335–1364.

Weaver, R.L. & Lobkis, O.I., 2001. Ultrasonics without a source: thermal fluctuation correlations at MHz frequencies, Phys. Rev. Lett., 87(13), 134301.