2. Mathematical and Numerical Procedures

We will describe the procedure for local spectra extraction and on the determination of the two proxy data, the local quantities E-max(k) and k-max(k) in the wavelet wavenumber k space. We now describe briefly the mathematical and numerical procedures used for transforming with wavelets the 3-D seismic data that is taken at one time instant, namely, today.

Wavelets can be viewed as a mathematical operation involving a 3-D convolution integral, which can be set to zoom in both physical space at a given location and in spectral space at a given scale. Here we employ the isotropic continuous wavelet transform in 3-D, as described by Murenzi and Antoine (1996). The scaling is isotropic and the wavelet is dilated equally in each direction with the parameter a. We also imposed that the L2-norm, a mathematical measure of the magnitude, of the wavelet is scale independent. The wavelet can then be translated in the 3-D physical space along each direction independently. The transformation of a signal f in Cartesian space is given explicitly by the 3-D integral

$$\tilde{f}(a,\vec{b})= {1\over {{a}^{3/2}}} \int_{0}^{Lx}\int... ...c{x}) {\psi\left({{\vec{x}-\vec{b}}\over a}\right)}d^3\vec{x},$$ (1)

Where Lx, Ly and Lz are the lengths of the periodic box in the Cartesian space.

The best way to calculate the wavelet transform is to go to the Fourier space, since the convolution involved in eqn.  1 becomes then a multiplicative operation. The limit of integration in spectral space is thus given by the biggest wavenumber. The convolution kernel used here is the second derivative of the Gaussian function, also known as the Mexican hat (Daubechies, 1992). We have employed higher-order Gaussian derivatives up to the eighth-order and we have found that the effects are not noticeable. The expression of the Mexican-hat function is given in Fourier space, where $\vec{k}$ is the global wavenumber, as

$$\hat{\psi}(\vec{k})=\left(2\pi\right)^{3/2}~ {\left\vert{\vec... ...ght\vert}^2~ \exp{(-{{\left\vert{\vec{k}}\right\vert}^2/2})}.$$ (2)

We compare here qualitatively the concept of local and global spectra. In the former case, the information concerning the location in physical space is lost. This phenomenon takes place with Fourier spectrum obtained in Cartesian space as well as spectrum obtained from spherical harmonic decomposition on the sphere (Fig. 2 left). On the other hand, local spectra analysis keeps an explicit account of the spatial location. In Fig. 2, this point is illustrated by conducting the analysis under Africa (the position parameter b) using the characteristic scales that go from northern Europe to the center of Africa (the scale parameter a).

  

Figure 2: Comparison between global spectrum obtained from spherical harmonic decomposition over the sphere (left) and local spectrum obtained using wavelet analysis in Cartesian space (right).

Quantitatively, the local power spectrum of the signal is given by the L2-norm of the wavelet transform, that is, the square of the modulus of eqn. 1. The sign of the signal is important. In the case of SVA, it determines whether the local anomaly is slow (negative) relative to a 1-D background Earth model such as the PREM (Dziewonski and Anderson, 1981) or fast (positive). Obviously this information is lost by using the L2-norm alone. Consequently, we have represented the local sign energy by using E_k':

$$E_{k'}= sign(f(\vec{b})) \times {\left\vert \tilde{f}(a,\vec{b}) \right\vert}^2,$$ (3)

where k' is the norm of the local wavenumber (the inverse of the scale) averaged over a small ellipsoidal shell of thickness $\vert d\vec{k}\vert=2\pi / L$ in the local wavenumber.

Following Hudgins et al. (1993), we can compare the global spectrum obtained from the usual Fourier method and the wavelet transformation. We averaged the ensemble of the local wavelet spectra over every position in the wavelet space. It follows that the wavelet power spectrum is the Fourier power spectrum of the signal convolved with the power spectrum of the bandpass filter ${\widehat \psi}$:

$$\int \vec{db} \left\vert \tilde{f}(a,\vec{b}) \right\vert^2... ...at \psi}(a\vec{k})\vert^2 ~\vert{\widehat f}(\vec{k})\vert^2.$$ (4)

Moreover, this choice of the continuous wavelet transform enables one to interpolate the wavelet spectrum of the signal between the usual Fourier modes, since the scale parameter can be varied arbitrarily. we note that this analysis is non-orthonormal, which increases the computer storage.

The type of information produced by the local wavelet spectra analysis is very difficult to visualize, since at each location we need to render a whole spectrum. It also means that for a  N³ grid size, we would need at least N⁴ memory storage. This approach is prohibitive for a large database such as the ones to be expected with much higher resolution seismic network: For the small 10 Mbytes dataset used in Bergeron et al. (1999) this would already represent a respectable 1 Gbytes data set, and the sensitivity of the wavelet method to capture plume-like structures has been demonstrated in Bergeron et al. (2000a). In order to synthesize and assimulate this information more efficiently, we perform data-compression and extract two proxy quantities: The maximum of the local energy, E-max, and the related local wavenumber, k-max. The local wavenumber k' is related to the inverse of the scale  a by k'=1/a. These two proxy quantities are displayed in Fig. 3 where the local energy versus k' shows the locations of E-max and k-max.

  

Figure 3: Schematic diagram showing the disposition of the two proxy quantities E-max and k-max with the local wavenumber.

The total number of points used in the integration of equation 4 varies from 128x128x64 for the analytical case to 360x180x24 in the case of the seismic velocity anomalies obtained from the P1200 model of Zhou (1996). In the latter, each scale a is equal to one degree in longitude (about 110km) and the corresponding wavenumber is K=1/a. The algorithm itself is given in appendix B and can be skipped at the first reading.