Source code for fatiando.seismic.conv

Zero-offset convolutional seismic modeling

Give a depth model and obtain a seismic zero-offset convolutional gather. You
can give the wavelet, if you already have, or use one of the existing, from
which we advise ricker wavelet (rickerwave function).

* :func:`~fatiando.seismic.conv.convolutional_model`: given the reflectivity
  series and wavelet, it returns the convolutional seismic gather.
* :func:`~fatiando.seismic.conv.reflectivity`: calculates the reflectivity
  series from the velocity model (and density model if present).
* :func:`~fatiando.seismic.conv.depth_2_time`: convert depth property model to
  the model in time.
* :func:`~fatiando.seismic.conv.rickerwave`: calculates a ricker wavelet.


Yilmaz, Oz,
Ch.2 Deconvolution. In: YILMAZ, Oz. Seismic Data Analysis: Processing,
Inversion, and Interpretation of Seismic Data. Tulsa: Seg, 2001. Cap. 2.
p. 159-270. Available at: <>

from __future__ import division
import numpy as np
from scipy import interpolate  # linear interpolation of velocity/density

[docs]def convolutional_model(rc, f, wavelet, dt): """ Calculate convolutional seismogram for a geological model Calculate the synthetic convolutional seismogram of a geological model, Vp is mandatory while density is optional. The given model in a matrix form is considered a mesh of square cells. .. warning:: Since the relative difference between the model is the important, being consistent with the units chosen for the parameters is the only requirement, whatever the units. Parameters: * rc : 2D-array Reflectivity values in the time domain. * f : float Dominant frequency of the ricker wavelet. * wavelet : float The function to consider as source in the seismic modelling. * dt: float Sample time of the ricker wavelet and of the resulting seismogram, in general a value of 2.e-3 is used. Returns: * synth_l : 2D-array Resulting seismogram. """ w = wavelet(f, dt) # convolution synth_l = np.zeros(np.shape(rc)) for j in range(0, rc.shape[1]): if np.shape(rc)[0] >= len(w): synth_l[:, j] = np.convolve(rc[:, j], w, mode='same') else: aux = np.floor(len(w)/2.) synth_l[:, j] = np.convolve(rc[:, j], w, mode='full')[aux:-aux] return synth_l
[docs]def reflectivity(model_t, rho): """ Calculate reflectivity series in the time domain, so it is necessary to use the function depth_2_time first if the model is in depth domain. Note this this function can also be used to one dimensional array. Parameters: * model_t : 2D-array Velocity values in time domain. * rho : 2D-array (optional) Density values for all the model, in time domain. Returns: * rc : 2D-array Calculated reflectivity series for all the model given. """ err_message = "Velocity and density matrix must have the same dimension." assert model_t.shape == rho.shape, err_message rc = np.zeros(np.shape(model_t)) rc[1:, :] = ((model_t[1:]*rho[1:]-model_t[:-1]*rho[:-1]) / (model_t[1:]*rho[1:]+model_t[:-1]*rho[:-1])) return rc
[docs]def depth_2_time(vel, model, dt, dz): """ Convert depth property model to time model. Parameters: * vel : 2D-array Velocity values in the depth domain. * model : 2D-array Model values in the depth domain. * dt: float Sample time of the ricker wavelet and of the resulting seismogram, in general a value of 2.e-3 is used. * dz : float Length of square grid cells. * rho : 2D-array (optional) Density values for all the model, in depth domain. Returns: * model_t : 2D-array Property model in time domain. """ err_message = "Velocity and model matrix must have the same dimension." assert vel.shape == model.shape, err_message # downsampled time rate to make a better interpolation n_samples, n_traces = [vel.shape[0], vel.shape[1]] dt_dwn = dt/10. if dt_dwn > dz/np.max(vel): dt_dwn = (dz/np.max(vel))/10. TWT = np.zeros((n_samples, n_traces)) TWT[0, :] = 2*dz/vel[0, :] for j in range(1, n_samples): TWT[j, :] = TWT[j-1]+2*dz/vel[j, :] TMAX = max(TWT[-1, :]) TMIN = min(TWT[0, :]) TWT_rs = np.zeros(np.ceil(TMAX/dt_dwn)) for j in range(1, len(TWT_rs)): TWT_rs[j] = TWT_rs[j-1]+dt_dwn resmpl = int(dt/dt_dwn) model_t = _resampling(model, TMAX, TWT, TWT_rs, dt, dt_dwn, n_traces) return model_t
def _resampling(model, TMAX, TWT, TWT_rs, dt, dt_dwn, n_traces): """ Resamples the input data to adjust it after time conversion with the chosen time sample rate, dt. Returns: * vel_l : 2D-array Resampled input data. """ vel = np.ones((np.ceil(TMAX/dt_dwn), n_traces)) for j in range(0, n_traces): kk = np.ceil(TWT[0, j]/dt_dwn) lim = np.ceil(TWT[-1, j]/dt_dwn)-1 # necessary do before resampling to have values in all points of time model tck = interpolate.interp1d(TWT[:, j], model[:, j]) vel[kk:lim, j] = tck(TWT_rs[kk:lim]) # the model is extended in time because of depth time conversion vel[lim:, j] = vel[lim-1, j] # because of time conversion, the values between 0 e kk need to be filed vel[0:kk, j] = model[0, j] # resampling from dt_dwn to dt vel_l = np.zeros((np.ceil(TMAX/dt), n_traces)) resmpl = int(dt/dt_dwn) vel_l[0, :] = vel[0, :] for j in range(0, n_traces): for jj in range(1, int(np.ceil(TMAX/dt))): vel_l[jj, j] = vel[resmpl*jj, j] return vel_l
[docs]def rickerwave(f, dt): r""" Given a frequency and time sampling rate, outputs ricker function. The length of the function varies according to f and dt, in order for the ricker function starts and ends as zero. It is also considered that the functions is causal, what means it starts at time zero. To satisfy sampling and stability: .. math:: f << \frac{1}{2 dt}. Here, we consider this as: .. math:: f < 0.2 \frac{1}{2 dt}. Parameters: * f : float dominant frequency value in Hz * dt : float time sampling rate in seconds (usually it is in the order of ms) Returns: * ricker : float Ricker function for the given parameters. """ assert f < 0.2*(1./(2.*dt)), "Frequency too high for the dt chosen." nw = 2.2/f/dt nw = 2*np.floor(nw/2)+1 nc = np.floor(nw/2) ricker = np.zeros(nw) k = np.arange(1, nw+1) alpha = (nc-k+1)*f*dt*np.pi beta = alpha**2 ricker = (1.-beta*2)*np.exp(-beta) return ricker