Nessa aula, vamos começar a estudar a sísmica de reflexão. Veremos como identificar a reflexão na simulação de ondas P e no gráfico de tempo x distância. Vamos simular também uma seção Common Mid Point, na qual movemos a fonte e o receptor mantendo o ponto médio fixo.
Utilizaremos as simulações de ondas da biblioteca Fatiando a Terra. Essas simulações utilizam o método de diferenças finitas para calcular soluções da equação da onda.
Esse documento é um Jupyter notebook, um documento interativo que mistura texto (como esse), código (como abaixo), e o resultado de executar o código (números, texto, figuras, videos, etc).
O notebook te fornecerá exemplos interativos que trabalham os temas abordados no questionário. Utilize esses exemplos para responder as perguntas.
As células com números ao lado, como In [1]:
, são código Python. Algumas dessas células não produzem resultado e servem de preparação para os exemplos interativos. Outras, produzem gráficos interativos. Você deve executar todas as células, uma de cada vez, mesmo as que não produzem gráficos.
Para executar uma célula, clique em cima dela e aperte Shift + Enter
. O foco (contorno verde ou cinza em torno da célula) deverá passar para a célula abaixo. Para rodá-la, aperte Shift + Enter
novamente e assim por diante. Você pode executar células de texto que não acontecerá nada.
Rode as células abaixo para carregar os módulos necessários para essa prática.
%matplotlib inline
from __future__ import division, print_function
from multiprocessing import Pool, cpu_count
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import ipywidgets as ipw
from fatiando.seismic import RickerWavelet, FDAcoustic2D
from fatiando.vis import anim_to_html
from fatiando.vis.mpl import seismic_image
/home/leo/miniconda3/envs/geofisica2/lib/python2.7/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`. from ._conv import register_converters as _register_converters
Vamos utilizar as simulações de onda P para gerar dados sintéticos. Dessa vez também vamos ignorar as ondas S para não complicar a vida mais do que o necessário.
O nosso modelo será uma interface plana separando dois meios.
shape = (300, 1000)
spacing = 10
extent = [0, shape[1]*spacing, shape[0]*spacing, 0]
interface = 100
density = np.zeros(shape, dtype='float32') + 2000
density[interface:,:] = 2500
velocity = np.zeros(shape, dtype='float32') + 4000
velocity[interface:,:] = 5000
shot = FDAcoustic2D(velocity, density, spacing=spacing, taper=0.004, padding=60)
fonte = (0, 100)
shot.add_point_source(fonte, RickerWavelet(1, 60))
shot.run(3000)
[##################################################] | 100% Completed | 1min 55.6s
Rode a célula abaixo para gerar uma animação da aquisição.
every = 20
frames = shot.simsize//every
receptores = np.arange(fonte[1] + 50, shape[1], 15, dtype='int')
x = receptores*spacing
dados = shot[:, 0, receptores]
times = np.linspace(0, shot.dt*shot.simsize, shot.simsize)
dados_dummy = np.zeros_like(dados)
fig, axes = plt.subplots(2, 1, sharex=True, sharey=False, figsize=(12, 8), facecolor='white')
ax1 = axes[0]
ax1.set_title('iteration: 0')
ax1.set_ylabel('time (s)')
cutoff = 0.4
section = ax1.imshow(dados_dummy, extent=[x.min(), x.max(), times.max(), times.min()],
aspect=1400, cmap='Greys', vmin=-cutoff, vmax=cutoff,
interpolation='nearest')
ax1.set_ylim(0, times.max())
ax2 = axes[1]
ax2.set_xlabel('x (m)')
ax2.set_ylabel('depth (m)')
cutoff = 1
wavefield = ax2.imshow(shot[0, :, :], extent=extent, vmin=-cutoff, vmax=cutoff, cmap='Greys')
ax2.plot(fonte[1]*spacing, 0, '*y', markersize=15)
ax2.plot(x, np.zeros_like(receptores), 'vb', markersize=10)
ax2.hlines(interface*spacing, 0, shape[1]*spacing)
fig.tight_layout(h_pad=0)
def anim_shot(frame):
ax1.set_title('iteration: {:d}'.format(frame*every))
u = shot[frame*every]
wavefield.set_array(u)
dados_dummy[:frame*every, :] = dados[:frame*every, :]
section.set_array(dados_dummy)
return wavefield
anim = FuncAnimation(fig, anim_shot, frames=frames)
anim_to_html(anim, fps=6, dpi=60)
Na animação acima:
Agora vamos simular a aquisição de uma seção Common Mid Point. A diferença é que dessa vez utilizaremos vários pares fonte-receptor a distâncias diferentes. Na simulação anterior, utilizamos apenas uma fonte e vários receptores. Cada fonte será uma simulação da qual vamos extrair as medições de um único receptor. A animação abaixo deixará isso mais claro.
shape = (100, 200)
spacing = 10
extent = [0, shape[1]*spacing, shape[0]*spacing, 0]
interface = shape[0]//2
density = np.zeros(shape, dtype='float32') + 2000
density[interface:,:] = 2500
velocity = np.zeros(shape, dtype='float32') + 4000
velocity[interface:,:] = 5000
step = 3
fontes = np.array(list(reversed(range(55, shape[1]//2 - step//2, step))))
recep = np.array([shape[1] - s for s in fontes])
offsets = (recep - fontes)*spacing
print("Utilizando {} fontes e {} receptores.".format(len(fontes), len(recep)))
print('Fontes (m): {}'.format(fontes))
print('Receptores (m): {}'.format(recep*spacing))
print('Offsets (m): {}'.format(offsets))
Utilizando 15 fontes e 15 receptores. Fontes (m): [97 94 91 88 85 82 79 76 73 70 67 64 61 58 55] Receptores (m): [1030 1060 1090 1120 1150 1180 1210 1240 1270 1300 1330 1360 1390 1420 1450] Offsets (m): [ 60 120 180 240 300 360 420 480 540 600 660 720 780 840 900]
Rode as células abaixo para rodar uma simulação para cada fonte acima. Não aparecerá a barrinha de progresso dessa vez pois vamos rodar as simulações em paralelo para agilizar o processo.
def simul_shot(fonte, its=800, verbose=False):
shot = FDAcoustic2D(velocity, density, spacing=spacing, taper=0.005, padding=50, verbose=verbose)
shot.add_point_source((0, fonte), RickerWavelet(1, 100))
shot.run(its)
return shot
%%time
print('Simulando...')
pool = Pool(processes=cpu_count())
shots = pool.map(simul_shot, fontes)
pool.close()
print('Terminado.')
Simulando... Terminado. CPU times: user 20.5 ms, sys: 56.8 ms, total: 77.4 ms Wall time: 11.8 s
Rode a célula abaixo para gerar uma animação da aquisição do CMP.
every = 30
frames = shots[0].simsize//every
dt = shots[0].dt
times = np.linspace(0, dt*shots[0].simsize, shots[0].simsize)
CMP = np.empty((shots[0].simsize, len(recep)))
for i, s in enumerate(shots):
CMP[:, i] = s[:, 0, recep[i]]
CMP_dummy = np.zeros_like(CMP)
fig, axes = plt.subplots(2, 1, sharex=False, sharey=False, figsize=(7, 8), facecolor='white')
ax1 = axes[1]
ax1.set_title('shot: 1')
ax1.set_xlabel('x (m)')
ax1.set_ylabel('depth (m)')
cutoff = 1
wavefield = ax1.imshow(shots[0][500, :, :], extent=extent, vmin=-cutoff, vmax=cutoff, cmap='Greys')
src, = ax1.plot(fontes[0]*spacing, 0, '*y', markersize=15)
rec, = ax1.plot(recep[0]*spacing, 0, 'vb', markersize=10)
ax1.hlines(interface*spacing, 0, shape[1]*spacing)
ax2 = axes[0]
ax2.set_title('CMP: 1')
ax2.set_xlabel('offset (m)')
ax2.set_ylabel('times (s)')
cutoff = 1
section = ax2.imshow(CMP_dummy, extent=[offsets.min(), offsets.max(), times.min(), times.max()],
aspect=690, cmap='Greys', vmin=-cutoff, vmax=cutoff, interpolation='nearest')
fig.tight_layout()
def anim_shot(i):
shot = i//frames
frame = i%frames
ax1.set_title('shot: {:d}'.format(shot + 1))
u = shots[shot][frame*every]
wavefield.set_array(u)
src.set_xdata(fontes[shot]*spacing)
rec.set_xdata(recep[shot]*spacing)
ax2.set_title('CMP: {:d}'.format(shot + 1))
CMP_dummy[:, :shot] = CMP[:, :shot]
CMP_dummy[:frame*every, shot] = CMP[:frame*every, shot]
section.set_array(CMP_dummy)
return wavefield, CMP_dummy
anim = FuncAnimation(fig, anim_shot, frames=frames*len(shots))
anim_to_html(anim, fps=6, dpi=60)
Na animação acima:
Course website: https://github.com/leouieda/geofisica2
Note: This notebook is part of the course "Geofísica 2" of Geology program of the Universidade do Estado do Rio de Janeiro. All content can be freely used and adapted under the terms of the Creative Commons Attribution 4.0 International License.