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

Authors: Joost van der Neut; Maria Tatanova; Jan Thorbecke; Evert Slob; Kees Wapenaar
Journal: International Journal of Geophysics
 (2011)

---

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


---

## Body

## 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 xS. Receivers are located at xA and xB 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 xA and a receiver at xB (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 xA 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 xA and xB 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.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;xS denote the source locations, xA denote the receiver locations that are turned into virtual sources, and xB denote other receiver locations.
(a)(b)(c)(d)(e)(f)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 Figure2 we give a more schematic illustration of the problem formulated in Figure 1(a).  The aim is to redatum the source locations xS to a receiver location xA “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  r→∞; xS denotes a source location at ∂𝕊 and nS is the unit normal vector of this surface;  xO denotes a source at ∂𝕆 and  nO 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)|Ŝ(ω)|2[Ĝ(xB,xA,ω)+Ĝ*(xA,xB,ω)]=ÎS(xB,xA,ω)+ÎO(xB,xA,ω).
On the left hand side we find the Green’s function Ĝ(xB,xA,ω)  as if there was an injection rate point source at xA and a receiver for acoustic pressure at  xB. It is given in the frequency domain, indicated by the hat and angular frequency ω. We also find its acausal counterpart Ĝ*(xA,xB,ω), where superscript * denotes complex conjugation. Note that the retrieved response is bandlimited by  |Ŝ(ω)|2, where Ŝ(ω) is the spectrum of the source wavelet. On the right-hand side, we find two integrals:  ÎS and ÎO. ÎS is referred to as the “known integral,” obtained by cross-correlations of wavefields from existing source locations xS at ∂𝕊. ÎO is the “missing integral,” obtained by cross-correlation of wavefields from nonexisting source locations xO at  ∂𝕆. Even though the source locations xO are not present in a realistic survey, we keep ÎO in the representation, allowing us to quantify its contribution.First, let us consider the known integralÎS, as derived by Wapenaar and Fokkema [38]:(2)ÎS(xB,xA,ω)=-∫∂S1jωρ(xS)×[{nS⋅∇Sp̂(xB,xS,ω)}{p̂(xA,xS,ω)}*-{p̂(xB,xS,ω)}{nS⋅∇Sp̂(xA,xS,ω)}*]dxS.
Here p̂(xA,xS,ω) represents the pressure recording at  xA due to a monopole source at xS. These recordings are assumed to be Green’s functions (or impulse responses) Ĝ(xA,xS,ω) convolved with source wavelet Ŝ(ω). Normal vector nS points perpendicular (outward) to the source array ∂𝕊. Further, ∇S=(∂/∂x1,S,∂/∂x2,S,∂/∂x3,S)T, where xS=(x1,S,x2,S,x3,S)T and superscript T denotes the transposed. Hence, nS·∇Sp̂(xA,xS,ω) is interpreted as the response at xA due to a dipole source at xS.  ρ(xS) is the mass density at the source array.The representation of missing integralÎO is very similar. In this case responses of nonexisting source locations xO needs to be cross-correlated and nS needs to be replaced by the normal vector  nO along ∂𝕆, yielding
(3)ÎO(xB,xA,ω)=-∫∂O1jωρ(xO)×[{nO⋅∇Op̂(xB,xO,ω)}{p̂(xA,xO,ω)}*-{p̂(xB,xO,ω)}{nO⋅∇Op̂(xA,xO,ω)}*]dxO,
where ∇O is the spatial gradient at source location xO. Before addressing the implications of not evaluating integral ÎO, we focus our attention on ÎS. 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 p̂in(xA,xS,ω) and p̂out(xA,xS,ω), where superscripts in and out denote ingoing and outgoing fields at xS. We substitute p̂=p̂out+p̂in 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
(4)ÎS(xB,xA,ω)=∫∂S2jωρ(xS)×[{p̂in(xB,xS,ω)}{ns⋅∇Sp̂in(xA,xS,ω)}*+{p̂out(xB,xS,ω)}{ns⋅∇Sp̂out(xA,xS,ω)}*{ns⋅∇Sp̂in(xA,xS,ω)}*]dxS.
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 nS·∇Sp̂in≈-(jω/c(xS))cos(ϕ)p̂in, where c(xS) is the acoustic wave velocity along the source array and ϕ is the incidence angle of the dominant wavefield with respect to nS [38]. Similarly, we find for outgoing constituents that  nS·∇Sp̂out≈+(jω/c(xS))cos(ϕ)p̂out. We assume that ϕ is close to zero such that cos(ϕ)≈1. Substituting these approximations into (4) yields(5)ÎS(xB,xA,ω)≈∫∂S2ρ(xS)c(xS)[{p̂in(xB,xS,ω)}{p̂in(xA,xS,ω)}*-{p̂out(xB,xS,ω)}{p̂out(xA,xS,ω)}*{ns⋅∇Sp̂in(xA,xS,ω)}*]dxS.
Before proceeding, it is helpful to provide a similar analysis for the missing integral ÎO. We assume that both density and wave velocity are constant at r→∞, such that all wavefields that would be recorded at receivers due to missing source locations xO are ingoing at ∂𝕆 (such that p̂=p̂in and p̂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(6)ÎO(xB,xA,ω)=∫∂O2jωρ(xO){p̂(xB,xO,ω)}{nO⋅∇Op̂(xA,xO,ω)}*dxO.
We can further simplify this equation by substituting the far-field approximation for ingoing fields nO·∇Op̂≈-(jω/c(xO))p̂, yielding(7)ÎO(xB,xA,ω)≈∫∂O2ρ(xO)c(xO){p̂(xB,xO,ω)}{p̂(xA,xO,ω)}*dxO.
So far we have shown that a Green’s function can be retrieved by evaluation of integrals ÎS and ÎO. In practice, we are generally not that fortunate. First, we miss the source locations xO to compute ÎO. Second, we cannot discriminate between ingoing and outgoing wavefields to evaluate ÎS. 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 Ĉ:(8)Ĉ(xB,xA,ω)=∫∂S2W(xS)ρ(xS)c(xS){p̂(xB,xS,ω)}{p̂(xA,xS,ω)}*dxS.
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 W(xS) 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 ϕ(xA,xS) of the dominant contribution of the incident wavefield p̂(xA,xS,ω), for instance by ray tracing. In such cases, it can be beneficial to replace weighting factor W(xS) with W(xA,xS)=cosϕ(xA,xS) 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   p̂=p̂out+p̂in into (8) and rewrite the result as(9)Ĉ(xB,xA,ω)=∫∂S2W(xS)ρ(xS)c(xS)[{p̂in(xB,xS,ω)}{p̂in(xA,xS,ω)}*-{p̂out(xB,xS,ω)}{p̂out(xA,xS,ω)}*{ns⋅∇Sp̂in(xA,xS,ω)}*]dxS+∫∂S2W(xS)ρ(xS)c(xS)[{p̂(xB,xS,ω)}{p̂out(xA,xS,ω)}*+{p̂out(xB,xS,ω)}{p̂(xA,xS,ω)}*{ns⋅∇Sp̂in(xA,xS,ω)}*]dxS.
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 ĝS:(10)ĝS(xB,xA,ω)=∫∂S2W(xS)ρ(xS)c(xS)[{p̂(xB,xS,ω)}{p̂out(xA,xS,ω)}*+{p̂out(xB,xS,ω)}{p̂(xA,xS,ω)}*{p̂out(xA,xS,ω)}*]dxS.
In a similar fashion we identify a second ghost, due to the missing source locations xO. This is done by redefining the integral ÎO in (7) as(11)ĝO(xB,xA,ω)=-∫∂O2ρ(xO)c(xO){p̂(xB,xO,ω)}{p̂(xA,xO,ω)}*dxO.
Next, we substitute (11), (7), (5), (1) and (10) into (9) and rewrite the result as(12)Ĉ(xB,xA,ω)≈|Ŝ(ω)|2[Ĝ(xB,xA,ω)+Ĝ*(xA,xB,ω)]+ĝS(xB,xA,ω)+ĝO(xB,xA,ω).
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 Ĝ(xB,xA,ω) and its acausal counterpart Ĝ*(xA,xB,ω) plus two additional ghost terms. The first ghost ĝS 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 p̂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 ĝO, 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 r→∞. However: these conditions do not apply for integrals of cross-correlation products like the one in (11). The integrand in (11) decays with order 𝒪(r-2), whereas the integration surface ∂𝕆 grows with order 𝒪(r2), 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 𝒪(r-2) 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 Figure3. 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  xA  (green star) and a receiver at the same location xB=xA. The real source is always located at the earth’s surface location xS and additional sources xO are introduced to evaluate ghost gO. Location xE will play a role only later in this paper.Impedance (inkg·m-2·s-1) versus depth (in m) for synthetic models A–D. The red star corresponds to source location  xS, the green star to receiver/virtual source location xA=xB, the blue star to nonexisting source location xO, and the magenta star to nonexisting source location xE; models A–C have no free surface incorporated; model D has a free surface at 0 m depth.
(a)(b)(c)(d)In Figure4(a) we computed the desired reference response with wavelet |Ŝ(ω)|2 for model A (Figure 3(a)) by placing an actual source at the virtual source location xA. We can clearly see the direct arrival at t=0s and an event at t=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 xS, and hence all wavefields that arrive at xA are ingoing at xS, such that gS=0 (10). We computed the second ghost gO, using the additional source xO below all medium heterogeneities. Since no reflection from this source is registered at xA, the only contribution to gO stems from correlations of the direct field, mapping at t=0s in this zero-offset case. In Figure 4(b) we show the ghosts gS+gO. Note that indeed no other contribution is found outside t=0s. Therefore the correlation function C (8) is very similar to the true reflection response, except around t=0s. In Figure 4(c) we show C (black) and the sum of reference response and ghosts (dashed green), matching well.Cross-correlation of full fields for model A: (a) reference response by placing an active source at the virtual source locationxA; (b) computed ghost terms gS+gO; (c) causal part of the correlation function C (black) and the sum of reference response + ghosts (dashed green).
(a)(b)(c)In the previous example, the ghost termsgS+gO did not give a significant contribution outside  t=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 t=0.4s but also the reflections of the overburden at t=0.1s and t=0.3s and their multiples; see Figure 5(a). The ghosts gS+gO and the correlation function C are shown in Figures 5(b) and 5(c), respectively. Since no reflectors are present above the source, gS=0. However, the second ghost gO does give a significant contribution. The events at  t=0.1s and t=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 t=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  t=0.5s and t=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.Cross-correlation of full fields for model B: (a) reference response by placing an active source at the virtual source locationxA; (b) computed ghost terms gS+gO; (c) causal part of the correlation function C (black) and the sum of reference response + ghosts (dashed green).
(a)(b)(c)It has been demonstrated by Wapenaar [42] that ghost gO 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 gS+gO and the correlation function C are shown in Figures 6(b) and 6(c), respectively. Indeed, the contributions of ghost gO are minor and randomly distributed (Figure 6(b)). It can be shown analytically that placing infinitely many contrasts even completely eliminates gO. 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 (t=0.4s) and overburden (t=0.1s and t=0.3s), including their multiples (t=0.5s and t=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 gS≠0 for such scenario.Cross-correlation of full fields for model C: (a) reference response by placing an active source at the virtual source locationxA; (b) computed ghost terms gS+gO; (c) causal part of the correlation function C (black) and the sum of reference response + ghosts (dashed green).
(a)(b)(c)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 Figure3(d). As a result, the reference response does not only contain the desired reflection at  t=0.4s but also the primaries of the overburden contrasts at t=0.6s and  t=0.7s, the free surface reflection at t=0.8s; and free-surface multiples at t=0.9s and  t=1.0s; see Figure 7(a). Ghosts gS+gO and correlation function C are shown in Figures 7(b) and 7(c), respectively. We see strong spurious events at t=0.1s and t=0.2s  in the retrieved response (and a weaker event at t=0.3s), stemming from interactions of the first reflectors with the free surface. We also observe that the reflections at t=0.6s and t=0.7s, their multiples, and the free surface reflection at t=0.8s are hardly retrieved due to the missing source at  xO.Cross-correlation of full fields for model D: (a) reference response by placing an active source at the virtual source locationxA; (b) computed ghost terms gS+gO; (c) causal part of the correlation function C (black) and the sum of reference response + ghosts (dashed green).
(a)(b)(c)


## 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 Figure8 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 ρ0(x) and velocity c0(x). In this reference state all heterogeneities outside volume 𝕍 (namely the target) have been removed. These heterogeneities are referred to as the perturbations in density Δρ(x) and velocity Δc(x), where ρ(x)=ρ0(x)+Δρ(x) and c(x)=c0(x)+Δc(x) represent the density and velocity of the true physical medium. Fields that propagate in the reference state are referred to as incident fields p̂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  p̂sc, that is p̂=p̂inc+p̂sc. Vasconcelos et al. [11] have derived a representation for Green’s function retrieval of the scattered field between virtual source xA and receiver xB from cross-correlations of the scattered field at xB and the incident field at xA:(13)|Ŝ(ω)|2Ĝsc(xB,xA,ω)=ÎptbS(xB,xA,ω)+ÎptbE(xB,xA,ω),
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 ÎptbS, stemming from the cross-correlations of incident and scattered fields from the actual sources at ∂𝕊:(14)ÎptbS(xB,xA,ω)=-∫∂S1jωρ(xS)×[{nS⋅∇Sp̂sc(xB,xS,ω)}{p̂inc(xA,xS,ω)}*-{p̂sc(xB,xS,ω)}{nS⋅∇Sp̂inc(xA,xS,ω)}*]dxS.
The second integral, ÎptbE, stems from similar cross-correlations of nonexisting sources xE at  ∂𝔼:(15)ÎptbE(xB,xA,ω)=-∫∂E1jωρ(xE)×[{nE⋅∇Ep̂sc(xB,xE,ω)}{p̂inc(xA,xE,ω)}*-{p̂sc(xB,xE,ω)}{nE⋅∇Ep̂inc(xA,xE,ω)}*]dxE,
where ∇E is the spatial gradient at xE, nE is the outward-pointing unit normal vector to ∂𝔼,  and ρE 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(16)ÎptbS(xB,xA,ω)=∫∂S2jωρ(xS)×[{p̂scin(xB,xS,ω)}{nS⋅∇Sp̂incin(xA,xS,ω)}*+{p̂scout(xB,xS,ω)}{nS⋅∇Sp̂incout(xA,xS,ω)}*{nS⋅∇p̂incin(xA,xS,ω)}*]dxS.
Assuming that the relevant wavefields propagate near normal incidence, we can substitute far-field approximations nS·∇Sp̂incout≈+(jω/c(xS))p̂incout and nS·∇Sp̂incin≈-(jω/c(xS))p̂incin into (16), yielding(17)ÎptbS(xB,xA,ω)≈∫∂S2ρ(xS)c(xS)[{p̂scin(xB,xS,ω)}{p̂incin(xA,xS,ω)}*-{p̂scout(xB,xS,ω)}{p̂incout(xA,xS,ω)}*{p̂incin(xA,xS,ω)}*]dxS.
For integral ÎptbE 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 p̂inc=p̂incin and p̂incout=0). Cross-correlations between ingoing and outgoing fields cancel and the remaining terms can be merged [11], such that(18)ÎptbE(xB,xA,ω)=∫∂E2jωρ(xE)×[{p̂scin(xB,xE,ω)}{nE⋅∇Ep̂inc(xA,xE,ω)}*]dxE.
The spatial derivative can be approximated with far-field approximation nE·∇Ep̂inc≈-(jω/c(xE))p̂inc, such that(19)ÎptbE(xB,xA,ω)≈∫∂E2ρ(xE)c(xE)[{p̂scin(xB,xE,ω)}{p̂inc(xA,xE,ω)}*]dxE.
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 Ĉptb:(20)Ĉptb(xB,xA,ω)=∫∂S2W(xS)ρ(xS)c(xS){p̂sc(xB,xS,ω)}{p̂inc(xA,xS,ω)}*dxS.
To evaluate the consequence of this choice, we separate the scattered and incident wavefields in up- and downgoing constituents according to p̂inc=p̂incin+p̂incout and p̂sc=p̂scin+p̂scout and substitute these representations into (20). The result can be written as(21)Ĉptb(xB,xA,ω)=∫∂S2W(xS)ρ(xS)c(xS)[{p̂scin(xB,xS,ω)}{p̂incin(xA,xS,ω)}*-{p̂scout(xB,xS,ω)}{p̂incout(xA,xS,ω)}*{p̂incin(xA,xS,ω)}*]dxS+∫∂S2W(xS)ρ(xS)c(xS)[{p̂sc(xB,xS,ω)}{p̂incout(xA,xS,ω)}*+{p̂scout(xB,xS,ω)}{p̂inc(xA,xS,ω)}*]dxS.
The first integral in (21) can be identified as the desired integral ÎptbS (see (17)). The second integral is a ghost term that we redefine as(22)ĝptbS(xB,xA,ω)=∫∂S2W(xS)ρ(xS)c(xS)[{p̂sc(xB,xS,ω)}{p̂incout(xA,xS,ω)}*+{p̂scout(xB,xS,ω)}{p̂inc(xA,xS,ω)}*]dxS.
The missing sources at xE make up for a second ghost term that we define by rewriting the integral in (19) in the following way:(23)ĝptbE(xB,xA,ω)=-∫∂E2ρ(xE)c(xE)[{p̂scin(xB,xE,ω)}{p̂inc(xA,xE,ω)}*]dxE.
Now, by substituting (23), (19), (17), (13), and (20) into (21) we arrive at(24)Ĉptb(xB,xA,ω)≈|Ŝ(ω)|2Ĝsc(xB,xA,ω)+ĝptbS(xB,xA,ω)+ĝptbE(xB,xA,ω).
Equation (24) is useful for identifying ghosts in interferometry by CC of perturbed fields. Computation of Cptb (see (20)) yields a scaled bandlimited version of the desired scattered Green’s function plus two ghost terms gptbS and gptbE. The first ghost gptbS stems from reflectors outside the integration volume and behaves very similar to the ghost gS that we found for CC of full fields. The second ghost gptbE is relatively small. It consists of ingoing waves at xE 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; xS denotes a source location at ∂𝕊 and nS is the unit normal vector of this surface;  xE denotes an additional source location at ∂𝔼 and  nE 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 (Figure3(b)), with the additional source xE located at 1000 m depth (instead of the source at xO). Incident fields are computed in a reference medium with all heterogeneities below xE 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 xA and computing the scattered response at the same receiver; see Figure 9(a). We can clearly see not only the desired reflector at t=0.4s but also the multiple reflections from the overburden at t=0.5s and t=0.7s as well as higher-order multiples. In Figures 9(b) and 9(c) we show the ghosts  gptbS+gptbE and the correlation function Cptb.  Note which the spurious event at t=0.2s, that was visible in Figure 5(c), is significantly weakened and hardly visible. The multiples at  t=0.5s and t=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  xE.Cross-correlation of perturbed fields for model B: (a) reference response by placing an active source at the virtual source locationxA; (b) computed ghost terms gptbS+gptbE; (c) correlation function Cptb (black) and the sum of reference response + ghosts (dashed green).
(a)(b)(c)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 (Figure3(d)) as well. In Figure 10(a) we show the reference response. We observe the target event at t=0.4s and a multiple at t=1.0s. Ghosts gptbS+gptbE and Cptb are shown in Figures 10(b) and 10(c), respectively. Note that we can still observe several spurious events and the undesired overburden reflections.Cross-correlation of perturbed fields for model D: (a) reference response by placing an active source at the virtual source locationxA; (b) computed ghost terms gptbS+gptbE; (c) correlation function Cptb (black) and the sum of reference response + ghosts (dashed green).
(a)(b)(c)


## 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 Figure2. We define p̂+(xA,xS,ω) and p̂-(xB,xS,ω) as the downgoing and upgoing fields at receivers xA and xB, respectively. Instead of evaluating the correlation function of full fields (see (8)), we define Ĉdcp as the correlation function between p̂+  and p̂-; that is(25)Ĉdcp(xB,xA,ω)=∫∂S2W(xS)ρ(xS)c(xS)p̂-(xB,xS,ω){p̂+(xA,xS,ω)}*dxS,
where subscript dcp denotes “decomposed.” As a consequence of this choice, we will reconstruct a Green’s function Ĝ-,+(xB,xA,ω) as if there was a downgoing field emitted at xA and an upgoing field registered at xB plus an acausal Green’s function {Ĝ+,-(xA,xB,ω)}* as if there was an upgoing field emitted at xB and a downgoing field registered at  xA [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  xA and upgoing at xB; that is(26)Ĉdcp(xB,xA,ω)≈|Ŝ(ω)|2[Ĝ-,+(xB,xA,ω)+{Ĝ+,-(xA,xB,ω)}*]+ĝdcpS(xB,xA,ω)+ĝdcpO(xB,xA,ω),
where the decomposed ghosts are given by(27)ĝdcpS(xB,xA,ω)=∫∂S2W(xS)ρ(xS)c(xS)[{p̂-(xB,xS,ω)}{p̂+,out(xA,xS,ω)}*+{p̂-,out(xB,xS,ω)}{p̂+(xA,xS,ω)}*]dxS,ĝdcpO(xB,xA,ω)=-∫∂O2ρ(xO)c(xO){p̂-(xB,xO,ω)}{p̂+(xA,xO,ω)}*dxO.
Here p̂+,out(xA,xS,ω) refers to the field that is outgoing at source xS and downgoing at receiver xA. p̂-,out(xB,xS,ω) refers to the field that is outgoing at source xS and upgoing at receiver xB.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 G-,+(xB,xA,t)  plus two additional ghost terms gdcpS and gdcpO. 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 t=0.4s and a weak multiple at t=0.9s; see Figure 11(a). The ghosts gdcpS+gdcpO and correlation function Cdcp 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 t=0.2s and t=0.6s. Compared to interferometry by cross-correlation of perturbed fields (Figure 9(c)), the multiples at t=0.5s and t=0.7s have been effectively removed.Cross-correlation of decomposed fields for model B: (a) reference response by placing an active source at the virtual source locationxA; (b) computed ghost terms gdcpS+gdcpO; (c) correlation function Cdcp (black) and the sum of reference response + ghosts (dashed green).
(a)(b)(c)Wavefield decomposition does still not offer a solution for free-surface-related artifacts. To illustrate this, we apply the method to model D (Figure3(d)). In Figure 12(a) we show the reference response, containing only the desired reflection at t=0.4s. Ghosts gdcpS+gdcpO and the retrieved correlation function  Cdcp 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)).Cross-correlation of decomposed fields for model D: (a) reference response by placing an active source at the virtual source locationxA; (b) computed ghost terms gdcpS+gdcpO; (c) correlation function Cdcp (black) and the sum of reference response + ghosts (dashed green).
(a)(b)(c)


## 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 p̂=p̂inc+p̂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 xS and a receiver at xB, where xS is inside volume 𝕍 and xB is outside this volume, that is(28)p̂sc(xB,xS,ω)=ĴptbA(xB,xS,ω)+ĴptbF(xB,xS,ω).
On the right-hand side we find two integrals.  ĴptbA involves integration over the receiver array ∂𝔸:(29)ĴptbA(xB,xS,ω)=∫∂A1jωρ(xA)×[{nA⋅∇AĜsc(xB,xA,ω)}{p̂inc(xA,xS,ω)}-{Ĝsc(xB,xA,ω)}{nA⋅∇Ap̂inc(xA,xS,ω)}]dxA.∇A is the spatial gradient at receiver location xA and nA is the outward-pointing unit normal vector of ∂𝔸. Integral ĴptbF is a similar integral over ∂𝔽, which vanishes since all pressure recordings are zero at this interface. We rewrite the wavefield p̂ at ∂𝔸 in terms of ingoing components p̂in and outgoing components p̂out. Similarly, Ĝsc=Ĝscout+Ĝscin. Since the reference medium is homogeneous outside 𝕍, p̂inc is purely outgoing at xA (so p̂inc=p̂incout and p̂incin=0). It can be shown that the convolutions between ingoing fields at the virtual source xA and outgoing fields p̂inc(xA,xS,ω) cancel and that the remaining terms can be merged [22]. Consequently, (28) and (29) can be rewritten as(30)p̂sc(xB,xS,ω)=∫∂AĜptb(xB,xA,ω)p̂inc(xA,xS,ω)dxA,
where we have introduced(31)Ĝptb(xB,xA,ω)=2jωρ(xA)nA⋅∇AĜscout(xB,xA,ω).
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 Ĝinc(xA′,xS,ω) and summing over xS (where xA′ is at ∂𝔸). Hence(32)Ĉptb(xB,xA′,ω)=∫∂AĜptb(xB,xA,ω)Γ̂ptb(xA,xA′,ω)dxA.
On the left-hand side of (32) we have the correlation function of the incident field at xA′ and the scattered field at  xB:(33)Ĉptb(xB,xA′,ω)=∑iW(xS)p̂sc(xB,xS(i),ω){p̂inc(xA′,xS(i),ω)}*,
Note that (33) is identical to a discrete scaled representation of (20), if the medium parameters are constant at ∂𝕊 and xA′ is replaced by xA. On the right-hand side of (32) we have the so-called point-spread-function (PSF) for perturbed fields, defined as(34)Γ̂ptb(xA,xA′,ω)=∑iW(xS(i))p̂inc(xA,xS(i),ω){p̂inc(xA′,xS(i),ω)}*.
In both expressions, W(xS) is an optional source-dependent weighting function. If sufficient source and receiver locations are present, a multidimensional inverse Γ̂ptbinv of the PSF can be introduced, according to (35)∫∂AΓ̂ptb(xA,xA′,ω)Γ̂ptbinv(xA′,xA′′,ω)dxA′=δ(xA-xA′′).
where xA′′ is at ∂𝔸. The desired dipole response can now be found by filtering the correlation function with Γ̂inv, according to(36)Ĝptb(xB,xA′′,ω)=∫∂AĈptb(xB,xA′,ω)Γ̂ptbinv(xA′,xA′′,ω)dxA′,
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  ∂𝔸; xS denotes a source location;  xA  and  xB denote receiver locations at ∂𝔸 and  nA is the unit normal vector of the receiver array; nF 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 (Figure3(b)). The result is shown in Figure 14(a) (black) and compared with the reference response Gptb (dashed green). Note that the MDD response is similar to the correlation function Cptb (Figure 9(c)), apart from a very weak spurious event at t=0.2s, a scaling factor and subtle amplitude variations. We discussed that Cptb can be interpreted as the desired response Gptb, 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 t=±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 Gptb and Cptb are indeed very similar apart from a scaling factor.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  Cptb (dashed green).
(a)(b)(c)This is not the case for model D (Figure3(d)). The response of MDD of perturbed fields (black) is compared with the reference response Gptb (dashed green) in Figure 15(a). Compared to the correlation function Cptb (Figure 10(c)), the MDD response is very simple, containing only the target reflection at t=0.4s and a multiple at t=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 Cptb.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  Cptb (dashed green).
(a)(b)(c)


## 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 r→∞. 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 G̅̂. We formulate a convolution based representation for the field of source xS at receiver xB, where xS is outside volume 𝕍 and xB is inside this volume, reading(37)p̂(xB,xS,ω)=ĴdcpA(xB,xS,ω)+ĴdcpO(xB,xS,ω).
On the right-hand side we find two integrals. ĴdcpA involves integration over the receiver array ∂𝔸:(38)ĴdcpA(xB,xS,ω)=-∫∂A1jωρ(xA)×[{nA⋅∇AG̅̂(xB,xA,ω)}{p̂(xA,xS,ω)}-{G̅̂(xB,xA,ω)}{nA⋅∇Ap̂(xA,xS,ω)}]dxA.
Integral ĴdcpO is a similar integral over the interface ∂𝕆. Since this integral contains convolutions and its radius r→∞, it will vanish due to Sommerfeld’s radiation conditions [41]. Since the reference state is homogeneous outside 𝕍, G̅̂(xB,xA,ω) is purely ingoing at virtual source xA (G̅̂=G̅̂in and  G̅̂out=0). We substitute p̂=p̂in+p̂out. It can be shown that the convolutions of ingoing fields at the virtual source xA with outgoing fields p̂out(xA,xS,ω) cancel and that the remaining terms can be merged, yielding (39)p̂(xB,xS,ω)=∫∂AG̅̂dcp(xB,xA,ω)p̂in(xA,xS,ω)dxA,
where we have introduced(40)G̅̂dcp(xB,xA,ω)=-2jωρ(xA)nA⋅∇AG̅̂(xB,xA,ω).
In a similar way as for MDD of perturbed fields, a normal equation can be derived:(41)Ĉdcp(xB,xA′,ω)=∫∂AG̅̂dcp(xB,xA,ω)Γ̂dcp(xA,xA′,ω)dxA,
where xA′ is at ∂𝔸. On the left-hand side we have the correlation function of the ingoing field at xA′ and the outgoing field at  xB:(42)Ĉdcp(xB,xA′,ω)=∑iW(xS)p̂out(xB,xS(i),ω){p̂in(xA′,xS(i),ω)}*.
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 xA′ is replaced by xA. On the right-hand side we have the Point-Spread-Function (PSF) for decomposed fields, given by(43)Γ̂dcp(xA,xA′,ω)=∑iW(xS)p̂in(xA,xS(i),ω){p̂in(xA′,xS(i),ω)}*.
If acquisition conditions allow, we can take a multidimensional inverse of the PSF Γ̂dcpinv and convolve it with the correlation function according to(44)G̅̂dcp(xB,xA′′,ω)=∫∂AĈdcp(xB,xA′,ω)Γ̂dcpinv(xA′,xA′′,ω)dxA′,
where xA′′ is at ∂𝔸 and(45)∫∂AΓ̂dcp(xA,xA′,ω)Γ̂dcpinv(xA′,xA′′,ω)dxA′=δ(xA-xA′′).
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 ∂𝕆; xS denotes a source location;  xA  and  xB denote receiver locations at ∂𝔸 and  nA is the unit normal vector of the receiver array; nO 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 (Figure3(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 t=±0.2s and t=±0.5s. We convolve the reference response with the PSF to show that indeed the correlation function emerges, see Figure 17(c).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  Cdcp (dashed green).
(a)(b)(c)The effects of MDD of decomposed fields are more dramatically exposed by model D (Figure3(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 t=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.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  Cdcp (dashed green).
(a)(b)(c)


## 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 Figure19 we show a salt flank model, defined as a function of coordinates x1 (horizontal distance) and x3 (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 x3∈[1000m,5000m] at x1=15000m. We place 401 sources at the surface with 30 m spacing along the interval x1∈[0m,12000m] at x3=0m. No free surface is incorporated. In the following we generate a virtual shot record as if there was a source at receiver xA=(15000m,4000m) in the well using CC of perturbed fields. We generate a reference response by placing an active source at xA 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 x1∈[0m,12000m] at x3=0m; receivers are located at boundary ∂𝔸, covering x3∈[1000m,5000m] at x1=15000m; a reference virtual source is located at xA=(15000m,4000m).Snapshots at 2.0 s, 2.7 s and 3.4 s of the wavefield emitted by a source atxA=(15000m,4000m); reflections 1, 2, and 3 are indicated.
(a)(b)(c)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 sourcexA 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 atxA=(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 atxA=(15000m,4000m) with the scattered field at xB=(15000m,2000m) for different source locations xS, 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 atxA=(15000m,4000m) and a receiver at xB=(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 xA=(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 xA is correlated with the incident field at xA′ (“integrand” of (34)). If  xA=xA′, this corresponds to auto-correlation, providing a significant contribution at t=0, see Figure 24. These contributions will interfere constructively, to generate the desired spike in the PSF. If  xA≠xA′, 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 atxA=(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 atxA=(15000m,4000m) and xA′=xA=(15000m,4000m).Figure 25
PSF correlation gather for receivers atxA=(15000m,4000m) and xA′=(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 Figure26. 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 (x=xA and t=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 atxA=(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 Figure28. 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 (x=xA and t=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 atxA=(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 Figure30 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 atxA=(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 atxA=(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 Figure32 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 aroundt=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.



---

*Source: 870819-2011-06-08.xml*