https://orcid.org/0000-0002-1620-8282
SUMMARY
Passive seismic interferometry enables the estimation of the reflection response of the subsurface using passive receiver recordings at the surface 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 distributed, 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 dominant 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 individual 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 illumination 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.
1 INTRODUCTION
Seismic interferometry (SI) aims to reconstruct the impulsive response between receivers, as if one of the receivers were a source (Schuster 2001; Weaver & Lobkis 2001; Campillo & Paul 2003; Wapenaar 2003; Snieder 2004; Schuster 2009; Galetti & Curtis 2012). 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; Nishida 2013). 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 & Claerbout 1999; Abe et al.2007; Draganov et al.2007; Poli et al.2012; Boué 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 reflection response as an intermediate result. The migration of correlated data from individual sources in the subsurface has been referred to as interferometric imaging (II, Schuster et al.2004). Nowack et al. (2006) showed another example of migration of correlated data from individual passive sources, carried out in this case by using slant-stack windows of the data and migrating the autocorrelated data by means of Gaussian beams.
The use of a single passive source does not result in destructive interference of correlated artefacts, which may cause errors during the migration process. Therefore, when migrating these inaccurate correlation responses, the goal is to minimize the appearance of features produced by the migration of the correlated artefacts which are not properly suppressed. This study builds on the work on interferometric imaging for passive seismics and illustrates how the data that are migrated can be limited to the events in stationary phase within the acquisition array. This is achieved by applying directional constraints during the migration process. In this paper, we apply the adapted migration on synthetic passive data and compare it to conventional SI imaging. We also apply this alternative migration scheme on passive field data to perform lithospheric imaging of the Moho.
2 CORRELATION FUNCTION
where ℜ stands for the real part, {}* denotes complex conjugation, ω is the angular frequency and |$\,\,{\hat{}}\,\,$| indicates the wavefield is in the frequency domain. In eq. (1), the observed wavefield |$\hat{v}_{3}^{obs}({\bf x}_{0},{\bf x}_{\mathcal {B}},\omega )=\hat{G}^{}_{3,3}({\bf x}_{0},{\bf x}_{\mathcal {B}},\omega )\hat{S}(\omega )$|, is the vertical particle-velocity Green’s function |$\hat{G}^{}_{3,3}({\bf x}_{0},{\bf x}_{\mathcal {B}},\omega )$| due to a vertical point-force at |${\bf x}_{\mathcal {B}}$|, multiplied by the Fourier transform of the source function |$\hat{S}(\omega )$|. The product on the right hand side of eq. (1) corresponds to a cross-correlation in the time domain. The integration over |${{\bf x}_{\mathcal {B}}}$| is defined by the distribution of passive point-force sources in the medium, and ρ and cP stand 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 |$\hat{R}_{3}^{\circlearrowright }({\bf x}_{\mathcal {A}},{\bf x}_{0},\omega )$| (this is the representation for a vertical point-force source at x0 and vertical particle-velocity wavefield at |${\bf x}_{\mathcal {A}}$|), multiplied by the power spectrum of the source function |$\vert \hat{S}(\omega )\vert ^{2}=\hat{S}(\omega )\lbrace \hat{S}(\omega )\rbrace ^{*}$|.
Eqs (2) and (3) ignore scaling factors ρ and cP at the source locations |${\bf x}_{\mathcal {B}}$| because they are not known in practice.
Fig. 1 shows an example of applying eq. (2) in an acoustic model (Fig. 1a) with three different source distributions. Figs 1(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 scenario 1(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 |$\hat{C}_{{\bf x}_{\mathcal {B}}}$| that is in stationary phase, and in minimizing correlated artefacts during the imaging process.
Three different interferometric results using SI by cross-correlation. (a) Acoustic velocity model with receivers at the free surface (|${\bf x}_{\mathcal {A}}$|, indicated by yellow triangles) and passive sources (|${\bf x}_{\mathcal {B}}$|, indicated by red dots). (b) Reference reflection response |${R}_{3}^{\circlearrowright}$| 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 DIRECTIONAL CONSTRAINTS
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 Figs 2(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.
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 |${{\bf x}_{\mathcal {A}}}$| records a field originating from a subsurface source (|${{\bf x}_{\mathcal {B}}}$|) after being scattered by a reflector. A receiver at x0 records 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 |${{\bf x}_{\mathcal {A}}}$| with the one at x0 retrieves the reflection response at |${{\bf x}_{\mathcal {A}}}$| as 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 (|${{\bf p}^{{\bf x}_{0}}_{{\bf x}_{\mathcal {B}}}}$|). (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 |${{\bf x}_{\mathcal {B}}}$|, is needed.
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 |${C}_{{\bf x}_{\mathcal {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).
Results obtained from illumination diagnosis of the virtual source at x0 = 2500 m from three different sources located at 3000 m depth, at |${x}_{\mathcal {B}}=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 |${C}_{{\bf x}_{\mathcal {B}}}({\bf x}_{\mathcal {A}},{\bf x}_{0},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 |$\widetilde{C}_{{\bf x}_{\mathcal {B}}}({\bf x}_{0},{\bf p})$| with the dominant direction represented by |${\bf p}^{{\bf x}_{0}}_{{\bf x}_{\mathcal {B}}}$| (blue square).
4 MIGRATION SCHEME
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 |$D^{}_{{\bf x}_{\mathcal {B}}}$| and |$U^{}_{{\bf x}_{\mathcal {B}}}$| at instant t = 0.78 s.
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.
5 SYNTHETIC RESULTS
We use the 2-D acoustic scenarios depicted in Figs 3(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}_{\mathcal {A}}$| and x0 between 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. 5 displays 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.
(a), (b) and (c) Partial image results |$I_{{\bf x}_{\mathcal {B}}}({\bf x})$| obtained from sources |${\bf x}_{\mathcal {B}_{1}}$|, |${\bf x}_{\mathcal {B}_{2}}$| and |${\bf x}_{\mathcal {B}_{3}}$|, following scenarios in Figs 3a, 3b and 3c. (d), (e) and (f) Same as in (a), (b) and (c) using pre-stack depth migration (WEM) on conventional passive SI imaging.
Figs 5(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 Figs 5(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 Figs 5(b) and (e), or the area above the shallowest reflector in Figs 5(a) and (d), and Figs 5(c) and (f).
6 IMAGE INTERFEROMETRY
This procedure of interchanging the integral order has previously been applied in Artman (2006), where he combined the observed wavefields |$v_{3}^{obs}$| 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 |${I}_{{\bf x}_{\mathcal {B}}}$| complement 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}_{{\bf x}_{\mathcal {B}}}^{}({\bf x},{\bf x}_{0},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 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 x0 and passive source |${\bf x}_{\mathcal {B}}$|, thus reducing their imprint in the final image result.
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 |${\bf x}_{\mathcal {B}}$| (red star) arrives at the surface receiver at x0 (yellow triangle) with different dominant ray parameters: |${\bf p}^{{\bf x}_{0}}_{{\bf x}_{\mathcal {B}}}$|, |${\bf p}^{{\bf x}_{0}}_{{\bf M}_{1}}$|, |${\bf p}^{{\bf x}_{0}}_{{\bf M}_{2}}$|, etc. The corresponding specular-reflection points to these events P, M1 and M2, are located along their corresponding specular rays and dominant ray parameters. (b) The direction constraint limits the migration to a single ray parameter (|${\bf p}^{{\bf x}_{0}}_{{\bf x}_{\mathcal {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’1 and M’2), thus reducing the total imprint of the free-surface multiples in the partial image result.
where the weights |$W_{{\bf x}_{\mathcal {B}}}$| serve 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 & Snieder 2008) 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 Figs 5(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 eqs 14 and 15). The latter result is the same as stacking the partial image results of Figs 5(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.
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.
7 FIELD-DATA EXAMPLE
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 & Clayton 2012). 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érez-Campos 2011).
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 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 & Engdahl 1991). 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.
The earthquake settings mapped in Figs 8(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 Figs 9(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 Figs 8(d), (e), (f), and 9(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 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.
(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.
Same as in Fig. 8, but for different earthquake locations and their respective recordings.
The migration results are shown in Figs 10 and 11. At 15 km depth, results Figs 10(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.
(a), (b) and (c) Image results obtained using directionally constrained migration from the earthquakes displayed in Figs 8d, 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).
Displays in Fig. 11, show the image results of earthquakes that featured a higher SNR and permitted a clearer migration result. In all three results Figs 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.
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 (Figs 9a, d and 11a). 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érez-Campos (2011). We use exactly the same velocity values of the lithosphere as from the standard velocity model iasp91 (Kennett & Engdahl 1991). 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.
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 Sp converted-wave artefact.
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 Figs 11(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 DISCUSSION
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 complete image of the subsurface, as the responses to more passive sources get recorded and migrated with advancing time. In the case of two passive sources occurring simultaneously at two different source locations in the subsurface, the proposed methodology (in particular eq. 7) identifies the dominant ray parameters of each individual passive source.
Another interesting aspect to point out concerns the number of virtual sources required in the migration process. Conventional migration needs all virtual sources of the array in order to obtain sufficient cancellation of migration artefacts. On the other hand, with the directionally constrained migration scheme the artefacts are limited. This enables to speed up the imaging process by increasing the sparsity of virtual sources employed during migration.
In elastic media, wavefield separation enables independent migration results with complementary illumination, since for non-vertical incidence, P and Swaves follow different ray paths. Depending on the P- or S-virtual-source function analyzed, the forward-propagated source wavefield would employ the corresponding velocity model for its construction. Likewise, the backprojected receiver wavefield would utilize either velocity model, depending on the P- or Sfield 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 Vidal 2013). 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 CONCLUSIONS
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.
ACKNOWLEDGEMENTS
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.
