Seismic interferometry versus spatial auto-correlation method on the regional coda of the NPE

Seismic Interferometry and the Spatial Auto-Correlation Method

Seismic interferometry refers to the principle of generating new seismic responses through cross-correlations of recorded seismic wavefields at receivers (Schuster, 2001). The first to derive this principle for deterministic wavefields in 1D media was Claerbout (1968), who showed that the reflection response of a 1D medium could be synthesized from the auto-correlation of the transmission response. Later derivations include many different approaches based upon: diffusivity of the wavefields (Weaver and Lobkis, 2001; Roux et al., 2005; Sánchez-Sesma et al., 2006; Sánchez-Sesma and Campillo, 2006), stationary phase analysis (Schuster et al., 2004; Snieder, 2004) and propagation invariants and reciprocity theorems (van Manen et al., 2005; Weaver and Lobkis, 2004; Wapenaar and Fokkema, 2006; Claerbout, 1976; Wapenaar, 2004).

Wapenaar and Fokkema (2006) derives an interferometric representation of the Green's function from energy principles. Under certain conditions, this representation can be simplified to a direct cross-correlation between two receiver stations. Consider a domain $\mathbf{D}$ in an arbitrary, inhomogeneous medium enclosing two points $\mathbf{x}_A$ and $\mathbf{x}_B$, bounded by an arbitrarily shaped surface $\partial\mathbf{D}$ with outward pointing normal vector $\mathbf{n}$. The interferometric representation of the Green's function for the vertical component of particle velocity measured at $\mathbf{x}_A$ in response to a vertical force-impulse point-source acting at $\mathbf{x}_B$ is given in the frequency-domain as follows (Wapenaar and Fokkema, 2006):

$$2\Re\{\hat{G}_{3,3}^{\upsilon,f}(\mathbf{x}_A,\mathbf{x}_B,\omega)\} = \hspace{9cm}$$ (1)

$$- \oint_{\partial\mathbf{D}} \left[ \hat{G}_{3,ij}^{\upsilon,h}(\... ...bf{x}_A,\mathbf{x},\omega) \right\}^*\right] n_j \mathrm{d}^2\mathbf{x} \notag,$$ (2)

where the asterisk denotes complex conjugation, and $\omega$ denotes angular frequency. The notation convention for Green's functions is that superscripts denote the receiver (first) and source type (second), and subscripts denote the components of the source (first) and receiver (second) fields. The fields and sources in the elastodynamic system are particle velocity $\boldsymbol{\upsilon}$, stress tensor $\boldsymbol{\tau}$ (used below), external volume force density $\mathbf{f}$, and external deformation rate density $\mathbf{h}$. Einstein's summation convention is applied on all repeated subscripts.

The interferometric integral in equation 1 represents the real part of the elastodynamic Green's function between two receiver stations located at A and B, as a summation of cross-correlations of independent measurements at the two receiver stations. (Independent measurements of responses of various source components and types located on a surface enclosing both receivers are required to evaluate this integral.) The integral can be modified to reflect the field configuration of the NPE, where the receivers are located just below the traction-free surface, that has $\mathbf{n} = (0,0,1)$. The domain integral is split into two segments, $\partial \mathbf{D}_0$ and $\partial \mathbf{D}_1$, which are the parts of the domain boundary that coincide with the traction-free surface and the remainder, respectively. Thus the interferometric representation, equation 1, can be split into two parts:

$$2\Re\{\hat{G}_{3,3}^{\upsilon,f}(\mathbf{x}_A,\mathbf{x}_B,\omega)\} = \hspace{9cm}$$ (3)

$$- \oint_{\partial\mathbf{D}_1} \left[ \hat{G}_{3,ij}^{\upsilon,h}... ...hbf{x}_A,\mathbf{x},\omega) \right\}^*\right] n_j \mathrm{d}^2\mathbf{x} \notag$$ (4)

$$- \oint_{\partial\mathbf{D}_0} \left[ \hat{G}_{3,i3}^{\upsilon,h}... ...athbf{x},\omega) \right\}^*\right] \mathrm{d}^2\mathbf{x} \hspace{.25cm}\notag,$$ (5)

where $\mathbf{n} = (0,0,1)$ has been substituted into the integral segment over the traction-free surface, $\partial \mathbf{D}_0$. According to source-receiver reciprocity, the required response $\hat{G}^{\upsilon,h}_{3,i3}$ at the traction-free surface satisfies $\hat{G}^{\upsilon,h}_{3,i3}(\mathbf{x}_b,\mathbf{x},\omega) = \hat{G}_{i3,3}^{\tau,f}(\mathbf{x},\mathbf{x}_B,\omega) = 0$. Thus the second integral on the right hand side of equation 2 is equal to zero.

Following Wapenaar and Fokkema (2006), consider the situation when the medium outside $\mathbf{D}$ is homogeneous and when the wavefield is generated by many mutually uncorrelated sources located on $\partial \mathbf{D}_1$, acting simultaneously with a weighted power spectrum, $w(\mathbf{x})\hat{S}(\omega)$. The integral over $\partial \mathbf{D}_1$ in equation 2, can be evaluated by a direct cross-correlation between recordings of particle velocity at receiver stations A and B:

$$2\Re\{\hat{G}_{3,3}^{\upsilon,f}(\mathbf{x}_A,\mathbf{x}_B,\omega... ...hbf{x}_A,\omega) \hat{\upsilon}_3^*(\mathbf{x}_B,\omega) \frac{}{} \right>_{x},$$ (6)

where $c_p$ is the P-wave velocity. Note that the weighting factor $w(\omega)$ has disappeared. The weighting factor depends on the local medium parameters and source types, as is discussed at length by Wapenaar and Fokkema (2006). The wavefield required to make this simplification is named equipartitioned or diffuse (Hennino et al., 2001; Sánchez-Sesma and Campillo, 2006). The spatial ensemble average $\left< \cdot \right>_x$ over sources is usually evaluated using a sufficiently long recording. Secondary scattering can render the coda of the NPE sufficiently equipartitioned. But in the coda of the NPE, most of the coherent energy between the receivers resides in the surface-wave mode (de Ridder, 2008) and does not provide significant energy for imaging.

A similar situation occurs for earthquake tremor. Aki Aki (1957) developed a technique named the spatial auto-correlation method. The close relationship between SI and SPAC was reported by Yokoi and Margaryan (2008). Their steps are briefly repeated here to derive from equation 3 a relationship used in the SPAC method. For a wavefield dominated by the surface modes, the frequency-domain Green's function for the vertical component of particle velocity measured at $\mathbf{x}_A$ in response to a vertical force-impulse point-source acting at $\mathbf{x}_B$ is

$$\hat{G}^{\upsilon,f}_{3,3}(\mathbf{x}_A,\mathbf{x}_B,\omega) \app... ...at{m}_2(k_n,x_{3,B}) J_0(k_n\left\vert \mathbf{x}_A - \mathbf{x}_B\right\vert),$$ (7)

where $\hat{m}_2$ is a normalized eigenfunction, and $J_0$ is a zero-order Bessel function of the first kind (Yokoi and Margaryan (2008) from Aki and Richards (2002)). Substituting this Green's function into the left side of equation 3 expands it to

$$2\Re\left\{ -\omega \sum_{n=0}^{\infty} \hat{m}_2(k_n,x_{3,A})\ha... ...hbf{x}_A,\omega) \hat{\upsilon}_3^*(\mathbf{x}_B,\omega) \frac{}{} \right>_{x}.$$ (8)

Normalizing by the auto-correlation of the recording at station A gives

$$\frac{ 2\Re\left\{-\omega \sum_{n=0}^{\infty} \hat{m}_2(k_n,x_{3,... ..._3(\mathbf{x}_A,\omega) \hat{\upsilon}_3^*(\mathbf{x}_A,\omega) \right>_{x}}.$$ (9)

When the fundamental surface-wave dominates, and both receivers are located at equal depth ( $x_{3,A}=x_{3,B}$), the higher-order terms can be neglected, and equation 6 simplifies to

$$J_0(k_0\left\vert \mathbf{x}_A - \mathbf{x}_B\right\vert) \approx... ...hbf{x}_A,\omega) \right>_{x}} \approx \phi(\mathbf{x}_A,\mathbf{x}_B,\omega),$$ (10)

where $\phi(\mathbf{x}_A,\mathbf{x}_B,\omega)$ is defined as the azimuthally averaged auto-correlation coefficient. The wavenumber of the fundamental surface-wave mode is given by a specific dispersion curve, $k_0 = \frac{\omega}{c(\omega)}$, where $c(\omega)$ is phase velocity.

Notice from equation 7 how the cross-spectra in frequency and space are predicted to obey Bessel functions, with oscillations determined by the phase velocity. The Bessel function of the first kind is real-valued. The cross-spectrum on the right hand-side is complex-valued, but if the conditions that lead to equation 7 are fulfilled, the imaginary component vanishes (Asten, 2006), ( $\hat{\upsilon}_3(\mathbf{x}_A,\omega) \hat{\upsilon}_3^*(\mathbf{x}_A,\omega)$ is always real). The real part of the cross-spectrum is retrieved as the zero-lag temporal cross-correlation, i.e., a spatial auto-correlation coefficient. It should be noted that the close relationship between SI and SPAC seems to hold only for surface-waves in horizontally stratified media. The SPAC method as commonly applied involves fitting Bessel functions to the computed auto-correlation coefficient with frequency, with explicit directional averaging of the wavefield in all directions. In the case of isotropic wavefields, this averaging is unnecessary (Aki, 1957; Okada, 2003). The coda of the NPE quickly becomes isotropic after the first break, thus this relationship seems suitable for the cross-spectra calculated from NPE data. To introduce an estimation of attenuation, the Green's function, equation 4, is supplemented with an exponential attenuation factor $Q(\omega)$, (Aki and Richards, 2002). The final model for the frequency-domain spatial auto-correlation coefficient $\phi(\mathbf{x}_A,\mathbf{x}_B,\omega)$ becomes

$$\phi(\mathbf{x}_A,\mathbf{x}_B,\omega) = J_0\left(\frac{\omega}{c... ...frac{1}{2Q(\omega)} \left\vert \mathbf{x}_A - \mathbf{x}_B \right\vert\right\}.$$ (11)

Seismic interferometry versus spatial auto-correlation method on the regional coda of the NPE