Seismic Interferometry: reconstructing the Earth’s reflection response

Can the reflection response of the subsurface be obtained without the need to place an active source at the surface? Using numerical modeling results based on one-way wavefields, we show that the reflection response of a 3-D inhomogeneous lossless medium can be reconstructed from the crosscorrelation of its transmission response at the surface. Depending on the type of subsurface noise sources, we can do this using two different procedures. When we can record distinct transmission response from each subsurface source, we first crosscorrelate the traces in each shot gather with a master trace from the same gather, then we resort the traces to common receiver gather and sum the traces in each receiver gather. When we can not make separate transmission recordings from the noise sources, but rather record their superposition, then assuming the sources to be uncorrelated in time we directly correlate the transmission panel with a master trace. The simulated reflection response reconstructs not only the kinematics, but also the amplitudes of the events. The reconstructed reflection response strongly depends on the number of the present subsurface sources. In the extreme of only a few sources, the simulated reflection events are of very bad quality. But if these reconstructed reflection responses are further migrated, the subsurface reflectors are imaged correctly even for the extreme of a few subsurface sources. Alternatively, the transmission responses at the surface can be migrated directly. The results of the two migration paths are the same. To reconstruct the reflection response in the elastic case, we crosscorrelate the decomposed one-way transmission responses. This implies that at the surface all the three component of the wavefield need to be measured. INTRODUCTION The idea to use ambient seismic noise in the exploration practice was first proposed by Claerbout (1968). There he showed that the reflection response of a 1-D medium can be synthesized from the autocorrelation of the transmission response measured at the surface. He named this method ”acoustic daylight imaging”. Later, he stated the conjecture that for 3-D media the crosscorrelation of the transmission responses measured at surface points A and B 2 from ambient noise sources in the subsurface will reconstruct the reflection response at point B as if from a source at A. Using numerical modeling, Rickett and Claerbout (1996); Rickett (1996) showed that to obtain a good reconstruction after the crosscorrelation, one need to have long time recordings and a lot of spatially uncorrelated white-noise sources. Schuster (2001) proposed to use crosscorrelation of existing reflection recordings at A and B to create a new reflection response at B as if from a source at A. If the simulated reflection responses are then migrated, an image of the subsurface is obtained. Schuster named this procedure ”seismic interferometric imaging”. In (Wapenaar et al., 2002), Wapenaar proved Claerbout’s conjecture for 3-D inhomogeneous anisotropic acoustic medium using one-way wavefield reciprocity theorem of the correlation type. Later, he extended this relation also for the elastic case (Wapenaar et al., 2004b). In the same article he also derived other relations between the reflection and the transmission response with applications in multiple removal, using multiples to create new primary response and reconstruction of the transmission response from the reflection response. FROM RECORDED ACOUSTIC TRANSMISSION TO SIMULATED REFLECTION Let us have a lossless, source-free 3-D inhomogeneous domain D (see Figure 1), embedded between the plan-parallel surfaces ∂D0 and ∂Dm. Just above ∂D0 we have a free surface and below ∂Dm the half space is homogeneous. For this configuration, we can write the following relation (Wapenaar et al., 2002) R (xA,xB, t) +R + (xA,xB,−t) = δ (xH,B − xH,A) δ (t) − ∫ ∂Dm T (xA,x,−t) ∗ T − (xB,x, t) dx. (1) In equation 1, R (xA,xB, t) is the upcoming reflection response of the subsurface measured at the surface point A from a source at surface point B; T (xB,x, t) is the upcoming transmission response of the subsurface measured at surface point B to a subsurface source at point with coordinate vector x; δ() stands for Dirac delta function; and ∗ stands for convolution operation. Both the reflection response and the transmission response include all internal and free-surface multiples. As we are using here decomposed wavefields, the sign ”+” after the reflection response R means that the source emitted down-going waves and the sign ”-” after T means that the source emitted up-going waves. Equation 1 states that the crosscorrelation of the transmission response measured at two surface points A and B in the presence of sources 3 along the surface ∂Dm will reconstruct, with a minus sign, the reflection response and its time reversal at A as if from an impulsive source at B. Even though the relation is derived for subsurface sources along a flat surface, the sources can be randomly distributed in depth as the crosscorrelation process eliminates the extra travel times. We will show in the following example how to apply formula (1) in practice. Figure 2 shows a 2-D subsurface velocity model used for acoustic modeling in this article. It consists of four layers separated by one syncline and two dipping line boundaries, respectively. There are 561 receivers at the free surface (x3 = 0) at a distance of 10 m from each other starting from x1 = 1200 m. The sources are situated between depth levels x3 = 700 m and x3 = 850 m with a random x3 coordinate. We place the subsurface sources in the horizontal direction every 25 m starting from x1 = 1200 m (in total 225 sources). For this configuration, we model a transmission shot gather for each source position using a finite-difference modeling code. Figure 3 shows an example transmission shot gather for a subsurface source at horizontal position x1 = 4000 m. We extract one ”master” trace from this panel (see Figure 4) and correlate it with the whole transmission panel. The result of this operation represents the term T (xA,x,−t)∗T − (xB,x, t) in equation 1. Now we perform the same procedure for all modeled transmission shot gathers. As a next step, we take all the traces with the same horizontal receiver position x1, for example 3000 m, from all correlated panels and sum them obtaining a single trace that is placed at position x1 = 3000 m in the resulting crosscorrelation panel. We repeat this operation for all receiver positions. The result represents the integral in equation 1. The output of the above procedure contains events at positive as well as at negative times. As the reflection response is a causal function of time, to obtain the final result, we can simply mute the negative times. The result is shown in Figure 5. Comparing the simulated reflection response from Figure 5 with the directly modeled reflection response as shown in Figure 6, we see that kinematically the two pictures are the same. Primaries as well as multiple arrivals are reconstructed correctly. Only the horizontally propagating modes in the simulated reflection are not reconstructed correctly as they are not included in the theory. Equation 1 was derived with the assumption that the medium below the sources is homogeneous. ? shows theoretically that the presence of inhomogeneities below the sources results in an extra error term on the left-hand side of relation 1. In (Draganov et al., 2004) we showed with numerical modeling results, that the presence of a reflector below 4 the sources causes ghost events to appear in the reconstructed reflection image. These ghost events, however, disappear when the subsurface sources are randomly distributed in depth. On the other hand, the real reflection arrivals from the reflector below the sources are correctly reconstructed after crosscorrelation (see the events indicated by the arrows in the Figures 5 and 6). Formula 1 does not limit itself only to kinematics. As it is exact in the light of the made assumptions, it also tells us that the amplitudes must be reconstructed correctly. So let us also compare the amplitudes of the corresponding traces from the simulated and the directly modeled reflection responses. Because equation 1 was derived using flux-normalized wavefields, to make a good comparison we need to flux-normalize both responses. Figure 7 shows the overlay of the directly modeled reflection response trace (green solid line) and the reconstructed reflection response trace (red dotted line) at the position of the surface impulsive source. We can see that also the amplitudes are reconstructed correctly. Comparing the corresponding traces away from the source position (Figures 8 and 9), we see that a discrepancy arises between the amplitudes of the two traces. This is explained with the fact that away from the surface source more energy propagates horizontally. As mentioned above, this energy is not included in the theory and the result is this amplitude difference between the reconstructed and the directly modeled traces. According to equation 1, to reconstruct the reflection response at the surface we need to have a continuous distribution of the subsurface sources (i.e., very closely spaced sources) and that we have to have distinct recordings of the transmission responses from each of the subsurface sources. In practice, it is not easy to find sources exploding at distinct times allowing separate transmission measurements. That is why, writing the integral over the source positions as a discrete sum and assuming the subsurface sources to be uncorrelated in time, we can rewrite equation 1 into R (xA,xB, t) +R + (xA,xB,−t) = δ (xH,B − xH,A) δ (t) − T obs (xA,x,−t) ∗ T − obs (xB,x, t) . (2) In equation 2, T obs (xA,x,−t) and T − obs (xB,x, t) represent the transmission responses recorded at the surface points A and B as if all the discretely present in the subsurface white-noise sources have exploded at the same time. Let us take again the subsurface model from Figure 2. In the presence of many subsurface white-noise sources exploding at the same time, the recorded at the surface transmission 5 response T obs (xB,x, t) will look like the one in Figure 10. In the figure we show the first 3 seconds of the complete transmission recording, which is 23 minutes long. As said above, each trace from this recording already contains in itself the sum over all present subsurface sources as they have exploded simultaneously. That is, according to equation 2, to reconstruct the reflection response we only need to extract one master trace from the transmission panel (like the example in Figure 11) and correlate it with the whole panel. The simulated reflection response is shown in Figure 12. Formulas 1 and 2 provide us with the means to reconstruct the reflection response in two different ways depending on the practical situation. When we can record at the surface distinct transmission responses from separate subsurface sources, for example different earthquakes, we can make use of equation 1. But when we do not know exactly when the subsurface source are active, we can just record at the surface the coming waves for some period of time and then make use of equation 2. Of course, in this case the quality of the obtained reconstructed reflection response will depend on the whiteness of the noise sources. FROM SIMULATED ACOUSTIC REFLECTION TO DEPTH IMAGE Once we have the reconstructed reflection responses we would like to see how the subsurface looks like when we migrate these reflections. We can do this in the standard way with a poststack or with a prestack migration scheme. Artman (Artman et al., 2004) proposed an alternative, where a depth image of the subsurface can be constructed by direct migration (which is also a crosscorrelation process naar proved mathematically this idea. Starting with a formula for downward extrapolation of the recorded at the surface reflection response he showed that the reflection response R (ξA, ξB, ω) can be written in the frequency domain as R (ξA, ξB, ω) = ∫ ∂D0 { W (ξA,xA, ω) }∗ T obs (xA, ω) dxA