Wavefield summation

To calculate the Fourier transform of the reflection response of the subsurface, $R({\bf x}_r,{\bf x}_s, \omega)$, Wapenaar et al. (2004), proves  

$$2\Re [R({\bf x}_r,{\bf x}_s,\omega)]= \delta({\bf x}_s-{\bf... ... T^*({\bf x}_s,{\bf \xi},\omega)\; \mbox{d}^2{\bf \xi} \;.$$ (1)

The vector ${\bf x}$ will correspond herein to horizontal coordinates, where subscripts r and s indicate different station locations from a transmission wavefield. After correlation they acquire the meaning of receiver and source locations, respectively, associated with an active survey. The RHS represents summing correlations of windows of passive data around the occurrence of individual sources, at locations $\xi$, on the domain boundary $\partial\! D_m$ that surrounds the subsurface region of interest. The transmission wavefields $T(\xi)$ contain only one subsurface source.

Equation 1 dictates that the correlations of transmission wavefields must only be from individual transmission wavefields recorded over an interval, t, during which a single source is actively probing the subsurface. In this case, the zero-time of the correlations are correctly shared by the output of each correlation operation since each result is zero-phase. If more than one source function is active during a time window, or it is impractical to window the raw data around individual sources, the result of correlating the raw data will not yield the reflection wavefield R.

In practice, raw data are collected over a long time and sources are weak and/or overlapping. For transmission wavefields, the time axis and the shot axis are naturally combined. If we assume that n_s individual sources, and the reflections that occur t seconds afterward, are distributed at intervals within the total recording time $\tau$, field data can only be parameterized $T_f({\bf x}_r,\tau)$. Both t and $\tau$ represent the real time axis, though I will parameterize different wavefields with them with the understanding that $\max(t)$ is the two-way time to the deepest reflector of interest and $\tau$ is the real time axis from the beginning to the end of the total recording time such that

$$\max(\tau)=\max(t)n_s+\sum_j^{n_s}{\tt wait}\!\mbox{-}\tt{time}_j\;.$$ (2)

I will further adopt the conventions $\omega$ for the frequency domain dual variable of t, and $\varpi$ for the frequency domain dual variable of $\tau$.

Without knowing when sources happen, and acknowledging that the wait-time between shots can also be negative, it is impossible to separate field data into individual wavefields when attempting to image with the ambient noisefield. In this case, equation 1 can only be implemented with a single time function of length $\tau$  

$$\tilde{R}({\bf x}_r,{\bf x}_s,\varpi)=T_f({\bf x}_r,\varpi) T_f^*({\bf x}_s,\varpi) .$$ (3)

The immediate ramification of this formulation can be considered by Fourier analysis. Correlation of more than a few hundred samples is more efficiently performed in the Fourier domain, $C(\omega)=B(\omega)A^*(\omega)$, so I will first consider the general definition of the discreet Fourier transform (DFT). The DFT of a signal $f(\tau)$ can be evaluated for a particular frequency $\varpi$,

$$F\vert _\varpi=\mbox{DFT}[f(\tau)]\vert _{\varpi}= \frac{1}{\sqrt{n_\tau}}\sum_{\tau} f(\tau)\mbox{e}^{-i\varpi \tau} \;.$$ (4)

If the long function $f(\tau)$ is broken into N short sections, g_n(t), of the same length, the amplitude of a particular frequency $\varpi$ can also be calculated  

$$\mbox{DFT}[f(\tau)]\vert _{\varpi}= N^{-3/2}\sum_{n=1}^N ... ...a}= N^{-3/2}\mbox{DFT}[\sum_{n=1}^J g_n(t)]\vert _{\omega} \;$$ (5)

by simply changing the order of summation for convenience. The only requirement for the equation above is for both of the two different length transforms to contain the particular frequency being calculated (ie. a particular frequency where $\varpi=\omega$ is possible as dictated by the Fourier sampling theorem). The longer transform has many more frequencies at a smaller sampling interval between the shared values calculated by the short transform of the stacked data.

For passive field data, $f(\tau)$ is $T_f(\tau)$, and the g_n(t) are $T({\bf x}_r,{\bf \xi},t)\vert _{\xi=n}$ as well as background noise between sources. Thus, correlating a single long recording from many sources implicitly stacks the wavefields from individual sources at each frequency. For passive data, the time axis and the shot axis are combined in nature (by sources refusing to wait in turn) and in processing (by seismologists incapable of or refusing to process individual time windows). However, this transform only supports a time signal of length t and is aliasing the long field record.

What of the intermediate frequencies $\varpi$ that would be lost by stacking the time windows? It is necessary to remove these. Fine sampling in frequency carries information about the late time samples of the signal's dual representation. After correlation, shot records in the time domain must be windowed to remove late lag correlations which are superpositions of correlations of the different sources convolved with the earth model. These are completely uninterpretable in terms of the desired product $R({\bf x}_r,{\bf x}_s,t)$ and will be noise in further processing. If the result needs truncation after inverse transform, it is more efficient to only transform the part of the result desired.

Time windowing has a Fourier dual operation. The Fourier sampling theorem, solved for $\Delta t$ is

$$\Delta t= 1/(N\Delta f) \; .$$

Subsampling the frequency axis increases $\Delta f$ by a, and reduces the number of samples to N/a. The new time domain trace length is $\Delta t \;N/a$. Removing every other frequency, a=2, halves the length of the trace in the time domain. This process is the Fourier dual of reducing the Nyquist frequency by subsampling the time axis. Importantly, the ability to rearrange the order of summation in equation 5 or subsample the Fourier domain representation of signal to alias it means that the spectra of the two different length transforms are not related by smoothing. Instead, the short one is a subsampled version of the long one where all frequency samples shared by the two are identical.

The left column of Figure  shows a processing flow of a simple time domain signal with a zoomed in view of each trace (the first 32^nd of the axis) on the right. The top trace is the input signal. The middle trace is its autocorrelation. The bottom trace maintains a part of the autocorrelation result deemed important. To compute the bottom trace, the input was subsampled by 8 in the Fourier domain, multiplied by its conjugate, and inverse transformed. To facilitate plotting, the trace was padded with zeros. Identical results are obtained by time-domain stacking and Fourier domain subsampling as long as the level of decimation does not alias the acausal lags of the correlation into the positive lags.

 

freq Figure 1 Right column is 32x zoom of left beginning at t=0. (top) Idealized signal of three identical subsurface sources. (middle) Autocorrelation. (bottom) Autocorrelation performed with every 8th frequency. Zero values are padded on the bottom trace to facilitate plotting.

Subsampling in one domain is not an identical operation to windowing in the dual domain. The periodicity of the Fourier transform dictates that subsampling leads to aliasing rather than true truncation. Aliasing the time domain is more efficiently performed by summing short time windows before making the correlations and thereby greatly reducing the amount of computation required. Knowing that the late time lags of correlation are aphysical for passive imaging, the above analysis shows that only the frequency samples at intervals $\omega$, associated with records of length t, need to be inverse transformed after correlation. Since correlation is linear, we only need to calculate the frequency domain representation of the field data at this reduced sampling interval. The definition of the DFT shows that this is equivalent to first stacking the long time axis of the data. The random background and instrument noise between sources is diminished by the stacking of the time axis which also decimates a potentially enormous data volume.

Under the assumption that all the source functions recorded in the data are white and uncorrelated, the summation of the source wavefields may not be to harmful and $\tilde{R}\approx R$. Further, if the sources are all continuously ringing, and thus zero-phase over the recording, the correlations will not have residual phase. In my previous reports, both of these assumptions were made (sometimes not intentionally), which I now believe highly improbably for a real earth experiment. To utilize bursts of local seismicity for imaging, cross-correlating traces to make shot gathers makes $\tilde{R} \not=R$. The situation is directly analogous to summing two or more shot-records. While it may be useful for some applications, this sum can not be treated as a single record with an impulsive source. Cross-talk is introduced due to the inability to separate energy from the distinct experiments.

 

• Stacking wavefields