449 lines
15 KiB
Plaintext
449 lines
15 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "17de1a34-2933-47ed-b1c5-ae9dd0466edd",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Assignment 7"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "681b47c6-f9af-4877-b745-febd443e1877",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Seismometer transfer function and correlation"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "85cb9941-eb24-4eb2-9c46-8610844e4063",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"%matplotlib inline\n",
|
|
"import numpy as np\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"from setupFigure import SetupFigure\n",
|
|
"from dftFast import dft_fast_coeff, dft_fast_synthesis\n",
|
|
"from seismicTrace import SeismicTrace"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "a627c30b-eb67-4319-930e-7f1acb8e3ee1",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Task 1: Write a function that calculates the transfer function of a seismometer for different damping $h=\\sin\\phi$."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "8cbc4a71-fc18-4c3e-a26a-0c02ae89314d",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def transfer_seismometer_velocity(fc, df, nf, phi, sigma):\n",
|
|
" \"\"\"\n",
|
|
" Calculates samples of the ground velocity seismometer transfer function (positive frequencies only)\n",
|
|
" :param fc: eigenfrequency of seismometer in Hz\n",
|
|
" :param df: frequency stepping\n",
|
|
" :param nf: number of frequencies\n",
|
|
" :param phi: damping angle in degrees (h=sin(phi))\n",
|
|
" :param sigma: generator constant\n",
|
|
" :return array with frequencies, array of complex values of transfer function\n",
|
|
" \"\"\"\n",
|
|
" f = df*np.arange(0, nf)\n",
|
|
" alfa = np.deg2rad(phi)\n",
|
|
" return f, -sigma*(f/fc)/(f/fc-np.exp(1j*alfa)) * (f/fc)/(f/fc-np.exp(1j*(np.pi-alfa)))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "7f7d857c-1bfd-4703-96a7-61c7af95bbdc",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### You may use this function for the Butterworth lowpass filter (or your own from Assignment 6)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "dd1b00e9-a953-4d64-b91f-7d7ad6949a3f",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def butterworth_lowpass(fc, df, nf, order):\n",
|
|
" \"\"\"\n",
|
|
" Calculates samples of the lowpass Butterworth filter transfer function (positive frequencies only)\n",
|
|
" :param fc: corner frequency in Hz\n",
|
|
" :param df: frequency stepping in Hz\n",
|
|
" :param nf: number of frequencies\n",
|
|
" :param order: filter order\n",
|
|
" \"\"\"\n",
|
|
" f = df*np.arange(0, nf)\n",
|
|
" h = np.ones(nf, dtype=np.complex64)\n",
|
|
" for k in range(order):\n",
|
|
" arg = 1j*(2*k+1)*np.pi/(2*order)\n",
|
|
" h = h*(-1j)/(f/fc-np.exp(arg))\n",
|
|
" return h"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "38d6be7b-6579-4726-abee-4666e6910f39",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### This is a useful function to plot traces next to each other with a vertical offset"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "48a7e17e-eaf2-4162-89f2-0d81eb9b7631",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def plotTraceGather(tracelist, ax):\n",
|
|
" \"\"\"\n",
|
|
" Plots a gather of traces with vertical offset\n",
|
|
" :param tracelist: list of SeismicTrace objects\n",
|
|
" :param ax: Matplotlib axis object\n",
|
|
" \"\"\"\n",
|
|
" for j, tr in enumerate(tracelist):\n",
|
|
" t = tr.tanf+tr.dt*np.arange(0, tr.nsamp)\n",
|
|
" norm = np.max(np.abs(tr.vals))\n",
|
|
" ax.plot(t, 0.5*tr.vals/norm+j, color='b', ls='-') # Normalization to 0.5*absmax and shift by multiples of 1\n",
|
|
" ax.text(tr.tanf, j+0.1, tr.staname, fontsize=14) # Add station name on top of trace"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "68bf389c-4464-43dc-853d-ff28e5791f99",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Task 2.1: Write a function that computes the cross correlation function for a given range of time delays"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "47a92e2f-da8c-45ce-b58c-7b8d153e38f5",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def crossCorrelation(x, y, dt, neglim, poslim):\n",
|
|
" \"\"\"\n",
|
|
" Correlate x with y according to c_n = Delta t*sum_(k=0)^(N-1) x_(n+k) y_k.\n",
|
|
" Assumes that both traces have same DT!\n",
|
|
" Consider x is a Gaussian with max at t2, y a Gaussian with max at t1 < t2.\n",
|
|
" To bring x and y into a match, x neeeds to be shifted to the left implying a positive n=n_s.\n",
|
|
" Thus, ns*dt is the delay of x relative to y or the advance of y relative to x.\n",
|
|
" If n+k < 0 or n+k >= x.size, x_(n+k) is considered as zero.\n",
|
|
" :param x: values of first trace\n",
|
|
" :param y: values of second trace\n",
|
|
" :param dt: sampling interval\n",
|
|
" :param neglim: most negative value for n\n",
|
|
" :param poslim: most positive value for n (poslim not included)\n",
|
|
" :return: array of delay times associated with the array elements of the cross correlation\n",
|
|
" :return: array of values of the cross correlation function\n",
|
|
" \"\"\"\n",
|
|
" nsamp = poslim-neglim\n",
|
|
" cc = np.zeros(nsamp)\n",
|
|
" for n in range(neglim, poslim):\n",
|
|
" nc = n-neglim # adjust index associating index 0 to neglim\n",
|
|
" for k in range(y.size):\n",
|
|
" if n+k >= 0 and n+k < x.size:\n",
|
|
" cc[nc] = cc[nc]+x[n+k]*y[k]\n",
|
|
" return dt*np.arange(neglim, poslim), cc*dt"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "79f44179-4943-4e87-9505-4c326c4237c2",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Task 2.2: Evaluate the amplitude spectrum of the seismometer transfer function for different angles\n",
|
|
"Use fc = 10 Hz, nf = 256, df = 0.2, sigma = 1. Choose angles 15, 25, 45, 60, 85 degrees. Plot the\n",
|
|
"amplitude spectra into one graph."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "56a730c7-ba3c-4404-a8e4-6d8c83d5fab3",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"fc = 10\n",
|
|
"nf = 256\n",
|
|
"df = 0.2\n",
|
|
"sigma = 1\n",
|
|
"phi = [15, 25, 45, 60, 85]\n",
|
|
"col = ['black', 'red', 'blue', 'orange', 'magenta']"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "dcfdc609-996d-425a-9750-f10cc453ef89",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"fig1, ax1 = SetupFigure(12, 5, \"Frequency [Hz]\", \"$|H(f)|$\", \"Amplitude of seismometer transfer function\", 14)\n",
|
|
"for j, p in enumerate(phi):\n",
|
|
" f, hseis = transfer_seismometer_velocity(fc, df, nf, p, sigma)\n",
|
|
" ax1.plot(f, np.absolute(hseis), color=col[j], ls='-', label=\"Phi = {:4.0f}\".format(p))\n",
|
|
"ax1.legend()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "b2576de7-c170-4829-a5a9-a107e6ad09bb",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Task 2.3: Evaluate the phase spectrum of the seismometer transfer function for different angles\n",
|
|
"Use fc = 10 Hz, nf = 256, df = 0.2, sigma = 1. Choose angles 15, 25, 45, 60, 85 degrees. Plot the\n",
|
|
"amplitude spectra into one graph."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "a796b315-f5e4-48cc-ae2e-d98323cc4996",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"fig2, ax2 = SetupFigure(12, 5, \"Frequency [Hz]\", \"Phase / PI\", \"Phase of seismometer transfer function\", 14)\n",
|
|
"for j, p in enumerate(phi):\n",
|
|
" f, hseis = transfer_seismometer_velocity(fc, df, nf, p, sigma)\n",
|
|
" ax2.plot(f, np.angle(hseis)/np.pi, color=col[j], ls='-', label=\"Phi = {:4.0f}\".format(p))\n",
|
|
"ax2.legend()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "d220be85-9530-418d-9069-6ad21d87872a",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Task 2.4: Setup two Gaussians $e^{-(t-\\mu)^2/\\sigma^2}$ of lengths 100 s and 80 s with dt=1 s, $\\mu$=45 and 55 and $\\sigma$ = 5. Plot the two into one graph. Then, compute the cross correlation between -30 and +30 s using your function from Task 2 and plot it into a separate graph. "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "f7537491-54c9-4f9c-ac22-e02d9cefae05",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"tx = np.arange(0,100)\n",
|
|
"ty = np.arange(0,80)\n",
|
|
"x = np.exp(-((tx-45)/5)**2)\n",
|
|
"y = np.exp(-((ty-55)/5)**2)\n",
|
|
"fig7, ax7 = SetupFigure(12, 5, \"Time [s]\", \"Amplitude\", \"Two Gaussians\", 14)\n",
|
|
"ax7.plot(tx, x, color='blue', ls='-', label='Time series x(t)')\n",
|
|
"ax7.plot(ty, y, color='red', ls='-', label='Time series y(t)')\n",
|
|
"ax7.legend()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "624805e5-2feb-4a42-bd91-5b48e798b01a",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"tc, cc = crossCorrelation(x, y, 1.0, -30, 30)\n",
|
|
"fig8, ax8 = SetupFigure(12, 5, \"Time delay [s]\", \"Correlation\", \"Cross correlation\", 14)\n",
|
|
"ax8.plot(tc, cc, color='blue', ls='-')"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "1529c72b-6ae3-416a-b4b2-ed1fd0bffba3",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Task 2.5: For comparison, use numpy's routine correlate to calculate the correlation. Try the different modes offered.\n",
|
|
"Change the lengths of x and y."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "da38f6a3-4d2e-4a69-a605-7de5a55f5cfd",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"fig9, ax9 = SetupFigure(12, 5, \"Time delay [s]\", \"Correlation\", \"Numpy cross correlation\", 14)\n",
|
|
"ccnp = np.correlate(x, y, mode='full')\n",
|
|
"tcnp = np.arange(-y.size+1, x.size) # This seems to be the correct relation between lag and index\n",
|
|
"print(x.size, y.size)\n",
|
|
"ax9.set_xlim(-30,30)\n",
|
|
"ax9.plot(tcnp, ccnp, color='red', ls='-')"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "5226e897-3874-45ba-86ab-883dd92fc13e",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Task 3.1: We want to prepare time delay measurements between two seismograms by cross-correlation. First, read the data of stations CH.PANIX and Z3.A261A.HHZ. Put the SeismicTrace objects into a list."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "0655529e-3212-4069-a60f-bef57644328c",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"stalist = [\"CH.PANIX.HHZ\", \"Z3.A261A.HHZ\"]\n",
|
|
"rawlist = []\n",
|
|
"for sta in stalist:\n",
|
|
" with open(sta, 'r') as fd:\n",
|
|
" tr = SeismicTrace.readFromAscii(fd)\n",
|
|
" tr.printInfo()\n",
|
|
" rawlist.append(tr)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "9f38af03-65f9-49d8-a28d-5fc0e2c4d965",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Task 3.2: Plot the traces using the function plotTraceGather."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "56551bfa-e2da-4773-9247-a33cae96a3c1",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"fig3, ax3 = SetupFigure(12, 5, \"Time [s]\", \"Displacement\", \"Raw traces\", 14)\n",
|
|
"plotTraceGather(rawlist, ax3)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "cec5430e-c602-4a99-a7c3-f3371e990d82",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Task 3.3: Filter the data using a Butterworth low pass with fc=0.5 Hz and order 6. Use either the provided function for the low-pass transfer function or take your own. First, Fourier transform the data, then multiply with the filter transfer function, and finally do the back transform. Create new SeismicTrace objects for the filtered traces, put them into a list and plot them."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "6ce5cc83-f5ff-4fb2-a9ac-d98595429c02",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"fc = 0.5; order = 6\n",
|
|
"filtlist = []\n",
|
|
"for tr in rawlist:\n",
|
|
" c = dft_fast_coeff(tr.vals)\n",
|
|
" df = 1./tr.tlen\n",
|
|
" hb = butterworth_lowpass(fc, df, c.size, order)\n",
|
|
" fts = c*hb\n",
|
|
" s = dft_fast_synthesis(fts, outnum='even')\n",
|
|
" trfil = SeismicTrace(tr.staname, tr.comp, tr.dt, tr.tanf, s)\n",
|
|
" filtlist.append(trfil)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "2fd73161-ade3-47e9-9880-774459625017",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"fig4, ax4 = SetupFigure(12, 5, \"Time [s]\", \"Displacement\", \"Low pass filtered traces, fc = {:5.2f} Hz, Order = {:d}\".format(fc, order), 14)\n",
|
|
"plotTraceGather(filtlist, ax4)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "ca21ecc8-0d60-46c4-8be8-c1e8fbe2c940",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"winlist = [[64285, 64300], [64280, 64295.1]]\n",
|
|
"cutlist = []\n",
|
|
"for j, tr in enumerate(filtlist):\n",
|
|
" win = winlist[j]\n",
|
|
" jfirst = int((win[0]-tr.tanf)/tr.dt)\n",
|
|
" jlast = int((win[1]-tr.tanf)/tr.dt)\n",
|
|
" nsamp = jlast-jfirst\n",
|
|
" tanew = jfirst*tr.dt+tr.tanf\n",
|
|
" print(tr.tanf, tanew, jfirst, jlast)\n",
|
|
" cuttr = SeismicTrace(tr.staname, tr.comp, tr.dt, tanew, tr.vals[jfirst:jlast])\n",
|
|
" cutlist.append(cuttr)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "f006344c-4728-4595-988e-d018bc88f994",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"fig5, ax5 = SetupFigure(12, 5, \"Time [s]\", \"Displacement\", \"Cut low pass filtered traces\".format(fc, order), 14)\n",
|
|
"plotTraceGather(cutlist, ax5)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "4993c95b-831d-490f-927c-cb38a2de1bfb",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"dt = cutlist[0].dt\n",
|
|
"neglim = -500; poslim = 500\n",
|
|
"t, cc = crossCorrelation(cutlist[0].vals, cutlist[1].vals, dt, neglim, poslim)\n",
|
|
"delay = (np.argmax(cc)+neglim)*dt+cutlist[0].tanf-cutlist[1].tanf\n",
|
|
"print(\"Tanf_0 = \", cutlist[0].tanf, \"Tanf_1 = \", cutlist[1].tanf, \"Index of max cc: \", np.argmax(cc)+neglim)\n",
|
|
"print(\"Delay of \", cutlist[0].staname,\" relative to \", cutlist[1].staname,\":\", delay)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "65cf3ace-f9ce-4da0-ac50-d348e7a129b3",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"fig6, ax6 = SetupFigure(12, 5, \"Time delay [s]\", \"Correlation\", \"Cross correlation\", 14)\n",
|
|
"ax6.plot(t, cc, color='blue', ls='-')"
|
|
]
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "Python 3 (ipykernel)",
|
|
"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.12.4"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 5
|
|
}
|