In [1]:
from math import exp, pi, sin, cos, tan, asin, acos, atan, atan2, sqrt, floor, ceil
import numpy as np
import scipy.linalg as la
from matplotlib import pyplot as plt
def RayleighMicroseisms( Dn, Nn, De, Ne, Dt, Nt, c, A, sigmadz, inie, test=0 ):
# by Bill Menke, March 2026
# makes synthetic seismograms of spatially-correlated microseismic noise
# caused by dispersive Rayleigh waves from all directions
# Cartesian Coodinate system is positive down, north, east (dz,dn,de)
# Dn, Nn: grid spacing and number of gridpoints in the n direction
# De, Nn: grid spacing and number of gridpoints in the e direction
# notes: Dn=De and Nn=Ne are best
# Nn, Ne, Nt must be even and powers of 2 are best
# Dt, Nt: grid spacing and number of gridpoints in position t
# c: Nw x 1 array of phase velocities, where Nw=Nt/2+1
# evenly spaced between 0 and wmax=2*pi*fmax with fmax=1/(2Dt)
# A: Nw x 1 array of source amplitudes, where Nw=Nt/2+1
# sigmadz: only shape of A is relevant, as a final normaization scales uz to
# standard deviation sigmadz
# inie: Nne x 2 integer array of (in,ie) grid nodes of seismograms
# that are to be returned; in in range 0 to Nn-1,
# ie in range 0 to Ne-1, Nne can be any length
# test=1 performs a test case, so ignore this
# returns three Nt-by-Nne arrays dz, dn, de
# Note: I am using numpy's fft routine np.fft.ifft2() to perform the
# (kn,ke,w) -> (n,e,w) transform and np.fft.irfft() to perform the
# (n,e,w) -> (n,e,t) transform. Both use the same sign of the exponent
# in the IFT, whereas most authors (e.g. Aki and Richards) switch the
# sign of (n,e,w) -> (n,e,t) transform so that waves of positive (kn,ke,w)
# move in the positive (n,e) direction with time. I compensate for this
# difference in the code, and have also examimed tests to verify that
# polarizations are correct.
Nne, i = np.shape(inie);
# The ellipticity is for a homogeneous halfspace of Poisson ration of 0.25
# Aki and Richards, 2ed , prob 7.7
# (r1,r2) at z=0
r1 = 1.0 - 0.5773;
r2 = 0.8475 - 1.4679;
# Aki and Richards, 2ed, eqn 7.99
# n-axis and kn-axis
Nno2 = int(Nn/2);
Nkn=Nno2+1; # number of non-negative frequencies
knmax = 2.0*pi/(2.0*Dn); # Nyquist frequency
Dkn = knmax/(Nn/2); # frequency increment
knp = np.zeros((Nkn,1));
knp[0:Nkn,0] = np.linspace(0.0,knmax,Nkn);
kn = np.zeros((Nn,1));
kn[0:Nn,0]=Dkn * np.concatenate(
(np.linspace(0,Nno2,Nkn), np.linspace(-Nno2+1,-1,Nno2-1) ), axis=0 );
# e-axis and ke-axis
Neo2 = int(Ne/2);
Nke=Neo2+1; # number of non-negative frequencies
kemax = 2.0*pi/(2.0*De); # Nyquist frequency
Dke = kemax/(Ne/2); # frequency increment
kep = np.zeros((Nke,1));
kep[0:Nke,0] = np.linspace(0.0,kemax,Nke);
ke = np.zeros((Ne,1));
ke[0:Ne,0]=Dke * np.concatenate(
(np.linspace(0,Neo2,Nke), np.linspace(-Neo2+1,-1,Neo2-1) ), axis=0 );
# t-axis and w-axis
tmax=Dt*(Nt-1);
t = np.zeros((Nt,1));
t[0:Nt,0] = Dt * np.linspace(0,Nt-1,Nt);
Nto2 = int(Nt/2);
Nw=Nto2+1; # number of non-negative frequencies
fmax = 1.0/(2.0*Dt); # Nyquist frequency
wmax = 2.0*pi*fmax;
Dw = wmax/(Nt/2); # frequency increment
wp = np.zeros((Nw,1));
wp[0:Nw,0] = np.linspace(0.0,wmax,Nw);
w = np.zeros((Nt,1));
w[0:Nt,0]=Dw * np.concatenate(
(np.linspace(0,Nto2,Nw), np.linspace(-Nto2+1,-1,Nto2-1) ), axis=0 );
# fourier transform of selected seismograms
dztne = np.zeros( (Nw,Nne), dtype=complex );
dntne = np.zeros( (Nw,Nne), dtype=complex );
detne = np.zeros( (Nw,Nne), dtype=complex );
# angular information
Nth = 2*(Nn+Ne);
thmax = 2.0*pi;
Dth = thmax/(Nth-1);
th = np.zeros((Nth,1));
th[0:Nth,0] = np.linspace(0.0,thmax,Nth);
sth = np.sin(th);
cth = np.cos(th);
# loop over frequencies
if( test==0 ):
mywrange = range(1,Nw);
else:
mywrange = range(1);
for iiw in mywrange:
if( test == 0 ):
iw = iiw;
else:
iw = 16;
# Fourier transform at all wavenumbers and constant fquency
dzt = np.zeros((Nn,Ne),dtype=complex);
dnt = np.zeros((Nn,Ne),dtype=complex);
det = np.zeros((Nn,Ne),dtype=complex);
wi = wp[iw,0];
ci = c[iw,0];
Ai = A[iw,0];
# random phases
nseR = np.random.normal( loc=0.0, scale=1.0, size=(Nth,1) );
nseI = np.random.normal( loc=0.0, scale=1.0, size=(Nth,1) );
# c = w/k so k = w/c
kr0 = wi/ci;
kn0 = kr0*cth;
ke0 = kr0*sth;
# populate circle in kn-ky plane with random numbers
if( test==0 ):
mythrange = range(Nth);
else:
mythrange=range(1);
for iith in mythrange:
if( test==0 ):
ith = iith;
elif( test==1 ):
ith = int(Nth/2);
else:
ith = int(5*Nth/8);
ikn0 = int(kn0[ith,0]/Dkn);
ike0 = int(ke0[ith,0]/Dke);
if( ikn0 < 0 ):
ikn0 = Nn+ikn0;
if( ikn0 < 0 ):
ikn0 = 0;
elif( ikn0 > (Nn-1) ):
ikn0 = (Nn-1);
if( ike0 < 0 ):
ike0 = Ne+ike0;
if( ike0 < 0 ):
ike0=0;
elif( ike0 > (Ne-1) ):
ike0 = (Ne-1);
dtR = nseR[ith,0];
dtI = nseR[ith,0];
dzt[ikn0, ike0] = r2*complex( 0.0, 1.0 )*complex(dtR,dtI);
dnt[ikn0, ike0] = r1*cth[ith,0]*complex( dtR, dtI );
det[ikn0, ike0] = r1*sth[ith,0]*complex( dtR, dtI );
sqrtpower = sqrt(np.sum(np.power(np.abs(dzt),2)));
dzt = Ai*dzt/sqrtpower;
dnt = Ai*dnt/sqrtpower;
det = Ai*det/sqrtpower;
dz = np.fft.ifft2(dzt); # inverse trasnform (kx,ky) to (x,y)
dn = np.fft.ifft2(dnt);
de = np.fft.ifft2(det);
# grab desired (x,y) positions
for ine in range(Nne):
dztne[ iw, ine ] = dz[ inie[ine,0], inie[ine,1] ];
dntne[ iw, ine ] = dn[ inie[ine,0], inie[ine,1] ];
detne[ iw, ine ] = de[ inie[ine,0], inie[ine,1] ];
# Fourier transform w to t
dz = np.zeros((Nt,Nne));
dn = np.zeros((Nt,Nne));
de = np.zeros((Nt,Nne));
for ine in range(Nne):
# z
uzt = dztne[0:Nw, ine:ine+1];
uz = np.fft.irfft(uzt,axis=0);
stduz = np.std(uz.ravel());
uz = sigmadz*uz/stduz
dz[0:Nt,ine:ine+1] = uz;
# n
unt = dntne[0:Nw, ine:ine+1];
un = np.fft.irfft(unt,axis=0);
un = sigmadz*un/stduz;
dn[0:Nt,ine:ine+1] = un;
# e
uet = detne[0:Nw, ine:ine+1];
ue = np.fft.irfft(uet,axis=0);
ue = sigmadz*ue/stduz;
de[0:Nt,ine:ine+1] = ue;
return(dz,dn,de);
def LoveMicroseisms( Dn, Nn, De, Ne, Dt, Nt, c, A, sigmad, inie, test=0 ):
# by Bill Menke, March 2026
# makes synthetic seismograms of spatially-correlated microseismic noise
# caused by dispersive Love waves from all directions
# Cartesian Coodinate system is positive down, north, east (dz,dn,de)
# Dn, Nn: grid spacing and number of gridpoints in the n direction
# De, Nn: grid spacing and number of gridpoints in the e direction
# notes: Dn=De and Nn=Ne are best
# Nn, Ne, Nt must be even and powers of 2 are best
# Dt, Nt: grid spacing and number of gridpoints in position t
# c: Nw x 1 array of phase velocities, where Nw=Nt/2+1
# evenly spaced between 0 and wmax=2*pi*fmax with fmax=1/(2Dt)
# A: Nw x 1 array of source amplitudes, where Nw=Nt/2+1
# sigmauz: only shape of A is relevant, as a final normaization scales dn and de to
# standard deviation sigmad
# inie: Nne x 2 integer array of (in,ie) grid nodes of seismograms
# that are to be returned; in in range 0 to Nn-1,
# ie in range 0 to Ne-1, Nne can be any length
# test=1 performs a test case, so ignore this
# returns two Nt-by-Nne arrays dn, de
# Note: I am using numpy's fft routine np.fft.ifft2() to perform the
# (kn,ke,w) -> (n,e,w) transform and np.fft.irfft() to perform the
# (n,e,w) -> (n,e,t) transform. Both use the same sign of the exponent
# in the IFT, whereas most authors (e.g. Aki and Richards) switch the
# sign of (n,e,w) -> (n,e,t) transform so that waves of positive (kn,ke,w)
# move in the positive (n,e) direction with time. I compensate for this
# difference in the code, and have also examimed tests to verify that
# polarizations are correct.
Nne, i = np.shape(inie);
# The ellipticity is for a homogeneous halfspace of Poisson ration of 0.25
# Aki and Richards, 2ed , prob 7.7
# (r1,r2) at z=0
r1 = 1.0 - 0.5773;
r2 = 0.8475 - 1.4679;
# Aki and Richards, 2ed, eqn 7.99
# n-axis and kn-axis
Nno2 = int(Nn/2);
Nkn=Nno2+1; # number of non-negative frequencies
knmax = 2.0*pi/(2.0*Dn); # Nyquist frequency
Dkn = knmax/(Nn/2); # frequency increment
knp = np.zeros((Nkn,1));
knp[0:Nkn,0] = np.linspace(0.0,knmax,Nkn);
kn = np.zeros((Nn,1));
kn[0:Nn,0]=Dkn * np.concatenate(
(np.linspace(0,Nno2,Nkn), np.linspace(-Nno2+1,-1,Nno2-1) ), axis=0 );
# e-axis and ke-axis
Neo2 = int(Ne/2);
Nke=Neo2+1; # number of non-negative frequencies
kemax = 2.0*pi/(2.0*De); # Nyquist frequency
Dke = kemax/(Ne/2); # frequency increment
kep = np.zeros((Nke,1));
kep[0:Nke,0] = np.linspace(0.0,kemax,Nke);
ke = np.zeros((Ne,1));
ke[0:Ne,0]=Dke * np.concatenate(
(np.linspace(0,Neo2,Nke), np.linspace(-Neo2+1,-1,Neo2-1) ), axis=0 );
# t-axis and w-axis
tmax=Dt*(Nt-1);
t = np.zeros((Nt,1));
t[0:Nt,0] = Dt * np.linspace(0,Nt-1,Nt);
Nto2 = int(Nt/2);
Nw=Nto2+1; # number of non-negative frequencies
fmax = 1.0/(2.0*Dt); # Nyquist frequency
wmax = 2.0*pi*fmax;
Dw = wmax/(Nt/2); # frequency increment
wp = np.zeros((Nw,1));
wp[0:Nw,0] = np.linspace(0.0,wmax,Nw);
w = np.zeros((Nt,1));
w[0:Nt,0]=Dw * np.concatenate(
(np.linspace(0,Nto2,Nw), np.linspace(-Nto2+1,-1,Nto2-1) ), axis=0 );
# fourier transform of selected seismograms
dntne = np.zeros( (Nw,Nne), dtype=complex );
detne = np.zeros( (Nw,Nne), dtype=complex );
# angular information
Nth = 2*(Nn+Ne);
thmax = 2.0*pi;
Dth = thmax/(Nth-1);
th = np.zeros((Nth,1));
th[0:Nth,0] = np.linspace(0.0,thmax,Nth);
sth = np.sin(th);
cth = np.cos(th);
# loop over frequencies
if( test==0 ):
mywrange = range(1,Nw);
else:
mywrange = range(1);
for iiw in mywrange:
if( test == 0 ):
iw = iiw;
else:
iw = 16;
# Fourier transform at all wavenumbers and constant fquency
dnt = np.zeros((Nn,Ne),dtype=complex);
det = np.zeros((Nn,Ne),dtype=complex);
wi = wp[iw,0];
ci = c[iw,0];
Ai = A[iw,0];
# random phases
nseR = np.random.normal( loc=0.0, scale=1.0, size=(Nth,1) );
nseI = np.random.normal( loc=0.0, scale=1.0, size=(Nth,1) );
# c = w/k so k = w/c
kr0 = wi/ci;
kn0 = kr0*cth;
ke0 = kr0*sth;
# populate circle in kn-ky plane with random numbers
if( test==0 ):
mythrange = range(Nth);
else:
mythrange=range(1);
for iith in mythrange:
if( test==0 ):
ith = iith;
elif( test==1 ):
ith = int(Nth/2);
else:
ith = int(5*Nth/8);
ikn0 = int(kn0[ith,0]/Dkn);
ike0 = int(ke0[ith,0]/Dke);
if( ikn0 < 0 ):
ikn0 = Nn+ikn0;
if( ikn0 < 0 ):
ikn0 = 0;
elif( ikn0 > (Nn-1) ):
ikn0 = (Nn-1);
if( ike0 < 0 ):
ike0 = Ne+ike0;
if( ike0 < 0 ):
ike0=0;
elif( ike0 > (Ne-1) ):
ike0 = (Ne-1);
dtR = nseR[ith,0];
dtI = nseR[ith,0];
dnt[ikn0, ike0] = sth[ith,0]*complex( dtR, dtI );
det[ikn0, ike0] = -cth[ith,0]*complex( dtR, dtI );
sqrtpower = sqrt(
np.sum(np.power(np.abs(dnt),2)) +
np.sum(np.power(np.abs(det),2)) );
dnt = Ai*dnt/sqrtpower;
det = Ai*det/sqrtpower;
dn = np.fft.ifft2(dnt);
de = np.fft.ifft2(det);
# grab desired (x,y) positions
for ine in range(Nne):
dntne[ iw, ine ] = dn[ inie[ine,0], inie[ine,1] ];
detne[ iw, ine ] = de[ inie[ine,0], inie[ine,1] ];
# Fourier transform w to t
dn = np.zeros((Nt,Nne));
de = np.zeros((Nt,Nne));
for ine in range(Nne):
# n comp
unt = dntne[0:Nw, ine:ine+1];
un = np.fft.irfft(unt,axis=0);
# e comp
uet = detne[0:Nw, ine:ine+1];
ue = np.fft.irfft(uet,axis=0);
# std dev
stdune = np.std( np.concatenate( (un.ravel(),ue.ravel()), axis=0) );
un = sigmad*un/stdune;
ue = sigmad*ue/stdune;
# save
dn[0:Nt,ine:ine+1] = un;
de[0:Nt,ine:ine+1] = ue;
return(dn,de);
def randomize_phase(d):
# randomizes the phase of Nx1 timeseries d. N must be even.
N,i = np.shape(d);
Nf = int(N/2)+1;
ph = np.random.uniform(low=0.0,high=2.0*pi,size=(Nf,1));
ph[0,0]=0.0;
ph[Nf-1,0]=0.0;
dt = np.fft.rfft(d, axis=0);
dt = np.multiply(dt, np.exp(complex(0.0,1.0)*ph));
d = np.real(np.fft.irfft(dt,axis=0));
return(d);
In [2]:
# demo code for RayleighMicroseisms
# for Aki-type noise, use RANDOMIZE=0
# to randomize phase but preserve power, use RANDOMIZE=1
RANDOMIZE=0;
# n-axis
Nn=512;
Dn = 1.0;
nmin=0.0;
n = np.zeros((Nn,1));
n[0:Nn,0] = Dn*np.linspace(0,Nn-1,Nn);
nmin = n[0,0];
nmax = n[Nn-1,0];
# e-axis
Ne=512;
De = 1.0;
emin=0.0;
e = np.zeros((Ne,1));
e[0:Ne,0] = De*np.linspace(0,Ne-1,Ne);
emin = e[0,0];
emax = e[Ne-1,0];
# time axis
N=1024;
Nt=N;
No2 = int(N/2);
Dt = 1.0;
tmin=0.0;
t = np.zeros((N,1));
t[0:N,0] = Dt*np.linspace(0,N-1,N);
tmin = t[0,0];
tmax = t[N-1,0];
# station positions
iaperature=10;
Nne = 3;
ine = np.zeros((Nne,2),dtype=int);
inbase=int(Nn/2);
iebase=int(Ne/2);
ine[0,0] = inbase; ine[0,1] = iebase;
ine[1,0] = inbase; ine[1,1] = iebase+iaperature;
ine[2,0] = inbase+iaperature; ine[2,1] = iebase;
n1 = n[ine[0,0],0]; e1 = e[ ine[0,1], 0];
n2 = n[ine[1,0],0]; e2 = e[ ine[1,1], 0];
n3 = n[ine[2,0],0]; e3 = e[ ine[2,1], 0];
print("stations at:");
print("%.2f %.2f" % (n1,e1) );
print("%.2f %.2f" % (n2,e2) );
print("%.2f %.2f" % (n3,e3) );
print(" ");
# generic frequency setup
fmax=1/(2.0*Dt); # nyquist frequency
Df=fmax/No2; # frequency sampling
Nf=No2+1; # number of non-negative frequencies
# f is full freqeuncy vector:
# f = [0, Df, 2Df ... fmax, -fmax+Df, ... -Df]
f = np.zeros((N,1));
f[0:N,0] = Df*np.concatenate(
(np.linspace(0,No2,Nf),np.linspace(-No2+1,-1,No2-1)),
axis=0 );
Dw=2*pi*Df; # angular frequency sampling
w=2*pi*f; # full angular frequency vector
Nw=Nf; # number of non-negative angular f's
fpos = np.zeros((Nf,1));
fpos[0:Nf,0] = Df * np.linspace(0,No2,Nf); # non-neg f's
wpos=2*pi*fpos; # non-negative angular frequencies
# phase velocity
f0 = fmax/10.0;
v = 3.0*np.ones((N,1))+1.0*np.exp( -np.abs(f)/f0 );
A = np.zeros((Nf,1));
c = np.zeros((Nf,1));
wc = 2.0*pi*0.05;
sw = 2.0*pi*0.05;
sw2 = sw**2;
for iw in range(Nf):
wi = w[iw,0];
c[iw,0] = v[iw,0];
A[iw,0] = exp(-0.5*(wi-wc)**2/sw2 );
fig=plt.figure();
ax=plt.subplot(1,1,1);
plt.axis( [0.0,0.1,0.0,1.0] );
plt.plot(fpos,A,'k-');
plt.xlabel('f');
plt.ylabel('A');
plt.show();
print('Caption: Amplitude curve');
print(" ");
print(" ");
fig=plt.figure();
ax=plt.subplot(1,1,1);
plt.axis( [0.0,0.1,0.0,5.0] );
plt.plot(fpos,c,'k-');
plt.xlabel('f');
plt.ylabel('c');
plt.show();
print('Caption: Dispersion curve');
print(" ");
print(" ");
# compute micro-seismic noise wavefield
sigmadz=0.0025;
dz, dn, de = RayleighMicroseisms( Dn, Nn, De, Ne, Dt, Nt, c, A, sigmadz, ine, test=0 );
if( RANDOMIZE ):
for i in range(Nne):
dz[0:Nt,i:i+1] = randomize_phase( dz[0:Nt,i:i+1] );
dn[0:Nt,i:i+1] = randomize_phase( dn[0:Nt,i:i+1] );
de[0:Nt,i:i+1] = randomize_phase( de[0:Nt,i:i+1] );
fig=plt.figure();
for i in range(Nne):
ax=plt.subplot(Nne,1,i+1);
plt.plot( t, dz[0:Nt,i:i+1], "g-" );
plt.plot( t, dn[0:Nt,i:i+1], "r-" );
plt.plot( t, de[0:Nt,i:i+1], "b-" );
plt.xlabel('t');
plt.ylabel('noise');
plt.show();
print('Caption: Microseismic noise uz (green) un (red) ue (blue)');
stations at: 256.00 256.00 256.00 266.00 266.00 256.00
Caption: Amplitude curve
Caption: Dispersion curve
Caption: Microseismic noise uz (green) un (red) ue (blue)
In [3]:
# demo code for LoveMicroseisms
# for Aki-type noise, use RANDOMIZE=0
# to randomize phase but preserve power, use RANDOMIZE=1
RANDOMIZE=0;
# n-axis
Nn=512;
Dn = 1.0;
nmin=0.0;
n = np.zeros((Nn,1));
n[0:Nn,0] = Dn*np.linspace(0,Nn-1,Nn);
nmin = n[0,0];
nmax = n[Nn-1,0];
# e-axis
Ne=512;
De = 1.0;
emin=0.0;
e = np.zeros((Ne,1));
e[0:Ne,0] = De*np.linspace(0,Ne-1,Ne);
emin = e[0,0];
emax = e[Ne-1,0];
# time axis
N=1024;
Nt=N;
No2 = int(N/2);
Dt = 1.0;
tmin=0.0;
t = np.zeros((N,1));
t[0:N,0] = Dt*np.linspace(0,N-1,N);
tmin = t[0,0];
tmax = t[N-1,0];
# station positions
iaperature=10;
Nne = 3;
ine = np.zeros((Nne,2),dtype=int);
inbase=int(Nn/2);
iebase=int(Ne/2);
ine[0,0] = inbase; ine[0,1] = iebase;
ine[1,0] = inbase; ine[1,1] = iebase+iaperature;
ine[2,0] = inbase+iaperature; ine[2,1] = iebase;
n1 = n[ine[0,0],0]; e1 = e[ ine[0,1], 0];
n2 = n[ine[1,0],0]; e2 = e[ ine[1,1], 0];
n3 = n[ine[2,0],0]; e3 = e[ ine[2,1], 0];
print("stations at:");
print("%.2f %.2f" % (n1,e1) );
print("%.2f %.2f" % (n2,e2) );
print("%.2f %.2f" % (n3,e3) );
print(" ");
# generic frequency setup
fmax=1/(2.0*Dt); # nyquist frequency
Df=fmax/No2; # frequency sampling
Nf=No2+1; # number of non-negative frequencies
# f is full freqeuncy vector:
# f = [0, Df, 2Df ... fmax, -fmax+Df, ... -Df]
f = np.zeros((N,1));
f[0:N,0] = Df*np.concatenate(
(np.linspace(0,No2,Nf),np.linspace(-No2+1,-1,No2-1)),
axis=0 );
Dw=2*pi*Df; # angular frequency sampling
w=2*pi*f; # full angular frequency vector
Nw=Nf; # number of non-negative angular f's
fpos = np.zeros((Nf,1));
fpos[0:Nf,0] = Df * np.linspace(0,No2,Nf); # non-neg f's
wpos=2*pi*fpos; # non-negative angular frequencies
# phase velocity
f0 = fmax/10.0;
v = 3.0*np.ones((N,1))+1.0*np.exp( -np.abs(f)/f0 );
A = np.zeros((Nf,1));
c = np.zeros((Nf,1));
wc = 2.0*pi*0.05;
sw = 2.0*pi*0.05;
sw2 = sw**2;
for iw in range(Nf):
wi = w[iw,0];
c[iw,0] = v[iw,0];
A[iw,0] = exp(-0.5*(wi-wc)**2/sw2 );
fig=plt.figure();
ax=plt.subplot(1,1,1);
plt.axis( [0.0,0.1,0.0,1.0] );
plt.plot(fpos,A,'k-');
plt.xlabel('f');
plt.ylabel('A');
plt.show();
print('Caption: Amplitude curve');
print(" ");
print(" ");
fig=plt.figure();
ax=plt.subplot(1,1,1);
plt.axis( [0.0,0.1,0.0,5.0] );
plt.plot(fpos,c,'k-');
plt.xlabel('f');
plt.ylabel('c');
plt.show();
print('Caption: Dispersion curve');
print(" ");
print(" ");
# compute micro-seismic noise wavefield
sigmadz=0.0025;
dn, de = LoveMicroseisms( Dn, Nn, De, Ne, Dt, Nt, c, A, sigmadz, ine, test=0 );
if( RANDOMIZE ):
for i in range(Nxy):
dn[0:Nt,i:i+1] = randomize_phase( dn[0:Nt,i:i+1] );
de[0:Nt,i:i+1] = randomize_phase( de[0:Nt,i:i+1] );
fig=plt.figure();
for i in range(Nne):
ax=plt.subplot(Nne,1,i+1);
plt.plot( t, dn[0:Nt,i:i+1], "r-" );
plt.plot( t, de[0:Nt,i:i+1], "g-" );
plt.xlabel('t');
plt.ylabel('noise');
plt.show();
print('Caption: Microseismic noise un (red) ue (green)');
stations at: 256.00 256.00 256.00 266.00 266.00 256.00
Caption: Amplitude curve
Caption: Dispersion curve
Caption: Microseismic noise un (red) ue (green)
In [4]:
# Rayleigh test code, checks polarities
# n-axis
Nn=512;
Dn = 1.0;
nmin=0.0;
n = np.zeros((Nn,1));
n[0:Nn,0] = Dn*np.linspace(0,Nn-1,Nn);
nmin = n[0,0];
nmax = n[Nn-1,0];
# e-axis
Ne=512;
De = 1.0;
emin=0.0;
e = np.zeros((Ne,1));
e[0:Ne,0] = De*np.linspace(0,Ne-1,Ne);
emin = e[0,0];
emax = e[Ne-1,0];
# time axis
N=1024;
Nt=N;
No2 = int(N/2);
Dt = 1.0;
tmin=0.0;
t = np.zeros((N,1));
t[0:N,0] = Dt*np.linspace(0,N-1,N);
tmin = t[0,0];
tmax = t[N-1,0];
# generic frequency setup
fmax=1/(2.0*Dt); # nyquist frequency
Df=fmax/No2; # frequency sampling
Nf=No2+1; # number of non-negative frequencies
# f is full freqeuncy vector:
# f = [0, Df, 2Df ... fmax, -fmax+Df, ... -Df]
f = np.zeros((N,1));
f[0:N,0] = Df*np.concatenate(
(np.linspace(0,No2,Nf),np.linspace(-No2+1,-1,No2-1)),
axis=0 );
Dw=2*pi*Df; # angular frequency sampling
w=2*pi*f; # full angular frequency vector
Nw=Nf; # number of non-negative angular f's
fpos = np.zeros((Nf,1));
fpos[0:Nf,0] = Df * np.linspace(0,No2,Nf); # non-neg f's
wpos=2*pi*fpos; # non-negative angular frequencies
# phase velocity
f0 = fmax/10.0;
v = 3.0*np.ones((N,1))+1.0*np.exp( -np.abs(f)/f0 );
A = np.zeros((Nf,1));
c = np.zeros((Nf,1));
wc = 2.0*pi*0.05;
sw = 2.0*pi*0.05;
sw2 = sw**2;
for iw in range(Nf):
wi = w[iw,0];
c[iw,0] = v[iw,0];
A[iw,0] = exp(-0.5*(wi-wc)**2/sw2 );
# Part 1 tests
# station positions
iaperature=20;
Nne = 2;
ine = np.zeros((Nne,2),dtype=int);
inbase=int(Nn/2);
iebase=int(Ne/2);
ine[0,0] = inbase; ine[0,1] = iebase;
ine[1,0] = inbase+iaperature; ine[1,1] = iebase;
n1 = n[ine[0,0],0]; e1 = e[ ine[0,1], 0];
n2 = n[ine[1,0],0]; e2 = e[ ine[1,1], 0];
print("stations at:");
print("%.2f %.2f" % (n1,e1) );
print("%.2f %.2f" % (n2,e2) );
print(" ");
# compute micro-seismic noise wavefield
sigmadz=1.0;
dz, dn, de = RayleighMicroseisms( Dn, Nn, De, Ne, Dt, Nt, c, A, sigmadz, ine, test=1 );
fig=plt.figure();
ax=plt.subplot(1,1,1);
plt.axis( [0, 200, -2, 2] );
plt.plot( t, dz[0:Nt,0:1], "r-" );
plt.plot( t, dz[0:Nt,1:2], "g-" );
plt.xlabel('t');
plt.ylabel('dz at two locations');
plt.show();
print('Caption: Green should be delayed with respect to red');
fig=plt.figure();
ax=plt.subplot(1,1,1);
plt.axis( [0, 200, -2, 2] );
plt.plot( t, dn[0:Nt,0:1], "r-" );
plt.plot( t, de[0:Nt,0:1], "g-" );
plt.xlabel('t');
plt.ylabel('dn and de at same location');
plt.show();
print('Caption: Green should be zero');
fig=plt.figure();
ax=plt.subplot(1,1,1);
plt.axis( [-2, 2, -2, 2] );
ax.invert_yaxis(); # so that z is down
plt.plot( dn[0:Nt,0:1], dz[0:Nt,0:1], "k-" );
plt.plot( dn[0:1,0:1], dz[0:1,0:1], "r." );
plt.plot( dn[1:2,0:1], dz[1:2,0:1], "g." );
plt.plot( dn[2:3,0:1], dz[2:3,0:1], "b." );
plt.xlabel('t');
plt.ylabel('un vs uz at same locations');
plt.show();
print('Caption: Check retrograde. order of dots red, green, blue');
# Part 2 tests
# station positions
iaperature=20;
Nne = 2;
ine = np.zeros((Nne,2),dtype=int);
inbase=int(Nn/2);
iebase=int(Ne/2);
ine[0,0] = inbase; ine[0,1] = iebase;
ine[1,0] = inbase+iaperature; ine[1,1] = iebase+iaperature;
n1 = n[ine[0,0],0]; e1 = e[ ine[0,1], 0];
n2 = n[ine[1,0],0]; e2 = e[ ine[1,1], 0];
print("stations at:");
print("%.2f %.2f" % (n1,e1) );
print("%.2f %.2f" % (n2,e2) );
print(" ");
# compute micro-seismic noise wavefield
sigmadz=1.0;
dz, dn, de = RayleighMicroseisms( Dn, Nn, De, Ne, Dt, Nt, c, A, sigmadz, ine, test=2 );
fig=plt.figure();
ax=plt.subplot(1,1,1);
plt.axis( [0, 200, -2, 2] );
plt.plot( t, dz[0:Nt,0:1], "r-" );
plt.plot( t, dz[0:Nt,1:2], "g-" );
plt.xlabel('t');
plt.ylabel('dz at two locations');
plt.show();
print('Caption: Green should be delayed with respect to red');
fig=plt.figure();
ax=plt.subplot(1,1,1);
plt.axis( [0, 200, -2, 2] );
plt.plot( t, dn[0:Nt,0:1], "r-", lw=3 );
plt.plot( t, de[0:Nt,0:1], "g-" );
plt.xlabel('t');
plt.ylabel('dn and de at same location');
plt.show();
print('Caption: Green and red should be the same');
fig=plt.figure();
ax=plt.subplot(1,1,1);
plt.axis( [-2, 2, -2, 2] );
ax.invert_yaxis(); # so that z is down
radial = (dn[0:Nt,0:1]+de[0:Nt,0:1])/sqrt(2.0);
plt.plot( radial, dz[0:Nt,0:1], "k-" );
plt.plot( radial[0:1,0:1], dz[0:1,0:1], "r." );
plt.plot( radial[1:2,0:1], dz[1:2,0:1], "g." );
plt.plot( radial[2:3,0:1], dz[2:3,0:1], "b." );
plt.xlabel('t');
plt.ylabel('ur vs uz at same locations');
plt.show();
print('Caption: Check retrograde. order of dots red, green, blue');
stations at: 256.00 256.00 276.00 256.00
Caption: Green should be delayed with respect to red
Caption: Green should be zero
Caption: Check retrograde. order of dots red, green, blue stations at: 256.00 256.00 276.00 276.00
Caption: Green should be delayed with respect to red
Caption: Green and red should be the same
Caption: Check retrograde. order of dots red, green, blue
In [5]:
# Love test code, checks polarities
# n-axis
Nn=512;
Dn = 1.0;
nmin=0.0;
n = np.zeros((Nn,1));
n[0:Nn,0] = Dn*np.linspace(0,Nn-1,Nn);
nmin = n[0,0];
nmax = n[Nn-1,0];
# e-axis
Ne=512;
De = 1.0;
emin=0.0;
e = np.zeros((Ne,1));
e[0:Ne,0] = De*np.linspace(0,Ne-1,Ne);
emin = e[0,0];
emax = e[Ne-1,0];
# time axis
N=1024;
Nt=N;
No2 = int(N/2);
Dt = 1.0;
tmin=0.0;
t = np.zeros((N,1));
t[0:N,0] = Dt*np.linspace(0,N-1,N);
tmin = t[0,0];
tmax = t[N-1,0];
# generic frequency setup
fmax=1/(2.0*Dt); # nyquist frequency
Df=fmax/No2; # frequency sampling
Nf=No2+1; # number of non-negative frequencies
# f is full freqeuncy vector:
# f = [0, Df, 2Df ... fmax, -fmax+Df, ... -Df]
f = np.zeros((N,1));
f[0:N,0] = Df*np.concatenate(
(np.linspace(0,No2,Nf),np.linspace(-No2+1,-1,No2-1)),
axis=0 );
Dw=2*pi*Df; # angular frequency sampling
w=2*pi*f; # full angular frequency vector
Nw=Nf; # number of non-negative angular f's
fpos = np.zeros((Nf,1));
fpos[0:Nf,0] = Df * np.linspace(0,No2,Nf); # non-neg f's
wpos=2*pi*fpos; # non-negative angular frequencies
# phase velocity
f0 = fmax/10.0;
v = 3.0*np.ones((N,1))+1.0*np.exp( -np.abs(f)/f0 );
A = np.zeros((Nf,1));
c = np.zeros((Nf,1));
wc = 2.0*pi*0.05;
sw = 2.0*pi*0.05;
sw2 = sw**2;
for iw in range(Nf):
wi = w[iw,0];
c[iw,0] = v[iw,0];
A[iw,0] = exp(-0.5*(wi-wc)**2/sw2 );
# Part 1 tests
# station positions
iaperature=20;
Nne = 2;
ine = np.zeros((Nne,2),dtype=int);
inbase=int(Nn/2);
iebase=int(Ne/2);
ine[0,0] = inbase; ine[0,1] = iebase;
ine[1,0] = inbase+iaperature; ine[1,1] = iebase;
n1 = n[ine[0,0],0]; e1 = e[ ine[0,1], 0];
n2 = n[ine[1,0],0]; e2 = e[ ine[1,1], 0];
print("stations at:");
print("%.2f %.2f" % (n1,e1) );
print("%.2f %.2f" % (n2,e2) );
print(" ");
# compute micro-seismic noise wavefield
sigmadz=1.0;
dn, de = LoveMicroseisms( Dn, Nn, De, Ne, Dt, Nt, c, A, sigmadz, ine, test=1 );
fig=plt.figure();
ax=plt.subplot(1,1,1);
plt.axis( [0, 200, -2, 2] );
plt.plot( t, de[0:Nt,0:1], "r-" );
plt.plot( t, de[0:Nt,1:2], "g-" );
plt.xlabel('t');
plt.ylabel('dn at two locations');
plt.show();
print('Caption: Green should be delayed with respect to red');
fig=plt.figure();
ax=plt.subplot(1,1,1);
plt.axis( [0, 200, -2, 2] );
plt.plot( t, dn[0:Nt,0:1], "r-" );
plt.plot( t, de[0:Nt,0:1], "g-" );
plt.xlabel('t');
plt.ylabel('dn and de at same location');
plt.show();
print('Caption: Red should be zero');
# Part 2 tests
# station positions
iaperature=20;
Nne = 2;
ine = np.zeros((Nne,2),dtype=int);
inbase=int(Nn/2);
iebase=int(Ne/2);
ine[0,0] = inbase; ine[0,1] = iebase;
ine[1,0] = inbase+iaperature; ine[1,1] = iebase+iaperature;
n1 = n[ine[0,0],0]; e1 = e[ ine[0,1], 0];
n2 = n[ine[1,0],0]; e2 = e[ ine[1,1], 0];
print("stations at:");
print("%.2f %.2f" % (n1,e1) );
print("%.2f %.2f" % (n2,e2) );
print(" ");
# compute micro-seismic noise wavefield
sigmadz=1.0;
dn, de = LoveMicroseisms( Dn, Nn, De, Ne, Dt, Nt, c, A, sigmadz, ine, test=2 );
fig=plt.figure();
ax=plt.subplot(1,1,1);
plt.axis( [0, 200, -2, 2] );
plt.plot( t, de[0:Nt,0:1], "r-" );
plt.plot( t, de[0:Nt,1:2], "g-" );
plt.xlabel('t');
plt.ylabel('de at two locations');
plt.show();
print('Caption: Green should be delayed with respect to red');
fig=plt.figure();
ax=plt.subplot(1,1,1);
plt.axis( [0, 200, -2, 2] );
plt.plot( t, dn[0:Nt,0:1], "r-", lw=3 );
plt.plot( t, de[0:Nt,0:1], "g-" );
plt.xlabel('t');
plt.ylabel('dn and de at same location');
plt.show();
print('Caption: Green and red should have opposite polarity');
stations at: 256.00 256.00 276.00 256.00
Caption: Green should be delayed with respect to red
Caption: Red should be zero stations at: 256.00 256.00 276.00 276.00
Caption: Green should be delayed with respect to red
Caption: Green and red should have opposite polarity
In [ ]: