## 1. Background

Ambient noise tomography has been a very popular imaging strategy of investigating the earth’s interior in past over two decades. The basic principle of this method is that stacked cross-correlations from noise recordings made by two receivers can be approximately the Green’s function from one receiver to another. We have the need of synthesizing ambient noise for validating this important idea, or exploring characteristics of ambient noise. Here we introduce a simple way of computing synthetic ambient noise.

## 2. Steps to doing synthetics

### 2.1 Configuring the receiver and source positions

Firstly we randomly generate source and receiver positions with given coordinate ranges using the following python codes:

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15  import numpy as np import matplotlib.pyplot as plt nr = 60 r = 250 th = np.linspace(0, 2*np.pi, nr) ra = np.linspace(0, r, nr) xy = np.zeros((nr, 2)) xy[:, 0] = ra * np.cos(th) xy[:, 1] = ra * np.sin(th) xy[:, 0] += (np.random.randn(nr) * 100) xy[:, 1] += (np.random.randn(nr) * 100) plt.figure(figsize=(6, 6)) plt.plot(xy[:, 0], xy[:, 1], 'bv', alpha=0.5) plt.show() np.savetxt('rec.txt', xy) 

We generated 60 receivers within the circle with radius of ~250 km, and coordinates are saved in rec.txt. Now we generate 2001 sources within a ring belt with radius limitations of 600 km and 1200km, and three-component forces. Note that source positions are stored in file src.txt, and force components are saved in file force.txt.

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18  import numpy as np import matplotlib.pyplot as plt rxy = np.loadtxt('rec.txt') r1 = 600; r2 = 1200 d1 = 1e-3; d2 = 1e-2 ns = 2001 th = np.linspace(0, 2*np.pi, ns) r = np.random.random(ns) dep = np.random.random(ns) r /= np.abs(r).max() r = r1 + r*(r2-r1) dep = dep / np.abs(dep).max() dep += 1 dep /= np.abs(dep).max() dep = d1 + dep*(d2-d1) force = np.random.random((ns, 3)) * 1.0e26 np.savetxt('force.txt', force) 

The receiver and source configuration 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 23 24 25 26 27 28 29 30 31 32 33 34  import numpy as np import matplotlib.pyplot as plt from obspy import read xy = np.loadtxt('rec.txt') sr = np.loadtxt('src.txt') f = np.loadtxt('force.txt') plt.figure(figsize=(14, 5.5)) ax1 = plt.subplot(121) ax1.scatter(xy[:, 0], xy[:, 1], s=50, marker='v', facecolor='#AAAAAA', edgecolor='k', lw=1, alpha=0.5) plt.scatter(sr[:, 0], sr[:, 1], s=20, c=sr[:, 2]*1e3, marker='o', cmap='coolwarm', facecolor='#AAAAAA', edgecolor='k', lw=0.5, alpha=0.75) cbar = plt.colorbar() cbar.ax.tick_params(labelsize=11) cbar.set_label('Depth (m)', fontsize=13) ax1.tick_params(labelsize=13) ax1.set_xlabel('Easting (km)', fontsize=15) ax1.set_ylabel('Northing (km)', fontsize=15) ax2 = plt.subplot(3, 2, 2) ax2.plot(f[:, 0], 'r.', alpha=0.50) ax2.set_ylabel('Fx (dyne)', fontsize=15) ax2.tick_params(labelsize=13) ax2.set_xticks([]) ax3 = plt.subplot(3, 2, 4) ax3.plot(f[:, 1], 'g.', alpha=0.30) ax3.set_ylabel('Fy (dyne)', fontsize=15) ax3.tick_params(labelsize=13) ax3.set_xticks([]) ax4 = plt.subplot(3, 2, 6) ax4.plot(f[:, 2], 'b.', alpha=0.15) ax4.set_xlabel('Source index', fontsize=15) ax4.set_ylabel('Fz (dyne)', fontsize=15) plt.show() 

### 2.2 Synthetic ambient noise

The ambient noise recordings can be considered as the wavefield excited by many events, like the primary and secondary microseisms from interactions between ocean waves and our solid earth, and maybe some weak seismic signals generated by small earthquakes. These events randomly occur in time and space. Therefore, we compute seismic waveforms triggered by many earthquakes, and randomly stack these earthquake-generated waveforms according to event occuring time and their positions. Finally we synthesize ambient seismic noise. The python codes of computing waveforms excited by earthquakes:

  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  import numpy as np import os # read in source and receiver position information. src = np.loadtxt('src.txt') rec = np.loadtxt('rec.txt') nr = rec.shape[0] evt_dir = 'event/' if not os.path.exists(evt_dir): os.mkdir(evt_dir) for i in range(nr): sta = evt_dir + 'station_%02d'%(i+1) if not os.path.exists(sta): os.mkdir(sta) nmode = 3 ns = src.shape[0] # read in three-component forces. force = np.loadtxt('force.txt') fx = force[:, 0] fy = force[:, 1] fz = force[:, 2] # compute synthetic waveforms for i in range(ns): sx, sy, dep = src[i] tag_s = i + sn1 + 1 for j in range(nr): dx = sx - rec[j, 0]; dy = sy - rec[j, 1] dist = ( dx**2 + dy**2 ) ** 0.5 az = np.pi/2 - np.arctan2(dy, dx) az = az * 180 / np.pi dfile = '%g 0.25 2048 0 0'%dist with open('dfile', 'w') as fout: fout.write(dfile) os.system('sprep96 -M model.md -NMOD %d -HS %g -HR 0 -d dfile -R'%(nmode, dep)) os.system('sdisp96') os.system('sregn96') os.system('spulse96 -d dfile -D -i -ALL -2 > file96') os.system('fmech96 -A %g -fx %f -fy %f -fz %f -ROT < file96 > tosac.file'%(az, fx[i], fy[i], fz[i])) os.system('f96tosac -B tosac.file') os.system('mv B00101Z00.sac %s/station_%02d/%04d.Z.%02d.SAC'%(evt_dir, j+1, tag_s, j+1)) os.system('rm *R00.sac') os.system('rm file96 sdisp96.dat sregn96.egn tosac.file sdisp96.ray') 

Now we synthesize ambient noise with the following codes:

  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  import numpy as np from obspy import read from obspy.core import Trace, AttribDict from glob import glob import os evt_dir = 'event/' noise_dir = 'noise/' if not os.path.exists(noise_dir): os.mkdir(noise_dir) print('Read in event data ...') st = [] nr = 60 ns = 2001 for i in range(1, nr+1): st.append(read('%sstation_%02d/*.Z.*.SAC'%(evt_dir, i))) print(len(st), len(st[0])) print('Synthesize noise data ...') dt = st[0][0].stats.delta ns = len(st[0]) nd = len(st[0][0].data) n = nd * 500 data = np.zeros((nr, n)) index = np.random.random(ns*50) index /= index.max() for i in range(len(index)): si = int(index[i]*ns) if si >= ns: si = ns - 1 ti = int(np.random.random()*n) if ti > n - 1 - nd: ti = n - 1 - nd for j in range(nr): data[j, ti: ti+nd] += st[j][si].data for i in range(nr): data[i] /= np.abs(data[i]).max() print('Write out noise data ...') t = np.arange(n) * dt sac = Trace(t) sac.stats.sac = AttribDict({}) sac.stats.delta = dt sac.stats.sac.b = 0 for i in range(nr): sac.data = data[i] sac.write('%snoise.Z.%02d.SAC'%(noise_dir, i+1), format='SAC') 

We randomly choose synthetic ambient noise waveforms from 5 receivers to visualize:

  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  import numpy as np import matplotlib.pyplot as plt from obspy import read nr = 60 n = 5 index = [] for i in range(200): if len(index) >= n: break ri = int(np.random.random() * nr) if ri not in index: index.append(ri) xy = np.loadtxt('rec.txt') plt.figure(figsize=(12, 5)) plt.subplot(121) plt.scatter(xy[:, 0], xy[:, 1], s=100, marker='v', facecolor='#AAAAAA', edgecolor='k', lw=1, alpha=0.5) plt.scatter(xy[index, 0], xy[index, 1], s=100, marker='v', facecolor='r', edgecolor='k', lw=1, alpha=0.5) plt.xlabel('Easting (km)') plt.ylabel('Northing (km)') plt.subplot(122) for i in range(n): tr = read('noise/noise.Z.%02d.SAC'%(index[i]+1))[0] dt = tr.stats.delta data = tr.data t = np.arange(len(data)) * dt plt.plot(t, data*0.6+i, lw=0.5, color='k') plt.xlabel('Time (s)') plt.show() 

## 3. Checking the relation between receiver-pair cross-correlation and the green’s function

Ideal source distribution (energy equipartion in all directions) can lead to the phase smilarity between stacked cross-correlations and the Green’s functions. For detail, the relation is given by $$G_{AB}(t) = -\frac{dC_{AB}(t)}{dt}=\frac{dC_{AB}(-t)}{dt}, \tag{1}$$ where, $G_{AB}(t)$ and $C_{AB}(t)$ are separately the Green’s function and stacked cross-correlation between receiver A and B. Now we compute the temporally average cross-correlation function between receivers.

  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  import numpy as np from obspy import read from obspy.core import Trace, AttribDict import os rxy = np.loadtxt('rec.txt') st = read('noise/noise*.SAC') nr = len(st) dt = st[0].stats.delta # ccf time window length and movement window length segt = 2000; pt = 1000 segn = int(segt/dt) pn = int(pt/dt) data = np.zeros((nr, st[0].stats.npts)) print(st[0].stats.npts, data.shape) for i in range(nr): data[i] = st[i].data ccf_dir = 'ncf_%s/'%(noise.split('_')[-1][:-1]) if not os.path.exists(ccf_dir): os.mkdir(ccf_dir) nfft = 2 * segn - 1 operator = np.ones(5) / 5 cfd = np.zeros(nfft, dtype=complex) ct = np.arange(nfft) * dt - nfft//2 * dt cor = Trace(ct) cor.stats.delta = dt cor.stats.sac = AttribDict({}) cor.stats.sac.b = ct[0] cor.stats.sac.lcalda = True # frequency band of ccf f1 = 0.01; f2 = 0.6 for i in range(nr-1): for j in range(i+1, nr): ccf = ccf_dir + 'ccf_%02d_%02d.SAC'%(i+1, j+1) print(ccf) dist = ( (rxy[i][0]-rxy[j][0])**2 + (rxy[i][1]-rxy[j][1])**2 ) ** 0.5 cor.stats.sac.dist = dist cfd = np.zeros(nfft, dtype=complex) nn = 0 while nn < st[0].stats.npts - segn: dd1 = data[i, nn: nn+segn].copy() dd2 = data[j, nn: nn+segn].copy() dd1 -= dd1.mean(); dd2 -= dd2.mean() fd1 = np.fft.fft(dd1, nfft) fd2 = np.fft.fft(dd2, nfft) fd1 /= np.convolve(np.abs(fd1), operator, 'same') fd2 /= np.convolve(np.abs(fd2), operator, 'same') fd1[np.isnan(np.abs(fd1))] = 0. fd2[np.isnan(np.abs(fd2))] = 0. cfd += (fd1 * np.conjugate(fd2)) nn += pn cor.data = np.fft.ifftshift(np.fft.ifft(cfd)).real cor.filter('bandpass', freqmin=f1, freqmax=f2, corners=4, zerophase=True) cor.write(ccf, format='SAC') 

Here we compare the empirical green’s function from cross-corrlation and corresponding Green’s function with distance of ~ 500 km on different frequency bands.

  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  import numpy as np import matplotlib.pyplot as plt from obspy import read tr1 = read('ncf_01/ccf_14_41.SAC')[0] tr1.differentiate() tr1.data *= (-1) tr2 = read('B00102R00.sac')[0] T = [[2, 5], [5, 10], [10, 20], [15, 50], [20, 100]] plt.figure(figsize=(10, 8)) for i in range(len(T)): T1 = T[i][0]; T2 = T[i][1] f1 = 1 / T2; f2 = 1 / T1 trc1 = tr1.copy() trc1.filter('bandpass', freqmin=f1, freqmax=f2, corners=4, zerophase=True) trc2 = tr2.copy() trc2.filter('bandpass', freqmin=f1, freqmax=f2, corners=4, zerophase=True) d1 = trc1.data; d1 /= np.abs(d1).max() d2 = trc2.data; d2 /= np.abs(d2).max() tt1 = np.arange(len(d1)) * tr1.stats.delta + tr1.stats.sac.b tt2 = np.arange(len(d2)) * tr2.stats.delta + tr2.stats.sac.b plt.subplot(len(T), 1, i+1) plt.plot(tt1, d1, 'r', lw=2, alpha=0.5, label='EGF') plt.plot(tt2, d2, 'b', lw=2, alpha=0.35, label='GF') plt.xlim(0, 300) if i < 4: plt.xticks([]) if i < 1: plt.title('Dist.: %.3f km Period: %d-%d s'%(tr1.stats.sac.dist, T1, T2), fontsize=13) plt.legend(fontsize=12, loc='upper right') else: plt.title('Period: %d-%d s'%(T1, T2), fontsize=13) plt.gca().tick_params(labelsize=12) plt.xlabel('Time (s)', fontsize=15) plt.show() 

The cross-correlation and Green’s function waveforms are very similar on different period ranges. Furthermore, we use the image transform technique (Yao et al., 2006) to estimate dispersion curve from the cross-correlation function, and compare it with the theoretical one.

  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 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101  import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import CubicSpline from scipy.signal import hilbert, hann from obspy import read tr = read('ncf_01/ccf_14_41.SAC')[0] d0 = tr.data.copy() tr.differentiate() tr.data *= (-1) T1 = 2 T2 = 120 gv1 = 2 gv2 = 5 dt = tr.stats.delta b = tr.stats.sac.b gcarc = tr.stats.sac.dist t1 = gcarc / gv2 t2 = gcarc / gv1 n1 = int((t1-b)/dt); n2 = int((t2-b)/dt) nt = n2 - n1 + 1 tr.data[:n1] = 0 tr.data[n2:] = 0 d = tr.data[n1: n2+1] t = np.arange(nt) * dt + t1 dT0 = 0.5 T = np.arange(T1, T2+dT0, dT0) nT = len(T) dT = 0.25 ng = 251 dv = 0.2 p = np.zeros((ng, nT)) vph = np.linspace(gv1, gv2, ng) dn = int(dv/(vph[1]-vph[0])) h = hann(2*dn) for i in range(nT): dTT = dT f1 = 1 / (T[i]+dTT) f2 = 1 / (T[i]-dTT) trc = tr.copy() trc.filter('bandpass', freqmin=f1, freqmax=f2, corners=4, zerophase=True) d = trc.data[n1: n2+1] d /= abs(d).max() tmp = gcarc/(t-T[i]/8) index = np.argsort(tmp) natural = CubicSpline(tmp[index], d[index], bc_type='natural') p[:, i] = natural(vph, nu=0) p[:dn, i] *= h[:dn]; p[-dn:, i] *= h[dn:] p[:, i] /= np.abs(p[:, i]).max() p[:dn, i] *= h[:dn]; p[-dn:, i] *= h[dn:] d0 /= np.abs(d0).max() d1 = tr.data; d1 /= np.abs(d1).max() fig = plt.figure(figsize=(15, 6)) ax1 = plt.subplot2grid((1, 4), (0, 0), rowspan=1, colspan=1) ax1.plot(d0[n1: n2+1], t, lw=1.25, color='b', alpha=0.5, label='CCF') ax1.plot(tr.data[n1: n2+1], t, lw=1.25, color='r', alpha=0.5, label='EGF') ax1.legend() ax1.set_ylabel('Lag time (s)') ax1.set_title('Inter-sta dist: %.2f km'%gcarc) ax1.set_ylim(t[0], t[-1]) ax1.invert_yaxis() ax2 = plt.subplot2grid((1, 4), (0, 1), rowspan=1, colspan=3) ax2.pcolormesh(T, vph, p, cmap='RdYlBu_r') ax2.plot(T, gcarc/T, 'k--', lw=2, label='One $\lambda$') ax2.plot(T, gcarc/T/2, 'g--', lw=2, label='Two $\lambda$') ax2.plot(T, gcarc/T/3, 'y--', lw=2, label='Three $\lambda$') ax2.set_xlabel('Period (s)') ax2.set_ylabel('Vph (km/s)') ax2.set_xlim(T1, T2) ax2.set_ylim(gv1, gv2) ax2.yaxis.set_label_position('right') ax2.yaxis.set_ticks_position('right') dv = 0.075; r = 1.0 tmp = [5.7, 3.253] if len(tmp) < 1: plt.close() sys.exit(1) plt.close() t0, v0 = tmp print(t0, v0) tn0 = int((t0-T.min())/(T[1]-T[0])+0.5) vn0 = int((v0-vph.min())/(vph[1]-vph[0])+0.5) vi = vn0 vm = [] for i in range(tn0, nT): dvn = int(dv/(vph[1]-vph[0])+0.5) pi = np.argmax(p[vi-dvn: vi+dvn, i]) vi = vi - dvn + pi vm.append([T[i], vph[vi]]) if gcarc/T[i]*r < vph[vi]: break vm = np.array(vm) ax2.plot(vm[:, 0], vm[:, 1], 'k.', lw=2, label='Measured') dc = np.loadtxt('disp.txt') ax2.plot(1/dc[:, 0], dc[:, 1], 'w--', label='Theoretical') plt.legend(loc='lower right', fontsize=15, ncol=2) plt.show() 

## 4. Noise source distribution from beamforming analysis

We utilize the beamforming array method to detect back-azimuth and slowness of noise sources to explore whether these sources come from all directions.

  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  import numpy as np import matplotlib.pyplot as plt from scipy import fft from obspy import read from glob import glob import re rxy = np.loadtxt('rec.txt') * 1e3 files = sorted(glob('ncf_01/*.SAC')) ccf = [] pos = [] t0 = 250 t1 = -t0; t2 = t0 for sac in files: tmp = re.split(r'[_|.]', sac.split('/')[-1]) ri1 = int(tmp[1]); ri2 = int(tmp[2]) tr = read(sac)[0] x1, y1 = rxy[ri1-1] x2, y2 = rxy[ri2-1] 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]; n = len(d) fd = fft.fft(d) ccf.append(fd); pos.append([x2-x1, y2-y1]) ccf = np.array(ccf); pos = np.array(pos) nf = ccf.shape[1]; fs = 1 / dt fre = np.arange(nf) * fs / nf s1 = 1e-4; s2 = 5e-4 sn = 81; bn = 181 slow = np.linspace(s1, s2, sn) baz = np.radians(np.linspace(0, 360, bn)) bazz = baz - np.pi T = [5, 10, 20, 30, 40, 50, 60, 70] dT = [0.1] * 8 plt.figure(figsize=(16, 9)) for ti in range(len(T)): print('Period of %d s ...'%T[ti]) f1 = 1 / (T[ti]+dT[ti]); f2 = 1 / (T[ti]-dT[ti]) fn1 = np.argmin(np.abs(fre-f1)) fn2 = np.argmin(np.abs(fre-f2)) p = np.zeros((sn, bn)) for fn in range(fn1, fn2+1): ff = fre[fn] cf = ccf[:, fn] cf /= np.abs(cf); cf[np.isnan(np.abs(cf))] = 0. tmp = -2j * np.pi * ff for i, ss in enumerate(slow): for j, bb in enumerate(bazz): shift = ss * (np.sin(bb)*pos[:, 0] + np.cos(bb)*pos[:, 1]) p[i, j] += np.abs(np.dot(cf, np.exp(tmp*shift))) p = p / cf.shape[0] / (fn2-fn1+1) ax = plt.subplot(2, 4, ti+1, projection='polar') plt.pcolormesh(baz, slow*1e3, p, cmap='magma') cbar = plt.colorbar(shrink=0.5, pad=0.1) cbar.set_label(r'Beam coherence', fontsize=13) cbar.ax.tick_params(labelsize=13) ax.grid(color='gray', ls=(10, (6, 5)), lw=0.5) ax.tick_params(axis='y', colors='c', labelsize=13) ax.set_theta_zero_location('N') ax.set_theta_direction(-1) ax.set_xlabel('Slowness (s/deg)', fontsize=13) plt.xticks(fontsize=13); plt.yticks(fontsize=13) plt.title('%d s'%(T[ti]), fontsize=13, pad=15) plt.show() 

## References

Yao, H., R. D. van der Hilst, and M. V. de Hoop (2006), Surface-wave array tomography in SE Tibet from ambient seismic noise and two-station analysis—I. Phase velocity maps, Geophys. J. Int., 166, 732 – 744, doi:710.1111/j.1365-1246X.2006.03028.x