288 lines
9.6 KiB
Plaintext
288 lines
9.6 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "17de1a34-2933-47ed-b1c5-ae9dd0466edd",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Assignment 8"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "681b47c6-f9af-4877-b745-febd443e1877",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Correlation picking"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 14,
|
|
"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": "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": 15,
|
|
"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": 16,
|
|
"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": [
|
|
"#### You may use this function that computes the cross correlation function for a given range of time delays"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 17,
|
|
"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": "3404b2f6-d296-4b3e-8234-a01c2bd061ea",
|
|
"metadata": {},
|
|
"source": [
|
|
"### We repeat the first parts of task 3 from the previous assignment to make this assignment self-contained."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "5226e897-3874-45ba-86ab-883dd92fc13e",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Task 1.1: We want to do 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": "markdown",
|
|
"id": "b8182997-34f6-4091-8f71-caa942e9d1f2",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Task 1.2: Plot the traces using the function plotTraceGather. "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "951456c6-a3d2-4737-a803-c88d813655e4",
|
|
"metadata": {},
|
|
"source": [
|
|
"This function plots the seismogram values versus time in seconds after midnight."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "db6ea78c-03ec-4e56-88c5-c7dc58efed9c",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Task 1.3: Filter the data using a Butterworth low pass with fc=0.5 Hz and order 6. "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "4f3c59ab-69b7-4ab0-bedc-e525134a687f",
|
|
"metadata": {},
|
|
"source": [
|
|
"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": "markdown",
|
|
"id": "3f1990b2-9684-4a1e-a678-448ebd0b4afe",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Task 1.4: Cut out the P-wave train from the recordings. Try different time windows for each trace. "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "e619d197-db9b-48ef-840d-aa123b5ff884",
|
|
"metadata": {},
|
|
"source": [
|
|
"As a start, use the times [64285, 64300] for CH.PANIX and the times [64280, 64295] for Z3.A261A. For cutting out samples, proceed as follows: (1) Calculate the sample indices for the beginning and end of the window, (2) calculate the new starting time, (3) create a new SeismicTrace object with same station name, component and DT but adjusted start time (tanf) and the data values within the\n",
|
|
"precalculated index range. (4) Again, put the windowed seismograms into a list."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "058f0df7-d40a-4a57-bec0-052fa2eb1c86",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Task 1.5: Plot the windowed traces using plotTraceGather."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "59ff1b7f-40f9-42d4-9354-10b14a3516bc",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Task 1.6: Cross-correlate the windowed seismograms and compute the time delay"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "175aabb2-a1b0-4c85-a84c-7f1e84ed2546",
|
|
"metadata": {},
|
|
"source": [
|
|
"Compute the correlation for delays from -500 to +500 samples. To get the time delay, search the index jmax of the maximum of the cross-correlation. Since sample zero of the cross-correlation belongs to a time of -500*dt, the time delay is given by $\\tau=(jmax-500)dt$. This is correct if the start times of the traces are identical. If not, we must correct for the start time difference. For example, if you correlate PANIX with A261A in this order, you obtain the delay of PANIX relative to A261A which is given by $\\tau$ + ta(PANIX)-ta(A261A).\n",
|
|
"\n",
|
|
"Print the index of the CC maximum, print the start times of the traces and print the time delay. Does the result meet your expectations from the figure of Task 1.5?"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "32263d3c-e959-4644-9fee-20c42c8b5fdb",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Task 1.7: Plot the cross-correlation function"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "0a9469e7-7d61-4d5e-8188-9259d982d60f",
|
|
"metadata": {},
|
|
"source": [
|
|
"Use the true delay time on the x-axis of your plot. "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "a87c81e1-5c21-4b27-83e7-73712d7b2be3",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Task 1.8: Compute the energy and rms of the windowed data of each trace"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "33e4046e-b168-4105-a861-3e1d9c7f903f",
|
|
"metadata": {},
|
|
"source": [
|
|
"The energy of the trace is given by E=$\\int s^2(t)dt$. The square root of the energy is the root-mean-square (rms)."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "a19552cd-14d4-4d60-a3d3-6ee13fe94f2b",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Task 1.9: Normalize the cross correlation function by the rms values of the correlated traces and plot it."
|
|
]
|
|
}
|
|
],
|
|
"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
|
|
}
|