{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib notebook" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\nStep and impulse responses\n==========================\n\nThese examples compare the analytical solution with `empymod` for time-domain\nstep and impulse responses for inline, x-directed source and receivers, for the\nfour different frequency-to-time methods **QWE**, **DLF**, **FFTLog**, and\n**FFT**. Which method is faster and which is more precise depends on the model\n(land or marine, source/receiver at air-interface or not) and the response\n(step or impulse).\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import empymod\nimport numpy as np\nfrom scipy.special import erf\nimport matplotlib.pyplot as plt\nfrom scipy.constants import mu_0 # Permeability of free space [H/m]\nplt.style.use('ggplot')\ncolors = [color['color'] for color in list(plt.rcParams['axes.prop_cycle'])]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Analytical solutions\n--------------------\n\nAnalytical solution for source and receiver at the interface between two\nhalf-spaces\n\nThe time-domain step and impulse responses for a source at the origin\n($x_s = y_s = z_s = 0$ m) and an in-line receiver at the surface\n($y_r = z_r = 0$ m), is given by the following equations, where\n$\\rho_h$ is horizontal resistivity ($\\Omega$ m),\n$\\lambda$ is anisotropy (-), with $\\lambda =\n\\sqrt{\\rho_v/\\rho_h}$, $r$ is offset (m), $t$ is time (s), and\n$\\tau_h = \\sqrt{\\mu_0 r^2/(\\rho_h t)}$; $\\mu_0$ is the magnetic\npermeability of free space (H/m).\n\nTime Domain: Step Response $\\mathbf{\\mathcal{H}(t)}$\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n\\begin{align}E_x(\\rho_h,\\lambda,r,t) = \\frac{\\rho_h}{2 \\pi r^3} \\left[ 2\\lambda +\n \\rm{erf}\\left(\\frac{\\tau_h}{2}\\right) - 2\\lambda\n \\rm{erf}\\left(\\frac{\\tau_h}{2\\lambda}\\right) + \\frac{\\tau_h}{\\sqrt{\\pi}}\n \\exp\\left(- \\frac{\\tau_h^2}{4\\lambda^2}\\right) \\right]\\end{align}\n\nTime Domain: Impulse Response $\\mathbf{\\delta(t)}$\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n\\begin{align}\\dot{E}_x(\\rho_h,\\lambda,r,t) =\n \\frac{\\rho_h}{2 \\pi r^3} \\left[ \\delta(t) + \\frac{\\tau_h}{2t\\sqrt{\\pi}}\n \\left\\{ - \\exp\\left(-\\frac{\\tau_h^2}{4}\\right) +\n \\left( \\frac{\\tau_h^2}{2 \\lambda^2} + 1 \\right) \\exp\\left(-\n \\frac{\\tau_h^2}{4 \\lambda^2}\\right) \\right\\} \\right]\\end{align}\n\nReference\n~~~~~~~~~\nEquations 3.2 and 3.3 in Werthm\u00fcller, D., 2009, Inversion of multi-transient\nEM data from anisotropic media: M.S. thesis, TU Delft, ETH Z\u00fcrich, RWTH\nAachen; UUID: `f4b071c1-8e55-4ec5-86c6-a2d54c3eda5a\n`_.\n\nAnalytical functions\n~~~~~~~~~~~~~~~~~~~~\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def ee_xx_impulse(res, aniso, off, time):\n \"\"\"VTI-Halfspace impulse response, xx, inline.\n\n res : horizontal resistivity [Ohm.m]\n aniso : anisotropy [-]\n off : offset [m]\n time : time(s) [s]\n \"\"\"\n tau_h = np.sqrt(mu_0*off**2/(res*time))\n t0 = tau_h/(2*time*np.sqrt(np.pi))\n t1 = np.exp(-tau_h**2/4)\n t2 = tau_h**2/(2*aniso**2) + 1\n t3 = np.exp(-tau_h**2/(4*aniso**2))\n Exx = res/(2*np.pi*off**3)*t0*(-t1 + t2*t3)\n Exx[time == 0] = res/(2*np.pi*off**3) # Delta dirac part\n return Exx\n\n\ndef ee_xx_step(res, aniso, off, time):\n \"\"\"VTI-Halfspace step response, xx, inline.\n\n res : horizontal resistivity [Ohm.m]\n aniso : anisotropy [-]\n off : offset [m]\n time : time(s) [s]\n \"\"\"\n tau_h = np.sqrt(mu_0*off**2/(res*time))\n t0 = erf(tau_h/2)\n t1 = 2*aniso*erf(tau_h/(2*aniso))\n t2 = tau_h/np.sqrt(np.pi)*np.exp(-tau_h**2/(4*aniso**2))\n Exx = res/(2*np.pi*off**3)*(2*aniso + t0 - t1 + t2)\n return Exx" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Example 1: Source and receiver at z=0m\n--------------------------------------\n\nComparison with analytical solution; put 1 mm below the interface, as they\nwould be regarded as in the air by `emmod` otherwise.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "src = [0, 0, 0.001] # Source at origin, slightly below interface\nrec = [6000, 0, 0.001] # Receivers in-line, 0.5m below interface\nres = [2e14, 10] # Resistivity: [air, half-space]\naniso = [1, 2] # Anisotropy: [air, half-space]\neperm = [0, 1] # Set el. perm. of air to 0 because of num. noise\nt = np.logspace(-2, 1, 301) # Desired times (s)\n\n# Collect parameters\ninparg = {'src': src, 'rec': rec, 'depth': 0, 'freqtime': t, 'res': res,\n 'aniso': aniso, 'epermH': eperm, 'ht': 'dlf', 'verb': 2}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Impulse response\n~~~~~~~~~~~~~~~~\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ex = ee_xx_impulse(res[1], aniso[1], rec[0], t)\n\ninparg['signal'] = 0 # signal 0 = impulse\nprint('QWE')\nqwe = empymod.dipole(**inparg, ft='qwe')\nprint('DLF (Sine)')\nsin = empymod.dipole(**inparg, ft='dlf', ftarg={'dlf': 'key_81_CosSin_2009'})\nprint('FFTLog')\nftl = empymod.dipole(**inparg, ft='fftlog')\nprint('FFT')\nfft = empymod.dipole(\n **inparg, ft='fft',\n ftarg={'dfreq': .0005, 'nfreq': 2**20, 'pts_per_dec': 10})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "=> `FFTLog` is the fastest by quite a margin, followed by the `Sine`-filter.\nWhat cannot see from the output (set `verb` to something bigger than 2 to see\nit) is how many frequencies each method used:\n\n- QWE: 159 (0.000794328 - 63095.7 Hz)\n- Sine: 116 (5.33905E-06 - 52028 Hz)\n- FFTLog: 60 (0.000178575 - 141.847 Hz)\n- FFT: 61 (0.0005 - 524.288 Hz)\n\nNote that for the actual transform, `FFT` used 2^20 = 1'048'576 frequencies!\nIt only computed 60 frequencies, and then interpolated the rest, as it\nrequires regularly spaced data.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure()\nplt.title(r'Impulse response for HS-model, $r=$' +\n str(int(rec[0]/1000)) + ' km.')\nplt.xlabel('Time (s)')\nplt.ylabel(r'Amplitude (V/m)')\nplt.semilogx(t, ex, 'k-', label='Analytical')\nplt.semilogx(t, qwe, 'C0-', label='QWE')\nplt.semilogx(t, sin, 'C1--', label='Sine Filter')\nplt.semilogx(t, ftl, 'C2-.', label='FFTLog')\nplt.semilogx(t, fft, 'C3:', label='FFT')\nplt.legend(loc='best')\nplt.ylim([-.1*np.max(ex), 1.1*np.max(ex)])\nplt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure()\nplt.title('Error')\nplt.xlabel('Time (s)')\nplt.ylabel('Relative error (-)')\nplt.loglog(t, abs(qwe-ex)/ex, 'C0-', label='QWE')\nplt.plot(t, abs(sin-ex)/ex, 'C1--', label='Sine Filter')\nplt.plot(t, abs(ftl-ex)/ex, 'C2-.', label='FFTLog')\nplt.plot(t, abs(fft-ex)/ex, 'C3:', label='FFT')\nplt.legend(loc='best')\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "=> The error is comparable in all cases. `FFT` is not too good at later\ntimes. This could be improved by computing lower frequencies. But because FFT\nneeds regularly spaced data, our vector would soon explode (and you would\nneed a lot of memory). In the current case we are already using 2^20 samples!\n\nStep response\n~~~~~~~~~~~~~\n\nStep responses are almost impossible with `FFT`. We can either try to model\nlate times with lots of low frequencies, or the step with lots of high\nfrequencies. I do not use `FFT` in the step-response examples.\n\nSwitch-on\n'''''''''\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ex = ee_xx_step(res[1], aniso[1], rec[0], t)\n\ninparg['signal'] = 1 # signal 1 = switch-on\nprint('QWE')\nqwe = empymod.dipole(**inparg, ft='qwe')\nprint('DLF (Sine)')\nsin = empymod.dipole(**inparg, ft='dlf', ftarg={'dlf': 'key_81_CosSin_2009'})\nprint('FFTLog')\nftl = empymod.dipole(**inparg, ft='fftlog', ftarg={'q': -0.6})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Used number of frequencies:\n\n- QWE: 159 (0.000794328 - 63095.7 Hz)\n- Sine: 116 (5.33905E-06 - 52028 Hz)\n- FFTLog: 60 (0.000178575 - 141.847 Hz)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure()\nplt.title(r'Switch-on response for HS-model, $r=$' +\n str(int(rec[0]/1000)) + ' km.')\nplt.xlabel('Time (s)')\nplt.ylabel('Amplitude (V/m)')\nplt.semilogx(t, ex, 'k-', label='Analytical')\nplt.semilogx(t, qwe, 'C0-', label='QWE')\nplt.semilogx(t, sin, 'C1--', label='Sine Filter')\nplt.semilogx(t, ftl, 'C2-.', label='FFTLog')\nplt.legend(loc='best')\nplt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure()\nplt.title('Error')\nplt.xlabel('Time (s)')\nplt.ylabel('Relative error (-)')\nplt.loglog(t, abs(qwe-ex)/ex, 'C0-', label='QWE')\nplt.plot(t, abs(sin-ex)/ex, 'C1--', label='Sine Filter')\nplt.plot(t, abs(ftl-ex)/ex, 'C2-.', label='FFTLog')\nplt.legend(loc='best')\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Switch-off\n''''''''''\nFor switch-off to work properly you need `empymod`-version bigger than 1.3.0!\nYou can do it with previous releases too, but you will have to do the\nDC-computation and subtraction manually, as is done here for `ee_xx_step`.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "exDC = ee_xx_step(res[1], aniso[1], rec[0], 60*60)\nex = exDC - ee_xx_step(res[1], aniso[1], rec[0], t)\n\ninparg['signal'] = -1 # signal -1 = switch-off\nprint('QWE')\nqwe = empymod.dipole(**inparg, ft='qwe')\nprint('DLF (Cosine/Sine)')\nsin = empymod.dipole(**inparg, ft='dlf', ftarg={'dlf': 'key_81_CosSin_2009'})\nprint('FFTLog')\nftl = empymod.dipole(**inparg, ft='fftlog', ftarg={'add_dec': [-5, 3]})" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure()\nplt.title(r'Switch-off response for HS-model, $r=$' +\n str(int(rec[0]/1000)) + ' km.')\nplt.xlabel('Time (s)')\nplt.ylabel('Amplitude (V/m)')\nplt.semilogx(t, ex, 'k-', label='Analytical')\nplt.semilogx(t, qwe, 'C0-', label='QWE')\nplt.semilogx(t, sin, 'C1--', label='Cosine/Sine Filter')\nplt.semilogx(t, ftl, 'C2-.', label='FFTLog')\nplt.legend(loc='best')\nplt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure()\nplt.title('Error')\nplt.xlabel('Time (s)')\nplt.ylabel('Relative error (-)')\nplt.loglog(t, abs(qwe-ex)/ex, 'C0-', label='QWE')\nplt.plot(t, abs(sin-ex)/ex, 'C1--', label='Sine Filter')\nplt.plot(t, abs(ftl-ex)/ex, 'C2-.', label='FFTLog')\nplt.legend(loc='best')\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Example 2: Air-seawater-halfspace\n---------------------------------\n\nIn seawater the transformation is generally much easier, as we do not have\nthe step or the impules at zero time.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "src = [0, 0, 950] # Source 50 m above seabottom\nrec = [6000, 0, 1000] # Receivers in-line, at seabottom\nres = [1e23, 1/3, 10] # Resistivity: [air, water, half-space]\naniso = [1, 1, 2] # Anisotropy: [air, water, half-space]\nt = np.logspace(-2, 1, 301) # Desired times (s)\n\n# Collect parameters\ninparg = {'src': src, 'rec': rec, 'depth': [0, 1000], 'freqtime': t,\n 'res': res, 'aniso': aniso, 'ht': 'dlf', 'verb': 2}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Impulse response\n~~~~~~~~~~~~~~~~\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "inparg['signal'] = 0 # signal 0 = impulse\nprint('QWE')\nqwe = empymod.dipole(**inparg, ft='qwe', ftarg={'maxint': 500})\nprint('DLF (Sine)')\nsin = empymod.dipole(**inparg, ft='dlf', ftarg={'dlf': 'key_81_CosSin_2009'})\nprint('FFTLog')\nftl = empymod.dipole(**inparg, ft='fftlog')\nprint('FFT')\nfft = empymod.dipole(\n **inparg, ft='fft',\n ftarg={'dfreq': .001, 'nfreq': 2**15, 'ntot': 2**16, 'pts_per_dec': 10}\n )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Used number of frequencies:\n\n- QWE: 167 (0.000794328 - 158489 Hz)\n- Sine: 116 (5.33905E-06 - 52028 Hz)\n- FFTLog: 60 (0.000178575 - 141.847 Hz)\n- FFT: 46 (0.001 - 32.768 Hz)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure()\nplt.title(r'Impulse response for HS-model, $r=$' +\n str(int(rec[0]/1000)) + ' km.')\nplt.xlabel('Time (s)')\nplt.ylabel(r'Amplitude (V/m)')\nplt.semilogx(t, qwe, 'C0-', label='QWE')\nplt.semilogx(t, sin, 'C1--', label='Sine Filter')\nplt.semilogx(t, ftl, 'C2-.', label='FFTLog')\nplt.semilogx(t, fft, 'C3:', label='FFT')\nplt.legend(loc='best')\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step response\n~~~~~~~~~~~~~\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "inparg['signal'] = 1 # signal 1 = step\nprint('QWE')\nqwe = empymod.dipole(**inparg, ft='qwe', ftarg={'nquad': 31, 'maxint': 500})\nprint('DLF (Sine)')\nsin = empymod.dipole(**inparg, ft='dlf', ftarg={'dlf': 'key_81_CosSin_2009'})\nprint('FFTLog')\nftl = empymod.dipole(**inparg, ft='fftlog', ftarg={'add_dec': [-2, 4]})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Used number of frequencies:\n\n- QWE: 173 (0.000398107 - 158489 Hz)\n- Sine: 116 (5.33905E-06 - 52028 Hz)\n- FFTLog: 90 (0.000178575 - 141847 Hz)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure()\nplt.title(r'Step response for HS-model, $r=$' + str(int(rec[0]/1000)) + ' km.')\nplt.xlabel('Time (s)')\nplt.ylabel('Amplitude (V/m)')\nplt.semilogx(t, qwe, 'C0-', label='QWE')\nplt.semilogx(t, sin, 'C1--', label='Sine Filter')\nplt.semilogx(t, ftl, 'C2-.', label='FFTLog')\nplt.ylim([-.1e-12, 1.5*qwe.max()])\nplt.legend(loc='best')\nplt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "empymod.Report()" ] } ], "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.7.3" } }, "nbformat": 4, "nbformat_minor": 0 }