Passive body-wave interferometric imaging with directionally constrained migration

Passive seismic interferometry enables the estimation of the reﬂection 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 reﬂection response in order to study the subsurface. Successful interferometric imaging relies on the availability of passive recordings from sufﬁcient 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 reﬂection 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 waveﬁeld to the dominant radiation direction of the virtual source. This constraint enables the construction of accurate images from individual sources with a signiﬁcantly 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 reﬂection-response domain, thus eliminating the need for the retrieval of the reﬂection response as an intermediate step.


I N T RO 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 & 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 Passive body-wave interferometric imaging with directionally constrained migration 1023 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.

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Ĝ(x A , x 0 , ω) between a receiver at x A and a virtual source at x 0 from recordings at these two positions of a continuous distribution of passive sources (individually located at x B ). 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,R 3 : 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 wavefieldv obs , is the vertical particle-velocity Green's function G 3,3 (x 0 , x B , ω) due to a vertical point-force at x B , multiplied by the Fourier transform of the source functionŜ(ω). The product on the right hand side of eq. (1) corresponds to a cross-correlation in the time domain. The integration over x B is defined by the distribution of passive point-force sources in the medium, and ρ and c P 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 responseR 3 (x A , x 0 , ω) (this is the representation for a vertical point-force source at x 0 and vertical particle-velocity wavefield at x A ), multiplied by the power spectrum of the source function |Ŝ(ω)| 2 =Ŝ(ω){Ŝ(ω)} * .
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Ŝ(ω), 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): whereĈ x B (x A , x 0 , ω) stands for the correlation function of a single passive source at x B : Eqs (2) and (3) ignore scaling factors ρ and c P at the source locations x 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Ĉ x B that is in stationary phase, and in minimizing correlated artefacts during the imaging process.

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 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 x 0 ) 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 C x 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).
Downloaded from https://academic.oup.com/gji/article-abstract/215/2/1022/5060754 by Technical University Delft user on 07 September 2018 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 x A records a field originating from a subsurface source (x B ) after being scattered by a reflector. A receiver at x 0 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 x A with the one at x 0 retrieves the reflection response at x A as if a source were located at receiver x 0 (red triangle). The locations of the passive source and the virtual source define a unique ray parameter (p 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 C where p is the ray parameter vector, x H,A correspond to the horizontal coordinates of x A and C x B is the illumination distribution of the virtual-source function for each virtual source x 0 . 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: A display of the illumination distribution of a virtual-source function C x B is shown in Figs 3(g), (h) and (i) (with their respective dominant ray parameter p x 0 x B ), corresponding to the parabolic slant-stack applied on the virtual-source functions in Figs 3(d), (e) and (f), respectively. All results correspond to the model scenario described in Figs 3(a), (b) and (c), for the same virtual-source location.

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 (Červený et al. 1982). In Popov et al. (2010), the Green's function between location x 0 and any point x in the 3-D medium is represented as the integration of individual Gaussian beams (û G B ) over different directions (described by azimuthal and polar angles θ and φ). This is expressed as: 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 (Popov 1982). 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). 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 C x B (x A , x 0 , 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 c P (x 0 ) at the virtual-source location, we convert the coordinates of the virtual-source function from ray parameters into angular directions: p), using the horizontal-slowness coordinates of the ray parameter p(θ, φ) = ( cos(θ) sin(φ) ). The approximated Green's function due to a directionally constrained virtual-source located at x 0 evaluated at x is weighted according to the radiation pattern This equation can be simplified by constraining it to the direction in which the ray parameter distribution attains a maximum, p x 0 x B : is now constructed by a single Gaussian beam in the direction of the direct arrival from the passive source, Fig. 4a) at the instant t is generated for the virtual-source position x 0 by using the Green's function approximation of eq. (8): whereŜ x 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 C x B and 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 : 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 C G B x B in all directions: The upgoing wavefield U x B contains the autocorrelation of the source signal provided by the correlation function C 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 B and U x B set the imaging condition to the image result (I x B , Fig. 4c): where I x B (x, x 0 ) is the partial image produced by the passive source x B and illuminated by the virtual source at x 0 . The contribution of every virtual source completes the final image: The result obtained in I x B identifies that part of the medium that can be reliably imaged for the limited ray parameter provided by the single passive source x B .

S Y N T H E T I C R E S U LT S
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 A and x 0 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).

I M A G E I N T E R F E RO 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: where ⊗ symbolises the cross-correlation product, H(t) is the Heaviside function and S ac (t) is an average of the auto-correlation of the respective passive sources S ac (t) = S x B (t) ⊗ S x 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 x A and virtual sources x 0 ): 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' 1 and 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 B we rewrite eq. (15): where we image first and subsequently integrate over the passive sources: This procedure of interchanging the integral order has previously been applied in Artman (2006), where he combined the observed wavefields v obs 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 x 0 and passive source x B , thus reducing their imprint in the final image result. In II we integrate eq. (17) by summing together the individual partial images I x B with weight factors: where the weights W x 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 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. 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). a large aperture, we focus on seismic tremors in the region with magnitudes larger than 4 M w . 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.

F I E L D -DATA E X A M P L E
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.
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 M x B is the result of implementing the following multidimensional deconvolution formula (Hartstra et al. 2017): In this equation, the recorded responses of the receivers in the arrayv obs 3 are organized as a column vectorv 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 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.
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.
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.

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.
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 (S p ). 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. 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 S waves 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 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 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.

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 nonphysical 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.