Skip to content

Instantly share code, notes, and snippets.

@eusoubrasileiro
Last active December 15, 2021 13:09
Show Gist options
  • Save eusoubrasileiro/8165804 to your computer and use it in GitHub Desktop.
Save eusoubrasileiro/8165804 to your computer and use it in GitHub Desktop.
Seismic: 2D finite difference simulation of elastic P and SV wave propagation Simulates Parana Basin geometry of basalt cover with 2 horizontal layers, acquires just z displacement (default seismic acquisition).
# This Python file uses the following encoding: utf-8
"""
Seismic: 2D finite difference simulation of elastic P and SV wave propagation
Simulates Parana Basin geometry of basalt cover with horizontal layers.
Acquires just z displacement (default seismic acquisition).
4 layers :
bauru ~ 100 meters ~ 6x scaled to 600 m
serra geral flow 1 ~ 20 meters ~ 6x scaled to 120 m
serra geral flow 2 ~ 20 meters ~ 6x scaled to 120 m
botucatu ~ till end
Example based on pg 281/283 - Análise do Sinal Sísmico Rosa, André L.
plus personal insight of vibroseis data of Parana Basin acquisition 0314 . scaled x 6
rock properties
bauru : vp = 2181 vs = 1166 rho = 2.12
serra geral flow 1 : vp = 4800 vs = 2611 rho = 2.75
serra geral flow 2 : vp = 3540 vs = 1906 rho = 2.54 (vugular)
botucatu : vp = 3200 vs = 1600 rho = 2.65 (~botucatu sandstone)
"""
import numpy as np
from matplotlib import animation
from fatiando.seismic import wavefd
from fatiando.vis import mpl
from fatiando import utils
# Set the parameters of the finite difference grid
shape = (200, 200) # z and x shape
area = [0, 8000, 0, 2000] # x min and xmax etc... dx = 40 ; dz = 10
dx = area[1]/shape[1]
dz = area[3]/shape[0]
# Geologic model make a density and S wave velocity model
density = 2650*np.ones(shape) # botucatu
density[0:int(600/dz), :] = 2120 # bauru
density[int(600/dz):int(600/dz)+int(120/dz), :] = 2750 # basalt flow 1
density[int(600/dz)+int(120/dz):int(600/dz)+int(240/dz), :] = 2540 # basalt flow 2
pvel = 3200*np.ones(shape) # botucatu
pvel[0:int(600/dz), :] = 2181 # bauru
pvel[int(600/dz):int(600/dz)+int(120/dz), :] = 4800 # basalt flow 1
pvel[int(600/dz)+int(120/dz):int(600/dz)+int(240/dz), :] = 3540 # basalt flow 2
svel = 1600*np.ones(shape) # botucatu
svel[0:int(600/dz), :] = 1166 # bauru
svel[int(600/dz):int(600/dz)+int(120/dz), :] = 2611 # basalt flow 1
svel[int(600/dz)+int(120/dz):int(600/dz)+int(240/dz), :] = 1900 # basalt flow 2
mu = wavefd.lame_mu(svel, density)
lamb = wavefd.lame_lamb(pvel, svel, density)
# Make a wave source from a gauss derivative that vibrates in the z direction only
# (removes the shear component amplitude equals zero, simulating a default seismic acquisition)
sources = [[wavefd.GaussSource(2000, 0, area, shape, 0.0, 5.)], # x source or shear source
[wavefd.GaussSource(2000, 0, area, shape, 100.0, 25.)]] # z source or z compressional source
# Get the iterator for the simulation
dt = wavefd.maxdt(area, shape, np.max(pvel))
duration = 3.0
maxit = int(duration/dt)
stations = [[2200+i*dx, 0] for i in range(130)] # x, z coordinate of the seismometers
snapshot = int(0.004/dt) # Plot a snapshot of the simulation every 4 miliseconds
print "dt for simulation is ", dt
print "max iteration for simulation is ", maxit
print "duration for simulation is ", duration
simulation = wavefd.elastic_psv(lamb, mu, density, area, dt, maxit, sources,
stations, snapshot, padding=50, taper=0.01)
# This part makes an animation using matplotlibs animation API
fig = mpl.figure(figsize=(12, 5))
# Start with everything zero and grab the plot so that it can be updated later
vp = pvel*10**-4
vpmean = np.mean(vp)
mpl.imshow(vp[::-1], extent=area, alpha=0.25)
wavefield = mpl.imshow(np.zeros(shape), extent=area, vmin=-10**-7, vmax=10**-7, alpha=0.75, cmap=mpl.cm.gray_r)
mpl.points(stations, '.k')
mpl.ylim(area[2:][::-1])
mpl.xlabel('x (km)')
mpl.ylabel('z (km)')
mpl.m2km()
times = np.linspace(0, maxit*dt, maxit)
# the only way for getting the feedback response of the seimograms
global seismograms
# This function updates the plot every few timesteps
def animate(i):
# Simulation will yield panels corresponding to P and S waves because
# xz2ps=True
t, ux, uz, xcomp, zcomp = simulation.next()
mpl.title('time: %0.3f s' % (times[t]))
# wavefield.set_array((p + s)[::-1])
wavefield.set_data(uz[::-1])
global seismograms
seismograms = zcomp
return wavefield, seismograms
anim = animation.FuncAnimation(fig, animate, frames=maxit/snapshot, interval=1)
mpl.show()
# turn the seismogram data in a numpy array
traces = np.array(seismograms)
parana_shot = utils.matrix2stream(traces)
mpl.seismic_wiggle(parana_shot, ranges=[200, 130*dx], scale=2, color='b', normalize=True)
mpl.xlabel('offset (m)')
mpl.ylabel('time (ms)')
mpl.show()
@vasconcelos82
Copy link

Olá André, como vai? Sou professor do curso de geofísica da UFBA.Achei interessante o código, e gostaria de tirar algumas dúvidas com você. Posso contacta-lo por email? Se sim, qual seria? Abs, Marcos

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment