Frequency-wavenumber-analysis

Seismic arrays are able to measure the properties of propagating waves. The arrays were built to detect the origin of seismic signals. This can be done by determining the vector velocity of a seismic wave. The vector velocity can be separated into the two components, horizontal slowness $u$ (i.e. the inverse of the apparent horizontal velocity) and backazimuth $\Theta$. Several methods to measure slowness and backazimuth with array data have been developed (Davies et al. 1971; Capon, 1973; Posmentier and Herrmann, 1971). In this study, the frequency-wavenumber analysis (fk-analysis) is used. The fk-analysis is a standard array technique which simultaneously calculates the power distributed among different slownesses and directions of approach (Capon, 1973; Harjes and Henger, 1973; Aki and Richards, 1980).
Most array methods use time delays of the signals recorded at different array stations and summation of the coherent data stream from each station to improve the signal-to-noise ratio (SNR) for a specific phase. The time delays required to bring the signals into phase provide a direct estimate of the backazimuth and the slowness of the signal. If the slowness and backazimuth of a signal are unknown, a grid search for all u and $\Theta$ combinations can be performed to find the best parameter combination, producing the highest amplitudes of the summed signal. This computation is performed in the spectral domain due to lower computational time.
The following derivation of the fk-analysis follows Kelly (1967) and Harjes and Henger (1973).
A signal arriving at a reference point within the array with a horizontal velocity, $v_{s}$, and a backazimuth, $\Theta$, can be written as $s(t)$. The n-th seismometer with the location vector, ${\underline r_{n}}$, relative to the array reference point records the signal, $x_n(t)$:

$$x_{n}(t) = s(t-\underline{u_0} \cdot \underline{r_{n}})$$ (1)

where $\underline{u_0}$ is the slowness vector with

$$\underline{u_0} = \frac {1}{v_{s}} \cdot (\cos(\Theta),\sin(\Theta)).$$ (2)

The maximum amplitude of the sum of all array seismometers is reached if the signals of all stations are in phase, i.e. if the time shifts $\underline{u_0} \cdot \underline{r_{n}}$ disappear (beamforming).
The output of the array can be computed by

$$y(t) = \frac{1}{N} \sum_{n=1}^N x_{n}(t + \underline{u_0} \cdot \underline{r_{n}})$$ (3)

for an array of N elements. For a signal with a different slowness vector $u$ the beam trace is computed using equation (5.1):

$$y(t) = \frac{1}{N} \sum_{n=1}^N s(t+((\underline{u_0} - \underline{u}) \cdot \underline{r_n}).$$ (4)

The total energy recorded at the array can be calculated by the integration of the squared summed amplitudes over time

$$E(k - k_0) = \int\limits_{-\infty}^{\infty} y^{2}(t) dt$$

$$= \frac{1}{2\pi} \int\limits_{-\infty}^{\infty} \left\vert S(\omega... ...\underline{k} - \underline{k_0}) \cdot \underline{r_n}} \right\vert^{2} d\omega$$ (5)

using Parseval's theorem. In equation (5.5), S$(\omega$) is the Fourier transformation of $s(t)$ and k is the wavenumber vector with

$$\underline{k} = (k_x, k_y) = \omega \cdot \underline{u}$$

$$= \frac{\omega}{v_s} \cdot (\cos(\Theta), \sin(\Theta))$$ (6)

and $\underline{k_0}$ the wavenumber vector for $\underline{u_{0}}$. The backazimuth determines the direction of $\underline{k}$ and the slowness the magnitude of $\underline{k}$. Equation (5.5) can be written as:

$$E(\underline{k} - \underline{k_0}) = \frac{1}{2\pi}\int\limi... ...vert A( \underline{k} - \underline{k_0} \right)\vert^2 d\omega$$ (7)

with

$$\left\vert A(\underline{k} - \underline{k_0})\right\vert^2 =... ...ine{k} - \underline{k_0}) \cdot \underline{r_n})}\right\vert^2$$ (8)

as array response function (ARF).
As can be seen from equation (5.7), the total energy recorded at the array is defined by the power spectral density $\left\vert S(\omega)\right\vert^2$ and the ARF $\left\vert A(\underline{k} - \underline{k_0})\right\vert^2$. The ARF is controlled by the design (aperture, configuration and interstation spacing) of the array.
The result of the fk-analysis is a power spectral density as a function of slowness and backazimuth. The slowness can be calculated from the wavenumber vector $\underline{k} = (k_x, k_y)$

$$\vert\underline{k}\vert = (k_x^2 + k_y^2)^{1/2}= \frac{2\pi}{\lambda_s} = \frac{\omega}{v_s}$$ (9)

with $\lambda_s$ as the apparent horizontal slowness.
The backazimuth $\Theta$ can be calculated by

$$\Theta = tan^{-1}\left(\frac{k_x}{k_y}\right).$$ (10)

The power spectral density is displayed in a polar-coordinate system called the fk-diagram. In the fk-diagram, the backazimuth is plotted on the azimuthal axis and the slowness on the radial axis.
Figure 5.2a) shows the array response of Yellowknife. The ARF was computed for a monochromatic wave with a frequency of 1 Hz arriving at an array with YKA configuration. It has been assumed that the power spectral density $S(\omega)$ is normalized

$$\frac{1}{2\pi} \int\limits_{-\infty}^{\infty} \left\vert S(\omega))\right\vert^{2} d\omega = 1$$ (11)

which leads to

$$E(\underline{k} - \underline{k_0}) = \left\vert A(\underline{k} - \underline{k_0})\right\vert^2$$ (12)

where the wave arrives with a slowness u = 0, i.e. $\underline{k_0}$ = 0. The logarithmic power is colour coded and the maximum is normalized to 0 dB. The lines are -1 dB contour isolines. The slowness and backazimuth of the maximum power is marked by the white circle. The maximum slowness displayed is 12 s/$^{\circ}$ with 2 s/$^{\circ}$ tick intervals. The azimuth describes the full circle from 0$^{\circ}$ to 360$^{\circ}$ clockwise. The cross-shaped form of the ARF resembles the geometry of YKA and is important for the detection of P$^d$P as discussed in chapter 5.2.1. Figure 5.2b) shows the logarithmic colourscale valid for all fk-diagrams displayed in this thesis.

Figure 5.2: a) Array response function of YKA computed for 1 Hz data. The logarithmic power spectral density is displayed as colour and -1 dB isolines are added. The maximum power density is normalized to 0 dB (red). The slowness with 2 s/$^{\circ}$ per tick is displayed on the radial axis, the backazimuth clockwise from 0$^{\circ}$ to 360$^{\circ}$. The maximum power is marked by the white circle.
b) Colour scale for all fk-diagrams in this thesis.
c) Same as a) but for a synthetic monochromatic (1 Hz) sinusoidal wave arriving at YKA with a slowness 7.5 s/$^{\circ}$ and backazimuth 225$^{\circ}$.
d) Same as a) but for a P arrival of real data (event: 04-jun-1993 10:49) at YKA. Theoretical slowness and backazimuth values are marked by the white diamonds. The lower slowness marks the P slowness (u$_P$ = 4.43). The other diamond marks the PP slowness (u$_{PP}$ = 7.65). The theoretical bazimuth is 295$^{\circ}$.

For coherent waves with 1 s period arriving from different backazimuths and with slownesses different from u = 0 s/$^{\circ}$, the maximum of the power spectral density is shifted to the corresponding values, without changing the form of the ARF. This is shown in Figure 5.2 c) for a signal arriving with a slowness of u = 7.5 s/$^{\circ}$ and a backazimuth of 225$^{\circ}$. For simplicity, less isolines are displayed in this figure than in Figure 5.2 a). For a signal with different frequency f$_{1}$ $\neq$ f$_{0}$ = 1 Hz the fk-diagram is scaled by the factor f$_{0}$/f$_{1}$ without a change of the cross shaped form of the ARF.
Real seismograms do not consist of monochromatic waves, but of wavelets showing a broad frequency spectrum. For these signals the fk-diagram is computed as an integration over all frequencies in a chosen frequency band. The cut-off frequencies for the fk-analysis coincide with the edges of the used band-pass filter discussed in chapter 5.1.
The fk-diagram of a real P arrival recorded at YKA (13-jul-1994 10:49) is displayed in Figure 5.2 d). The phase shows a slowness comparable to the theoretical slowness. The theoretical slowness and backazimuth for P (u = 4.43 s/$^{\circ}$, $\Theta$ = 295.5$^{\circ}$) and PP (u = 7.65 s/$^{\circ}$, $\Theta$ = 295.5$^{\circ}$) are marked by the white diamonds. The deviation of the backazimuth from $\Theta_{theo}$ is within the azimuthal resolution of the fk-analysis at YKA. The form of the fk-diagram is similar to the ARF. This resemblance indicates a large coherency of the phase in the fk time window. The additional information obtained by the fk-analysis is used to identify coherent arrivals and to distinguish them from incoherent noise in the P$^d$P study.