About this Journal Submit a Manuscript Table of Contents
International Journal of Geophysics
Volume 2011 (2011), Article ID 870819, 28 pages
http://dx.doi.org/10.1155/2011/870819
Research Article

Deghosting, Demultiple, and Deblurring in Controlled-Source Seismic Interferometry

1Department of Geotechnology, Faculty of Civil Engineering and Geosciences, Delft University of Technology, Stevingweg 1, 2628 CN, Delft, The Netherlands
2Earth's Physics Department, Saint Petersburg State University, Ulyanovskaya Street 1, Saint Petersburg 198504, Russia

Received 21 September 2010; Accepted 29 March 2011

Academic Editor: Sergey Fomel

Copyright © 2011 Joost van der Neut et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

With controlled-source seismic interferometry we aim to redatum sources to downhole receiver locations without requiring a velocity model. Interferometry is generally based on a source integral over cross-correlation (CC) pairs of full, perturbed (time-gated), or decomposed wavefields. We provide an overview of ghosts, multiples, and spatial blurring effects that can occur for different types of interferometry. We show that replacing cross-correlation by multidimensional deconvolution (MDD) can deghost, demultiple, and deblur retrieved data. We derive and analyze MDD for perturbed and decomposed wavefields. An interferometric point spread function (PSF) is introduced that can be obtained directly from downhole data. Ghosts, multiples, and blurring effects that may populate the retrieved gathers can be locally diagnosed with the PSF. MDD of perturbed fields can remove ghosts and deblur retrieved data, but it leaves particular multiples in place. To remove all overburden-related effects, MDD of decomposed fields should be applied.

1. Introduction

Seismic interferometry is an effective tool to redatum sources to receiver locations, without the need of a velocity model. Recently, we have seen an increase of various applications; see Curtis et al. [1] and Wapenaar et al. [2]. In this paper we restrict ourselves to controlled-source interferometry for data-driven redatuming. In a recent publication of Schuster [3], numerous applications in this field can be found. Among them is the well-known virtual-source method of Bakulin and Calvert [4].

Traditionally, the theory of interferometry has been derived from a reciprocity theorem of the correlation type or from time-reversal arguments [5, 6]. A few special applications are based on wavefield convolutions [7, 8]. For controlled-source applications, the theory is generally applied with one-sided illumination, meaning that sources are located at the earth’s surface only and are not—as often assumed in theory—enclosing a volume. Moreover, interactions with the free surface and intrinsic losses are generally not taken into account. Because of these factors, spurious events can enter the retrieved gathers [9] and true amplitudes are generally not preserved [10].

To mitigate some of these artifacts, several methods have been proposed. In perturbation-based interferometry [11], incident and scattered wavefields are separated prior to cross-correlation. In the virtual source method [4, 12], a similar separation is achieved by time-gating the direct arrival prior to cross-correlation. Mehta et al. [13] showed that separation of up- and downgoing waves with multicomponent sensors can yield even further improvements. Vasconcelos et al. [14] demonstrate a variety of these methods in complex synthetic subsalt environments.

A different group of methods is based on deconvolution instead of cross-correlation (CC). Replacing cross-correlation by deconvolution can remove undesired multiples from the overburden [15], a concept that has also been referred to as Noah deconvolution [16] or Einstein deconvolution [17, 18]. An additional advantage of deconvolution is that the source wavelet is deconvolved before stacking, which can be beneficial if the source has a complicated signature [19, 20]. Various authors have suggested to redatum data by summing over single-station deconvolution traces [4, 21]. However, to retrieve an exact Green’s function by deconvolution in 3D heterogeneous media, single-station deconvolution should be replaced by multidimensional deconvolution (MDD), as shown by Wapenaar et al. [22]. MDD is based on the inversion of a forward problem that is generally derived for decomposed wavefields. The method has a lot in common with Betti deconvolution, as implemented by Amundsen et al. [23] and Holvik and Amundsen [24] to remove free-surface multiples from ocean bottom cable (OBC) data. For various applications of MDD, see Wapenaar et al. [22]. Van der Neut et al. [25] showed that MDD can correct for attenuation and improve interferometric imaging below complex overburden. Minato et al. [26] applied MDD to virtual cross-well data. MDD can also be applied to ground penetrating radar [27], controlled-source electromagnetic exploration [28, 29], and lithospheric-scale imaging [30].

A typical application of controlled-source seismic interferometry is to redatum sources to a downhole receiver array below a complex overburden. Bakulin and Calvert [4] were pioneering in this field using the so-called virtual source method. A typical configuration is shown in Figure 1(a). Sources are situated at the earth’s surface locations . Receivers are located at and in a well that can be horizontal, deviated, or vertical. The aim is to transform the data obtained with the configuration shown in Figure 1(a) into virtual data as if there was a source at and a receiver at (Figure 1(b)). Like Bakulin and Calvert [4], we will do so without requiring a velocity model, thus bypassing all complexities of the overburden. Schuster’s group considered a range of other configurations [3], one of them being shown in Figure 1(c). Here the aim is to create a virtual source at location by exploiting scattered or dived waves that illuminate the target (e.g., a salt flank) under angles that are unseen in conventional processing (Figure 1(d)); see Xiao et al. [31], Hornby and Yu [32], Lu et al. [33] and Ferrandis et al. [34] for applications. Another application is virtual cross-well acquisition, where and are located in separate wells that can be vertical [26, 35, 36], horizontal [37], or deviated; see Figures 1(e) and 1(f). Many of the formulations that appear in this paper require spatial integrals not only over source locations but also over receiver locations. For a 3D heterogeneous medium, this means that 2D arrays of both sources and receivers should be deployed. Since we assume the presence of downhole receivers, this is generally not feasible. That is why we restrict ourselves to wave propagation in a 2D plane, ignoring out-of-plane reflections.

fig1
Figure 1: Typical redatuming problems: (a) input data for redatuming through complex overburden; (b) output data after redatuming through complex overburden (with some types of interferometry, the overburden is replaced by a homogeneous halfspace); (c) input data for flank imaging; (d) output data for flank imaging after redatuming; (e) input data for virtual cross-well; (f) output data after virtual cross-well; denote the source locations, denote the receiver locations that are turned into virtual sources, and denote other receiver locations.

In this paper we distinguish between ghosts and multiples. With ghosts we refer to spurious events that populate our retrieved gathers, because initial assumptions were not properly fulfilled. With multiples we refer to physical events stemming from multiple reflections. Blurring effects can occur if illumination conditions are imperfect. We analyze the ghosts and multiples that can occur in interferometry by CC of full, perturbed, and decomposed fields. Next we introduce MDD of perturbed and decomposed fields, which can be applied with single-component sensors. We analyze to what extent these methods can be used to remove ghosts, multiples, and blurring artifacts on single- and multicomponent data.

2. Cross-Correlation of Full Fields

Various authors have shown that cross-correlation of wavefields at two receivers followed by summation over a closed boundary of sources can provide a Green’s function as if there was a virtual source at one of the receiver locations and a receiver at the other. In Figure 2 we give a more schematic illustration of the problem formulated in Figure 1(a). The aim is to redatum the source locations to a receiver location “below” the overburden but “above” the target of imaging, without requiring a velocity model of the medium. Note that the terminology of “below” and “above” can be interchanged with “left” and “right” for the situation of salt flank imaging (Figure 1(c)) or the virtual cross-well (Figure 1(e)).

870819.fig.002
Figure 2: Configuration for controlled-source interferometry by cross-correlation of full fields; volume is enclosed by the source array and an additional half sphere whose radius ; denotes a source location at and is the unit normal vector of this surface; denotes a source at and is the unit normal vector of this surface; heterogeneities between sources and receivers are referred to as “overburden,” whereas heterogeneities below the receivers are referred to as “target.”

To retrieve an exact Green’s function, both monopole and dipole sources are required along the closed boundary spanned by and . The medium should be free of intrinsic losses inside . Under these conditions, the following representation can be derived for Green’s function retrieval [38]: On the left hand side we find the Green’s function as if there was an injection rate point source at and a receiver for acoustic pressure at . It is given in the frequency domain, indicated by the hat and angular frequency . We also find its acausal counterpart , where superscript denotes complex conjugation. Note that the retrieved response is bandlimited by , where is the spectrum of the source wavelet. On the right-hand side, we find two integrals: and . is referred to as the “known integral,” obtained by cross-correlations of wavefields from existing source locations at . is the “missing integral,” obtained by cross-correlation of wavefields from nonexisting source locations at . Even though the source locations are not present in a realistic survey, we keep in the representation, allowing us to quantify its contribution.

First, let us consider the known integral , as derived by Wapenaar and Fokkema [38]: Here represents the pressure recording at due to a monopole source at . These recordings are assumed to be Green’s functions (or impulse responses) convolved with source wavelet . Normal vector points perpendicular (outward) to the source array . Further, , where and superscript denotes the transposed. Hence, is interpreted as the response at due to a dipole source at .   is the mass density at the source array.

The representation of missing integral is very similar. In this case responses of nonexisting source locations needs to be cross-correlated and needs to be replaced by the normal vector along , yielding where is the spatial gradient at source location . Before addressing the implications of not evaluating integral , we focus our attention on . Evaluation in its present form would require both monopole and dipole sources at . In practice, interferometry is generally applied with monopole sources only. To overcome this limitation, one often introduces a so-called far-field approximation [3, 38]. This approximation can only be made if the direction of wave propagation with respect to the source array is known. Therefore, we separate ingoing and outgoing wavefields with respect to the volume . We introduce and , where superscripts and denote ingoing and outgoing fields at . We substitute into (2). It can be shown that the cross-correlations of ingoing and outgoing waves cancel and that the remaining terms can be merged [38], such that Next, the far-field high-frequency approximation can be introduced. It is assumed that the medium is smooth in a small region around . We find for ingoing constituents that , where is the acoustic wave velocity along the source array and is the incidence angle of the dominant wavefield with respect to [38]. Similarly, we find for outgoing constituents that . We assume that is close to zero such that . Substituting these approximations into (4) yields Before proceeding, it is helpful to provide a similar analysis for the missing integral . We assume that both density and wave velocity are constant at , such that all wavefields that would be recorded at receivers due to missing source locations are ingoing at (such that and ). Given these considerations, it can be shown that both terms in the integrand of (3) give equal contributions to the stationary points, and therefore this equation can be rewritten as We can further simplify this equation by substituting the far-field approximation for ingoing fields , yielding So far we have shown that a Green’s function can be retrieved by evaluation of integrals and . In practice, we are generally not that fortunate. First, we miss the source locations to compute . Second, we cannot discriminate between ingoing and outgoing wavefields to evaluate . Instead, we cross-correlate the full fields as emitted by the sources and integrate over . We refer to the obtained function as the correlation function : We assume that the density and wave velocity are known at the source array. In most applications of interferometry, however, they are assumed constant and not evaluated inside the integrand. We have introduced an additional weighting factor that can be used to rebalance the sources before integration or to taper the edges of a finite source array in practical applications [39]. We can also define a weighting function that takes the wave propagation angle into account, reducing the artifacts introduced by the far-field approximation, where we neglected a -term. Sometimes we can approximate the incidence angle of the dominant contribution of the incident wavefield , for instance by ray tracing. In such cases, it can be beneficial to replace weighting factor with to improve the retrieved amplitude variations with angle; see Schuster [3] for examples.

The fields in (8) consist of both ingoing and outgoing constituents at . We substitute into (8) and rewrite the result as Note that the first integral at the right-hand side is identical to the desired integral in (5). The second integral is undesired and referred to as ghost : In a similar fashion we identify a second ghost, due to the missing source locations . This is done by redefining the integral in (7) as Next, we substitute (11), (7), (5), (1) and (10) into (9) and rewrite the result as Equation (12) is useful for identifying ghosts in interferometry by CC of full fields. We have shown that implementation of the correlation function (8) yields a bandlimited version of the desired Green’s function and its acausal counterpart plus two additional ghost terms. The first ghost is described by (10), which is due to the presence of any undesired reflectors above the source array . Note that in typical controlled source applications, the free surface acts as such a reflector, scattering waves back into the volume , such that . Consequently, spurious events can be expected to populate retrieved gathers due to free surface interactions if free surface ghosts and multiples are not eliminated prior to applying interferometry. Similar artifacts have also been found in passive seismic interferometry; see Draganov et al. [40]. The second ghost , described by (11), stems from the missing integral at . For convolution-based reciprocity theorems, one often applies Sommerfeld’s radiation condition [41] to show that boundary integrals over convolution products vanish when . However: these conditions do not apply for integrals of cross-correlation products like the one in (11). The integrand in (11) decays with order , whereas the integration surface grows with order , which is insufficient to cancel the integral. However, Wapenaar [42] showed that if sufficient scattering takes place in the volume , the decay of the integrand is stronger than and the integral can indeed be neglected. This condition has been referred to as the antiradiation condition [3]. Not obeying this condition can lead not only to incorrect amplitudes but also to the emergence of spurious events in the retrieved data [9].

To illustrate the ghosts in interferometry by CC of full fields, we define four synthetic 1D models in Figure 3. For simplicity, the velocity is kept constant at 2000 m/s with density contrasts introduced. In each model the aim is to generate a (zero-offset) response as if there was a virtual source at receiver location (green star) and a receiver at the same location . The real source is always located at the earth’s surface location and additional sources are introduced to evaluate ghost . Location will play a role only later in this paper.

fig3
Figure 3: Impedance (in ) versus depth (in m) for synthetic models A–D. The red star corresponds to source location , the green star to receiver/virtual source location , the blue star to nonexisting source location , and the magenta star to nonexisting source location ; models A–C have no free surface incorporated; model D has a free surface at 0 m depth.

In Figure 4(a) we computed the desired reference response with wavelet for model A (Figure 3(a)) by placing an actual source at the virtual source location . We can clearly see the direct arrival at and an event at . The latter event is our target reflection, relating to the impedance contrast at 1200 m depth. Since no free surface is incorporated, no heterogeneities exist above the true source location , and hence all wavefields that arrive at are ingoing at , such that (10). We computed the second ghost , using the additional source below all medium heterogeneities. Since no reflection from this source is registered at , the only contribution to stems from correlations of the direct field, mapping at in this zero-offset case. In Figure 4(b) we show the ghosts . Note that indeed no other contribution is found outside . Therefore the correlation function (8) is very similar to the true reflection response, except around . In Figure 4(c) we show (black) and the sum of reference response and ghosts (dashed green), matching well.

fig4
Figure 4: Cross-correlation of full fields for model A: (a) reference response by placing an active source at the virtual source location ; (b) computed ghost terms ; (c) causal part of the correlation function (black) and the sum of reference response + ghosts (dashed green).

In the previous example, the ghost terms did not give a significant contribution outside . This is not the case for model B (Figure 3(b)), which is the same as model A, except for two additional contrasts at 500 m and 700 m depth. The reference response reveals not only the desired reflection at but also the reflections of the overburden at and and their multiples; see Figure 5(a). The ghosts and the correlation function are shown in Figures 5(b) and 5(c), respectively. Since no reflectors are present above the source, . However, the second ghost does give a significant contribution. The events at and appear with opposite polarity (Figure 5(b)) compared to the reference response (Figure 5(a)). Therefore, these events have incorrect amplitudes in the correlation function and are hardly visible (Figure 5(c)). More importantly, there is a ghost at (Figure 5(c)) that is not visible in the reference response (Figure 5(a)). This ghost originates from an internal multiple between the contrasts at 500 m and 700 m. Finally we note that the multiples at and also appear as ghost terms, with different amplitudes and polarities than the reference response. Therefore these responses are retrieved with incorrect amplitudes (Figure 5(c)). For practical applications this can be seen as advantageous, since these multiple reflections are generally not desired for further processing.

fig5
Figure 5: Cross-correlation of full fields for model B: (a) reference response by placing an active source at the virtual source location ; (b) computed ghost terms ; (c) causal part of the correlation function (black) and the sum of reference response + ghosts (dashed green).

It has been demonstrated by Wapenaar [42] that ghost vanishes if sufficient scattering occurs below the receiver array. To demonstrate this effect, we introduce model C, being similar to model B except for a series of additional random contrasts deeper in the subsurface superimposed by a trend of acoustic impedance increases with depth, see Figure 3(c) (note the differently scaled axes in Figures 3(b) and 3(c)). All contrasts lay sufficiently deep, such that the reference response is no different from that of model B within the display window. The ghosts and the correlation function are shown in Figures 6(b) and 6(c), respectively. Indeed, the contributions of ghost are minor and randomly distributed (Figure 6(b)). It can be shown analytically that placing infinitely many contrasts even completely eliminates . Note that the reference response (Figure 6(a)) and the correlation function (Figure 6(c)) agree relatively well. The so-called antiradiation condition has thus been successful in eliminating the effects of one-sided illumination. Reflections of the target () and overburden ( and ), including their multiples ( and ), have all been retrieved with true amplitudes. Note that recording times need to be sufficiently long for the antiradiation condition to hold [42]. In this example, the total recording time is . Moreover, it should be mentioned that the antiradiation condition can not eliminate the effects of reflectors placed outside the volume , as for such scenario.

fig6
Figure 6: Cross-correlation of full fields for model C: (a) reference response by placing an active source at the virtual source location ; (b) computed ghost terms ; (c) causal part of the correlation function (black) and the sum of reference response + ghosts (dashed green).

In model D we demonstrate the effect of free surface interactions, by placing a free surface at 0 m and two contrasts at 100 m and 200 m depth, see Figure 3(d). As a result, the reference response does not only contain the desired reflection at but also the primaries of the overburden contrasts at and , the free surface reflection at ; and free-surface multiples at and ; see Figure 7(a). Ghosts and correlation function are shown in Figures 7(b) and 7(c), respectively. We see strong spurious events at and in the retrieved response (and a weaker event at ), stemming from interactions of the first reflectors with the free surface. We also observe that the reflections at and , their multiples, and the free surface reflection at are hardly retrieved due to the missing source at .

fig7
Figure 7: Cross-correlation of full fields for model D: (a) reference response by placing an active source at the virtual source location ; (b) computed ghost terms ; (c) causal part of the correlation function (black) and the sum of reference response + ghosts (dashed green).

3. Cross-Correlation of Perturbed Fields

In many cases our aim is not to retrieve a full Green’s function but to retrieve only a part of it. For controlled-source applications, for instance, a full Green’s function would contain not only reflections from the target area but also reflections from the overburden. In practice we often wish to eliminate the latter by restricting a virtual source to radiate downwards only. In the virtual source method [4, 12], this is effectively achieved by time-gating the first arrival prior to cross-correlation. In perturbation-based methodology [11, 14], a similar discrimination is made between the incident field and the scattered field. These fields are usually obtained by time-gating the full fields.

In Figure 8 we show a typical configuration for interferometry by CC of perturbed wavefields. Note the similarities with Figure 2. The only difference is that boundary has been replaced with a boundary , located between the overburden and the target area. Volume is now enclosed by and . We define a reference state by density and velocity . In this reference state all heterogeneities outside volume (namely the target) have been removed. These heterogeneities are referred to as the perturbations in density and velocity , where and represent the density and velocity of the true physical medium. Fields that propagate in the reference state are referred to as incident fields . Fields that propagate in the physical medium can now be expressed as a superposition of the incident field and the so-called scattered field , that is . Vasconcelos et al. [11] have derived a representation for Green’s function retrieval of the scattered field between virtual source and receiver from cross-correlations of the scattered field at and the incident field at : where subscript stands for “perturbed.” On the left-hand side we find the desired scattered Green’s function, imprinted by the squared amplitude spectrum of the source wavelet. Note that no acausal Green’s function is retrieved. On the right-hand side we find integral , stemming from the cross-correlations of incident and scattered fields from the actual sources at : The second integral, , stems from similar cross-correlations of nonexisting sources at : where is the spatial gradient at , is the outward-pointing unit normal vector to , and is the density along this surface. We write the wavefields in (14) in terms of ingoing and outgoing constituents at . It can be shown that the cross-correlations between ingoing and outgoing wavefields cancel each other and the remaining terms can be merged [11], such that Assuming that the relevant wavefields propagate near normal incidence, we can substitute far-field approximations and into (16), yielding For integral the situation is slightly different. Since we have chosen the reference state to have no heterogeneities below , all wavefields arriving at the receivers are ingoing at (such that and ). Cross-correlations between ingoing and outgoing fields cancel and the remaining terms can be merged [11], such that The spatial derivative can be approximated with far-field approximation , such that Note that evaluation of (17) still requires separation of ingoing and outgoing waves. In practice, we generally cross-correlate complete scattered and incident fields at the receivers to obtain the correlation function of perturbed wave fields : To evaluate the consequence of this choice, we separate the scattered and incident wavefields in up- and downgoing constituents according to and and substitute these representations into (20). The result can be written as The first integral in (21) can be identified as the desired integral (see (17)). The second integral is a ghost term that we redefine as The missing sources at make up for a second ghost term that we define by rewriting the integral in (19) in the following way: Now, by substituting (23), (19), (17), (13), and (20) into (21) we arrive at Equation (24) is useful for identifying ghosts in interferometry by CC of perturbed fields. Computation of (see (20)) yields a scaled bandlimited version of the desired scattered Green’s function plus two ghost terms and . The first ghost stems from reflectors outside the integration volume and behaves very similar to the ghost that we found for CC of full fields. The second ghost is relatively small. It consists of ingoing waves at that have scattered in the target area before arriving at the receivers. However, to reach the target area, these waves should have also scattered in the overburden. This means that such fields have at least scattered twice. It is reasoned by Vasconcelos et al. [11] that such contributions can generally be neglected.

870819.fig.008
Figure 8: Configuration for controlled-source interferometry by cross-correlation of perturbed fields; volume is enclosed by the source boundary and an additional boundary between overburden and target; denotes a source location at and is the unit normal vector of this surface; denotes an additional source location at and is the unit normal vector of this surface; heterogeneities between sources and receivers are referred to as “overburden,” whereas heterogeneities below the receivers are referred to as “target.”

To illustrate the effectiveness of cross-correlation of perturbed fields we apply this methodology to model B (Figure 3(b)), with the additional source located at 1000 m depth (instead of the source at ). Incident fields are computed in a reference medium with all heterogeneities below removed. Scattered fields are computed by subtraction of these incident fields from full fields. We generate the reference response by placing an active source at and computing the scattered response at the same receiver; see Figure 9(a). We can clearly see not only the desired reflector at but also the multiple reflections from the overburden at and as well as higher-order multiples. In Figures 9(b) and 9(c) we show the ghosts and the correlation function . Note which the spurious event at , that was visible in Figure 5(c), is significantly weakened and hardly visible. The multiples at and have smaller amplitudes in the correlation function (Figure 9(c)) than in the reference response (Figure 9(a)). This is a consequence of not having the source at .

fig9
Figure 9: Cross-correlation of perturbed fields for model B: (a) reference response by placing an active source at the virtual source location ; (b) computed ghost terms ; (c) correlation function (black) and the sum of reference response + ghosts (dashed green).

To show that CC of perturbed fields does not offer a solution for reflectors above the source array, we apply the methodology to model D (Figure 3(d)) as well. In Figure 10(a) we show the reference response. We observe the target event at and a multiple at . Ghosts and are shown in Figures 10(b) and 10(c), respectively. Note that we can still observe several spurious events and the undesired overburden reflections.

fig10
Figure 10: Cross-correlation of perturbed fields for model D: (a) reference response by placing an active source at the virtual source location ; (b) computed ghost terms ; (c) correlation function (black) and the sum of reference response + ghosts (dashed green).

4. Cross-Correlation of Decomposed Fields

Another way to eliminate ghosts in interferometry is to separate up- and downgoing waves prior to cross-correlation, as proposed by Mehta et al. [13]. If receivers are installed in horizontal wells, decomposition can be implemented by combining pressure and particle velocity recordings along the receiver array [43, 44]. If wave propagation is close to normal incidence, these methods can be approximated by a weighted summation of pressure and particle velocity at a single receiver [13, 36, 45]. If receivers are installed in vertical wells, decomposition can be implemented by FK-filtering or taking vertical derivatives along the borehole [46]. Note that similar decomposition schemes can be applied for salt flank imaging in vertical wells (Figures 1(c) and 1(d)), where up- and downgoing have to be replaced by left and right going.

We use the same configuration as in Figure 2. We define and as the downgoing and upgoing fields at receivers and , respectively. Instead of evaluating the correlation function of full fields (see (8)), we define as the correlation function between and ; that is where subscript denotes “decomposed.” As a consequence of this choice, we will reconstruct a Green’s function as if there was a downgoing field emitted at and an upgoing field registered at plus an acausal Green’s function as if there was an upgoing field emitted at and a downgoing field registered at [25]. We can find a representation for this case by rewriting (12) in terms of up- and downgoing constituents and removing all but those that are downgoing at and upgoing at ; that is where the decomposed ghosts are given by Here refers to the field that is outgoing at source and downgoing at receiver . refers to the field that is outgoing at source and upgoing at receiver .

Equation (26) is useful for identifying ghosts in interferometry by CC of decomposed fields. It shows that cross-correlation of up- and downgoing constituents (25) yields as a causal part the desired Green’s function plus two additional ghost terms and . Note that the acausal part can be used to generate virtual sources that radiate upwards, as demonstrated by van der Neut et al. [25].

To illustrate the consequence of decomposition, we demonstrate (25) on model B (Figure 3(b)). The reference response consists of a single reflection at and a weak multiple at ; see Figure 11(a). The ghosts and correlation function are shown in Figures 11(b) and 11(c), respectively. Note that the ghosts are close to zero, meaning that the reflection response is almost perfectly reconstructed, apart from two weak spurious events at and . Compared to interferometry by cross-correlation of perturbed fields (Figure 9(c)), the multiples at and have been effectively removed.

fig11
Figure 11: Cross-correlation of decomposed fields for model B: (a) reference response by placing an active source at the virtual source location ; (b) computed ghost terms ; (c) correlation function (black) and the sum of reference response + ghosts (dashed green).

Wavefield decomposition does still not offer a solution for free-surface-related artifacts. To illustrate this, we apply the method to model D (Figure 3(d)). In Figure 12(a) we show the reference response, containing only the desired reflection at . Ghosts and the retrieved correlation function are shown in Figures 12(b) and 12(c), respectively. Note that artifacts are almost similar to those of perturbation-based interferometry (Figure 10(c)).

fig12
Figure 12: Cross-correlation of decomposed fields for model D: (a) reference response by placing an active source at the virtual source location ; (b) computed ghost terms ; (c) correlation function (black) and the sum of reference response + ghosts (dashed green).

5. Multidimensional Deconvolution of Perturbed Fields

Another strategy to eliminate ghosts in interferometry is to replace cross-correlation by deconvolution or, more exact, multidimensional deconvolution (MDD) [22]. In MDD, a Green’s function is retrieved by solving an inverse problem. This inverse problem is generally derived for decomposed wavefields and requires the installation of multicomponent receivers or two parallel receiver arrays close to each other. Before addressing MDD of decomposed fields, we present an alternative approach for perturbed fields, which can be applied with single-component sensors and time-gating.

To derive MDD for perturbed fields, we define a new volume , bounded by free surface and receiver array , see Figure 13. We define a reference medium, which is similar to the physical medium within the integration volume, but homogeneous below . All heterogeneities in the physical medium outside are referred to as perturbations. Subscript indicates those fields that propagate in the reference medium. Subscript is used for scattered fields, where denotes the full field in the physical medium. Vasconcelos et al. [11] derived a convolution-based representation for the scattered field with a source at and a receiver at , where is inside volume and is outside this volume, that is On the right-hand side we find two integrals. involves integration over the receiver array : is the spatial gradient at receiver location and is the outward-pointing unit normal vector of . Integral is a similar integral over , which vanishes since all pressure recordings are zero at this interface. We rewrite the wavefield at in terms of ingoing components and outgoing components . Similarly, . Since the reference medium is homogeneous outside , is purely outgoing at (so and ). It can be shown that the convolutions between ingoing fields at the virtual source and outgoing fields cancel and that the remaining terms can be merged [22]. Consequently, (28) and (29) can be rewritten as where we have introduced Equation (31) represents the exact scaled dipole Green’s function that can be solved by inversion of (30). Since an exact solution of (30) is generally not feasible, we apply least-squares inversion. We show in Appendix A that finding a least-squares solution of this problem is equivalent to solving the so-called normal equation, which is obtained by correlating both sides of (30) with and summing over (where is at ). Hence On the left-hand side of (32) we have the correlation function of the incident field at and the scattered field at : Note that (33) is identical to a discrete scaled representation of (20), if the medium parameters are constant at and is replaced by . On the right-hand side of (32) we have the so-called point-spread-function (PSF) for perturbed fields, defined as In both expressions, is an optional source-dependent weighting function. If sufficient source and receiver locations are present, a multidimensional inverse of the PSF can be introduced, according to where is at . The desired dipole response can now be found by filtering the correlation function with , according to Implementation of (36) is referred to as MDD of perturbed fields. This method allows us to deghost and deblur the correlation function of perturbed fields. We retrieve a Green’s function of an outward-radiating virtual source. However, this retrieved Green’s function lives in the physical medium and multiples from the overburden can still populate the retrieved gathers.

870819.fig.0013
Figure 13: Configuration for controlled-source interferometry by multidimensional deconvolution of perturbed fields; integration volume is enclosed by the free surface and receiver array ; denotes a source location; and denote receiver locations at and is the unit normal vector of the receiver array; is the unit normal vector of the free surface; heterogeneities between sources and receivers are referred to as “overburden,” whereas heterogeneities below the receivers are referred to as “target.”

We apply MDD of perturbed fields to model B (Figure 3(b)). The result is shown in Figure 14(a) (black) and compared with the reference response (dashed green). Note that the MDD response is similar to the correlation function (Figure 9(c)), apart from a very weak spurious event at , a scaling factor and subtle amplitude variations. We discussed that can be interpreted as the desired response , convolved with the PSF (see (32)). In Figure 14(b) we show that the PSF is close to a scaled band-limited delta function, with two additional weak events observed at . If we convolve the reference response with the PSF, we find indeed the correlation function, see Figure 14(c). From the simple structure of the PSF, it could have been directly concluded that and are indeed very similar apart from a scaling factor.

fig14
Figure 14: Deconvolution of perturbed fields for model B: (a) retrieved Green’s function (black) and the reference response (dashed green); (b) point spread function ; (c) convolution of the reflection response with the PSF (black) and correlation function (dashed green).

This is not the case for model D (Figure 3(d)). The response of MDD of perturbed fields (black) is compared with the reference response (dashed green) in Figure 15(a). Compared to the correlation function (Figure 10(c)), the MDD response is very simple, containing only the target reflection at and a multiple at . The strong disagreement of the MDD result and the correlation function is reflected by the PSF; see Figure 15(b). In Figure 15(c) we show that convolution of the PSF and reference response yields indeed the correlation function .

fig15
Figure 15: Deconvolution of perturbed fields for model D: (a) retrieved Green’s function (black) and the reference response (dashed green); (b) point spread function ; (c) convolution of the reflection response with the PSF (black) and the correlation function (dashed green).

6. Multidimensional Deconvolution of Decomposed Fields

The response retrieved by MDD of perturbed fields can still contain undesired reflections from the overburden. If ingoing and outgoing waves are separated at the receiver level prior to MDD, also these multiples can be removed [22]. For this purpose we redefine our volume once more, see Figure 16. Instead of enclosing the volume above the receivers, we now define as the volume enclosed by receiver array and a halfsphere below the receivers with radius . We define a reference state, in which all heterogeneities above are removed. Fields that propagate in this reference state are indicated with a bar; that is . We formulate a convolution based representation for the field of source at receiver , where is outside volume and is inside this volume, reading On the right-hand side we find two integrals. involves integration over the receiver array : Integral is a similar integral over the interface . Since this integral contains convolutions and its radius , it will vanish due to Sommerfeld’s radiation conditions [41]. Since the reference state is homogeneous outside , is purely ingoing at virtual source ( and ). We substitute . It can be shown that the convolutions of ingoing fields at the virtual source with outgoing fields cancel and that the remaining terms can be merged, yielding where we have introduced In a similar way as for MDD of perturbed fields, a normal equation can be derived: where is at . On the left-hand side we have the correlation function of the ingoing field at and the outgoing field at : Note that (42) is similar to a discrete scaled version of (25), if “ingoing” and “outgoing” are interchanged with “downgoing” and “upgoing,” medium parameters are constant at the source array and is replaced by . On the right-hand side we have the Point-Spread-Function (PSF) for decomposed fields, given by If acquisition conditions allow, we can take a multidimensional inverse of the PSF and convolve it with the correlation function according to where is at and Implementation of (44) is referred to as MDD of decomposed fields. This method allows us to deghost and deblur the correlation function of decomposed fields. Moreover, we retrieve a response that lives in the reference medium, where all multiples from the overburden have been removed.

870819.fig.0016
Figure 16: Configuration for controlled-source interferometry by multidimensional deconvolution of decomposed fields; integration volume is enclosed by receiver array and an additional halfsphere ; denotes a source location; and denote receiver locations at and is the unit normal vector of the receiver array; is the unit normal vector of ; heterogeneities between sources and receivers are referred to as “overburden,” whereas heterogeneities below the receivers are referred to as “target.”

We return to model B (Figure 3(b)). In Figure 17(a) we show the result of MDD of decomposed fields (black) and the reference response (dashed green). The response does not contain any of the multiples that have been retrieved by MDD of perturbed fields (Figure 14(a)). Note that indeed the weak ghosts and multiples that we observed in the correlation function (Figure 11(c)) have been eliminated, as predicted by theory. The weakness of these events is reflected in the PSF, showing a scaled band-limited delta function with weak events at and . We convolve the reference response with the PSF to show that indeed the correlation function emerges, see Figure 17(c).

fig17
Figure 17: Deconvolution of decomposed fields for model B: (a) retrieved Green’s function (black) and the reference response (dashed green); (b) point spread function ; (c) convolution of the reflection response with the PSF (black) and the correlation function (dashed green).

The effects of MDD of decomposed fields are more dramatically exposed by model D (Figure 3(d)). In Figure 18(a) we show that also for this model we can retrieve a response that is free of ghosts and multiples. Compared to MDD of perturbed fields (Figure 15(a)), we observe that the multiple at has been eliminated. The complex character of the PSF (Figure 18(b)) exposes the difference between the MDD response and the correlation function (Figure 12(c)). In Figure 18(c) we show that the correlation function can indeed be found by convolving the reference response with the PSF.

fig18
Figure 18: Deconvolution of decomposed fields for model D: (a) retrieved Green’s function (black) and the reference response (dashed green); (b) point spread function ; (c) convolution of the reflection response with the PSF (black) and the correlation function (dashed green).

7. Spatial Aspects

So far we have mostly concentrated on temporal aspects (ghosts and multiples) of interferometry. There are spatial aspects as well. In the representations (32) and (41) we have shown that the correlation functions of perturbed and decomposed fields can be interpreted as the desired reflection response, blurred in time and space with the PSF. This means that we can only retrieve an accurate response by cross-correlation if the PSF is focused at the virtual source location. However, due to unbalanced source distributions, intrinsic losses, or heterogeneities in the overburden, the PSF can be spatially defocused. As a result, the retrieved data by cross-correlation will be blurred. We illustrate this with a salt flank example.

In Figure 19 we show a salt flank model, defined as a function of coordinates (horizontal distance) and (depth). Note the velocity gradient in the medium, producing diving waves that are useful for salt flank imaging. We place 201 receivers in a vertical well with 20 m vertical spacing along the interval at . We place 401 sources at the surface with 30 m spacing along the interval at . No free surface is incorporated. In the following we generate a virtual shot record as if there was a source at receiver in the well using CC of perturbed fields. We generate a reference response by placing an active source at and modeling the wavefield. In Figure 20 we show three snapshots of the emitted wavefield that was used for the reference response. We indicate three reflections by numbers 1, 2, and 3. These are the reflections we aim to retrieve.

870819.fig.0019
Figure 19: Velocity model for the salt flank imaging example. Sources are located at boundary , covering at ; receivers are located at boundary , covering at ; a reference virtual source is located at .
fig20
Figure 20: Snapshots at 2.0 s, 2.7 s and 3.4 s of the wavefield emitted by a source at ; reflections 1, 2, and 3 are indicated.

We time-gate the incident fields of the observed data at the receiver array and subtract it from the full fields to extract the scattered fields. Next, we obtain the correlation function at virtual source by stacking cross-correlations of incident and scattered fields over all 401 sources (see (33)). A Hanning taper is applied to the first 20 and last 20 sources. In Figure 21 we show the retrieved response (red) and the reference response (black). For display purposes, only every 20th trace is shown. Note that the match is perfect and the reflections 1, 2, and 3 that were indicated in Figure 20 can easily be recognized.

870819.fig.0021
Figure 21: Reference response (black) and the correlation function (red) for a virtual source at ; reflections 1, 2, and 3 are indicated.

In seismic interferometry, the retrieval of a reflection response is often explained by summations over correlation gathers. To briefly illustrate this concept, we cross-correlated the incident field at with the scattered field at for different source locations , see Figure 22. For display purposes, only every 20th trace is shown. To retrieve the trace at 2000 m of the virtual shot record (Figure 21), a stack over all traces in Figure 22 is required. It can be shown that constructive interference takes place at the stationary points, being the maxima in time of each reflector in Figure 22, obeying Fermat’s principle. Destructive interference takes place outside these stationary points.

870819.fig.0022
Figure 22: Correlation gather for a virtual source at and a receiver at ; reflections 1, 2, and 3 are indicated.

Alternatively, we could interpret interferometry as the process of focusing a virtual source, meaning that the PSF in (34) is converging towards a spatial and temporal delta function, indicating that the correlation function does indeed represent the desired reflection response, free of blurring. In Figure 23 we show the PSF for a virtual source at . As for a virtual shot gather, the traces of the PSF can be interpreted as summations over correlation gathers. In the so-called PSF correlation gathers, the incident field at is correlated with the incident field at (“integrand” of (34)). If , this corresponds to auto-correlation, providing a significant contribution at , see Figure 24. These contributions will interfere constructively, to generate the desired spike in the PSF. If , cross-correlations of different source locations will map at different times; see Figure 25. If the virtual source is uniformly illuminated in each spatial direction, all such contributions will interfere destructively and the PSF will converge to the desired delta function. If this is the case, the virtual source will be perfectly focused, as in Figure 23. The two “legs” that can be observed in the PSF are caused by the finite source aperture. Finite source-aperture artifacts can be reduced by tapering the edges of the source array [39], as has already been done in this example. It was shown by van der Neut and Thorbecke [47] that the legs of the PSF can be used to diagnose the need for and the effects of such tapering.

870819.fig.0023
Figure 23: PSF for a virtual source at ; the red-blue plot is clipped at 20% of the maximum amplitude; the black and green traces represent true-amplitude temporal and spatial traces, respectively.
870819.fig.0024
Figure 24: PSF correlation gather for receivers at and .
870819.fig.0025
Figure 25: PSF correlation gather for receivers at and .

If illumination conditions vary with incidence angle, spatial blurring can occur. Illumination variations can be introduced at the source side, for instance if sources are nonuniformly distributed or source characteristics are spatially varying due to source coupling conditions. We mimic such a situation by assigning location- and frequency-dependent spectra to the sources, see Figure 26. In this particular example, the peak frequency and source strength vary randomly with the source location, superimposed by an additional trend in source strength along the array. Note that the low source numbers (corresponding to the left side of the array in Figure 19) are overilluminated with respect to the high source numbers (corresponding to the right side of the array in Figure 19). As a consequence of these variations, the PSF of the virtual source is no longer optimally focused, as shown in Figure 27. The events in the PSF that intersect the focus point ( and ) (apart from the legs), as seen in Figure 27, stem from incomplete destructive interference and are typical for spatial defocusing of the virtual source, leading to a distorted virtual source radiation pattern [10].

870819.fig.0026
Figure 26: Source spectra to mimic varying source-side conditions (dB scale).
870819.fig.0027
Figure 27: PSF for a virtual source at for a model with varying source spectra; the red-blue plot is clipped at 20% of the maximum amplitude; the black and green traces represent true-amplitude temporal and spatial traces, respectively.

Virtual source defocusing can also be a consequence of velocity or density variations in the overburden. We illustrate this by introducing a gas cloud in the model, see Figure 28. In Figure 29 we show the PSF in the medium with gas cloud, where the variations in source spectra of Figure 26 have also been incorporated. Note that besides spatial defocusing, we also observe temporal defocusing stemming from gas cloud scattering. Events in the PSF that do not intersect the focus point ( and ), as seen in Figure 29, are typical indicators of temporal defocusing, which can be related to the ghosts and multiples that we discussed earlier in this paper.

870819.fig.0028
Figure 28: Velocity model for the salt flank imaging example with gas cloud.
870819.fig.0029
Figure 29: PSF for a virtual source at for a model with varying source spectra and a gas cloud; the red-blue plot is clipped at 20% of the maximum amplitude; the black and green traces represent true-amplitude temporal and spatial traces, respectively.

In Figure 30 we show the correlation function for the model with varying source spectra and a gas cloud. Note that the correlation function is no longer perfectly matching the reference response, especially not in terms of amplitudes. This mismatch can be effectively removed from the retrieved data by MDD, meaning that the correlation function is filtered with the inverse of the PSF (see (36)). The retrieved gather after MDD is shown in Figure 31. Note that the match between retrieved and reference response has improved considerably. We do point out, however, that inversion artifacts can be created due to limited illumination conditions, as highlighted in the top of Figure 31.

870819.fig.0030
Figure 30: Reference response (black) and the correlation function (red) for a virtual source at in a model with varying source spectra and a gas cloud; reflections 1, 2, and 3 are indicated.
870819.fig.0031
Figure 31: Reference response (black) and the retrieved response by MDD (red) for a virtual source at in a model with varying source spectra and a gas cloud; reflections 1, 2, and 3 and an inversion artifact are indicated.

In Figure 32 we show an image of the salt flank in the model with source spectra variations and a gas cloud. The image is obtained by one-way shot profile migration of the correlation functions at all possible virtual source locations. In dashed red we show the location of the salt flank, as taken from the velocity model. Note that interferometry has enabled us to image the top of the flank “from below” by smart utilization of the diving waves. Also note that artifacts can be observed, which might be interpreted as additional ghost reflectors or diffractors. As MDD allows us to refocus the virtual sources before migration, several of these artifacts can be removed by migrating data after MDD, as shown in Figure 33. We observe that several potential “ghost reflectors” have been eliminated throughout the gather. The continuity of the salt flank amplitude marked by “A” is improved. The strong artifact marked by “B” as well as various other artifacts in this part of the gather have been suppressed. Artifacts such as indicated by “C” have been eliminated, but new artifacts such as “D” have emerged. We point with special attention to events “E” and “F” that might well be mistakenly interpreted as the continuation of the lower part of the salt flank. Note that MDD has completely eliminated these spurious arrivals. Finally, we note that inversion artifacts have hardened the interpretation of the lowest part of the salt flank in “G.” We conclude that MDD can indeed improve the image and remove defocusing effects, but care should be taken for potential inversion artifacts that can deteriorate parts of the image and mislead interpretation.

870819.fig.0032
Figure 32: Migrated section of the salt flank from the CC data in a model with varying source spectra and a gas cloud; the location of the salt flank as taken from the velocity model is given by the red dashed line; markers A–G are discussed in the main text.
870819.fig.0033
Figure 33: Migrated section of the salt flank from the MDD data in a model with varying source spectra and a gas cloud; the location of the salt flank as taken from the velocity model is given by the red dashed line; markers A–G are discussed in the main text.

8. Discussion

An overview has been given of the ghosts that appear in correlation interferometry of full, perturbed and decomposed fields. Equations (12), (24), and (26) describe these ghosts in an additive way. Alternatively, the correlation functions of perturbed or decomposed fields can be interpreted as the desired reflection responses, convolved in space and time with the PSF, see (32) and (41). Analysis of the PSF’s allows us to diagnose the quality of virtual source focusing in time and space. Along the temporal axis, the PSF gives information on possible ghosts and undesired multiples that may hamper the retrieved data. Along the spatial axis, the PSF gives information on focusing, that can be blurred due to unbalanced acquisition, intrinsic losses or complexities in the velocity model. The correlation function can be deblurred by filtering with the inverse of the PSF. This process is multidimensional deconvolution, allowing us to deghost, demultiple, and deblur the retrieved data.

The PSF can also provide insight in the effects of time-gating, which is often applied for separation of incident and scattered wavefields. In the Virtual Source method of Bakulin and Calvert [4], it is advocated to cross-correlate only the direct field instead of the full incident field. In Appendix B we use the PSF to show that narrow time-gates can improve the quality of virtual source focusing.

To compute the PSF for perturbed fields we require separation of the incident and scattered fields for each source. However, such time-gating can sometimes be problematic. Instead, an approximation of the PSF could also be obtained by time-gating the contributions around in the correlation function of full fields. This PSF can provide valuable insights in spatial virtual source focusing in various types of applications. Wapenaar et al. [22], for instance, showed how an estimate of the PSF could be obtained from cross-correlations of ambient seismic noise records. These PSFs could then be used to correct Green’s functions as retrieved by seismic interferometry for nonuniform passive source distributions.

Similarities exist between the derived methodology and model-driven redatuming [48]. Correlation-based interferometry can be related to correlation-based redatuming schemes such as those derived by Berryhill [49]. Multidimensional deconvolution of perturbed fields can be compared with rigorous redatuming [50]. Parallels can also be found with seismic migration and inversion [51]). The PSF that we defined for perturbed wavefields has close similarities with the resolution function in seismic inversion [5254]. Similarly, inversion of the PSF can be compared with migration deconvolution [55] or refocusing migrated images [56]. However, in all these cases it is important to realize that having actual subsurface receivers allows us to redatum wavefields, including multiple scattered events, much more effectively than with any model-driven method.

9. Conclusion

Controlled-source seismic interferometry is generally explained from cross-correlation based theory. Although this theory is exact, the required assumptions are often not fulfilled in practice. Because of one-sided illumination, complex subsurface structures, intrinsic losses, usage of single source types and free surface interactions, virtual sources can defocus and unphysical ghosts can enter the retrieved gathers. Even if all assumptions are fulfilled, particular undesired reflections from the overburden can still be retrieved by cross-correlation. Separation of incident and scattered fields or wavefield decomposition prior to cross-correlation can remove particular ghosts and multiples. Multidimensional deconvolution of perturbed (time-gated) fields allows us to refocus defocused virtual sources and remove additional ghosts and multiples. However, the method leaves particular multiples in place. To remove all multiples, multidimensional deconvolution should be applied to decomposed fields. It can be hard to stabilize the required inversion and artifacts can easily be introduced, especially if illumination conditions are limited. Through the interferometric Point Spread Function, we can diagnose illumination variations, ghosts and multiples. As this function can be obtained directly from the data, it can be a useful tool for analyzing virtual source focusing and, consequently, the quality of the retrieved data.

Appendices

A. Least-Squares Inversion

Since (30) does generally not have a unique solution, we aim to minimize a cost function instead. In least squares theory, this cost function is generally defined as [57] where is introduced as an additional weighting factor and is the misfit between the left- and right-hand side of (30), that is In least-squares inversion, our goal is to minimize the cost function at each receiver and frequency-component . However, as such inversion is generally unstable, we pose an additional constraint on minimizing the solution length : Instead of minimizing , we minimize , where determines the balance between minimizing the misfit and minimizing the solution length. Next we start to search for the solution , obeying After some algebra, (A.4) can be rewritten as [57] where is at . Quantities and are defined in (33) and (34). Inversion of (A.5) is equal to finding its least-squares inverse. By setting , (A.5) is similar to (32). This result is often referred to as the normal equation. For more details, see Menke [57].

B. Time Gating

To separate incident and scattered fields, we generally rely on time-gating. Incident fields generally contain not only primaries but also multiples from the overburden. In the virtual source method of Bakulin and Calvert [4] it is advocated to cross-correlate only the direct field (instead of the full incident field) at the virtual source location with the scattered fields at the other receivers. In this appendix we study the advantage of such strategy, using the point spread function.

First, let us introduce the virtual source correlation function of the direct field at receiver with the scattered field at receiver : where subscripts and stand for “virtual source method” and “direct field,” respectively. Obviously, the direct field does not contain the full incident field . A particular section is not captured by the time-gate. We may substitute into (30), to show that with Equation (B.2) may be solved by MDD. With similar reasoning as expressed in Appendix A we can show that this yields the following normal equation: where is the correlation function of the virtual source method (B.1) evaluated at receiver instead of . is the PSF of the virtual source method: is a ghost term associated with this strategy: With (B.4) we have shown that the response as retrieved by the virtual source method can be interpreted as the desired Green’s function blurred by plus an additional ghost . In (32) we derived that cross-correlation of the full incident field instead of the direct field yields the same Green’s function blurred by without a ghost term. By setting the time-gate, we have thus introduced an additional ghost. However, since does only contain the cross-correlations of the direct field, this function behaves generally much more like the desired delta function than . In other words: time-gating the direct field tends to focus the virtual source, which can also eliminate ghosts, multiples and blurring. A more detailed discussion on the aspects of time-gating is beyond the scope of this paper, but we refer to van der Neut and Bakulin [10] for an analysis in layered media.

Acknowledgments

This work was supported by the Dutch Technology Foundation STW, applied science division of NWO and the Technology Program of the Ministry of Economic Affairs. Part of this research was conducted at St Petersburg State University during a two-months stay of the first author. The authors thank Professor Boris Kashtan for hosting and for stimulating discussions at this university.

References

  1. A. Curtis, P. Gerstoft, H. Sato, R. Snieder, and K. Wapenaar, “Seismic interferometry—turning noise into signal,” Leading Edge, vol. 25, no. 9, pp. 1082–1092, 2006. View at Publisher · View at Google Scholar · View at Scopus
  2. K. Wapenaar, E. Slob, R. Snieder, and A. Curtis, “Tutorial on seismic interferometry—part 2: underlying theory and new advances,” Geophysics, vol. 75, no. 5, pp. 75A211–75A227, 2010.
  3. G. T. Schuster, Seismic Interferometry, Cambridge University Press, Cambridge, UK, 2009.
  4. A. Bakulin and R. Calvert, “The virtual source method: theory and case study,” Geophysics, vol. 71, no. 4, pp. SI139–SI150, 2006. View at Publisher · View at Google Scholar · View at Scopus
  5. A. Derode, E. Larose, M. Tanter et al., “Recovering the Green's function from field-field correlations in an open scattering medium (L),” Journal of the Acoustical Society of America, vol. 113, no. 6, pp. 2973–2976, 2003. View at Publisher · View at Google Scholar · View at Scopus
  6. K. Wapenaar, J. Fokkema, and R. Snieder, “Retrieving the Green's function in an open system by cross correlation: a comparison of approaches (L),” Journal of the Acoustical Society of America, vol. 118, no. 5, pp. 2783–2786, 2005. View at Publisher · View at Google Scholar · View at Scopus
  7. E. C. Slob and K. Wapenaar, “Electromagnetic Green's functions retrieval by cross-correlation and cross-convolution in media with losses,” Geophysical Research Letters, vol. 34, no. 5, Article ID L05307, 2007. View at Publisher · View at Google Scholar · View at Scopus
  8. F. Poletto and K. Wapenaar, “Virtual reflector representation theorem (acoustic medium),” Journal of the Acoustical Society of America, vol. 125, no. 4, pp. EL111–EL116, 2009. View at Publisher · View at Google Scholar · View at Scopus
  9. R. Snieder, K. Wapenaar, and K. Larner, “Spurious multiples in seismic interferometry of primaries,” Geophysics, vol. 71, no. 4, pp. SI111–SI124, 2006. View at Publisher · View at Google Scholar · View at Scopus
  10. J. van der Neut and A. Bakulin, “Estimating and correcting the amplitude radiation pattern of a virtual source,” Geophysics, vol. 74, no. 2, pp. SI27–SI36, 2009. View at Publisher · View at Google Scholar · View at Scopus
  11. I. Vasconcelos, R. Snieder, and H. Douma, “Representation theorems and Green's function retrieval for scattering in acoustic media,” Physical Review E, vol. 80, no. 3, Article ID 036605, 2009. View at Publisher · View at Google Scholar · View at Scopus
  12. A. Bakulin, A. Mateeva, K. Mehta et al., “Virtual source applications to imaging and reservoir monitoring,” Leading Edge, vol. 26, no. 6, pp. 732–740, 2007. View at Publisher · View at Google Scholar · View at Scopus
  13. K. Mehta, A. Bakulin, J. Sheiman, R. Calvert, and R. Snieder, “Improving the virtual source method by wavefield separation,” Geophysics, vol. 72, no. 4, pp. V79–V86, 2007. View at Publisher · View at Google Scholar · View at Scopus
  14. I. Vasconcelos, R. Snieder, and B. Hornby, “Imaging internal multiples from subsalt VSP data—examples of target-oriented interferometry,” Geophysics, vol. 73, no. 4, pp. S157–S168, 2008. View at Publisher · View at Google Scholar · View at Scopus
  15. R. Snieder, J. Sheiman, and R. Calvert, “Equivalence of the virtual-source method and wave-field deconvolution in seismic interferometry,” Physical Review E, vol. 73, no. 6, Article ID 066620, 2006. View at Publisher · View at Google Scholar · View at Scopus
  16. D. C. Riley and J. F. Claerbout, “2-D multiple reflections,” Geophysics, vol. 41, pp. 592–620, 1976.
  17. D. Loewenthal and E. A. Robinson, “On unified dual fields and Einstein deconvolution,” Geophysics, vol. 65, no. 1, pp. 293–303, 2000. View at Scopus
  18. F. Poletto, M. Malusa, F. Miranda, and U. Tinivella, “Seismic-while-drilling by using dual sensors in drill strings,” Geophysics, vol. 69, no. 5, pp. 1261–1271, 2004. View at Scopus
  19. F. Poletto and L. Petronio, “Seismic interferometry with a TBM source of transmitted and reflected waves,” Geophysics, vol. 71, no. 4, pp. SI85–SI93, 2006. View at Publisher · View at Google Scholar · View at Scopus
  20. I. Vasconcelos and R. Snieder, “Interferometry by deconvolution—part 2: theory for elastic waves and application to drill-bit seismic imaging,” Geophysics, vol. 73, no. 3, pp. S129–S141, 2008. View at Publisher · View at Google Scholar · View at Scopus
  21. I. Vasconcelos and R. Snieder, “Interferometry by deconvolution—part 1: theory for acoustic waves and numerical examples,” Geophysics, vol. 73, no. 3, pp. S115–S128, 2008. View at Publisher · View at Google Scholar · View at Scopus
  22. K. Wapenaar, J. van der Neut, E. Ruigrok, et al., “Seismic interferometry by crosscorrelation and by multi-dimensional deconvolution: a systematic comparison,” Geophysical Journal International, vol. 185, pp. 1335–1364, 2001.
  23. L. Amundsen, L. T. Ikelle, and L. E. Berg, “Multidimensional signature deconvolution and free-surface multiple elimination of marine multicomponent ocean-bottom seismic data,” Geophysics, vol. 66, no. 5, pp. 1594–1604, 2001. View at Scopus
  24. E. Holvik and L. Amundsen, “Elimination of the overburden response from multicomponent source and receiver seismic data, with source designature and decomposition into PP-, PS-, and SS-wave responses,” Geophysics, vol. 70, no. 2, pp. S43–S59, 2005. View at Publisher · View at Google Scholar · View at Scopus
  25. J. van der Neut, K. Mehta, J. Thorbecke, E. Slob, and K. Wapenaar, “Controlled-source interferometric redatuming by cross-correlation and multi-dimensional deconvolution in elastic media,” Geophysics, vol. 76, no. 4, pp. SA63–SA76, 2011.
  26. S. Minato, T. Matsuoka, T. Tsuji, D. Draganov, J. Hunziker, and K. Wapenaar, “Seismic interferometry using multidimensional deconvolution and crosscorrelation for crosswell seismic reflection data without borehole sources,” Geophysics, vol. 76, no. 1, pp. SA19–SA34, 2011.
  27. E. Slob, “Interferometry by deconvolution of multicomponent multioffset GPR data,” IEEE Transactions on Geoscience and Remote Sensing, vol. 47, no. 3, Article ID 4773210, pp. 828–838, 2009. View at Publisher · View at Google Scholar · View at Scopus
  28. J. Hunziker, E. C. Slob, and K. Wapenaar, “Controlled source electromagnetic interferometry by multidimensional deconvolution—spatial sampling aspects,” in Proceedings of the 71st Annual International Conference & Technical Exhibition, vol. P074, EAGE, Extended Abstracts, 2009.
  29. Y. Fan, R. Snieder, and J. Singer, “3-D controlled source electromagnetic (CSEM) interferometry by multi-dimensional deconvolution,” in Proceedings of the 79th Annual International Meeting, vol. 28, pp. 779–783, SEG Expanded Abstracts, 2009.
  30. E. Ruigrok, J. van der Neut, H. Djikpesse, C. Chen, and K. Wapenaar, “A feasibility study for the application of seismic interferometry by multi-dimensional deconvolution for lithosperic-scale imaging,” EGU General Assembly, Geophysical Research Abstracts, vol. 12, p. EGU2010-9370, 2010.
  31. X. Xiao, M. Zhou, and G. T. Schuster, “Salt-flank delineation by interferometric imaging of transmitted P- to S-waves,” Geophysics, vol. 71, no. 4, pp. SI197–SI207, 2006. View at Publisher · View at Google Scholar · View at Scopus
  32. B. E. Hornby and J. Yu, “Interferometric imaging of a salt flank using walkaway VSP data,” Leading Edge, vol. 26, no. 6, pp. 760–763, 2007. View at Publisher · View at Google Scholar · View at Scopus
  33. R. Lu, M. Willis, X. Campman, J. Ajo-Franklin, and M. N. Toksöz, “Redatuming through a salt canopy and target-oriented salt-flank imaging,” Geophysics, vol. 73, no. 3, pp. S63–S71, 2008. View at Publisher · View at Google Scholar · View at Scopus
  34. J. Ferrandis, A. Mateeva, P. Jorgensen, J. Lopez, and H. Dijkerman, “Application of virtual-source technology to the Zuidwending gas storage project,” The Leading Edge, vol. 28, pp. 296–301, 2009.
  35. F. Poletto, L. Petronio, and F. Miranda, “Seismic interferometry in a crosshole-RVSP experiment,” in 69th European Association of Geoscientists and Engineers Conference and Exhibition 2007 - "Securing The Future", pp. 169–173, gbr, June 2007. View at Scopus
  36. K. Mehta, D. Kiyashchenko, P. Jorgensen, J. Lopez, J. Ferrandis, and M. Costello, “Virtual source method applied to crosswell and horizontal well geometries,” Leading Edge, vol. 29, no. 6, pp. 712–723, 2010. View at Publisher · View at Google Scholar · View at Scopus
  37. J. Byun, J. Yu, and S. J. Seol, “Crosswell monitoring using virtual sources and horizontal wells,” Geophysics, vol. 75, no. 3, pp. SA37–SA43, 2010.
  38. K. Wapenaar and J. Fokkema, “Green's function representations for seismic interferometry,” Geophysics, vol. 71, no. 4, pp. SI33–SI46, 2006. View at Publisher · View at Google Scholar · View at Scopus
  39. K. Mehta, R. Snieder, R. Calvert, and J. Sheiman, “Acquisition geometry requirements for generating virtual-source data,” Leading Edge, vol. 27, no. 5, pp. 620–629, 2008. View at Publisher · View at Google Scholar · View at Scopus
  40. D. Draganov, K. Wapenaar, and J. Thorbecke, “Passive seismic imaging in the presence of white noise sources,” Leading Edge, vol. 23, no. 9, pp. 889–892, 2004. View at Scopus
  41. J. T. Fokkema and P. M. van den Berg, Seismic Applications of Acoustic Reciprocity, Elsevier, Amsterdam, The Netherlands, 1993.
  42. K. Wapenaar, “Green's function retrieval by cross-correlation in case of one-sided illumination,” Geophysical Research Letters, vol. 33, no. 19, Article ID L19304, 2006. View at Publisher · View at Google Scholar · View at Scopus
  43. C. P. A. Wapenaar, P. Herrmann, D. J. Verschuur, and A. J. Berkhout, “Decomposition of multicomponent seismic data into primary P- and S- wave responses,” Geophysical Prospecting, vol. 38, no. 6, pp. 633–661, 1990. View at Scopus
  44. L. Amundsen and A. Reitan, “Decomposition of multicomponent sea-floor data into upgoing and downgoing P- and S-waves,” Geophysics, vol. 60, no. 2, pp. 563–572, 1995. View at Scopus
  45. F. J. Barr, “Dual-sensor OBC technology,” Leading Edge, vol. 16, no. 1, pp. 45–51, 1997. View at Scopus
  46. L. Petronio and F. Poletto, “Dual signals separation in shallow borehole hydrophone data,” in Proceedings of the 72nd EAGE Conference & Exhibition, vol. F033, Extended Abstracts, 2010.
  47. J. van der Neut and J. Thorbecke, “Resolution function for controlled-source seismic interferometry: a data-driven diagnosis,” in Proceedings of the 79th International Meeting Society of Exploration Geophysicists, vol. 28, pp. 4090–4094, Expanded Abstracts, 2009.
  48. G. T. Schuster and M. Zhou, “A theoretical overview of model-based and correlation-based redatuming methods,” Geophysics, vol. 71, no. 4, pp. SI103–SI110, 2006. View at Publisher · View at Google Scholar · View at Scopus
  49. J. R. Berryhill, “Wave-equation datuming,” Geophysics, vol. 44, no. 8, pp. 1329–1344, 1979. View at Scopus
  50. W. A. Mulder, “Rigorous redatuming,” Geophysical Journal International, vol. 161, no. 2, pp. 401–415, 2005. View at Publisher · View at Google Scholar · View at Scopus
  51. J. Thorbecke and K. Wapenaar, “On the relation between seismic interferometry and the migration resolution function,” Geophysics, vol. 72, no. 6, pp. T61–T66, 2007. View at Publisher · View at Google Scholar · View at Scopus
  52. G. T. Schuster and J. Hu, “Green's function for migration: continuous recording geometry,” Geophysics, vol. 65, no. 1, pp. 167–175, 2000. View at Scopus
  53. G. Toxopeus, J. Thorbecke, K. Wapenaar, S. Petersen, E. Slob, and J. Fokkema, “Simulating migrated and inverted seismic data by filtering a geologic model,” Geophysics, vol. 73, no. 2, pp. T1–T10, 2008. View at Publisher · View at Google Scholar · View at Scopus
  54. E. J. van Veldhuizen, G. Blacquière, and A. J. Berkhout, “Acquisition geometry analysis in complex 3D media,” Geophysics, vol. 73, no. 5, pp. Q43–Q58, 2008. View at Publisher · View at Google Scholar · View at Scopus
  55. J. Yu, J. Hu, G. T. Schuster, and R. Estill, “Prestack migration deconvolution,” Geophysics, vol. 71, no. 2, pp. S53–S62, 2006. View at Publisher · View at Google Scholar · View at Scopus
  56. C. Zhang and T. J. Ulrych, “Refocusing migrated seismic images in absorptive media,” Geophysics, vol. 75, no. 3, pp. S103–S110, 2010.
  57. W. Menke, Geophysical Data Analysis, Academic Press, San Diego, Calif, USA, 1989.