SurfaceWaveDispersion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main focus of this notebook are surface waves and their dispersion. \n",
    "Using synthetic seismograms generated with Instaseis (http://instaseis.net/), we try to recover the dispersion curves of Love and Rayleigh waves in the PREM model.\n",
    "\n",
    "Author: Thomas Ulrich"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import instaseis\n",
    "import obspy\n",
    "from obspy.taup.taup_geo import calc_dist\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import pycwt as wavelet\n",
    "import pandas as pd\n",
    "from scipy import interpolate\n",
    "\n",
    "# db = instaseis.open_db(\"syngine://prem_a_10s\")\n",
    "db = instaseis.open_db(\"syngine://prem_i_2s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first generate synthetic seismograms with Instaseis by specifying moment tensor, source and receiver locations.\n",
    "The source and receivers locations are fully arbitrary (we align source and receivers along a meridian for commodity)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Source\n",
    "lats = 89.91\n",
    "lons = 74.4940\n",
    "# receiver\n",
    "latr = 42.6390\n",
    "lonr = 74.4940"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "receiver = instaseis.Receiver(latitude=latr, longitude=lonr, network=\"AB\", station=\"CED\")\n",
    "t0 = obspy.UTCDateTime(2011, 1, 2, 3, 4, 5)\n",
    "source = instaseis.Source(\n",
    "    latitude=lats,\n",
    "    longitude=lons,\n",
    "    depth_in_m=12000,\n",
    "    m_rr=4.710000e24 / 1e7,\n",
    "    m_tt=3.810000e22 / 1e7,\n",
    "    m_pp=-4.740000e24 / 1e7,\n",
    "    m_rt=3.990000e23 / 1e7,\n",
    "    m_rp=-8.050000e23 / 1e7,\n",
    "    m_tp=-1.230000e24 / 1e7,\n",
    "    origin_time=t0,\n",
    ")\n",
    "st = db.get_seismograms(\n",
    "    source=source, receiver=receiver, components=[\"Z\", \"R\", \"T\"], kind=\"displacement\"\n",
    ")\n",
    "# convert to micro_m\n",
    "for i in range(3):\n",
    "    st[i].data *= 1e6\n",
    "st.plot(show=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Group Velocity\n",
    "\n",
    "If dispersion occurs, waves of different wavelengths travel at different speeds.  \n",
    "In the following, we band-pass filter the seismograms around various central periods, ranging from 5 to 240s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "iZRT = 1\n",
    "dt = st[iZRT].stats.delta\n",
    "duration = st[iZRT].stats.npts * dt\n",
    "myPeriods = np.logspace(np.log2(10.0 * dt), np.log2(duration / 15), 10, base=2)\n",
    "\n",
    "nper = len(myPeriods)\n",
    "\n",
    "figall, axarr = plt.subplots(nper + 1, 1, figsize=(14, 9), dpi=160, sharex=True, sharey=False)\n",
    "\n",
    "\n",
    "ti = st[iZRT].times(reftime=t0)\n",
    "axarr[0].plot(ti, st[iZRT].data, label=\"no filter\")\n",
    "axarr[0].legend(loc=\"upper left\")\n",
    "\n",
    "pickedArrival = np.zeros((nper,))\n",
    "for i, T in enumerate(myPeriods):\n",
    "    st_temp = st.copy()\n",
    "    st_temp.filter(\"bandpass\", freqmin=0.85 / T, freqmax=1.15 / T, corners=4, zerophase=True)\n",
    "    axarr[i + 1].plot(\n",
    "        ti,\n",
    "        st_temp[iZRT].data,\n",
    "        label=f\"T={T:.1f}s\",\n",
    "    )\n",
    "    # Envelope of filtered data\n",
    "    data_envelope = obspy.signal.filter.envelope(st_temp[iZRT].data)\n",
    "    axarr[i + 1].plot(ti, data_envelope.data)\n",
    "    idmax = np.argmax(data_envelope.data)\n",
    "    pickedArrival[i] = ti[idmax]\n",
    "    axarr[i + 1].legend(loc=\"upper left\")\n",
    "\n",
    "\n",
    "axarr[0].set_xlim([0, 2200])\n",
    "\n",
    "axarr[nper].set_xlabel(\"time (s)\")\n",
    "axarr[nper // 2].set_ylabel(\"displacement (um)\")\n",
    "plt.show()\n",
    "\n",
    "d = {\"filter_Tc\": myPeriods, \"picked_arrival\": pickedArrival}\n",
    "df = pd.DataFrame(data=d)\n",
    "pd.options.display.float_format = \"{:.0f}\".format\n",
    "print(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The idea of using a narrow filter to pick the arrival is imprecise and challenging with real data. Instead, a wavelet transform approach allows tracking more accurately the central frequency with time (\"A wavelet is a wave-like oscillation with an amplitude that begins at zero, increases or decreases, and then returns to zero one or more times\", https://en.wikipedia.org/wiki/Wavelet). Such a solution is implemented below (inspired from http://regeirk.github.io/pycwt/tutorial.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_wavelet_power_spectrum(mydata, dt):\n",
    "    # we normalize the data\n",
    "    std = mydata.std()\n",
    "    dat_norm = mydata / std\n",
    "\n",
    "    mother = wavelet.Morlet(6)\n",
    "\n",
    "    s0 = 2 * dt  # Starting scale\n",
    "    dj = 1 / 12.0  # Twelve sub-octaves per octaves\n",
    "    J = 8.0 / dj  # Height powers of two with dj sub-octaves\n",
    "\n",
    "    # wave is shaped (nfreq, nti), with nfreq = J + 1 = 8 * 12 +1,\n",
    "    # and nti the number of samples of the signal\n",
    "    wave, scales, freqs, coi, fft, fftfreqs = wavelet.cwt(dat_norm, dt, dj, s0, J, mother)\n",
    "    iwave = wavelet.icwt(wave, scales, dt, dj, mother) * std\n",
    "    period = 1.0 / freqs\n",
    "\n",
    "    # compute power spectrum\n",
    "    power = (np.abs(wave)) ** 2\n",
    "    # rectify the power spectrum according to the suggestions proposed by Liu et al. (2007)\n",
    "    power /= scales[:, None]\n",
    "    return period, power, wave, coi\n",
    "\n",
    "\n",
    "periodsWT, power, wave, coi = compute_wavelet_power_spectrum(st[iZRT].data, dt)\n",
    "\n",
    "# Plot wavelet power spectrum figure\n",
    "figall, ax = plt.subplots()\n",
    "levels = [16 / 2**i for i in range(10, -1, -1)]\n",
    "ax.contourf(\n",
    "    ti, np.log2(periodsWT), np.log2(power), np.log2(levels), extend=\"both\", cmap=plt.cm.viridis\n",
    ")\n",
    "\n",
    "if False:\n",
    "    # change to True to also plot the \"cone of influence\" as a hatched area\n",
    "    ax.fill(\n",
    "        np.concatenate([ti, ti[-1:] + dt, ti[-1:] + dt, ti[:1] - dt, ti[:1] - dt]),\n",
    "        np.concatenate(\n",
    "            [np.log2(coi), [1e-9], np.log2(periodsWT[-1:]), np.log2(periodsWT[-1:]), [1e-9]]\n",
    "        ),\n",
    "        \"k\",\n",
    "        alpha=0.3,\n",
    "        hatch=\"x\",\n",
    "    )\n",
    "\n",
    "\n",
    "ax.set_title(f\"Wavelet Power Spectrum\")\n",
    "ax.set_xlabel(\"time (s)\")\n",
    "ax.set_ylabel(\"period (s)\")\n",
    "ax.grid()\n",
    "Yticks = 2 ** np.arange(np.ceil(np.log2(periodsWT.min())), np.ceil(np.log2(periodsWT.max())))\n",
    "\n",
    "ax.set_yticks(np.log2(Yticks))\n",
    "ax.set_yticklabels(Yticks)\n",
    "\n",
    "plt.xlim(0, 2200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pickArrivalWavelet(ti, power, period, tmin=1050, tmax=3000):\n",
    "    # pick arrival time from the wavelet Power Spectrum data\n",
    "    maxs = np.amax(power, axis=1)\n",
    "    # enforce that the arrival time should be greater than tmin\n",
    "    idtmin = np.argmin(np.abs(ti - tmin))\n",
    "    idtmax = np.argmin(np.abs(ti - tmax))\n",
    "    idmaxs = np.argmax(power[:, idtmin:idtmax], axis=1) + idtmin\n",
    "\n",
    "    pickedArrivalWavelet = np.zeros_like(period)\n",
    "    for i in range(power.shape[0]):\n",
    "        pickedArrivalWavelet[i] = ti[idmaxs[i]]\n",
    "    return pickedArrivalWavelet\n",
    "\n",
    "\n",
    "pickedArrivalWavelet = pickArrivalWavelet(ti, power, periodsWT)\n",
    "plt.plot(periodsWT, pickedArrivalWavelet)\n",
    "plt.xlabel(\"period (s)\")\n",
    "plt.ylabel(\"picked time (s)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that for period in the range 30 to 60s, there is energy arriving very early (between 500s and 1000s following the earthquake onset). We need to manually enforce that the picked time is greater than 1050s, to avoid picking wrong arrivals.   \n",
    "Next, we compute the distance between source and receivers using obspy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_distance_source_receiver(lons, lats, lonr, latr):\n",
    "    dist_degree_sr = calc_dist(\n",
    "        source_latitude_in_deg=lats,\n",
    "        source_longitude_in_deg=lons,\n",
    "        receiver_latitude_in_deg=latr,\n",
    "        receiver_longitude_in_deg=lonr,\n",
    "        radius_of_planet_in_km=6371,\n",
    "        flattening_of_planet=0,\n",
    "    )\n",
    "    dist_km_sr = 6371 * dist_degree_sr * 2 * np.pi / 360.0\n",
    "    return dist_degree_sr, dist_km_sr\n",
    "\n",
    "\n",
    "dist_degree_sr, dist_km_sr = compute_distance_source_receiver(lons, lats, lonr, latr)\n",
    "print(f\"distance source-receivers: {dist_km_sr:.1f} km ({dist_degree_sr:.1f} deg)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we plot the inferred dispersion curves and compare them with theoretical estimates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "group_velocity_from_wavelet = dist_km_sr / pickedArrivalWavelet\n",
    "plt.plot(myPeriods, dist_km_sr / pickedArrival, \"x\", label=\"this notebook\")\n",
    "plt.plot(periodsWT, group_velocity_from_wavelet, \"o\", label=\"this notebook, wavelet\")\n",
    "plt.xlabel(\"period (s)\")\n",
    "plt.ylabel(\"group velocity (km/s)\")\n",
    "\n",
    "surfaceType = \"Love\" if iZRT == 2 else \"Rayleigh\"\n",
    "PREM_groupVel = np.loadtxt(f\"PREM_groupVel{surfaceType}.dat\", skiprows=1, delimiter=\",\")\n",
    "\n",
    "plt.plot(PREM_groupVel[:, 0], PREM_groupVel[:, 1], label=\"PREM (G. Laske)\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Phase velocity\n",
    "\n",
    "The frequency-dependent phase of a seismogram can be obtained by taking the Fourrier transform of it. It can be written as follows: \n",
    "\n",
    "\\begin{equation*}\n",
    "\\Phi(\\omega) = \\omega t - k(\\omega) x + \\Phi_i(\\omega) + 2n\\pi = \\omega t - \\omega x / c(\\omega)+ \\Phi_i(\\omega) + 2n\\pi\n",
    "\\end{equation*}\n",
    "\n",
    "Where $\\omega t - \\omega x / c(\\omega)$ is the phase due to the propagation in space and time,  \n",
    "$\\Phi_i(\\omega)$ depends on the source and  \n",
    "$2n\\pi$ reflects the periodicity of the phase.  \n",
    "To compute the phase velocity we need to track the same wave at two receivers.  \n",
    "Assuming that the receivers are at distance $x_1$ and $x_2$ from the source, and the wave passes at time $t_1$ and $t_2$ at each receivers, then the phase of both receivers would be:  \n",
    "\n",
    "\\begin{equation*}\n",
    "\\Phi_1(\\omega) =   \\omega t_1 - \\omega x_1 / c(\\omega)+ \\Phi_i(\\omega) + 2n\\pi\n",
    "\\end{equation*}\n",
    "and  \n",
    "\\begin{equation*}\n",
    "\\Phi_2(\\omega) = \\omega t_2 - \\omega x_2 / c(\\omega)+ \\Phi_i(\\omega) + 2m\\pi\n",
    "\\end{equation*}\n",
    "\n",
    "Then c can be obtained from the difference $\\Phi_{21}(\\omega)$ between $\\Phi_2(\\omega)$ and $\\Phi_1(\\omega)$:\n",
    "\\begin{equation*}\n",
    "c(\\omega) = \\omega (x_2-x_1)/(\\omega (t_2-t_1)+ 2(m-n)\\pi- \\Phi_{21}(\\omega))\n",
    "\\end{equation*}\n",
    "the term m-n is found empirically ensuring that the phase velocity is reasonable for long periods.  \n",
    "source: Stein and Wysession."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "receiver2 = instaseis.Receiver(latitude=latr - 6, longitude=lonr)\n",
    "st2 = db.get_seismograms(\n",
    "    source=source, receiver=receiver2, components=[\"Z\", \"R\", \"T\"], kind=\"displacement\"\n",
    ")\n",
    "# convert to micro_m\n",
    "for i in range(3):\n",
    "    st2[i].data *= 1e6\n",
    "st.plot()\n",
    "st2.plot(show=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To facilitate comparison of the phase of the signals at the 2 stations, we align both signals by compensating for the across stations travel time of the wavepacket. We estimate it period-wise by interpolating the group slowness."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_degree_sr_2, dist_km_sr_2 = compute_distance_source_receiver(lons, lats, lonr, latr - 6)\n",
    "dist_km_src12 = dist_km_sr_2 - dist_km_sr\n",
    "print(dist_km_src12)\n",
    "f_group_slowness = interpolate.interp1d(periodsWT, 1.0 / group_velocity_from_wavelet)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The phase differences (and associated time shifts) between the aligned signals is associated with the difference between phase and group velocity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "figall, axarr = plt.subplots(nper + 1, 1, figsize=(14, 18), dpi=160, sharex=False, sharey=False)\n",
    "\n",
    "axarr[0].plot(ti, st[iZRT].data, label=\"no filter, r1\")\n",
    "axarr[0].plot(ti, st2[iZRT].data, label=\"no filter, r2\")\n",
    "axarr[0].legend()\n",
    "\n",
    "pickedArrival = np.zeros((nper, 3))\n",
    "pickedArrival[:, 0] = myPeriods\n",
    "\n",
    "for i, T in enumerate(myPeriods):\n",
    "    for k, sti in enumerate([st, st2]):\n",
    "        st_temp = sti.copy()\n",
    "        st_temp.filter(\"bandpass\", freqmin=0.85 / T, freqmax=1.15 / T, corners=4, zerophase=True)\n",
    "        # align signals given interpolated group velocity\n",
    "        shift = 0 if k == 0 else f_group_slowness(T) * dist_km_src12\n",
    "        axarr[i + 1].plot(\n",
    "            ti - shift,\n",
    "            st_temp[iZRT].data,\n",
    "            label=f\"T={T:.1f}s\",\n",
    "        )\n",
    "        # Envelope of filtered data\n",
    "        data_envelope = obspy.signal.filter.envelope(st_temp[iZRT].data)\n",
    "        axarr[i + 1].plot(ti - shift, data_envelope.data)\n",
    "        idmax = np.argmax(data_envelope.data)\n",
    "        pickedArrival[i, k + 1] = st_temp[iZRT].times(reftime=t0)[idmax]\n",
    "        axarr[i + 1].legend(loc=\"upper left\")\n",
    "        # center the plot on the max energy of the filtered signal\n",
    "        tc = f_group_slowness(T) * dist_km_sr\n",
    "        axarr[i + 1].set_xlim(tc - 8 * T, tc + 8 * T)\n",
    "\n",
    "axarr[-1].set_xlabel(\"time (s)\")\n",
    "axarr[nper // 2].set_ylabel(\"displacement (um)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Phase differences can be evaluated with cross-correlation (and is sensitive to phase shift). Below we provide a simple code estimating phase velocity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from obspy.signal.cross_correlation import correlate\n",
    "\n",
    "myPeriods2 = np.logspace(np.log2(10.0 * dt), np.log2(duration / 15), 40, base=2)\n",
    "phase_vel = np.zeros_like(myPeriods2)\n",
    "\n",
    "nsamples = 1000\n",
    "atshift = np.linspace(-nsamples, nsamples, 2 * nsamples + 1) * dt\n",
    "cc_data = np.zeros((myPeriods2.shape[0], 2 * nsamples + 1))\n",
    "\n",
    "# filter and cross-correlate both signals\n",
    "x_filter = 0.25\n",
    "for i, T in enumerate(myPeriods2):\n",
    "    st_temp = st.copy()\n",
    "    st_temp2 = st2.copy()\n",
    "    for sti in [st_temp, st_temp2]:\n",
    "        sti.filter(\n",
    "            \"bandpass\",\n",
    "            freqmin=(1.0 - x_filter) / T,\n",
    "            freqmax=(1.0 + x_filter) / T,\n",
    "            corners=4,\n",
    "            zerophase=True,\n",
    "        )\n",
    "    cc = correlate(st_temp2[iZRT], st_temp[iZRT], nsamples)\n",
    "    cc_data[i, :] = cc\n",
    "\n",
    "\n",
    "figall, ax = plt.subplots()\n",
    "Tc = 80\n",
    "idT = np.argmin(np.abs(myPeriods2 - Tc))\n",
    "plt.plot(atshift, cc_data[idT, :])\n",
    "plt.ylabel(\"cross correlation amplitude\")\n",
    "plt.xlabel(\" time shift(s)\")\n",
    "plt.title(f\"cross-correlogram for central period {Tc}s\")\n",
    "plt.show()\n",
    "\n",
    "\n",
    "# plot cross correlation array\n",
    "figall, ax = plt.subplots()\n",
    "X, Y = np.meshgrid(myPeriods2, atshift)\n",
    "plt.pcolormesh(X, Y, cc_data.T, cmap=\"seismic\")\n",
    "plt.xlabel(\"period (s)\")\n",
    "plt.ylabel(\" time shift(s)\")\n",
    "plt.title(\"cross-correlation array\")\n",
    "figall, ax = plt.subplots()\n",
    "\n",
    "# Plot only possible values\n",
    "Vmax = 7.0\n",
    "id_realistic = np.where(atshift > dist_km_src12 / Vmax)[0]\n",
    "atshift = atshift[id_realistic]\n",
    "cc_data = cc_data[:, id_realistic]\n",
    "V_phase = dist_km_src12 / atshift\n",
    "\n",
    "# generate 2D plot\n",
    "X, Y = np.meshgrid(myPeriods2, V_phase)\n",
    "plt.pcolormesh(X, Y, cc_data.T, cmap=\"seismic\")\n",
    "# Plot reference data\n",
    "surfaceType = \"Love\" if iZRT == 2 else \"Rayleigh\"\n",
    "PREM_phaseVel = np.loadtxt(f\"PREM_phaseVel{surfaceType}.dat\", skiprows=1, delimiter=\",\")\n",
    "plt.plot(PREM_phaseVel[:, 0], PREM_phaseVel[:, 1], label=\"PREM (G. Laske)\")\n",
    "\n",
    "\n",
    "plt.xlabel(\"period (s)\")\n",
    "plt.ylabel(\"phase velocity (km/s)\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}