Seismic array

With a seismic array, the signal-to-noise ratio (SNR) of a seismic signal can be improved by summing the coherent signals from the individual array sites. The most important point during the beamforming process is to find the best delay times by which the single traces must be shifted before summation in order to get the largest amplitudes due to coherent interference of the signals.

A wavefront coming from north-east and crossing a seismic array[1]

For distances from the source much larger than about 10 wavelengths, a seismic wave approaches an array as a wavefront that is close to planar. The directions of approach and propagation of the wavefront projected onto the horizontal plane are defined by the angles Φ and Θ.

• Φ Backazimuth (BAZ) = angle of wavefront approach, measured clockwise from the North to the direction towards the epicenter in degree.
• Θ Direction in which the wavefront propagates, measured in degree from the North, with Θ = Φ ±180°.
• d_j Horizontal distances between array site j and center site in [km].
• s Slowness vector with absolute value s = 1/ v_app
• v_app Apparent velocity vector with the absolute value v_app = 1/s . v_app = (v_app,x ,v_app,y ,v_app,z), where v_app,x ,v_app,y ,v_app,z are the single apparent velocity components in [km/s] of the wavefront crossing an array.
• v_app,h Absolute value of the horizontal component of the apparent velocity. ${\displaystyle v_{app,h}={\sqrt {v_{app,x}^{2}+v_{app,y}^{2}}}}$

In most cases, the elevation differences between single array sites are so small that travel-time differences due to elevation differences are negligible. In this case, we cannot measure the vertical component of the wavefront propagation. The time delay τ_j between the center site 0 and site j with the relative coordinates (x_j, y_j) is

${\displaystyle \tau _{j}={\frac {d_{j}}{v_{app,h}}}={\frac {-x_{j}sin\Phi -y_{j}cos\Phi }{v_{app,h}}}}$

In some cases, not all array sites are located on one horizontal plane. The time delays τ_j also depends on the local crustal velocities (v_c) below the given site j. The calculation of τ_j with coordinates (x_j, y_j, z_j) is

${\displaystyle \tau _{j}={\frac {-x_{j}sin\Phi -y_{j}cos\Phi }{v_{app,h}}}+{\frac {z_{j}cos\Phi }{v_{c}}}}$

In both the calculation can be written in vector syntax with position vector ${\displaystyle r_{j}}$ and slowness vector ${\displaystyle s_{j}}$:

${\displaystyle \tau _{j}=r_{j}\cdot s_{j}}$

Let w_j(t) be the digital sample of the seismometer from site j at time t, then the beam of the whole array is defined as

${\displaystyle b(t)={\frac {1}{M}}\sum _{j=1}^{M}w_{j}(t+\tau _{j})}$

If seismic waves are harmonic waves S(t) without noise, with identical site responses, and without attenuation, then the above operation would reproduce the signal S(t) accurately. Real data w(t) are the sum of background noise n(t) plus the signal of interest S(t), i.e. w(t) = S(t) + n(t). Assuming that the signal is coherent and not attenuated, calculating the sum of M observations and including noise we get

${\displaystyle B(t)=M\cdot S(t)+\sum _{j=1}^{M}n_{j}(t+\tau _{j})}$

Assuming that the noise n_j(t) has a normal amplitude distribution with zero mean and variance σ² at all sites, then the variance of the noise after summation is ${\displaystyle \sigma _{s}^{2}=M\sigma ^{2}}$ and the standard deviation is ${\displaystyle {\sqrt {M}}\sigma ^{2}}$. That means the standard deviation of the noise is multiplied by ${\displaystyle {\sqrt {M}}}$ while the coherent signal is multiplied by ${\displaystyle M}$. The theoretical improvement of the SNR by beamforming (aka array gain) will be ${\displaystyle G={\sqrt {M}}}$ for an array containing M sites.[1]

N-th root process is a non-linear method to enhance the SNR during beamforming. Before summing up the single seismic traces, the N-th root is calculated for each trace retaining the sign information. signum{w_j(t)} is a function defined as -1 or +1, depending on the sign of the actual sample w_j(t). N is an integer that has to be chosen by the analyst

${\displaystyle B_{N}(t)=\sum _{j=1}^{M}{\sqrt[{N}]{n_{j}(t+\tau _{j})}}\cdot signum\{w_{j}(t)\}}$

Here the value of the function ${\displaystyle signum\{w_{j}(t)\}}$ is defined as ±1 depending on the sign of the actual sample w_j(t). After this summation, the beam has to be raised to the power of N

${\displaystyle b(t)=|B_{N}(t)|^{N}\cdot signum\{w_{j}(t)\}}$

The N-th root process was first proposed by K. J. Muirhead and Ram Dattin in 1976.[2] With the N-th root process, the suppression of uncorrelated noise is better than with linear beamforming. However, it weighs the coherency of a signal higher than the amplitudes, which results in a distortion of the waveforms.

Schimmel and Paulssen introduced another non-linear stacking technique in 1997[3] to enhance signals through the reduction of incoherent noise, which shows a smaller waveform distortion than the N-th root process. Kennett proposed the use of the semblance of the signal as a weighting function in 2000[4] and achieved a similar resolution.

An easily implementable weighted stack method would be to weight the amplitudes of the single sites of an array with the SNR of the signal at this site before beamforming, but this does not directly exploit the coherency of the signals across the array. All weighted stack methods can increase the slowness resolution of velocity spectrum analysis.

A cluster of earthquakes can be used as a source array to analyze coherent signals in the seismic coda. This idea was consequently expanded by Krüger et al. in 1993 by analyzing seismic array data from well-known source locations with the so-called "double beam method".[5] The principle of reciprocity is used for source and receiver arrays to further increase the resolution and the SNR for small amplitude signals by combining both arrays in a single analysis.

The array transfer function describes sensitivity and resolution of an array for seismic signals with different frequency contents and slownesses. With an array, we are able to observe the wavenumber ${\displaystyle k=2\pi /\lambda =2\pi \cdot f\cdot s}$ of this wave defined by its frequency f and its slowness s. While time-domain analog-to-digital conversion may give aliasing effects in the time domain, the spatial sampling may give aliasing effects in the wavenumber domain. Thus the wavelength range of seismic signals and the sensitivity at different wavelengths must be estimated.[1]

The difference between a signal w at the reference site A and the signal w_n at any other sensor A_n is the travel time between the arrivals at the sensors. A plane wave is defined by its slowness vector s_o

${\displaystyle w_{n}(t)=w(t-r_{n}\cdot s_{0})}$, where ${\displaystyle r_{n}}$ is the position vector of site n

The best beam of an array with M sensors for a seismic signal for the slowness s_o is defined as

${\displaystyle b(t)={\frac {1}{M}}\sum _{j=1}^{M}w_{j}(t+r_{j}\cdot s_{0})}$

If we calculate all time shifts for a signal with the slowness s_o with respect to any other slowness s, the calculated beam becomes

${\displaystyle b(t)={\frac {1}{M}}\sum _{j=1}^{M}w_{j}(t+r_{j}\cdot (s_{0}-s))}$

The seismic energy of this beam can be calculated by integrating over the squared amplitudes

${\displaystyle E(t)=\int _{-\infty }^{\infty }b^{2}(t)dt=\int _{-\infty }^{\infty }[{\frac {1}{M}}\sum _{j=1}^{M}w_{j}(t+r_{j}\cdot (s_{0}-s))]^{2}dt}$

This equation can be written in the frequency domain with ${\displaystyle {\bar {w}}(\omega )}$ being the Fourier transform of the seismogram w(t), using the definition of the wavenumber vector k = ω⋅ s

${\displaystyle E(\omega ,k_{0}-k)={\frac {1}{2\pi }}\int _{-\infty }^{\infty }|{\bar {w}}(\omega )|^{2}\cdot |C(k_{0}-k)|^{2}d\omega }$, where ${\displaystyle C(k_{0}-k)={\frac {1}{M}}\sum _{j=1}^{M}e^{iwr_{j}(k_{0}-k)}}$

This equation is called the transfer function of an array. If the slowness difference is zero, the factor ${\displaystyle |C(k_{0}-k)|^{2}}$ becomes 1.0 and the array is optimally tuned for this slowness. All other energy propagating with a different slowness will be suppressed.[1]

Slowness estimation is a matter of forming beams with different slowness vectors and comparing the amplitudes or the power of the beams, and finding out the best beam by looking for the v_app and backazimuth combination with the highest energy on the beam.

Frequency-wavenumber analysis is used as a reference tool in array processing for estimating slowness. This method was proposed by Capon in 1969[6] and further developed to include wide-band analysis, maximum-likelihood estimation techniques, and three-component data in the 1980s.[7]

The methodology exploits the deterministic, non-periodic character of seismic wave propagation to calculate the frequency-wavenumber spectrum of the signals by applying the multidimensional Fourier transform. A monochromatic plane wave w(x,t) will propagate along the x direction according to equation

${\displaystyle w(x,t)=Ae^{i2\pi (f_{0}t-k_{0}x)}}$

It can be rewritten in frequency domain as

${\displaystyle W(k_{x},f)=A\delta (f-f_{0})\delta (k_{x}-k_{0})}$

which suggests the possibility to map a monochromatic plane wave in the frequency-wavenumber domain to a point with coordinates (f, k_x) = (f₀, k₀).

Practically, f-k analysis is performed in the frequency domain and represents in principle beamforming in the frequency domain for a number of different slowness values. At NORSAR slowness values between -0.4 and 0.4 s/km are used equally spaced over 51 by 51 points. For every one of these points the beam power is evaluated, giving an equally spaced grid of 2601 points with power information.[8]

A beampacking scheme was developed at NORSAR to apply f-k analysis of regional phases to data of large array.[8] This algorithm performs time-domain beamforming over a predefined grid of slowness points and measures the power of the beam.

In practice the beampacking process gives the same slowness estimate as for the f-k analysis in the frequency domain. Compared to the f-k process, the beampacking process results in a slightly (about 10%) narrower peak for the maximum power.

Another way of estimating slowness is to pick carefully times of the first onset or any other common distinguishable part of the same phase (same cycle) for all instruments in an array.[1] Let t_i be the arrival time picked at site i, and t_ref be the arrival time at the reference site, then τ_i = t_i − t_ref is the observed time delay at site i. We observe the plane wave at M sites. With M ≥ 3. The horizontal components (s_x, s_y) of the slowness vector s can be estimated by

${\displaystyle {\hat {s}}={\underset {s}{min}}\sum _{j=1}^{M}(\tau _{j}-r_{j}\cdot s)^{2}}$

Plane wave fitting requires interactive analyst's work. However, to obtain automatic time picks and thereby provide a slowness estimate automatically, techniques like cross-correlation or just picking of peak amplitude within a time window may be used.[9] Because of the amount of required computations, plane wave fitting is most effective for arrays with a smaller number of sites or for subarray configurations.

The Gräfenberg array is the first digital broadband array that has a continuous data history from 1976 until today. This array consists of 13 broadband stations in the Fränkische Alb. It extends approximately 100 kilometers north-south and approximately 40 kilometers east-west.

YKA or Yellowknife Seismological Array is a medium size seismic array established near Yellowknife in the Northwest Territories, Canada, in 1962, in cooperative agreement between the Department of Mines and Technical Surveys (now Natural Resources Canada) and the United Kingdom Atomic Energy Authority (UKAEA), to investigate the feasibility of teleseismic detection and identification of nuclear explosions. YKA currently consists of 19 short period seismic sensors in the form of a cross with an aperture of 2.5 km, plus 4 broadband seismograph sites with instruments able to detect a wide range of seismic wave frequencies.[10]

Configuration of large aperture array NORSAR and small aperture array NORES.[8]

LASA or Large Aperture Seismic Array is the first large seismic array. It was built in Montana, USA, in 1965.[11]

NORSAR or Norwegian Seismic Array was established at Kjeller, Norway in 1968 as part of the Norwegian-US agreement for the detection of earthquakes and nuclear explosions. It has been an independent, not-for-profit, research foundation within the field of geo-science since 1999. NORSAR was constructed as a large aperture array with a diameter of 100 km. It is the largest stand-alone array in the world.[8]

NORES was the first regional seismic array constructed in southern Norway in 1984. A sister array ARCES was established in northern Norway in 1987. NORES and ARCES are small aperture arrays with a diameter of only 3 km.[8]

GERES is a small aperture array built in the Bavarian Forest near the border triangle of Germany, Austria and Czech, in 1988. It consists of 25 individual seismic stations arranged in 4 concentric rings with radius of 200m, 430m, 925m and 1988m.[12]

SPITS is a very small aperture array at Spitsbergen, Norway. It was originally installed in 1992 and upgraded to IMS standard in 2007 by NORSAR.[13]