Background

Passive dense array surface wave tomography plays a significant role in probing the earth’s interior with the development of a large number of dense seismic arrays, attracting researchers’ growing interest in imaging applications in various scales. These array methods, like beamforming and MASW (multichannel analysis of surface wave) have been widely used to estimate dispersion curves in past several decades. A new array stacking technique, the frequency-Bessel (F-J) transform, was proposed about five years ago with some main advantages of extracting dispersion curves of high-resolution and broad frequency band (Wang et al., 2019). More importantly, this method facilitates the estimation of multimodal dispersion curves for investimating underground shear wave velocity structure (Wu et al., 2020; Li et al., 2021). Here we introduce the basic principle of this newly proposed array stacking method.

Principles of the F-J transform

The vertical-vertical (Z-Z) component cross-correlations in frequency domain $C_{ZZ}(\omega, r)$ from ambient seismic noise is approximately proportional to the imaginary part of the Green’s function $G_{ZZ}(\omega, r)$ between two receivers with distance of $r$, $$ C_{ZZ}(\omega,r) \approx A\cdot Im\{G_{ZZ}(\omega, r)\}, \tag{1} $$ in which, $\omega$ is angular frequency, and $Im\{\cdot\}$ represents imaginary part. Wang et al (2019) defined the spectrum $I(\omega, k)$ of the F-J transform as $$ I(\omega, k)=\int_0^{+\infty}C_{ZZ}(\omega, r)J_0(kr)rdr, \tag{2} $$ where, $J_0$ is the zeroth-order Bessel function of the first kind. For a flat multi-layered medium, $G_{ZZ}(\omega, r)$ due to an isotropic source is expressed as $$ G_{ZZ}(\omega, r)=\int_0^{+\infty}g_z(\omega, \kappa)J_0(\kappa r)\kappa d\kappa, \tag{3} $$ where, $g_z(\omega, k)$ is the kernel function of the medium, and it is inversely proportional to the surface wave dispersion equation, which means that $g_{z}(\omega, \kappa)$ reaches infinity at roots (dispersion curve) of dispersion equation. Now we use $(1)$, $(2)$ and $(3)$ to obtain $$ I(\omega, k)=\int_0^{+\infty}Im\{\int_0^{+\infty}g_z(\omega, \kappa)J_0(\kappa r) \kappa d\kappa\}J_0(kr)rdr. \tag{4} $$ Note the orthogonality of $J_0$, $\int_0^{+\infty}J_0(kr)rdr\int_0^{+\infty}J_0(\kappa r)\kappa d\kappa=\frac{\delta(k-\kappa)}{k}$, and thus we have $$ I(\omega, k) \approx A \cdot Im\{g_z(\omega, k)\}. \tag{5} $$ $(5)$ shows that $I(\omega, k)$ is infinity at the dispersion point. This is the basic idea that we can estimate dispersion curves utilizing the F-J transform.

Utilizing the CC-FJpy package to estimate dispersion curves

There is a convenient Python package named CC-FJpy for extracting dispersion curves developed by members from Xiaofei Chen’s team. We can find this package on this website (https://github.com/ColinLii/CC-FJpy). Here we show the dispersion diagram result computed using this package. Cross-corrrelation functions are from synthetic ambient noise based on the following source-receiver configuration. noise

We utilize 100 receivers and 2000 sources to synthesize seismic event. Synthetic ambient noise is generated according to steps described in this blog. Here we show the dispersion diagram using CC-FJpy.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69

import numpy as np
from obspy import read
import ccfj
import matplotlib.pyplot as plt
from glob import glob

# Match cross-correlation time function files
files = sorted(glob('ccf/*.SAC'))
# Frequency limitations of dispersion diagram
f1 = 0; f2 = 60
# Half time window width of cross-correlations
t0 = 10
t1 = -t0; t2 = t0
dist = []; data = []
ccfr = []; ccfi = []
index = 0
for sac in files:
    tr = read(sac)[0]
	# inster-receiver distance in km is stored in "dist" attribute
    dist.append(tr.stats.sac.dist)
    b = tr.stats.sac.b
    dt = tr.stats.delta
    n1 = int((t1-b)/dt); n2 = int((t2-b)/dt)
    d = tr.data[n1: n2+1]
    d = np.fft.fftshift(d)
    tmp = np.fft.fft(d)
    n = len(tmp)
    tmp = tmp[:n//2]
    ccfr.append(tmp.real)
    index += 1
ccfr = np.array(ccfr)
nd = len(d)
dist = np.array(dist) * 1e3
index = np.argsort(dist)

dist = dist[index]
ccfr = ccfr[index]
# Phase velocity testing range
c1 = 0.25; c2 = 0.75; nc = 501
c = np.linspace(c1*1e3, c2*1e3, nc)
f = np.arange(nd) / dt / (nd-1)
fn1 = int(f1*nd*dt); fn2 = int(f2*nd*dt)
nf = fn2 - fn1 + 1; nr = len(dist)

# Compute dispersion diagram using the F-J transform
im = ccfj.fj_noise(ccfr[:, fn1: fn2+1], dist, c,
                   f[fn1: fn2+1], fstride=1, itype=1,
                   func=0, num=1)
im = np.maximum(im, 0)
# Normalization along the frequancy axis
for i in range(len(im[0])):
    im[:, i] /= abs(im[:, i]).max()
f = f[fn1: fn2+1]
# Load theoretical dispersion curves
dc = np.loadtxt('syn-disp.txt')
dc = dc[::4]
# Visualize dispersion diagram
plt.figure(figsize=(12, 8))
plt.pcolormesh(f, c/1e3, im**2, cmap='jet')
plt.scatter(dc[:, 2], dc[:, 3], marker='o', s=0.5, facecolor='w', alpha=0.5)
plt.plot(f, f*dist.max()/1e3, 'w--')
plt.xlabel('Frequency (Hz)', fontsize=20)
plt.ylabel('Phase velocity (km/s)', fontsize=20)
plt.xlim(f1, f2)
plt.ylim(c1, c2)
plt.xticks(fontsize=17)
plt.yticks(fontsize=17)
plt.tight_layout()
plt.show()

disp White dots in the above figure are theoretical dispersion curves from the fundamental mode to the third higher mode.