Janis Heuel
a47d87bc9b
Reviewed-on: #13 Co-authored-by: Janis Heuel <janis.heuel@ruhr-uni-bochum.de> Co-committed-by: Janis Heuel <janis.heuel@ruhr-uni-bochum.de>
269 lines
7.9 KiB
Plaintext
269 lines
7.9 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Two-station method\n",
|
|
"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."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Preparation: load packages, set some basic options \n",
|
|
"%matplotlib inline\n",
|
|
"from obspy import UTCDateTime\n",
|
|
"from obspy.clients.fdsn import Client\n",
|
|
"from obspy.geodetics import locations2degrees\n",
|
|
"import numpy as np\n",
|
|
"import matplotlib.pylab as plt\n",
|
|
"from multipleFilter import multipleFilterAnalysis, normalizeMFT \n",
|
|
"plt.style.use('ggplot')\n",
|
|
"plt.rcParams['figure.figsize'] = 15, 4\n",
|
|
"plt.rcParams['lines.linewidth'] = 0.5"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Get data from earthquake (Gulf of Alaska, m=7.9, d = 24 km)\n",
|
|
"client = Client(\"BGR\")\n",
|
|
"t1 = UTCDateTime(\"2018-01-23T09:31:42.000\")\n",
|
|
"st1raw = client.get_waveforms(\"GR\",\"TNS\",\"\",\"LHZ\",t1,t1 + 2 * 3600, # data of station TNS of GRSN\n",
|
|
" attach_response = True)\n",
|
|
"inv1 = client.get_stations(network=\"GR\",station=\"TNS\",\n",
|
|
" channel='LHZ',level=\"channel\",starttime=t1) # inventory info about TNS\n",
|
|
"client = Client(\"GFZ\")\n",
|
|
"st2raw = client.get_waveforms(\"GE\",\"STU\",\"\",\"LHZ\",t1,t1 + 2 * 3600, # data of station STU of GEOFON\n",
|
|
" attach_response = True)\n",
|
|
"inv2 = client.get_stations(network=\"GE\",station=\"STU\",\n",
|
|
" channel='LHZ',level=\"channel\",starttime=t1) # inventory info about STU"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"st1 = st1raw.copy() # copy data\n",
|
|
"st2 = st2raw.copy()\n",
|
|
"st1.detrend('linear') # detrend\n",
|
|
"st2.detrend('linear')\n",
|
|
"st1.remove_response(output='VEL',zero_mean = True, # remove response (pre-filter is important)\n",
|
|
" pre_filt = (0.001, 0.005, 0.4, 0.5))\n",
|
|
"st2.remove_response(output='VEL',zero_mean = True, # remove response(pre-filter is important)\n",
|
|
" pre_filt = (0.001, 0.005, 0.4, 0.5))\n",
|
|
"st1.plot() # plot data\n",
|
|
"st2.plot()\n",
|
|
"print(st1)\n",
|
|
"print(st2)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"co = inv1.get_coordinates(\"GR.TNS..LHZ\") # extract coordinates from inventory\n",
|
|
"tns = (co[\"longitude\"],co[\"latitude\"])\n",
|
|
"co = inv2.get_coordinates(\"GE.STU..LHZ\")\n",
|
|
"stu = (co[\"longitude\"],co[\"latitude\"])\n",
|
|
"event = (-149.09,56.05) # coordinates of event\n",
|
|
"dis = locations2degrees(tns[1],tns[0],stu[1],stu[0])*np.pi/180.*6371. # epicentral distance\n",
|
|
"print(\"Distance in km: \",dis)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<img src=\"../images/phase_velocity_map.png\">"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"tw1 = UTCDateTime(\"2018-01-23T10:05:00.000\") # cut seismograms to surface wave train\n",
|
|
"tw2 = UTCDateTime(\"2018-01-23T10:16:00.000\")\n",
|
|
"st1.trim(tw1,tw2)\n",
|
|
"st2.trim(tw1,tw2)\n",
|
|
"st1.plot()\n",
|
|
"st2.plot()\n",
|
|
"npts = st1[0].stats.npts\n",
|
|
"dt = st1[0].stats.delta\n",
|
|
"print(\"Number of samples: \",npts,\"Sampling interval: \",dt)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"cc = np.correlate(st2[0],st1[0],mode=\"full\") # Cross correlation: cc(k)=sum_n a(n+k)b(n)\n",
|
|
"tm = np.arange(-(npts-1),npts)*dt # includes negative lags as well\n",
|
|
"plt.plot(tm,cc,'-k')\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
},
|
|
"scrolled": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"ns = 256 # limit length of cross correlation to first ns points\n",
|
|
"ccp = cc[npts-1:npts-1+ns] # cross correlation for positive lag\n",
|
|
"print(len(ccp))\n",
|
|
"tm = np.arange(0,ns)*dt # associated time values\n",
|
|
"plt.plot(tm,ccp,'-k')\n",
|
|
"plt.xlabel(\"Time\")\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"nd = int(pow(2, np.ceil(np.log(ns)/np.log(2)))) # Fourier transfrom of the cross correlation\n",
|
|
"fcp = np.fft.rfft(ccp,nd)\n",
|
|
"freq = np.fft.rfftfreq(nd,dt)\n",
|
|
"print(nd,dt,len(freq),freq[0:4])"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"fmax = 0.1 # compute phase up to maximum frequency fmax\n",
|
|
"fny = 1./(2.*dt)\n",
|
|
"jmax = int(fmax/fny*(nd//2-1))\n",
|
|
"print(jmax)\n",
|
|
"spec = fcp[0:jmax]\n",
|
|
"frs = freq[0:jmax]\n",
|
|
"phase = np.unwrap(np.angle(spec)) # unwrap the phase owing to 2*pi multiples"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"plt.subplot(311)\n",
|
|
"plt.plot(frs,np.abs(spec),\"-k\") # plot abs of spectrum\n",
|
|
"plt.subplot(312)\n",
|
|
"plt.plot(frs,np.angle(spec),\"-k\") # plot phase of spectrum\n",
|
|
"plt.subplot(313)\n",
|
|
"plt.plot(frs,phase,\"-k\") # plot unwrapped phase\n",
|
|
"plt.xlabel(\"Frequency [Hz]\")\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"jupyter": {
|
|
"outputs_hidden": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"j1 = int(jmax*0.1)\n",
|
|
"j2 = int(jmax*0.6)\n",
|
|
"phavel = -2.*np.pi*frs[j1:j2]*dis/phase[j1:j2] # calculate phase velocity and plot\n",
|
|
"plt.plot(frs[j1:j2],phavel,'-k')\n",
|
|
"plt.xlabel(\"Frequency [Hz]\")\n",
|
|
"plt.ylabel(\"Phase velocity [km/s]\")\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": []
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "Python 3",
|
|
"language": "python",
|
|
"name": "python3"
|
|
},
|
|
"language_info": {
|
|
"codemirror_mode": {
|
|
"name": "ipython",
|
|
"version": 3
|
|
},
|
|
"file_extension": ".py",
|
|
"mimetype": "text/x-python",
|
|
"name": "python",
|
|
"nbconvert_exporter": "python",
|
|
"pygments_lexer": "ipython3",
|
|
"version": "3.6.7"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 4
|
|
}
|