Short background
Ambient noise tomography of regional and continental scales is much benefited from
the system composed of the atmosphere, ocean and the solid earth. Ambient noise below
1Hz is largely due to oceans, and interactions between ocean waves and/or the sea floor
generate pressure impulse propagating into the solid earth. Land seismic stations are seizing
seismic waves which originating from microseisms, and these seismic waves include surface waves
and body waves. Individuals can easily retrieve surface waves from ambient noise cross-correlations.
Body wave extraction from ambient noise, however, remains difficult. We often explore body waves from
ambient noise using array stacking techniques for array stacking leads to high signal-to-noise ratios.
Here we perform the body wave source location utilizing beamforming and back-projection methods,
and we process one-day seismic recordings by a seismic array in China, which reveals that the source location
is in the North Atlantic Ocean, the south of Iceland.
Data and Methods
We collect one-day continuous seismic data of vertical component(here BHZ, 2008-03-10) recorded by a seismic array nearly in
the center part of China (see the following figure), and compute the cross-spectral matrix using one-hour length
segment with a overlap of 50% (1800 seconds):
$$
C(\omega) = \frac{1}{N}\sum_{n=1}^N\boldsymbol{u}_n(\omega) \cdot \boldsymbol{u}_n^H(\omega), \tag{1}
$$
in which, $\boldsymbol{u}_n(\omega)=[u^1_n(\omega), u^2_n(\omega), \cdots, u^M_n]^T$ is the Fourier spectra
of the $n_{th}$ segment data recorded by $M$ stations, $H$ denotes Hermitian transpose, and $N$ is the total number of data segments.
Then we create steering vector as follows,
$$
\boldsymbol{a}(\omega; \boldsymbol{s}, \boldsymbol{x}) = [e^{-i\omega \boldsymbol{s} \cdot \boldsymbol{x}_1}, e^{-i\omega \boldsymbol{s} \cdot \boldsymbol{x}_2}, \cdots, e^{-i\omega \boldsymbol{s} \cdot \boldsymbol{x}_M}]^T, \tag{2}
$$
where, $\boldsymbol{s}$ and $\boldsymbol{x}_m$ separately indicate horizontal slowness (ray parameter for body waves) and position of the $m_{th}$ station.
The beam power is given by
$$
B(\omega; \boldsymbol{s}) = \boldsymbol{a}^H \boldsymbol{C}(\omega) \boldsymbol{a}. \tag{3}
$$
For more details about beamforming, someone can refer to this blog.
An Example of Locating the distant P wave source
We first import required Python libraries.
1
2
3
4
5
6
7
8
|
import numpy as np
import matplotlib.pyplot as plt
from obspy import read
from glob import glob
from mpl_toolkits.basemap import Basemap
from scipy.signal.windows import hann
import matplotlib as mpl
from obspy.taup import TauPyModel
|
The 13 seismic stations are installed in the center part of China, and we collect the data of March 10, 2008 to perform this work.
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
|
files = sorted(glob('2008-03-10/*BHZ.SAC'))
sta = {}
with open('X4_IRIS_BHZ.txt', 'r') as fin:
for line in fin.readlines():
tmp = line.strip().split('|')
sta[tmp[1]] = [float(tmp[3]), float(tmp[2])]
xy = []
for f in files:
tmp = f.split('/')[-1].split('.')[1]
xy.append(sta[tmp])
xy = np.array(xy)
plt.figure(figsize=(12, 4.5))
plt.subplot2grid((1, 3), (0, 0), rowspan=1, colspan=2)
bmap1 = Basemap(projection='cyl', lon_0=0)
bmap1.shadedrelief(scale=0.2)
xx, yy = bmap1(106, 34)
bmap1.scatter(xx, yy, marker='v', s=200, facecolor='#87CEFA', edgecolor='k', lw=1, alpha=1)
bmap1.drawmeridians(meridians=np.linspace(-180, 180, 7),
labels=[0, 0, 0, 1], color='none', fontsize=17)
bmap1.drawparallels(circles=np.linspace(-90, 90, 7),
labels=[1, 0, 0, 0], color='none', fontsize=17)
plt.subplot2grid((1, 3), (0, 2), rowspan=1, colspan=1)
bmap2 = Basemap(projection='merc', llcrnrlat=30, llcrnrlon=100,
urcrnrlon=110, urcrnrlat=39)
bmap2.shadedrelief(scale=1)
xx, yy = bmap2(xy[:, 0], xy[:, 1])
bmap2.drawcoastlines(color='gray')
bmap2.scatter(xx, yy, marker='v', s=150, facecolor='#87CEFA', edgecolor='k', lw=1, alpha=1)
bmap2.drawmeridians(meridians=np.linspace(100, 110, 6),
labels=[0, 0, 0, 1], color='none', fontsize=13)
bmap2.drawparallels(circles=np.linspace(30, 40, 6),
labels=[1, 0, 0, 0], color='none', fontsize=13)
plt.show()
|
It is worth noting that we compute the beam power in Cartesian coordinate system while we only have the longitudes and latitudes of these stations.
Thus, we have to convert these longitudes and latitudes to Cartesian coordinates, so we define the following useful functions.
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
|
def distance(la1, lo1, LA2, LO2, R=6371):
d2r = np.pi / 180
dla = (LA2-la1) * d2r
dlo = (LO2-lo1) * d2r
a = np.sin(dla/2) ** 2 + np.cos(la1*d2r) * np.cos(LA2*d2r) * np.sin(dlo/2) ** 2
c = 2 * np.arctan2(a**0.5, (1-a)**0.5)
return R * c * 1e3
def azimuth(la1, lo1, LA2, LO2):
d2r = np.pi / 180
dla = (LA2 - la1) * d2r
dlo = (LO2 - lo1) * d2r
y = np.sin(dlo) * np.cos(LA2*d2r)
x = np.cos(la1*d2r) * np.sin(LA2*d2r) - np.sin(la1*d2r) * np.cos(LA2*d2r) * np.cos(dlo)
th = np.arctan2(y, x)
return (th / d2r + 360) % 360
def lola2xy(xy):
n = xy.shape[0]
lo0 = xy[:, 0].mean()
la0 = xy[:, 1].mean()
dist = distance(la0, lo0, xy[:, 1], xy[:, 0])
az = azimuth(la0, lo0, xy[:, 1], xy[:, 0]) * np.pi / 180
new_xy = np.zeros_like(xy)
new_xy[:, 0] = np.sin(az) * dist
new_xy[:, 1] = np.cos(az) * dist
return new_xy
def process_data(d, operator, p=0.01):
nd = len(d)
pn = int(nd*p)
h = hann(2*pn)
d -= d.mean()
d[:pn] *= h[:pn]
d[-pn:] *= h[-pn]
fd = np.fft.fft(d)
fd /= np.convolve(np.abs(fd), operator, 'same')
fd[np.isnan(fd)] = 0.
return fd
|
We need to estimate the average cross-spectral matrix of seismic data by this array before compute the beam power.
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
|
st = read('2008-03-10/*BHZ.SAC')
utcb = st[0].stats.starttime
utce = st[0].stats.endtime
dt = st[0].stats.delta; fs = 1 / dt
xy = []
for tr in st:
if utcb > tr.stats.starttime:
utcb = tr.stats.starttime
if utce < tr.stats.endtime:
utce = tr.stats.endtime
xy.append(sta[tr.stats.station])
xy = np.array(xy)
xy = lola2xy(xy)
n = xy.shape[0]
nd = int((utce-utcb)/dt)
data = np.zeros((n, nd+50))
for ti, tr in enumerate(st):
n0 = int((tr.stats.starttime-utcb)/dt)
data[ti, n0: n0+tr.stats.npts] = tr.data
pt = 3600; pn = int(pt/dt); on = pn // 2
hn = 15; operator = np.ones(hn*2+1) / (hn*2+1)
tmp_data = []
n1 = 0
while n1 < data.shape[1]-pn:
tmp = []
n2 = n1 + pn
for i in range(n):
tmp.append(process_data(data[i, n1:n2], operator))
tmp_data.append(tmp)
n1 += on
tmp_data = np.array(tmp_data)
f1 = 1 / 9; f2 = 1 / 5
fn1 = int(f1*pn/fs); fn2 = int(f2*pn/fs); nf = fn2 - fn1 +1
c = np.zeros((nf, n, n), dtype=complex)
u = np.zeros((n, 1), dtype=complex)
for i in range(fn1, fn2+1):
for j in range(tmp_data.shape[0]):
u[:, 0] = tmp_data[j, :, i]
c[i-fn1] += u @ np.conjugate(u.T)
|
Now we compute the beam power using formula $(1)$ to estimate the back-azmith and corresponding ray parameter of the potential body wave.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
|
s1 = 0; s2 = 10; ns = 101
g = 111194.92664455
s = np.linspace(s1, s2, ns) / g
b1 = 0; b2 = 360; nb = 181
b = np.linspace(b1, b2, nb) / 180 * np.pi
p1 = np.zeros((ns, nb), dtype=complex)
a = np.zeros((n, 1), dtype=complex)
for i in range(ns):
for j in range(nb):
shift = -s[i] * (np.sin(b[j])*xy[:, 0] + np.cos(b[j])*xy[:, 1])
for k in range(fn1, fn2+1):
w = 2 * np.pi * k * fs / pn
a[:, 0] = np.exp(-1j*w*shift)
cc = c[k-fn1].copy(); cc /= np.abs(cc); cc[np.isnan(cc)] = 0.
p1[i, j] += np.conjugate(a.T) @ cc @ a
p1 = np.abs(p1) / nf / n**2
|
Someone may find that we normalize the beam power using the number of frequency points and the number of stations, which results in
the coherence of signals caught by this array.
The maximum of backazimuth-slowness beam coherence panel suggests that the signal source comes from the direction ranging 315-360 degress with horizonral slowness ranging 4-6 sec/degree, and thus we identify this seismic phase as P wave. We calculate the ray parameter of P wave using the taup module, and the slowness-distance curve is shown by the following figure.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
|
model = TauPyModel(model='iasp91')
dd = 0.02
d = np.arange(0, 95+dd, dd)
pw = np.zeros((len(d), 2))
phase = ['P']
for i in range(len(d)):
ar = model.get_travel_times(source_depth_in_km=0,
distance_in_degree=d[i],
phase_list=phase)
pw[i] = np.array([d[i], ar[0].ray_param])
pw[:, 1] = pw[:, 1] * np.pi / 180
plt.figure(figsize=(8, 4))
plt.plot(pw[:, 0], pw[:, 1], lw=2.5, color='#00BFFF')
plt.xlabel('Distance (degree)', fontsize=18)
plt.ylabel('Ray parameter (sec/degree)', fontsize=18)
plt.gca().tick_params(labelsize=15)
plt.xlim(30, 95)
plt.ylim(4, 10)
plt.grid(ls=(10, (6, 5)))
plt.title('Ray parameter of seismic phase P', fontsize=17)
plt.show()
|
According to the beam coherence and the slowness-distance curve of P wave, we estimate that the distance between the source location the seismic array center is about 70 degrees. We utilize grid-searching rule to find the source location by projecting the backazimuth-slowness beam coherence to geographical coordinates.
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
|
def distance_matrix(la1, lo1, LA2, LO2):
d2r = np.pi / 180
dla = (LA2-la1) * d2r
dlo = (LO2-lo1) * d2r
a = np.sin(dla/2) ** 2 + np.cos(la1*d2r) * np.cos(LA2*d2r) * np.sin(dlo/2) ** 2
c = 2 * np.arctan2(a**0.5, (1-a)**0.5)
return c * 180 / np.pi
def azimuth_matrix(la1, lo1, LA2, LO2):
d2r = np.pi / 180
dla = (LA2 - la1) * d2r
dlo = (LO2 - lo1) * d2r
y = np.sin(dlo) * np.cos(LA2*d2r)
x = np.cos(la1*d2r) * np.sin(LA2*d2r) - np.sin(la1*d2r) * np.cos(LA2*d2r) * np.cos(dlo)
th = np.arctan2(y, x)
return (th / d2r + 360) % 360
xy = []
for tr in st:
xy.append(sta[tr.stats.station])
xy = np.array(xy)
lo0 = xy[:, 0].mean()
la0 = xy[:, 1].mean()
ns, nb = p1.shape
ss = np.linspace(0, 10, ns)
bb = np.linspace(0, 360, nb)
ds = ss[1] - ss[0]
db = bb[1] - bb[0]
lo1 = -180; lo2 = 180; nlo = 181
la1 = -90; la2 = 90; nla = 91
lo = np.linspace(lo1, lo2, nlo)
la = np.linspace(la1, la2, nla)
LO, LA = np.meshgrid(lo, la)
AZ = azimuth_matrix(la0, lo0, LA, LO)
DIST = distance_matrix(la0, lo0, LA, LO)
bp = np.zeros((nla, nlo))
for i in range(nla):
for j in range(nlo):
index = np.argmin(np.abs(DIST[i, j]-pw[:, 0]))
if np.abs(DIST[i, j]-pw[index, 0]) > 0.05 or pw[index, 1] > ss.max():
continue
si = int((pw[index, 1]-ss.min())/ds)
ai = int((AZ[i, j]-bb.min())/db)
bp[i, j] = p1[si, ai]
plt.figure(figsize=(12, 6))
bmap = Basemap(projection='cyl', lon_0=0)
plt.contourf(LO, LA, bp, cmap=cmap, levels=101)
cbar = plt.colorbar(shrink=0.8)
cbar.ax.tick_params(labelsize=15)
cbar.set_label('Cohenrence', fontsize=15)
bmap.fillcontinents(color='#3CB371')
bmap.drawcoastlines(linewidth=0.3, color='gray')
xx, yy = bmap(lo0, la0)
bmap.scatter(xx, yy, marker='v', s=150, facecolor='b', alpha=0.6, zorder=2)
bmap.drawparallels(circles=np.linspace(-90, 90, 7), labels=[1, 0, 0, 0], fontsize=17, color='none')
bmap.drawmeridians(meridians=np.linspace(-180, 180, 9), labels=[0, 0, 0, 1], fontsize=17, color='none')
plt.show()
|
This back-projection result shows that the possible source location is in the North Atlantic Ocean, which may relate to the irregular topography of the middle ocean rige in this place.
Perspective
We use this simple example to explain the process of using body waves from ambient noise cross-correlations to investigate the potential source locations of these indentified seismic phases. Indeed, it still holds many mysteries to reveal in the seismic ambient noise world. Thus, we will continue the journey of exploring the coupling process among the atmosphere, oceans and our solid planet Earth.
