#!/usr/bin/env python
# coding: utf-8
# # 8 - Sísmica de reflexão: Resolução e armadilhas na interpretação
#
# Nessa prática, vamos simular diversos cenários geológicos e prever como seria uma sessão sísmica já empilhada e em tempo. Você vai poder desenhar seus modelos em um programa estilo "MS Paint", designar velocidades para cada litologia e ver o resultado esperado da sísmica para cada situação.
#
# Exemplos de modelos e a inspiração para essa aula vieram das aulas de ["Pitfalls in Seismic Interpretation"](http://pages.geo.wvu.edu/~wilson/geol554.htm) do [Prof. Thomas H. Wilson](http://pages.geo.wvu.edu/~wilson/).
#
# Recomendo também a leitura do texto [Tuning Geology](https://www.agilegeoscience.com/blog/2011/7/8/tuning-geology.html) do Matt Hall.
# ## Objetivos
#
# * Familiarizar como o processo de modelagem direta para testar hipóteses
# * Observar o efeito de anomalias estáticas e de velocidade no dado sísmico
# * Aprender sobre situações que geralmente induzem a erros na interpretação
# ## Questão para entregar
#
#
#
# Explique os detalhes aos quais devemos nos atentar quando interpretando uma sísmica em tempo.
#
#
#
# **Dicas**:
#
# * Fale sobre o efeito da resolução quando interpretamos camadas finas.
# * Como as velocidades podem afetar a forma dos horizontes?
#
#
# ### Regras para a resposta
#
# * Coloque **nome, data e o número da prática** em sua resposta.
# * A resposta pode ter no **máximo 1 página** (não uma folha).
# * **Execute o notebook** antes de responder. As simulações abaixo foram feitas para te ajudar.
# * **Pense e organize** sua resposta andtes de começar a escrever.
# ## Instruções
#
# Esse documento é um [Jupyter notebook](http://jupyter.org/), 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](http://python.org/). 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.
# ## Setup
#
# Rode as células abaixo para carregar os módulos necessários para essa prática.
# In[ ]:
get_ipython().run_line_magic('matplotlib', 'inline')
from __future__ import division, print_function
import math
from PIL import Image
import numpy as np
import scipy.misc
import matplotlib.pyplot as plt
import ipywidgets as ipw
from fatiando.seismic import conv
from fatiando.vis.mpl import seismic_image, seismic_wiggle
# ## Resolução vertical e camadas finas
#
# Vimos em aula que a resolução vertical máxima do dado sísmico é aproximadamente $\lambda/4$ e
#
# $$\lambda = \frac{V}{f}$$
#
# Vamos explorar agora como isso afeta a nossa interpretação de uma camada fina. No exemplo abaixo, podemos variar a espessura e velocidade da camada e a frequência da onda sísmica.
# In[ ]:
def camada_fina(velocidade, espessura, frequencia):
n_samples, n_traces = [400, 100]
dt, dz = 1e-4, 1
velocity = 3000*np.ones((n_samples, n_traces))
top = 150
bottom = top + int(np.floor(espessura/dz))
velocity[top : bottom, :] = velocidade
# Gardner's relation: https://en.wikipedia.org/wiki/Gardner%27s_relation
rho = 0.31*(velocity**0.25)
vel_l = conv.depth_2_time(velocity, velocity, dt=dt, dz=dz)
rho_l = conv.depth_2_time(velocity, rho, dt=dt, dz=dz)
rc = conv.reflectivity(vel_l, rho_l)
data = conv.convolutional_model(rc, frequencia, conv.rickerwave, dt=dt)
fig, axes = plt.subplots(1, 2, figsize=(8, 5))
ax = axes[0]
ax.set_title("Modelo (em profundidade)")
tmp = ax.imshow(velocity, extent=[0, n_traces, n_samples*dz, 0],
cmap="viridis", aspect='auto', origin='upper',
vmin=1000, vmax=5000)
fig.colorbar(tmp, ax=ax, pad=0, aspect=50)
ax.set_ylabel('Profundidade (m)')
ax = axes[1]
ax.set_title("Sismograma")
seismic_image(data, dt=dt, cmap="RdBu_r", aspect='auto')
ax.set_ylabel('Tempo (s)')
ax.set_ylim(0.3, 0)
plt.tight_layout()
ipw.interactive(camada_fina,
velocidade=ipw.FloatSlider(min=2000, max=5000, step=200, value=1000),
espessura=ipw.FloatSlider(min=5, max=150, step=5, value=100),
frequencia=ipw.FloatSlider(min=10, max=60, step=5, value=30))
# ### Para pensar
#
# * O que acontece com o tempo de chegada da reflexão no topo da camada quando você diminui a espessura? Isso deveria acontecer?
# * Como a velocidade e a frequência influenciam sua capacidade de ver o topo e a base da camada separados?
# ## Cenários geológicos
#
# Agora, vamos modelar diferentes cenários geológicos utilizando um programa do tipo MS Paint.
#
# Rode as células abaixo para criar as funções que vamos precisar para nossa modelagem.
# In[ ]:
def img_template(fname, pixel_thresh=10, return_colors=False):
img = Image.open(fname)
count, colors = zip(*[[n, c] for n, c in img.getcolors() if n > pixel_thresh])
sort = np.argsort(count)
colors = np.array(colors)[sort][::-1]
data = scipy.misc.fromimage(img)
template = np.zeros(data.shape[:2], dtype=np.int)
for i, c in enumerate(colors):
template[np.all(data == c, axis=2)] = i
# For now, any pixel not of the top colors (eliminated by pixel_thresh)
# is assigned index 0. A better way would be interpolate or copy the neighbors.
if return_colors:
return template, colors
else:
return template
# In[ ]:
def painter_widget(template, vmin, vmax, step, cmap='copper_r', prop='V', plot_func=None):
colors = np.unique(template)
args = dict(min=vmin, max=vmax, step=step)
sliders = dict([['{}{}'.format(prop, i + 1), ipw.FloatSlider(**args)]
for i in range(len(colors))])
model = np.empty_like(template)
def plot_func(vp):
rho = 0.31*(vp**0.25)
dz = 2000/vp.shape[0]
dt = 1e-3
f = 15
vp_t = conv.depth_2_time(vp, vp, dt=dt, dz=dz)
rho_t = conv.depth_2_time(vp, rho, dt=dt, dz=dz)
rc = conv.reflectivity(vp_t, rho_t)
data = conv.convolutional_model(rc, f, conv.rickerwave, dt=dt)
times = np.arange(vp_t.shape[0])*dt
# Plots
fig = plt.figure(figsize=(12, 4))
ax = plt.subplot(121)
plt.imshow(vp, cmap='copper_r', extent=[0, vp.shape[1], vp.shape[0]*dz, 0],
vmin=vmin, vmax=vmax, aspect='auto')
plt.colorbar(pad=0.01, aspect=40).set_label('Velocidade (m/s)')
ax.set_title('Modelo')
ax.set_ylabel('profundidade')
ax = plt.subplot(122)
seismic_image(data, dt=dt, cmap='RdBu_r', aspect='auto')
ax.set_title(u'Sessão sísmica')
ax.set_ylabel('tempo (s)')
plt.tight_layout()
return fig
def callback(**kwargs):
for v, c in zip(kwargs, colors):
model[template == c] = kwargs[v]
plot_func(model)
return ipw.interactive(callback, **sliders)
# ### Desenhando modelos
#
# Utilize o Pinta para desenhar seus modelos:
#
# 1. Crie uma figura de **300 x 150 pixels** (clique em "New Figure").
# 2. Utilize as ferramentas de retângulo e pincel para desenhar.
# 3. Coloque **cada litologia** com uma cor diferente. Não importa qual.
# 4. Salve sua figura como **.png** na mesma pasta que está o notebook.
#
# ### Modelos
#
# Você deve criar 3 modelos:
#
# 1. Um domo de sal com as camadas sedimentares ao redor e camadas planas abaixo dele ([exemplo](../data/modelo-sal.png)).
# 2. Uma cada irregular na superfície com camadas plano-paralelas abaixo dela ([exemplo](../data/modelo-camada-intemperismo.png)).
# 3. Um carbonato (recife) inserido no meio de camadas plano-paralelas (faça o topo dele ser um pouco acima da camada na qual está inserido) ([exemplo](../data/modelo-carbonato.png)).
#
# Coloque o nome da sua figura nas células abaixo e execute as. Cada uma será uma figura interativa na qual você pode controlar as velocidades de cada camada em seu modelo.
#
# Utilize a tabela abaixo (retirada de Kearey et al, 2002) para escolher as velocidades para seu modelo.
#
# 
#
# In[ ]:
template_sal = img_template('NOME-DA-SUA-FIGURA.png')
painter_widget(template_sal, 1000, 8000, 100)
# In[ ]:
template_carbonato = img_template('NOME-DA-SUA-FIGURA.png')
painter_widget(template_carbonato, 1000, 8000, 100)
# In[ ]:
template_camada_irregular = img_template('NOME-DA-SUA-FIGURA.png')
painter_widget(template_camada_irregular, 1000, 8000, 100)
# ### Para pensar
#
# * O que acontece com as camadas abaixo do sal?
# * O que acontece com os cantos das camadas que foram deformadas pelo sal?
# * O que acontece com as camadas abaixo da camada irregular?
# * Como você interpretaria o dado do carbonado se não soubesse que ali tem um carbonato?
# ## Referências
#
# Kearey, P., M. Brooks, and I. Hill (2002), An Introduction to Geophysical Exploration, 3 edition., Wiley-Blackwell, Malden, MA.
# ## License and information
#
# **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](http://www.uerj.br/).
# All content can be freely used and adapted under the terms of the
# [Creative Commons Attribution 4.0 International License](http://creativecommons.org/licenses/by/4.0/).
#
# 