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

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.

(a)

(b)

(c)

(d)

(e)

(f)

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

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]: 𝑆(𝜔)²𝐺𝐱_𝐵,𝐱_𝐴,𝜔+𝐺^∗𝐱_𝐴,𝐱_𝐵=𝐼̂,𝜔_𝑆𝐱_𝐵,𝐱_𝐴,𝜔+𝐼^𝑂̂𝐱_𝐵,𝐱_𝐴,𝜔. (1) 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]:𝐼^𝑆̂𝐱_𝐵,𝐱_𝐴,𝜔=−_𝜕𝑆1𝑗𝜔𝜌𝐱_𝑆×𝐧_𝑆⋅∇_𝑆𝑝𝐱_𝐵,𝐱_𝑆,𝜔𝑝𝐱_𝐴,𝐱_𝑆,𝜔^∗−𝑝𝐱_𝐵,𝐱_𝑆,𝜔𝐧_𝑆⋅∇_𝑆𝑝𝐱_𝐴,𝐱_𝑆,𝜔^∗𝑑𝐱_𝑆. (2) 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, ∇_𝑆=(𝜕/𝜕𝑥_1,𝑆,𝜕/𝜕𝑥_2,𝑆,𝜕/𝜕𝑥_3,𝑆)^𝑇, where 𝐱_𝑆=(𝑥_1,𝑆,𝑥_2,𝑆,𝑥_3,𝑆)^𝑇 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 𝐼^𝑂̂𝐱_𝐵,𝐱_𝐴,𝜔=−_𝜕𝑂1𝑗𝜔𝜌𝐱_𝑂×𝐧_𝑂⋅∇_𝑂𝑝𝐱_𝐵,𝐱_𝑂,𝜔𝑝𝐱_𝐴,𝐱_𝑂,𝜔^∗−𝑝𝐱_𝐵,𝐱_𝑂,𝜔𝐧_𝑂⋅∇_𝑂𝑝𝐱_𝐴,𝐱_𝑂,𝜔^∗𝑑𝐱_𝑂, (3) 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 𝑝^out(𝐱_𝐴,𝐱_𝑆,𝜔), where superscripts in and out denote ingoing and outgoing fields at 𝐱_𝑆. We substitute 𝑝=𝑝^out+𝑝ⁱⁿ 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 𝐼^𝑆̂𝐱_𝐵,𝐱_𝐴=,𝜔_𝜕𝑆2𝑗𝜔𝜌𝐱_𝑆×𝑝ⁱⁿ𝐱_𝐵,𝐱_𝑆,𝜔𝐧_𝑠⋅∇_𝑆𝑝ⁱⁿ𝐱_𝐴,𝐱_𝑆,𝜔^∗+𝑝^out𝐱_𝐵,𝐱_𝑆,𝜔𝐧_𝑠⋅∇_𝑆𝑝^out𝐱_𝐴,𝐱_𝑆,𝜔^∗𝑑𝐱_𝑆. (4) 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 𝐧_𝑆⋅∇_𝑆𝑝ⁱⁿ≈−(𝑗𝜔/𝑐(𝐱_𝑆))cos(𝜙)𝑝ⁱⁿ, 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 𝐧_𝑆⋅∇_𝑆𝑝^out≈+(𝑗𝜔/𝑐(𝐱_𝑆))cos(𝜙)𝑝^out. We assume that 𝜙 is close to zero such that cos(𝜙)≈1. Substituting these approximations into (4) yields𝐼^𝑆̂𝐱_𝐵,𝐱_𝐴≈,𝜔_𝜕𝑆2𝜌𝐱_𝑆𝑐𝐱_𝑆𝑝ⁱⁿ𝐱_𝐵,𝐱_𝑆,𝜔𝑝ⁱⁿ𝐱_𝐴,𝐱_𝑆,𝜔^∗−𝑝^out𝐱_𝐵,𝐱_𝑆,𝜔𝑝^out𝐱_𝐴,𝐱_𝑆,𝜔^∗𝑑𝐱_𝑆. (5) 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 𝑝^out=0). 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𝐼^𝑂̂𝐱_𝐵,𝐱_𝐴=,𝜔_𝜕𝑂2𝑗𝜔𝜌𝐱_𝑂𝑝𝐱_𝐵,𝐱_𝑂,𝜔𝐧_𝑂⋅∇_𝑂𝑝𝐱_𝐴,𝐱_𝑂,𝜔^∗𝑑𝐱_𝑂. (6) We can further simplify this equation by substituting the far-field approximation for ingoing fields 𝐧_𝑂⋅∇_𝑂𝑝≈−(𝑗𝜔/𝑐(𝐱_𝑂))𝑝, yielding𝐼^𝑂̂𝐱_𝐵,𝐱_𝐴≈,𝜔_𝜕𝑂2𝜌𝐱_𝑂𝑐𝐱_𝑂𝑝𝐱_𝐵,𝐱_𝑂,𝜔𝑝𝐱_𝐴,𝐱_𝑂,𝜔^∗𝑑𝐱_𝑂. (7) 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 𝐶:𝐶𝐱_𝐵,𝐱_𝐴=,𝜔_𝜕𝑆2𝑊𝐱_𝑆𝜌𝐱_𝑆𝑐𝐱_𝑆𝑝𝐱_𝐵,𝐱_𝑆,𝜔𝑝𝐱_𝐴,𝐱_𝑆,𝜔^∗𝑑𝐱_𝑆. (8) 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 cos𝜙-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 𝑊(𝐱_𝐴,𝐱_𝑆)=cos𝜙(𝐱_𝐴,𝐱_𝑆) 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 𝑝=𝑝^out+𝑝ⁱⁿ into (8) and rewrite the result as𝐶𝐱_𝐵,𝐱_𝐴=,𝜔_𝜕𝑆2𝑊𝐱_𝑆𝜌𝐱_𝑆𝑐𝐱_𝑆𝑝ⁱⁿ𝐱_𝐵,𝐱_𝑆,𝜔𝑝ⁱⁿ𝐱_𝐴,𝐱_𝑆,𝜔^∗−𝑝^out𝐱_𝐵,𝐱_𝑆,𝜔𝑝^out𝐱_𝐴,𝐱_𝑆,𝜔^∗𝑑𝐱_𝑆+_𝜕𝑆2𝑊𝐱_𝑆𝜌𝐱_𝑆𝑐𝐱_𝑆𝑝𝐱_𝐵,𝐱_𝑆,𝜔𝑝^out𝐱_𝐴,𝐱_𝑆,𝜔^∗+𝑝^out𝐱_𝐵,𝐱_𝑆,𝜔𝑝𝐱_𝐴,𝐱_𝑆,𝜔^∗𝑑𝐱_𝑆. (9) 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 𝑔^𝑆:𝑔^𝑆𝐱_𝐵,𝐱_𝐴,𝜔=_𝜕𝑆2𝑊𝐱_𝑆𝜌𝐱_𝑆𝑐𝐱_𝑆𝑝𝐱_𝐵,𝐱_𝑆,𝜔𝑝^out𝐱_𝐴,𝐱_𝑆,𝜔^∗+𝑝^out𝐱_𝐵,𝐱_𝑆,𝜔𝑝𝐱_𝐴,𝐱_𝑆,𝜔^∗𝑑𝐱_𝑆. (10) 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𝑔^𝑂𝐱_𝐵,𝐱_𝐴,𝜔=−_𝜕𝑂2𝜌𝐱_𝑂𝑐𝐱_𝑂𝑝𝐱_𝐵,𝐱_𝑂,𝜔𝑝𝐱_𝐴,𝐱_𝑂,𝜔^∗𝑑𝐱_𝑂. (11) Next, we substitute (11), (7), (5), (1) and (10) into (9) and rewrite the result as𝐶𝐱_𝐵,𝐱_𝐴,𝜔≈𝑆(𝜔)²𝐺𝐱_𝐵,𝐱_𝐴,𝜔+𝐺^∗𝐱_𝐴,𝐱_𝐵+𝑔,𝜔_𝑆𝐱_𝐵,𝐱_𝐴,𝜔+𝑔^𝑂𝐱_𝐵,𝐱_𝐴,𝜔. (12) 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 𝑝^out≠0. 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.

(a)

(b)

(c)

(d)

Figure 3 

Impedance (in kg⋅m−2⋅s−1) 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 𝑡=0s and an event at 𝑡=0.4s. 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 𝑔^𝑆=0 (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 𝑡=0s in this zero-offset case. In Figure 4(b) we show the ghosts 𝑔^𝑆+𝑔^𝑂. Note that indeed no other contribution is found outside 𝑡=0s. Therefore the correlation function 𝐶 (8) is very similar to the true reflection response, except around 𝑡=0s. In Figure 4(c) we show 𝐶 (black) and the sum of reference response and ghosts (dashed green), matching well.

(a)

(b)

(c)

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 𝑡=0s. 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 𝑡=0.4s but also the reflections of the overburden at 𝑡=0.1s and 𝑡=0.3s 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, 𝑔^𝑆=0. However, the second ghost 𝑔^𝑂 does give a significant contribution. The events at 𝑡=0.1s and 𝑡=0.3s 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 𝑡=0.2s (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 𝑡=0.5s and 𝑡=0.7s 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.

(a)

(b)

(c)

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 (𝑡=0.4s) and overburden (𝑡=0.1s and 𝑡=0.3s), including their multiples (𝑡=0.5s and 𝑡=0.7s), 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 32s. Moreover, it should be mentioned that the antiradiation condition can not eliminate the effects of reflectors placed outside the volume 𝕍, as 𝑔^𝑆≠0 for such scenario.

(a)

(b)

(c)

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 𝑡=0.4s but also the primaries of the overburden contrasts at 𝑡=0.6s and 𝑡=0.7s, the free surface reflection at 𝑡=0.8s; and free-surface multiples at 𝑡=0.9s and 𝑡=1.0s; see Figure 7(a). Ghosts 𝑔^𝑆+𝑔^𝑂 and correlation function 𝐶 are shown in Figures 7(b) and 7(c), respectively. We see strong spurious events at 𝑡=0.1s and 𝑡=0.2s in the retrieved response (and a weaker event at 𝑡=0.3s), stemming from interactions of the first reflectors with the free surface. We also observe that the reflections at 𝑡=0.6s and 𝑡=0.7s, their multiples, and the free surface reflection at 𝑡=0.8s are hardly retrieved due to the missing source at 𝐱_𝑂.

(a)

(b)

(c)

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 𝑝_inc. Fields that propagate in the physical medium can now be expressed as a superposition of the incident field and the so-called scattered field 𝑝_sc, that is 𝑝=𝑝_inc+𝑝_sc. 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 𝐱_𝐴:𝑆(𝜔)²𝐺_sc𝐱_𝐵,𝐱_𝐴,𝜔=𝐼̂^𝑆_ptb𝐱_𝐵,𝐱_𝐴,𝜔+𝐼̂^𝐸_ptb𝐱_𝐵,𝐱_𝐴,𝜔, (13) where subscript ptb 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 𝐼̂^𝑆_ptb, stemming from the cross-correlations of incident and scattered fields from the actual sources at 𝜕𝕊:𝐼̂^𝑆_ptb𝐱_𝐵,𝐱_𝐴,𝜔=−_𝜕𝑆1𝑗𝜔𝜌𝐱_𝑆×𝐧_𝑆⋅∇_𝑆𝑝_sc𝐱_𝐵,𝐱_𝑆,𝜔𝑝_inc𝐱_𝐴,𝐱_𝑆,𝜔^∗−𝑝_sc𝐱_𝐵,𝐱_𝑆,𝜔𝐧_𝑆⋅∇_𝑆𝑝_inc𝐱_𝐴,𝐱_𝑆,𝜔^∗𝑑𝐱_𝑆. (14) The second integral, 𝐼̂^𝐸_ptb, stems from similar cross-correlations of nonexisting sources 𝐱_𝐸 at 𝜕𝔼:𝐼̂^𝐸_ptb𝐱_𝐵,𝐱_𝐴,𝜔=−_𝜕𝐸1𝑗𝜔𝜌𝐱_𝐸×𝐧_𝐸⋅∇_𝐸𝑝_sc𝐱_𝐵,𝐱_𝐸,𝜔𝑝_inc𝐱_𝐴,𝐱_𝐸,𝜔^∗−𝑝_sc𝐱_𝐵,𝐱_𝐸,𝜔𝐧_𝐸⋅∇_𝐸𝑝_inc𝐱_𝐴,𝐱_𝐸,𝜔^∗𝑑𝐱_𝐸, (15) 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𝐼̂^𝑆_ptb𝐱_𝐵,𝐱_𝐴,𝜔=_𝜕𝑆2𝑗𝜔𝜌𝐱_𝑆×𝑝ⁱⁿ_sc𝐱_𝐵,𝐱_𝑆,𝜔𝐧_𝑆⋅∇_𝑆𝑝ⁱⁿ_inc𝐱_𝐴,𝐱_𝑆,𝜔^∗+𝑝^out_sc𝐱_𝐵,𝐱_𝑆,𝜔𝐧_𝑆⋅∇_𝑆𝑝^out_inc𝐱_𝐴,𝐱_𝑆,𝜔^∗𝑑𝐱_𝑆. (16) Assuming that the relevant wavefields propagate near normal incidence, we can substitute far-field approximations 𝐧_𝑆⋅∇_𝑆𝑝^out_inc≈+(𝑗𝜔/𝑐(𝐱_𝑆))𝑝^out_inc and 𝐧_𝑆⋅∇_𝑆𝑝ⁱⁿ_inc≈−(𝑗𝜔/𝑐(𝐱_𝑆))𝑝ⁱⁿ_inc into (16), yielding𝐼̂^𝑆_ptb𝐱_𝐵,𝐱_𝐴,𝜔≈_𝜕𝑆2𝜌𝐱_𝑆𝑐𝐱_𝑆𝑝ⁱⁿ_sc𝐱_𝐵,𝐱_𝑆,𝜔𝑝ⁱⁿ_inc𝐱_𝐴,𝐱_𝑆,𝜔^∗−𝑝^out_sc𝐱_𝐵,𝐱_𝑆,𝜔𝑝^out_inc𝐱_𝐴,𝐱_𝑆,𝜔^∗𝑑𝐱_𝑆. (17) For integral 𝐼̂^𝐸_ptb 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 𝑝_inc=𝑝ⁱⁿ_inc and 𝑝^out_inc=0). Cross-correlations between ingoing and outgoing fields cancel and the remaining terms can be merged [11], such that𝐼̂^𝐸_ptb𝐱_𝐵,𝐱_𝐴=,𝜔_𝜕𝐸2𝑗𝜔𝜌𝐱_𝐸×𝑝ⁱⁿ_sc𝐱_𝐵,𝐱_𝐸,𝜔𝐧_𝐸⋅∇_𝐸𝑝_inc𝐱_𝐴,𝐱_𝐸,𝜔^∗𝑑𝐱_𝐸. (18) The spatial derivative can be approximated with far-field approximation 𝐧_𝐸⋅∇_𝐸𝑝_inc≈−(𝑗𝜔/𝑐(𝐱_𝐸))𝑝_inc, such that𝐼̂^𝐸_ptb𝐱_𝐵,𝐱_𝐴≈,𝜔_𝜕𝐸2𝜌𝐱_𝐸𝑐𝐱_𝐸𝑝ⁱⁿ_sc𝐱_𝐵,𝐱_𝐸,𝜔𝑝_inc𝐱_𝐴,𝐱_𝐸,𝜔^∗𝑑𝐱_𝐸. (19) 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 𝐶_ptb:𝐶_ptb𝐱_𝐵,𝐱_𝐴=,𝜔_𝜕𝑆2𝑊𝐱_𝑆𝜌𝐱_𝑆𝑐𝐱_𝑆𝑝_sc𝐱_𝐵,𝐱_𝑆,𝜔𝑝_inc𝐱_𝐴,𝐱_𝑆,𝜔^∗𝑑𝐱_𝑆. (20) To evaluate the consequence of this choice, we separate the scattered and incident wavefields in up- and downgoing constituents according to 𝑝_inc=𝑝ⁱⁿ_inc+𝑝^out_inc and 𝑝_sc=𝑝ⁱⁿ_sc+𝑝^out_sc and substitute these representations into (20). The result can be written as𝐶_ptb𝐱_𝐵,𝐱_𝐴=,𝜔_𝜕𝑆2𝑊𝐱_𝑆𝜌𝐱_𝑆𝑐𝐱_𝑆𝑝ⁱⁿ_sc𝐱_𝐵,𝐱_𝑆,𝜔𝑝ⁱⁿ_inc𝐱_𝐴,𝐱_𝑆,𝜔^∗−𝑝^out_sc𝐱_𝐵,𝐱_𝑆,𝜔𝑝^out_inc𝐱_𝐴,𝐱_𝑆,𝜔^∗𝑑𝐱_𝑆+_𝜕𝑆2𝑊𝐱_𝑆𝜌𝐱_𝑆𝑐𝐱_𝑆𝑝_sc𝐱_𝐵,𝐱_𝑆,𝜔𝑝^out_inc𝐱_𝐴,𝐱_𝑆,𝜔^∗+𝑝^out_sc𝐱_𝐵,𝐱_𝑆,𝜔𝑝_inc𝐱_𝐴,𝐱_𝑆,𝜔^∗𝑑𝐱_𝑆. (21) The first integral in (21) can be identified as the desired integral 𝐼̂^𝑆_ptb (see (17)). The second integral is a ghost term that we redefine as𝑔^𝑆_ptb𝐱_𝐵,𝐱_𝐴,𝜔=_𝜕𝑆2𝑊𝐱_𝑆𝜌𝐱_𝑆𝑐𝐱_𝑆𝑝_sc𝐱_𝐵,𝐱_𝑆,𝜔𝑝^out_inc𝐱_𝐴,𝐱_𝑆,𝜔^∗+𝑝^out_sc𝐱_𝐵,𝐱_𝑆,𝜔𝑝_inc𝐱_𝐴,𝐱_𝑆,𝜔^∗𝑑𝐱_𝑆. (22) The missing sources at 𝐱_𝐸 make up for a second ghost term that we define by rewriting the integral in (19) in the following way:𝑔^𝐸_ptb𝐱_𝐵,𝐱_𝐴,𝜔=−_𝜕𝐸2𝜌𝐱_𝐸𝑐𝐱_𝐸𝑝ⁱⁿ_sc𝐱_𝐵,𝐱_𝐸,𝜔𝑝_inc𝐱_𝐴,𝐱_𝐸,𝜔^∗𝑑𝐱_𝐸. (23) Now, by substituting (23), (19), (17), (13), and (20) into (21) we arrive at𝐶_ptb𝐱_𝐵,𝐱_𝐴,𝜔≈𝑆(𝜔)²𝐺_sc𝐱_𝐵,𝐱_𝐴+𝑔,𝜔_𝑆ptb𝐱_𝐵,𝐱_𝐴,𝜔+𝑔^𝐸_ptb𝐱_𝐵,𝐱_𝐴,𝜔. (24) Equation (24) is useful for identifying ghosts in interferometry by CC of perturbed fields. Computation of 𝐶_ptb (see (20)) yields a scaled bandlimited version of the desired scattered Green’s function plus two ghost terms 𝑔^𝑆_ptb and 𝑔^𝐸_ptb. The first ghost 𝑔^𝑆_ptb 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 𝑔^𝐸_ptb 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.

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 𝑡=0.4s but also the multiple reflections from the overburden at 𝑡=0.5s and 𝑡=0.7s as well as higher-order multiples. In Figures 9(b) and 9(c) we show the ghosts 𝑔^𝑆_ptb+𝑔^𝐸_ptb and the correlation function 𝐶_ptb. Note which the spurious event at 𝑡=0.2s, that was visible in Figure 5(c), is significantly weakened and hardly visible. The multiples at 𝑡=0.5s and 𝑡=0.7s 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 𝐱_𝐸.

(a)

(b)

(c)

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 𝑔𝑆ptb+𝑔𝐸ptb; (c) correlation function 𝐶ptb (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 𝑡=0.4s and a multiple at 𝑡=1.0s. Ghosts 𝑔^𝑆_ptb+𝑔^𝐸_ptb and 𝐶_ptb are shown in Figures 10(b) and 10(c), respectively. Note that we can still observe several spurious events and the undesired overburden reflections.

(a)

(b)

(c)

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 𝑔𝑆ptb+𝑔𝐸ptb; (c) correlation function 𝐶ptb (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 𝐶_dcp as the correlation function between 𝑝⁺ and 𝑝⁻; that is𝐶_dcp𝐱_𝐵,𝐱_𝐴=,𝜔_𝜕𝑆2𝑊𝐱_𝑆𝜌𝐱_𝑆𝑐𝐱_𝑆𝑝⁻𝐱_𝐵,𝐱_𝑆,𝜔𝑝⁺𝐱_𝐴,𝐱_𝑆,𝜔^∗𝑑𝐱_𝑆, (25) where subscript dcp 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𝐶_dcp𝐱_𝐵,𝐱_𝐴≈𝑆(𝜔,𝜔)²𝐺^−,+𝐱_𝐵,𝐱_𝐴,𝜔+𝐺^+,−𝐱_𝐴,𝐱_𝐵,𝜔^∗+𝑔^𝑆_dcp𝐱_𝐵,𝐱_𝐴,𝜔+𝑔^𝑂_dcp𝐱_𝐵,𝐱_𝐴,𝜔, (26) where the decomposed ghosts are given by𝑔^𝑆_dcp𝐱_𝐵,𝐱_𝐴,𝜔=_𝜕𝑆2𝑊𝐱_𝑆𝜌𝐱_𝑆𝑐𝐱_𝑆𝑝⁻𝐱_𝐵,𝐱_𝑆,𝜔𝑝^+,out𝐱_𝐴,𝐱_𝑆,𝜔^∗+𝑝^−,out𝐱_𝐵,𝐱_𝑆,𝜔𝑝⁺𝐱_𝐴,𝐱_𝑆,𝜔^∗𝑑𝐱_𝑆,𝑔^𝑂_dcp𝐱_𝐵,𝐱_𝐴,𝜔=−_𝜕𝑂2𝜌𝐱_𝑂𝑐𝐱_𝑂𝑝⁻𝐱_𝐵,𝐱_𝑂,𝜔𝑝⁺𝐱_𝐴,𝐱_𝑂,𝜔^∗𝑑𝐱_𝑂. (27) Here 𝑝^+,out(𝐱_𝐴,𝐱_𝑆,𝜔) refers to the field that is outgoing at source 𝐱_𝑆 and downgoing at receiver 𝐱_𝐴. 𝑝^−,out(𝐱_𝐵,𝐱_𝑆,𝜔) 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 𝑔^𝑆_dcp and 𝑔^𝑂_dcp. 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 𝑡=0.4s and a weak multiple at 𝑡=0.9s; see Figure 11(a). The ghosts 𝑔^𝑆_dcp+𝑔^𝑂_dcp and correlation function 𝐶_dcp 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 𝑡=0.2s and 𝑡=0.6s. Compared to interferometry by cross-correlation of perturbed fields (Figure 9(c)), the multiples at 𝑡=0.5s and 𝑡=0.7s have been effectively removed.

(a)

(b)

(c)

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 𝑔𝑆dcp+𝑔𝑂dcp; (c) correlation function 𝐶dcp (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 𝑡=0.4s. Ghosts 𝑔^𝑆_dcp+𝑔^𝑂_dcp and the retrieved correlation function 𝐶_dcp 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)).

(a)

(b)

(c)

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 𝑔𝑆dcp+𝑔𝑂dcp; (c) correlation function 𝐶dcp (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 inc indicates those fields that propagate in the reference medium. Subscript sc is used for scattered fields, where 𝑝=𝑝_inc+𝑝_sc 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𝑝_sc𝐱_𝐵,𝐱_𝑆,𝜔=𝐽̂^𝐴_ptb𝐱_𝐵,𝐱_𝑆,𝜔+𝐽̂^𝐹_ptb𝐱_𝐵,𝐱_𝑆,𝜔. (28) On the right-hand side we find two integrals. 𝐽̂^𝐴_ptb involves integration over the receiver array 𝜕𝔸:𝐽̂^𝐴_ptb𝐱_𝐵,𝐱_𝑆,𝜔=_𝜕𝐴1𝑗𝜔𝜌𝐱_𝐴×𝐧_𝐴⋅∇_𝐴𝐺_sc𝐱_𝐵,𝐱_𝐴,𝜔𝑝_inc𝐱_𝐴,𝐱_𝑆−𝐺,𝜔_sc𝐱_𝐵,𝐱_𝐴,𝜔𝐧_𝐴⋅∇_𝐴𝑝_inc𝐱_𝐴,𝐱_𝑆,𝜔𝑑𝐱_𝐴. (29)∇_𝐴 is the spatial gradient at receiver location 𝐱_𝐴 and 𝐧_𝐴 is the outward-pointing unit normal vector of 𝜕𝔸. Integral 𝐽̂^𝐹_ptb 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 𝑝^out. Similarly, 𝐺_sc=𝐺^out_sc+𝐺ⁱⁿ_sc. Since the reference medium is homogeneous outside 𝕍, 𝑝_inc is purely outgoing at 𝐱_𝐴 (so 𝑝_inc=𝑝^out_inc and 𝑝ⁱⁿ_inc=0). It can be shown that the convolutions between ingoing fields at the virtual source 𝐱_𝐴 and outgoing fields 𝑝_inc(𝐱_𝐴,𝐱_𝑆,𝜔) cancel and that the remaining terms can be merged [22]. Consequently, (28) and (29) can be rewritten as𝑝_sc𝐱_𝐵,𝐱_𝑆,𝜔=_𝜕𝐴𝐺_ptb𝐱_𝐵,𝐱_𝐴,𝜔𝑝_inc𝐱_𝐴,𝐱_𝑆,𝜔𝑑𝐱_𝐴, (30) where we have introduced𝐺_ptb𝐱_𝐵,𝐱_𝐴=2,𝜔𝑗𝜔𝜌𝐱_𝐴𝐧_𝐴⋅∇_𝐴𝐺^out_sc𝐱_𝐵,𝐱_𝐴,𝜔. (31) 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 𝐺_inc(𝐱^′_𝐴,𝐱_𝑆,𝜔) and summing over 𝐱_𝑆 (where 𝐱^′_𝐴 is at 𝜕𝔸). Hence𝐶_ptb𝐱_𝐵,𝐱^′_𝐴,𝜔=_𝜕𝐴𝐺_ptb𝐱_𝐵,𝐱_𝐴,𝜔Γ_ptb𝐱_𝐴,𝐱^′_𝐴,𝜔𝑑𝐱_𝐴. (32) On the left-hand side of (32) we have the correlation function of the incident field at 𝐱^′_𝐴 and the scattered field at 𝐱_𝐵:𝐶_ptb𝐱_𝐵,𝐱^′_𝐴,𝜔=_𝑖𝑊𝐱_𝑆𝑝_sc𝐱_𝐵,𝐱_𝑆^(𝑖),𝜔𝑝_inc𝐱^′_𝐴,𝐱_𝑆^(𝑖),𝜔^∗, (33) 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Γ_ptb𝐱_𝐴,𝐱^′_𝐴,𝜔=_𝑖𝑊𝐱_𝑆^(𝑖)𝑝_inc𝐱_𝐴,𝐱_𝑆^(𝑖),𝜔𝑝_inc𝐱^′_𝐴,𝐱_𝑆^(𝑖),𝜔^∗. (34) In both expressions, 𝑊(𝐱_𝑆) is an optional source-dependent weighting function. If sufficient source and receiver locations are present, a multidimensional inverse Γ^inv_ptb of the PSF can be introduced, according to _𝜕𝐴Γ_ptb𝐱_𝐴,𝐱^′_𝐴,𝜔Γ^inv_ptb𝐱^′_𝐴,𝐱_𝐴^′′,𝜔𝑑𝐱^′_𝐴=𝛿𝐱_𝐴−𝐱_𝐴^′′. (35) where 𝐱_𝐴^′′ is at 𝜕𝔸. The desired dipole response can now be found by filtering the correlation function with Γ^inv, according to𝐺_ptb𝐱_𝐵,𝐱_𝐴^′′,𝜔=_𝜕𝐴𝐶_ptb𝐱_𝐵,𝐱^′_𝐴,𝜔Γ^inv_ptb𝐱^′_𝐴,𝐱_𝐴^′′,𝜔𝑑𝐱^′_𝐴, (36) 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.

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 𝐺_ptb (dashed green). Note that the MDD response is similar to the correlation function 𝐶_ptb (Figure 9(c)), apart from a very weak spurious event at 𝑡=0.2s, a scaling factor and subtle amplitude variations. We discussed that 𝐶_ptb can be interpreted as the desired response 𝐺_ptb, 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 𝑡=±0.2s. 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 𝐺_ptb and 𝐶_ptb are indeed very similar apart from a scaling factor.

(a)

(b)

(c)

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 Γptb; (c) convolution of the reflection response with the PSF (black) and correlation function 𝐶ptb (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 𝐺_ptb (dashed green) in Figure 15(a). Compared to the correlation function 𝐶_ptb (Figure 10(c)), the MDD response is very simple, containing only the target reflection at 𝑡=0.4s and a multiple at 𝑡=1.0s. 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 𝐶_ptb.

(a)

(b)

(c)

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 Γptb; (c) convolution of the reflection response with the PSF (black) and the correlation function 𝐶ptb (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𝑝𝐱_𝐵,𝐱_𝑆,𝜔=𝐽̂^𝐴_dcp𝐱_𝐵,𝐱_𝑆,𝜔+𝐽̂^𝑂_dcp𝐱_𝐵,𝐱_𝑆,𝜔. (37) On the right-hand side we find two integrals. 𝐽̂^𝐴_dcp involves integration over the receiver array 𝜕𝔸:𝐽̂^𝐴_dcp𝐱_𝐵,𝐱_𝑆,𝜔=−_𝜕𝐴1𝑗𝜔𝜌𝐱_𝐴×𝐧_𝐴⋅∇_𝐴𝐺𝐱_𝐵,𝐱_𝐴,𝜔𝑝𝐱_𝐴,𝐱_𝑆−,𝜔𝐺𝐱_𝐵,𝐱_𝐴,𝜔𝐧_𝐴⋅∇_𝐴𝑝𝐱_𝐴,𝐱_𝑆,𝜔𝑑𝐱_𝐴. (38) Integral 𝐽̂^𝑂_dcp 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 𝐺^out=0). We substitute 𝑝=𝑝ⁱⁿ+𝑝^out. It can be shown that the convolutions of ingoing fields at the virtual source 𝐱_𝐴 with outgoing fields 𝑝^out(𝐱_𝐴,𝐱_𝑆,𝜔) cancel and that the remaining terms can be merged, yielding 𝑝𝐱_𝐵,𝐱_𝑆,𝜔=_𝜕𝐴𝐺_dcp𝐱_𝐵,𝐱_𝐴,𝜔𝑝ⁱⁿ𝐱_𝐴,𝐱_𝑆,𝜔𝑑𝐱_𝐴, (39) where we have introduced𝐺_dcp𝐱_𝐵,𝐱_𝐴2,𝜔=−𝑗𝜔𝜌𝐱_𝐴𝐧_𝐴⋅∇_𝐴𝐺𝐱_𝐵,𝐱_𝐴,𝜔. (40) In a similar way as for MDD of perturbed fields, a normal equation can be derived:𝐶_dcp𝐱_𝐵,𝐱^′_𝐴,𝜔=_𝜕𝐴𝐺_dcp𝐱_𝐵,𝐱_𝐴,𝜔Γ_dcp𝐱_𝐴,𝐱^′_𝐴,𝜔𝑑𝐱_𝐴, (41) where 𝐱^′_𝐴 is at 𝜕𝔸. On the left-hand side we have the correlation function of the ingoing field at 𝐱^′_𝐴 and the outgoing field at 𝐱_𝐵:𝐶_dcp𝐱_𝐵,𝐱^′_𝐴,𝜔=_𝑖𝑊𝐱_𝑆𝑝^out𝐱_𝐵,𝐱_𝑆^(𝑖),𝜔𝑝ⁱⁿ𝐱^′_𝐴,𝐱_𝑆^(𝑖),𝜔^∗. (42) 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Γ_dcp𝐱_𝐴,𝐱^′_𝐴,𝜔=_𝑖𝑊𝐱_𝑆𝑝ⁱⁿ𝐱_𝐴,𝐱_𝑆^(𝑖),𝜔𝑝ⁱⁿ𝐱^′_𝐴,𝐱_𝑆^(𝑖),𝜔^∗. (43) If acquisition conditions allow, we can take a multidimensional inverse of the PSF Γ^inv_dcp and convolve it with the correlation function according to𝐺_dcp𝐱_𝐵,𝐱_𝐴^′′,𝜔=_𝜕𝐴𝐶_dcp𝐱_𝐵,𝐱^′_𝐴,𝜔Γ^inv_dcp𝐱^′_𝐴,𝐱_𝐴^′′,𝜔𝑑𝐱^′_𝐴, (44) where 𝐱_𝐴^′′ is at 𝜕𝔸 and_𝜕𝐴Γ_dcp𝐱_𝐴,𝐱^′_𝐴,𝜔Γ^inv_dcp𝐱^′_𝐴,𝐱_𝐴^′′,𝜔𝑑𝐱^′_𝐴=𝛿𝐱_𝐴−𝐱_𝐴^′′. (45) 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.

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 𝑡=±0.2s and 𝑡=±0.5s. We convolve the reference response with the PSF to show that indeed the correlation function emerges, see Figure 17(c).

(a)

(b)

(c)

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 Γdcp; (c) convolution of the reflection response with the PSF (black) and the correlation function 𝐶dcp (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 𝑡=1.0s 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.

(a)

(b)

(c)

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 Γdcp; (c) convolution of the reflection response with the PSF (black) and the correlation function 𝐶dcp (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 𝑥₃∈[1000m,5000m] at 𝑥₁=15000m. We place 401 sources at the surface with 30 m spacing along the interval 𝑥₁∈[0m,12000m] at 𝑥₃=0m. No free surface is incorporated. In the following we generate a virtual shot record as if there was a source at receiver 𝐱_𝐴=(15000m,4000m) 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.

Figure 19 

Velocity model for the salt flank imaging example. Sources are located at boundary 𝜕𝕊, covering 𝑥1∈[0m,12000m] at 𝑥3=0m; receivers are located at boundary 𝜕𝔸, covering 𝑥3∈[1000m,5000m] at 𝑥1=15000m; a reference virtual source is located at 𝐱𝐴=(15000m,4000m).

(a)

(b)

(c)

Figure 20 

Snapshots at 2.0 s, 2.7 s and 3.4 s of the wavefield emitted by a source at 𝐱𝐴=(15000m,4000m); 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.

Figure 21 

Reference response (black) and the correlation function (red) for a virtual source at 𝐱𝐴=(15000m,4000m); 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 𝐱_𝐴=(15000m,4000m) with the scattered field at 𝐱_𝐵=(15000m,2000m) 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.

Figure 22 

Correlation gather for a virtual source at 𝐱𝐴=(15000m,4000m) and a receiver at 𝐱𝐵=(15000m,2000m); 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 𝐱_𝐴=(15000m,4000m). 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 𝑡=0, 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.

Figure 23 

PSF for a virtual source at 𝐱𝐴=(15000m,4000m); 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.

Figure 24 

PSF correlation gather for receivers at 𝐱𝐴=(15000m,4000m) and 𝐱′𝐴=𝐱𝐴=(15000m,4000m).

Figure 25 

PSF correlation gather for receivers at 𝐱𝐴=(15000m,4000m) and 𝐱′𝐴=(15000m,2000m).

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 𝑡=0) (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].

Figure 26 

Source spectra to mimic varying source-side conditions (dB scale).

Figure 27 

PSF for a virtual source at 𝐱𝐴=(15000m,4000m) 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 𝑡=0), 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.

Figure 28 

Velocity model for the salt flank imaging example with gas cloud.

Figure 29 

PSF for a virtual source at 𝐱𝐴=(15000m,4000m) 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.

Figure 30 

Reference response (black) and the correlation function (red) for a virtual source at 𝐱𝐴=(15000m,4000m) in a model with varying source spectra and a gas cloud; reflections 1, 2, and 3 are indicated.

Figure 31 

Reference response (black) and the retrieved response by MDD (red) for a virtual source at 𝐱𝐴=(15000m,4000m) 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.

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.

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 𝑡=0 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 [52–54]. 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]𝐸𝐱_𝐵,𝜔=_𝑖𝑊𝐱_𝑆^(𝑖)𝑒̂𝐱_𝐵,𝐱_𝑆^(𝑖),𝜔𝑒^∗̂𝐱_𝐵,𝐱_𝑆^(𝑖),𝜔, (A.1) where 𝑊 is introduced as an additional weighting factor and 𝑒̂ is the misfit between the left- and right-hand side of (30), that is𝑒̂𝐱_𝐵,𝐱_𝑆,𝜔=𝑝_sc𝐱_𝐵,𝐱_𝑆−,𝜔_𝜕𝐴𝐺_ptb𝐱_𝐵,𝐱_𝐴,𝜔𝑝_inc𝐱_𝐴,𝐱_𝑆,𝜔𝑑𝐱_𝐴. (A.2) 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 𝐿:𝐿𝐱_𝐵,𝜔=_𝜕𝐴𝐺_ptb𝐱_𝐵,𝐱_𝐴,𝜔𝐺_ptb𝐱_𝐵,𝐱_𝐴,𝜔^∗𝑑𝐱_𝐴. (A.3) 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 𝐺_ptb, obeying𝜕𝐸𝐱_𝐵,𝜔+𝜀²𝐿𝐱_𝐵,𝜔𝜕𝐺_ptb𝐱_𝐵,𝐱_𝐴,𝜔=0. (A.4) After some algebra, (A.4) can be rewritten as [57]𝐶_ptb𝐱_𝐵,𝐱^′_𝐴=,𝜔_𝜕𝐴𝐺_ptb𝐱_𝐴,𝐱_𝐵,𝜔Γ_ptb𝐱_𝐴,𝐱^′_𝐴,𝜔+𝜀²𝛿𝐱_𝐴−𝐱^′_𝐴𝑑𝐱_𝐴, (A.5) where 𝐱^′_𝐴 is at 𝜕𝔸. Quantities 𝐶_ptb and Γ_ptb are defined in (33) and (34). Inversion of (A.5) is equal to finding its least-squares inverse. By setting 𝜀=0, (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 𝐶_vsm of the direct field at receiver 𝐱_𝐴 with the scattered field at receiver 𝐱_𝐵: 𝐶_vsm𝐱_𝐵,𝐱_𝐴=,𝜔_𝑖𝑊𝐱_𝑆^(𝑖)𝑝_sc𝐱_𝐵,𝐱_𝑆^(𝑖),𝜔𝑝_dir𝐱_𝐴,𝐱_𝑆^(𝑖),𝜔^∗, (B.1) where subscripts vsm and dir stand for “virtual source method” and “direct field,” respectively. Obviously, the direct field does not contain the full incident field 𝑝_inc. A particular section Δ𝑝_inc is not captured by the time-gate. We may substitute 𝑝_inc=𝑝_dir+Δ𝑝_inc into (30), to show that𝑝_sc𝐱_𝐵,𝐱_𝑆,𝜔−Δ𝑝_sc𝐱_𝐵,𝐱_𝑆=,𝜔_𝜕𝐴𝐺_ptb𝐱_𝐵,𝐱_𝐴,𝜔𝑝_dir𝐱_𝐴,𝐱_𝑆,𝜔𝑑𝐱_𝐴, (B.2) withΔ𝑝_sc𝐱_𝐵,𝐱_𝑆,𝜔=_𝜕𝐴𝐺_ptb𝐱_𝐵,𝐱_𝐴,𝜔Δ𝑝_inc𝐱_𝐴,𝐱_𝑆,𝜔𝑑𝐱_𝐴. (B.3) 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:𝐶_vsm𝐱_𝐵,𝐱^′_𝐴=,𝜔_𝜕𝐴𝐺_ptb𝐱_𝐵,𝐱_𝐴,𝜔Γ_vsm𝐱_𝐴,𝐱^′_𝐴,𝜔𝑑𝐱_𝐴+𝑔_vsm𝐱_𝐵,𝐱^′_𝐴,𝜔, (B.4) where 𝐶_vsm is the correlation function of the virtual source method (B.1) evaluated at receiver 𝐱^′_𝐴 instead of 𝐱_𝐴. Γ_vsm is the PSF of the virtual source method:Γ_vsm𝐱^′_𝐴,𝐱_𝐴=,𝜔_𝜕𝐴𝑊𝐱_𝑆𝑝_dir𝐱^′_𝐴,𝐱_𝑆,𝜔𝑝_dir𝐱_𝐴,𝐱_𝑆,𝜔^∗𝑑𝐱_𝑆. (B.5)𝑔_vsm is a ghost term associated with this strategy:𝑔_vsm𝐱_𝐵,𝐱_𝐴,𝜔=−_𝜕𝐴𝑊𝐱_𝑆Δ𝑝_sc𝐱_𝐵,𝐱_𝑆,𝜔𝑝_dir𝐱_𝐴,𝐱_𝑆,𝜔^∗𝑑𝐱_𝑆. (B.6) 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 𝐺_ptb blurred by Γ_vsm plus an additional ghost 𝑔_vsm. 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 Γ_ptb without a ghost term. By setting the time-gate, we have thus introduced an additional ghost. However, since Γ_vsm does only contain the cross-correlations of the direct field, this function behaves generally much more like the desired delta function than Γ_ptb. 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.