# Step and impulse responses¶

These examples compare the analytical solution with empymod for time-domain step and impulse responses for inline, x-directed source and receivers, for the four different frequency-to-time methods QWE, DLF, FFTLog, and FFT. Which method is faster and which is more precise depends on the model (land or marine, source/receiver at air-interface or not) and the response (step or impulse).

import empymod
import numpy as np
from scipy.special import erf
import matplotlib.pyplot as plt
from scipy.constants import mu_0       # Permeability of free space  [H/m]
plt.style.use('ggplot')
colors = [color['color'] for color in list(plt.rcParams['axes.prop_cycle'])]


## Analytical solutions¶

Analytical solution for source and receiver at the interface between two half-spaces

The time-domain step and impulse responses for a source at the origin ($$x_s = y_s = z_s = 0$$ m) and an in-line receiver at the surface ($$y_r = z_r = 0$$ m), is given by the following equations, where $$\rho_h$$ is horizontal resistivity ($$\Omega$$ m), $$\lambda$$ is anisotropy (-), with $$\lambda = \sqrt{\rho_v/\rho_h}$$, $$r$$ is offset (m), $$t$$ is time (s), and $$\tau_h = \sqrt{\mu_0 r^2/(\rho_h t)}$$; $$\mu_0$$ is the magnetic permeability of free space (H/m).

### Time Domain: Step Response $$\mathbf{\mathcal{H}(t)}$$¶

$E_x(\rho_h,\lambda,r,t) = \frac{\rho_h}{2 \pi r^3} \left[ 2\lambda + \rm{erf}\left(\frac{\tau_h}{2}\right) - 2\lambda \rm{erf}\left(\frac{\tau_h}{2\lambda}\right) + \frac{\tau_h}{\sqrt{\pi}} \exp\left(- \frac{\tau_h^2}{4\lambda^2}\right) \right]$

### Time Domain: Impulse Response $$\mathbf{\delta(t)}$$¶

$\dot{E}_x(\rho_h,\lambda,r,t) = \frac{\rho_h}{2 \pi r^3} \left[ \delta(t) + \frac{\tau_h}{2t\sqrt{\pi}} \left\{ - \exp\left(-\frac{\tau_h^2}{4}\right) + \left( \frac{\tau_h^2}{2 \lambda^2} + 1 \right) \exp\left(- \frac{\tau_h^2}{4 \lambda^2}\right) \right\} \right]$

### Reference¶

Equations 3.2 and 3.3 in Werthmüller, D., 2009, Inversion of multi-transient EM data from anisotropic media: M.S. thesis, TU Delft, ETH Zürich, RWTH Aachen; UUID: f4b071c1-8e55-4ec5-86c6-a2d54c3eda5a.

### Analytical functions¶

def ee_xx_impulse(res, aniso, off, time):
"""VTI-Halfspace impulse response, xx, inline.

res   : horizontal resistivity [Ohm.m]
aniso : anisotropy [-]
off   : offset [m]
time  : time(s) [s]
"""
tau_h = np.sqrt(mu_0*off**2/(res*time))
t0 = tau_h/(2*time*np.sqrt(np.pi))
t1 = np.exp(-tau_h**2/4)
t2 = tau_h**2/(2*aniso**2) + 1
t3 = np.exp(-tau_h**2/(4*aniso**2))
Exx = res/(2*np.pi*off**3)*t0*(-t1 + t2*t3)
Exx[time == 0] = res/(2*np.pi*off**3)  # Delta dirac part
return Exx

def ee_xx_step(res, aniso, off, time):
"""VTI-Halfspace step response, xx, inline.

res   : horizontal resistivity [Ohm.m]
aniso : anisotropy [-]
off   : offset [m]
time  : time(s) [s]
"""
tau_h = np.sqrt(mu_0*off**2/(res*time))
t0 = erf(tau_h/2)
t1 = 2*aniso*erf(tau_h/(2*aniso))
t2 = tau_h/np.sqrt(np.pi)*np.exp(-tau_h**2/(4*aniso**2))
Exx = res/(2*np.pi*off**3)*(2*aniso + t0 - t1 + t2)
return Exx


## Example 1: Source and receiver at z=0m¶

Comparison with analytical solution; put 1 mm below the interface, as they would be regarded as in the air by emmod otherwise.

src = [0, 0, 0.001]          # Source at origin, slightly below interface
rec = [6000, 0, 0.001]       # Receivers in-line, 0.5m below interface
res = [2e14, 10]             # Resistivity: [air, half-space]
aniso = [1, 2]               # Anisotropy: [air, half-space]
eperm = [0, 1]               # Set el. perm. of air to 0 because of num. noise
t = np.logspace(-2, 1, 301)  # Desired times (s)

# Collect parameters
inparg = {'src': src, 'rec': rec, 'depth': 0, 'freqtime': t, 'res': res,
'aniso': aniso, 'epermH': eperm, 'ht': 'dlf', 'verb': 2}


### Impulse response¶

ex = ee_xx_impulse(res[1], aniso[1], rec[0], t)

inparg['signal'] = 0  # signal 0 = impulse
print('QWE')
qwe = empymod.dipole(**inparg, ft='qwe')
print('DLF (Sine)')
sin = empymod.dipole(**inparg, ft='dlf', ftarg={'dlf': 'key_81_CosSin_2009'})
print('FFTLog')
ftl = empymod.dipole(**inparg, ft='fftlog')
print('FFT')
fft = empymod.dipole(
**inparg, ft='fft',
ftarg={'dfreq': .0005, 'nfreq': 2**20, 'pts_per_dec': 10})


Out:

QWE
* WARNING :: Fourier-quadrature did not converge at least once;
=> desired atol and rtol might not be achieved.

:: empymod END; runtime = 0:00:03.487460 :: 1 kernel call(s)

DLF (Sine)

:: empymod END; runtime = 0:00:00.012493 :: 1 kernel call(s)

FFTLog

:: empymod END; runtime = 0:00:00.006674 :: 1 kernel call(s)

FFT

:: empymod END; runtime = 0:00:00.719433 :: 1 kernel call(s)


=> FFTLog is the fastest by quite a margin, followed by the Sine-filter. What cannot see from the output (set verb to something bigger than 2 to see it) is how many frequencies each method used:

• QWE: 159 (0.000794328 - 63095.7 Hz)

• Sine: 116 (5.33905E-06 - 52028 Hz)

• FFTLog: 60 (0.000178575 - 141.847 Hz)

• FFT: 61 (0.0005 - 524.288 Hz)

Note that for the actual transform, FFT used 2^20 = 1’048’576 frequencies! It only computed 60 frequencies, and then interpolated the rest, as it requires regularly spaced data.

plt.figure()
plt.title(r'Impulse response for HS-model, $r=$' +
str(int(rec[0]/1000)) + ' km.')
plt.xlabel('Time (s)')
plt.ylabel(r'Amplitude (V/m)')
plt.semilogx(t, ex, 'k-', label='Analytical')
plt.semilogx(t, qwe, 'C0-', label='QWE')
plt.semilogx(t, sin, 'C1--', label='Sine Filter')
plt.semilogx(t, ftl, 'C2-.', label='FFTLog')
plt.semilogx(t, fft, 'C3:', label='FFT')
plt.legend(loc='best')
plt.ylim([-.1*np.max(ex), 1.1*np.max(ex)])
plt.show()

plt.figure()
plt.title('Error')
plt.xlabel('Time (s)')
plt.ylabel('Relative error (-)')
plt.loglog(t, abs(qwe-ex)/ex, 'C0-', label='QWE')
plt.plot(t, abs(sin-ex)/ex, 'C1--', label='Sine Filter')
plt.plot(t, abs(ftl-ex)/ex, 'C2-.', label='FFTLog')
plt.plot(t, abs(fft-ex)/ex, 'C3:', label='FFT')
plt.legend(loc='best')
plt.show()


=> The error is comparable in all cases. FFT is not too good at later times. This could be improved by computing lower frequencies. But because FFT needs regularly spaced data, our vector would soon explode (and you would need a lot of memory). In the current case we are already using 2^20 samples!

### Step response¶

Step responses are almost impossible with FFT. We can either try to model late times with lots of low frequencies, or the step with lots of high frequencies. I do not use FFT in the step-response examples.

#### Switch-on¶

ex = ee_xx_step(res[1], aniso[1], rec[0], t)

inparg['signal'] = 1  # signal 1 = switch-on
print('QWE')
qwe = empymod.dipole(**inparg, ft='qwe')
print('DLF (Sine)')
sin = empymod.dipole(**inparg, ft='dlf', ftarg={'dlf': 'key_81_CosSin_2009'})
print('FFTLog')
ftl = empymod.dipole(**inparg, ft='fftlog', ftarg={'q': -0.6})


Out:

QWE
* WARNING :: Fourier-quadrature did not converge at least once;
=> desired atol and rtol might not be achieved.

:: empymod END; runtime = 0:00:14.095677 :: 1 kernel call(s)

DLF (Sine)

:: empymod END; runtime = 0:00:00.011704 :: 1 kernel call(s)

FFTLog

:: empymod END; runtime = 0:00:00.007237 :: 1 kernel call(s)


Used number of frequencies:

• QWE: 159 (0.000794328 - 63095.7 Hz)

• Sine: 116 (5.33905E-06 - 52028 Hz)

• FFTLog: 60 (0.000178575 - 141.847 Hz)

plt.figure()
plt.title(r'Switch-on response for HS-model, $r=$' +
str(int(rec[0]/1000)) + ' km.')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude (V/m)')
plt.semilogx(t, ex, 'k-', label='Analytical')
plt.semilogx(t, qwe, 'C0-', label='QWE')
plt.semilogx(t, sin, 'C1--', label='Sine Filter')
plt.semilogx(t, ftl, 'C2-.', label='FFTLog')
plt.legend(loc='best')
plt.show()

plt.figure()
plt.title('Error')
plt.xlabel('Time (s)')
plt.ylabel('Relative error (-)')
plt.loglog(t, abs(qwe-ex)/ex, 'C0-', label='QWE')
plt.plot(t, abs(sin-ex)/ex, 'C1--', label='Sine Filter')
plt.plot(t, abs(ftl-ex)/ex, 'C2-.', label='FFTLog')
plt.legend(loc='best')
plt.show()


#### Switch-off¶

For switch-off to work properly you need empymod-version bigger than 1.3.0! You can do it with previous releases too, but you will have to do the DC-computation and subtraction manually, as is done here for ee_xx_step.

exDC = ee_xx_step(res[1], aniso[1], rec[0], 60*60)
ex = exDC - ee_xx_step(res[1], aniso[1], rec[0], t)

inparg['signal'] = -1  # signal -1 = switch-off
print('QWE')
qwe = empymod.dipole(**inparg, ft='qwe')
print('DLF (Cosine/Sine)')
sin = empymod.dipole(**inparg, ft='dlf', ftarg={'dlf': 'key_81_CosSin_2009'})
print('FFTLog')
ftl = empymod.dipole(**inparg, ft='fftlog', ftarg={'add_dec': [-5, 3]})


Out:

QWE
* WARNING :: Fourier-quadrature did not converge at least once;
=> desired atol and rtol might not be achieved.

:: empymod END; runtime = 0:00:10.464879 :: 1 kernel call(s)

DLF (Cosine/Sine)

:: empymod END; runtime = 0:00:00.011193 :: 1 kernel call(s)

FFTLog

:: empymod END; runtime = 0:00:00.010142 :: 1 kernel call(s)

plt.figure()
plt.title(r'Switch-off response for HS-model, $r=$' +
str(int(rec[0]/1000)) + ' km.')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude (V/m)')
plt.semilogx(t, ex, 'k-', label='Analytical')
plt.semilogx(t, qwe, 'C0-', label='QWE')
plt.semilogx(t, sin, 'C1--', label='Cosine/Sine Filter')
plt.semilogx(t, ftl, 'C2-.', label='FFTLog')
plt.legend(loc='best')
plt.show()

plt.figure()
plt.title('Error')
plt.xlabel('Time (s)')
plt.ylabel('Relative error (-)')
plt.loglog(t, abs(qwe-ex)/ex, 'C0-', label='QWE')
plt.plot(t, abs(sin-ex)/ex, 'C1--', label='Sine Filter')
plt.plot(t, abs(ftl-ex)/ex, 'C2-.', label='FFTLog')
plt.legend(loc='best')
plt.show()


## Example 2: Air-seawater-halfspace¶

In seawater the transformation is generally much easier, as we do not have the step or the impules at zero time.

src = [0, 0, 950]            # Source 50 m above seabottom
rec = [6000, 0, 1000]        # Receivers in-line, at seabottom
res = [1e23, 1/3, 10]        # Resistivity: [air, water, half-space]
aniso = [1, 1, 2]            # Anisotropy: [air, water, half-space]
t = np.logspace(-2, 1, 301)  # Desired times (s)

# Collect parameters
inparg = {'src': src, 'rec': rec, 'depth': [0, 1000], 'freqtime': t,
'res': res, 'aniso': aniso, 'ht': 'dlf', 'verb': 2}


### Impulse response¶

inparg['signal'] = 0  # signal 0 = impulse
print('QWE')
qwe = empymod.dipole(**inparg, ft='qwe', ftarg={'maxint': 500})
print('DLF (Sine)')
sin = empymod.dipole(**inparg, ft='dlf', ftarg={'dlf': 'key_81_CosSin_2009'})
print('FFTLog')
ftl = empymod.dipole(**inparg, ft='fftlog')
print('FFT')
fft = empymod.dipole(
**inparg, ft='fft',
ftarg={'dfreq': .001, 'nfreq': 2**15, 'ntot': 2**16, 'pts_per_dec': 10}
)


Out:

QWE

:: empymod END; runtime = 0:00:15.203947 :: 1 kernel call(s)

DLF (Sine)

:: empymod END; runtime = 0:00:00.028618 :: 1 kernel call(s)

FFTLog

:: empymod END; runtime = 0:00:00.015045 :: 1 kernel call(s)

FFT

:: empymod END; runtime = 0:00:00.044509 :: 1 kernel call(s)


Used number of frequencies:

• QWE: 167 (0.000794328 - 158489 Hz)

• Sine: 116 (5.33905E-06 - 52028 Hz)

• FFTLog: 60 (0.000178575 - 141.847 Hz)

• FFT: 46 (0.001 - 32.768 Hz)

plt.figure()
plt.title(r'Impulse response for HS-model, $r=$' +
str(int(rec[0]/1000)) + ' km.')
plt.xlabel('Time (s)')
plt.ylabel(r'Amplitude (V/m)')
plt.semilogx(t, qwe, 'C0-', label='QWE')
plt.semilogx(t, sin, 'C1--', label='Sine Filter')
plt.semilogx(t, ftl, 'C2-.', label='FFTLog')
plt.semilogx(t, fft, 'C3:', label='FFT')
plt.legend(loc='best')
plt.show()


### Step response¶

inparg['signal'] = 1  # signal 1 = step
print('QWE')
qwe = empymod.dipole(**inparg, ft='qwe', ftarg={'nquad': 31, 'maxint': 500})
print('DLF (Sine)')
sin = empymod.dipole(**inparg, ft='dlf', ftarg={'dlf': 'key_81_CosSin_2009'})
print('FFTLog')
ftl = empymod.dipole(**inparg, ft='fftlog', ftarg={'add_dec': [-2, 4]})


Out:

QWE
* WARNING :: Fourier-quadrature did not converge at least once;
=> desired atol and rtol might not be achieved.

:: empymod END; runtime = 0:00:12.765028 :: 1 kernel call(s)

DLF (Sine)

:: empymod END; runtime = 0:00:00.028450 :: 1 kernel call(s)

FFTLog

:: empymod END; runtime = 0:00:00.022921 :: 1 kernel call(s)


Used number of frequencies:

• QWE: 173 (0.000398107 - 158489 Hz)

• Sine: 116 (5.33905E-06 - 52028 Hz)

• FFTLog: 90 (0.000178575 - 141847 Hz)

plt.figure()
plt.title(r'Step response for HS-model, $r=$' + str(int(rec[0]/1000)) + ' km.')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude (V/m)')
plt.semilogx(t, qwe, 'C0-', label='QWE')
plt.semilogx(t, sin, 'C1--', label='Sine Filter')
plt.semilogx(t, ftl, 'C2-.', label='FFTLog')
plt.ylim([-.1e-12, 1.5*qwe.max()])
plt.legend(loc='best')
plt.show()

empymod.Report()

 Wed Jun 23 07:09:49 2021 UTC OS Linux CPU(s) 2 Machine x86_64 Architecture 64bit RAM 7.5 GiB Environment Python Python 3.8.6 (default, Oct 19 2020, 15:10:29) [GCC 7.5.0] numpy 1.21.0 scipy 1.7.0 numba 0.53.1 empymod 2.1.1 IPython 7.24.1 matplotlib 3.4.2

Total running time of the script: ( 1 minutes 1.154 seconds)

Estimated memory usage: 315 MB

Gallery generated by Sphinx-Gallery