Creating a virtual source inside a medium from reflection data: heuristic derivation and stationary-phase analysis

Summary

1 Introduction

Broggini et al. (2011, 2012a) discuss a new approach to creating the response to a virtual source inside a medium that goes beyond seismic interferometry. They show that, given the reflection response of a 1-D layered medium, it is possible to obtain the response to a virtual source inside the medium, without the need to know the medium parameters. The method consists of an iterative scheme, akin to earlier work of Rose (2001). Interestingly, the response retrieved by this new method contains all scattering effects of the layered medium. Note that to obtain the same virtual-source response by seismic interferometry one would need a receiver at the position of the virtual source inside the medium, and real sources at the top and bottom of the medium. Hence, the advantage of the new approach over 1-D seismic interferometry is that no receivers are needed inside the medium and that the medium needs to be illuminated from one side only. Broggini et al. (2011) speculate that the 1-D method can be extended to three dimensions. This would imply that the 3-D response to a virtual source in the subsurface could be retrieved from 3-D reflection measurements at the surface, without knowing the parameters of the 3-D medium. Hence, unlike for controlled-source interferometric methods (Schuster et al. 2004; Bakulin & Calvert 2006), no receivers would be required in the subsurface, nor would the lack of sources illuminating the medium from below cause spurious multiples (Snieder et al. 2006).

Recently we made a first step towards generalizing the method of Broggini et al. (2011) to the 3-D situation (Wapenaar et al. 2011). Using physical arguments we proposed an iterative scheme that aims to transform the reflection response of a 3-D medium into the response to a virtual source inside the medium. The proposed scheme requires, apart from the reflection response, an estimate of the direct arrivals between the virtual source and the acquisition surface. It is, in fact, through the arrival time of direct arrivals that one specifies the location of the virtual source. Hence, the method is not fully model-independent. Note, however, that a model that relates direct arrivals to a source position can be much simpler than a model that explains all internal multiple scattering. In the proposed method the multiple-scattering part of the retrieved virtual-source response comes entirely from the reflection data.

The proposed method has not yet been proven mathematically (except for the 1-D situation), nor have the limitations been exhaustively investigated. Here, we present a heuristic derivation of the method. We follow a two-step procedure, analogous to Broggini et al. (2011, 2012a). In step 1 (Section 3) we design a downgoing wavefield at the surface which focuses at the virtual-source position. In step 2 (Section 4) we use this downgoing field and the reflected upgoing field to create the virtual-source response. We discuss these two steps at the hand of a simple 2-D configuration (introduced in Section 2), using the method of stationary phase.

2 The Medium Configuration and the Reflection Impulse Response

We consider a configuration of two parallel dipping reflectors in a lossless, constant velocity, variable density medium (Fig. 1a). The only reason for choosing a constant velocity is that all responses obey simple analytical expressions. The proposed scheme is, however, not restricted to constant velocity media. We denote spatial coordinates as x = (x, z). The acquisition surface is located at z = 0 and is transparent (i.e. the upper half-space has the same medium parameters as the first layer). The first dipping reflector is defined as z = z₁−ax, with z₁ = 1000 m and a = 1/4. The red dot in Fig. 1(a) denotes the position of the virtual source, with coordinates x_VS = (x_VS, z_VS) = (100, 1400). The second reflector is parallel to the first reflector, so that all mirror images of the virtual source lie on a line perpendicular to the reflectors. This line obeys the relation z = z₁+x/a. The second reflector intersects this line at x = (150, 1600). The velocity in the medium is set to c = 2000 ms⁻¹. The densities in the three layers are ρ₁ =ρ₃ = 1000 kgm⁻³, and ρ₂ = 5000 kgm⁻³. The reflection coefficients for downgoing waves at the two interfaces are r₁ = (ρ₂−ρ₁)/(ρ₂+ρ₁) = 2/3 and r₂ = (ρ₃−ρ₂)/(ρ₃+ρ₂) =−2/3, respectively. The reflection coefficients for upgoing waves are −r₁ and −r₂. The transmission coefficients for downgoing (+) and upgoing (−) waves are and .

Figure 1

(a) Configuration with two dipping reflectors. (b) Initial incident downgoing field and the reflected upgoing field , both at z = 0.

Open in new tabDownload slide

We introduce the Green's function G(x, x_S, t) as a solution of the wave equation , with . Defined in this way, the Green's function is the response to an impulsive point source of volume injection rate at x_S (de Hoop 1995). In the frequency domain, obeys the equation , with . Here j is the imaginary unit, ω denotes angular frequency and the circumflex denotes the frequency domain. We write , where superscripts d and s stand for direct and scattered, respectively. We define the reflection impulse response at the surface in terms of via

1

for z_R = z_S = 0. Note that (apart from a factor −2) represents the pressure at x_R, related to an impulsive vertical force source at x_S or, via reciprocity, the vertical particle velocity at x_S, related to an impulsive point source of volume injection rate at x_R. Since we assumed that the acquisition surface is transparent, no surface-related multiples are present. In the Supporting Information we give a high-frequency expression for the reflection response of the medium of Fig. 1(a). In practice, , or in the time domain R(x_R, x_S, t), is obtained from the measured reflection data by surface-related multiple elimination and deconvolution for the source wavelet (Verschuur et al. 1992).

3 Focusing A Wave Field At the Virtual Source Position

We discuss an iterative procedure to design a downgoing wave field at z = 0 which, when injected into the medium, focuses at t = 0 at the virtual source position x_VS (and at x_VS only). This procedure uses the reflection impulse response at the surface and an estimate of the direct arrivals between the virtual source at x_VS and the surface. These direct arrivals will be represented by the direct Green's function G^d(x, x_VS, t) for z = 0. For the constant velocity model of Fig. 1(a), the high-frequency approximation of the Fourier transform of G^d(x, x_VS, t) is given by , with

2

with t^d(x, x_VS) =|x−x_VS|/c and μ= sign(ω). Note that the scaled direct Green's function requires only information on the traveltime t^d(x, x_VS) of the direct arrivals and not on the position and properties of the interfaces (for more general configurations we would require a smooth subsurface model to define the traveltimes of the direct arrivals).

When is back-propagated through the medium it focuses, amongst other points, at x_VS. Equivalently, when the time-reversal of is injected into the medium, it focuses at the same positions. Through an iterative procedure we determine a downgoing wave that focuses at x_VS only. The iteration starts with the initial incident downgoing wavefield at z = 0, defined as the time-reversal of the scaled direct Green's function, convolved with a symmetric wavelet s(t). Hence, , with A⁽⁰⁾ = 1, see Fig. 1(b) for t < 0 (the asterisk denotes temporal convolution). When propagating through the medium, this field not only focuses at x_VS, but after reflection at the first interface it creates a “ghost focus” at t = 0 at , which is the mirror image of x_VS with respect to the first reflector, see Fig. 1(a). We want a field that focuses at x_VS only, so the question is how we can modify the downgoing wavefield at z = 0 so that this ghost focus is cancelled. In practice we do not have a model of the medium available, but only its reflection impulse response R(x_R, x_S, t) at the surface. We discuss how we can use this response to obtain the desired modification of the downgoing field. To this end we first derive the reflected upgoing wavefield at the surface by convolving the incident downgoing field with the reflection impulse response and integrating over the source positions, according to

3

for z_R = 0. We analyse this integral with the stationary-phase method in the Supporting Information. The result is , see Fig. 1(b) for t > 0. Here A⁽ⁿ⁾ is defined as the product of reflection and transmission coefficients contributing to the nth event of the reflection impulse response R(x_R, x_S, t), hence A⁽¹⁾ = r₁, , , etc. Furthermore, was defined above, is the mirror image of x_VS in the second reflector, is the mirror image of x_VS in the mirror image of the first reflector with respect to the second reflector, etc. Note that the first reflection event, (indicated by r₁ in Fig. 1b), seems to originate from . Hence, this event gives us indirect information about the location of the ghost focus. If we back-propagate this event through the medium it focuses at . Equivalently, if we forward propagate its time-reversal through the medium it also focuses at . Hence, if we want to cancel the ghost focus at resulting from the initial downgoing field , we merely need to modify this downgoing field by subtracting the time-reversal of the first event of the reflected upgoing field (indicated by r₁ in Fig. 1b).

In general, the modified downgoing field, when injected into the medium, may give rise to new ghost foci (all shallower than x_VS), which in turn need to be cancelled by again modifying the downgoing field, etc. We propose the following iterative scheme

4

5

for x and x_R at z = 0, where w(x, t) is a window function to be discussed below. Eq. (4) describes the subtraction of the time-reversal of a windowed version of the reflected upgoing field from the initial downgoing field and eq. (5) gives an update of the reflected upgoing field. We next discuss the desired action of the window function in eq. (4). For the considered example, for k− 1 = 0, this window function should pass only the first event of (indicated by r₁ in Fig. 1b), which seems to originate from the ghost focus at . For more general configurations and arbitrary k it should pass all events that seem to originate from all ghost foci. Since these ghost foci occur at shallower depths than x_VS, the traveltimes of the corresponding events in are smaller than t^d(x, x_VS) (except for large horizontal source–receiver distances). Hence, we define the window function as w(x, t) =H(t^d(x, x_VS) −ε−t), where H(t) is the Heaviside function. Here ε is half the duration of the symmetric wavelet s(t), hence, a band-limited event centred at t^d(x, x_VS) is not passed by the window.

Assuming the iterative scheme of eqs (4) and (5) converges, the final downgoing field is equal to minus the time-reversed upgoing field between the traveltime curve t^d(x, x_VS) −ε and its time-reversal −t^d(x, x_VS) +ε. Moreover, the final downgoing field focuses at t = 0 at x_VS, whereas all ghost foci are cancelled. As a result of the antisymmetry between and , this cancellation even takes place when these ghost foci are dispersed, for example, due to reflector curvature or model errors.

We now apply the iterative scheme to the configuration of Fig. 1(a). The upgoing field is shown in Fig. 1(b) (t > 0). The traveltimes t^d(x, x_VS) −ε and −t^d(x, x_VS) +ε are indicated by the red curves in this figure. The window function passes the part of above the lower curve (i.e. the event indicated by r₁ in Fig. 1b). Hence, for k = 1 we have . Subtracting its time-reversal from , according to eq. (4), gives the modified incident wavefield . This is shown in Fig. 2(a) (t < 0). Using eq. (5) we evaluate the response to this modified incident wavefield. This response is the superposition of , evaluated above, and the response to . Following the same reasoning as before, the first two terms of this additional response seem to originate from mirror images of in the first and second reflector, hence, from (the original virtual source) and , respectively. The amplitude factors of these two terms are and , respectively. Higher order terms are evaluated in a similar way. Fig. 2(a) (t > 0) shows the total response, that is, the total upgoing field , with the amplitude factors of the different events indicated in the right margin. Note that, between the red curves, is equal to minus the time-reversal of , which means that the scheme converged in one iteration. Hence, and . This is a consequence of the simple configuration of Fig. 1(a). For more complex configurations more iterations are required.

Figure 2

(a) Result of the iterative focusing process. The final incident downgoing field focuses at t = 0 at x_VS. Between the red curves the final incident downgoing field is minus the time-reversed final reflected upgoing field. (b) Superposition of total field and its time-reversed version. This field is proportional to G(x, x_VS, t) +G(x, x_VS, −t).

Open in new tabDownload slide

4 Creating the Virtual-Source Response

The final downgoing field focuses at t = 0 at x_VS. The focused wavefield propagates through the focus point (similar as in interferometry, van Manen et al. (2006)), diverges again, scatters at the interfaces, and a part of the field arrives at the surface as . Because it apparently originates from the focus at x_VS, this upgoing field resembles the response to a virtual source at x_VS, but closer inspection reveals that it is not identical to it. In our example, the final upgoing field in Fig. 2(a) contains the scattering coda of the virtual-source response (this will be confirmed below), but the event at the direct arrival time t^d(x, x_VS) (just below the lower red curve) has the wrong sign and amplitude (). Moreover, the event indicated by r₁ arrives prior to the direct arrival and thus does not belong to the virtual-source response.

We construct the virtual-source response as follows. We define p(x, t) as the superposition of the final downgoing and upgoing wavefields: . For this example, iteration k = 1 was the final iteration, hence . This total field is shown in Fig. 2(a) (all t). Note that p(x, t) obeys the wave equation in the inhomogeneous medium of Fig. 1(a); Fig. 2(a) is just the cross-section of this field at z = 0. Between the red curves the field is antisymmetric in time (this is a consequence of the iterative scheme). Hence, if we superpose the total field and its time-reversed version, that is, p(x, t) +p(x, −t), all events between these curves cancel each other, see Fig. 2(b). Because we consider a lossless medium, p(x, t) +p(x, −t) obeys the wave equation. The causal part of this superposition is equal to (Fig. 2b, t > 0) and the acausal part is equal to (Fig. 2b, t < 0). Since time-reversal changes the propagation direction, it follows that the causal part is upward propagating at z = 0 and the acausal part is downward propagating at z = 0. The first arrival of the causal part of Fig. 2(b) corresponds with the direct arrival of the response to the virtual source at x_VS. Given this last observation, combined with the fact that the causal part is upward propagating at z = 0, and that the total field obeys the wave equation in the inhomogeneous medium and is time-symmetric, it is plausible that the causal part is proportional to the virtual-source response G(x, x_VS, t) *s(t) and hence the total field is proportional to {G(x, x_VS, t) +G(x, x_VS, −t)}*s(t). This argumentation holds for more general situations, but we will check it here for the response in Fig. 2(b). From the procedure that led to this response, we find for the causal part (taking into account that )

6

with the virtual source position and its mirror images defined in Sections 2 and 3. For the configuration of Fig. 1(a), this is equal to . Hence, the total field is given by

7

where G_h(x, x_VS, t) =G(x, x_VS, t) +G(x, x_VS, −t). Because G obeys the wave equation (Section 2), G_h obeys the homogeneous equation LG_h = 0. This is in agreement with the fact that p(x, t) +p(x, −t) has been constructed without introducing a singularity at x_VS. G_h is called the homogeneous Green's function, after Porter (1970) and Oristaglio (1989) (but note that these authors take the difference instead of the sum of the Green's function and its time-reversal, because of a different definition of the source in the wave equation).

5 Concluding Remarks

The iterative scheme of eqs (4) and (5) is akin to an iterative solution of inverse scattering methods. As shown by Burridge (1980), the Gel'fand-Levitan equation, the Marchenko equation and the Gopinath-Sondhi equations of inverse scattering in 1-D can be written (in symbolic notation) as

8

where R denotes the reflection data, F a known function that may depend on the reflection data, and K the function that one solves for. The integral ∫_W denotes integration over a windowed time interval. The integral eq. (8) can be solved by iteration by writing it as K = F−∫_WKR and by inserting the left-hand side of iteration k− 1 into the right-hand side of the kth iteration (Ge 1987). This gives

9

The iterative system of eqs (4) and (5) has the same structure as eq. (9), with replacing K_k, and , as defined in eq. (3), replacing F. The connection between iteration, time-reversal and inverse scattering based on the Marchenko equation in 1-D has been clarified in detail by Rose (2001, 2002a,b). Broggini et al. (2011, 2012a) show how to use the 1-D scheme to create a virtual source in an unknown medium. Our present work aims at generalizing this to 2-D and 3-D media that are illuminated from above.

The heuristic derivation in this paper gives insight in the mechanism of the iterative scheme and the stationary-phase analysis confirms the creation of the virtual-source response for a simple 2-D configuration. Following the arguments in Sections 3 and 4, it is plausible that the proposed methodology will also hold for more complex environments. Of course the proposed method will also have its limitations. The scaling factor in eq. (7) is in more general situations replaced by the cumulative angle-dependent transmission effects, which imposes an apparent radiation pattern upon the virtual source. The effects of a finite acquisition aperture, strong velocity variations, triplications, head waves, diving waves, fine-layering, errors in the direct arrivals, etc. need further investigation. A first numerical test with a variable-velocity syncline model and non-exact direct arrivals (Broggini et al. 2012b) shows promising results with respect to the handling of triplications. Errors in the estimated direct arrivals cause defocusing and mispositioning of the virtual source (similar as in standard imaging algorithms), but these errors do not deteriorate the reconstruction of the internal multiples (which come from the response of the actual medium). Hence, a smooth subsurface model that approximately explains the direct arrivals suffices. This is contrary to other wavefield extrapolation methods that account for internal multiples (Wapenaar et al. 1987; Vasconcelos et al. 2010), which are sensitive to the positioning of the discontinuities in the subsurface model.

For those configurations for which the proposed methodology applies, the potential applications are fascinating. Since no actual receivers are needed inside the medium, virtual sources can be created anywhere. The virtual-source responses contain the correct internal multiples, which can be used to improve seismic imaging and suppress the internal multiple ghosts (Wapenaar et al. 2012). Because the created virtual sources are independent of each other, the prediction and removal of internal multiples will not suffer from error propagation, unlike other imaging methods that aim at internal multiple suppression.

Acknowledgments

This work is supported by the Netherlands Research Centre for Integrated Solid Earth Science (ISES). We thank David Halliday and Andrew Curtis for their constructive reviews and for challenging us to improve the explanation of the method.

Supporting Information

Additional Supporting Information may be found in the online version of this article:

Appendix. Stationary-phase analysis.

Please note: Wiley-Blackwell are not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

Supporting Information

Additional Supporting Information may be found in the online version of this article:

Appendix. Stationary-phase analysis.

Please note: Wiley-Blackwell are not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

References