DataAnalysis2021/06-Surface_Waves/phaseVelocityCrossCorr.ipynb

Two-station method

Here we compute the phase velocity form seismograms of two stations lying on a common great circle. We make use of the fact that the phase of the Fourier tranform of the cross-correlation of the two seismograms is equal to the phase difference between the surface wave trains from which we calclate the phase velocity.

In [ ]:

# Preparation: load packages, set some basic options  
%matplotlib inline
from obspy import UTCDateTime
from obspy.clients.fdsn import Client
from obspy.geodetics import locations2degrees
import numpy as np
import matplotlib.pylab as plt
from multipleFilter import multipleFilterAnalysis, normalizeMFT 
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = 15, 4
plt.rcParams['lines.linewidth'] = 0.5

In [ ]:

# Get data from earthquake (Gulf of Alaska, m=7.9, d = 24 km)
client = Client("BGR")
t1 = UTCDateTime("2018-01-23T09:31:42.000")
st1raw = client.get_waveforms("GR","TNS","","LHZ",t1,t1 + 2 * 3600,       # data of station TNS of GRSN
                          attach_response = True)
inv1 = client.get_stations(network="GR",station="TNS",
                           channel='LHZ',level="channel",starttime=t1)    # inventory info about TNS
client = Client("GFZ")
st2raw = client.get_waveforms("GE","STU","","LHZ",t1,t1 + 2 * 3600,       # data of station STU of GEOFON
                          attach_response = True)
inv2 = client.get_stations(network="GE",station="STU",
                           channel='LHZ',level="channel",starttime=t1)    # inventory info about STU

In [ ]:

st1 = st1raw.copy()                                                 # copy data
st2 = st2raw.copy()
st1.detrend('linear')                                               # detrend
st2.detrend('linear')
st1.remove_response(output='VEL',zero_mean = True,                  # remove response (pre-filter is important)
                    pre_filt = (0.001, 0.005, 0.4, 0.5))
st2.remove_response(output='VEL',zero_mean = True,                  # remove response(pre-filter is important)
                    pre_filt = (0.001, 0.005, 0.4, 0.5))
st1.plot()                                                          # plot data
st2.plot()
print(st1)
print(st2)

In [ ]:

co = inv1.get_coordinates("GR.TNS..LHZ")                                # extract coordinates from inventory
tns = (co["longitude"],co["latitude"])
co = inv2.get_coordinates("GE.STU..LHZ")
stu = (co["longitude"],co["latitude"])
event = (-149.09,56.05)                                                 # coordinates of event
dis = locations2degrees(tns[1],tns[0],stu[1],stu[0])*np.pi/180.*6371.   # epicentral distance
print("Distance in km: ",dis)

In [ ]:

tw1 = UTCDateTime("2018-01-23T10:05:00.000")          # cut seismograms to surface wave train
tw2 = UTCDateTime("2018-01-23T10:16:00.000")
st1.trim(tw1,tw2)
st2.trim(tw1,tw2)
st1.plot()
st2.plot()
npts = st1[0].stats.npts
dt = st1[0].stats.delta
print("Number of samples: ",npts,"Sampling interval: ",dt)

In [ ]:

cc = np.correlate(st2[0],st1[0],mode="full")          # Cross correlation: cc(k)=sum_n a(n+k)b(n)
tm = np.arange(-(npts-1),npts)*dt                     # includes negative lags as well
plt.plot(tm,cc,'-k')
plt.show()

In [ ]:

ns = 256                                       # limit length of cross correlation to first ns points
ccp = cc[npts-1:npts-1+ns]                     # cross correlation for positive lag
print(len(ccp))
tm = np.arange(0,ns)*dt                        # associated time values
plt.plot(tm,ccp,'-k')
plt.xlabel("Time")
plt.show()

In [ ]:

nd = int(pow(2, np.ceil(np.log(ns)/np.log(2))))          # Fourier transfrom of the cross correlation
fcp = np.fft.rfft(ccp,nd)
freq = np.fft.rfftfreq(nd,dt)
print(nd,dt,len(freq),freq[0:4])

In [ ]:

fmax = 0.1                                               # compute phase up to maximum frequency fmax
fny = 1./(2.*dt)
jmax = int(fmax/fny*(nd//2-1))
print(jmax)
spec = fcp[0:jmax]
frs = freq[0:jmax]
phase = np.unwrap(np.angle(spec))                        # unwrap the phase owing to 2*pi multiples

In [ ]:

plt.subplot(311)
plt.plot(frs,np.abs(spec),"-k")                          # plot abs of spectrum
plt.subplot(312)
plt.plot(frs,np.angle(spec),"-k")                        # plot phase of spectrum
plt.subplot(313)
plt.plot(frs,phase,"-k")                                 # plot unwrapped phase
plt.xlabel("Frequency [Hz]")
plt.show()

In [ ]:

j1 = int(jmax*0.1)
j2 = int(jmax*0.6)
phavel = -2.*np.pi*frs[j1:j2]*dis/phase[j1:j2]                # calculate phase velocity and plot
plt.plot(frs[j1:j2],phavel,'-k')
plt.xlabel("Frequency [Hz]")
plt.ylabel("Phase velocity [km/s]")
plt.show()

In [ ]: