Course website: http://www.leouieda.com/geofisica1
Note: This notebook is part of the course "Geofísica 1" 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.
Como na prática 5, esse documento que você está usando é um IPython notebook. Esta prática usará a parte de inversão 2D da biblioteca Fatiando a Terra de modelagem geofísica. Para isso, precisamos carregar (import
) os módulos que vamos usar e também o módulo numpy.
O notebook é divido em células (como esta). Para editar o conteúdo de uma célula, clique nela (clique nesta para editar esse texto). Para executar uma célula, aperte Shift + Enter
. Execute as duas células abaixo.
from IPython.core.pylabtools import print_figure
from IPython.html import widgets
from IPython.display import Image, display
from fatiando.gravmag import talwani
from fatiando.gravmag.basin2d import PolygonalBasinGravity
from fatiando.gravmag.interactive import Moulder
from fatiando.inversion.regularization import Damping, Smoothness1D, Smoothness2D
from fatiando.vis import mpl
from fatiando import utils
import fatiando
import numpy as np
print('Versão do Fatiando a Terra: {}'.format(fatiando.__version__))
Faça um modelo de uma bacia sedimentar e passe para o colega. Algumas instruções:
area = [0, 100e3, 0, 10e3]
x = np.linspace(area[0], area[1], 80)
z = -np.ones_like(x)
meu_modelo = Moulder(area, x, z, density_range=[-1000, 0])
meu_modelo.run()
Rode a célula abaixo para salvar o seu modelo. Mude o nome do arquivo e troque com um colega.
meu_modelo.save('modelo_bacia')
Primeiro precisamos pegar o modelo do colega e extrair os dados (sem ruído), o polígono que descreve a bacia e a densidade correta.
Mude o nome do arquivo 'modelo_bacia'
abaixo para o nome do arquivo do colega.
modelo = Moulder.load('modelo_bacia')
modelo.plot(figsize=(12, 6))
dado = modelo.predicted
bacia = modelo.model[0]
densidade = bacia.props.copy()
A célula abaixo prepara a inversão. pontos
define quantos pontos terá o polígono que descreverá a bacia. Esse será o número de parâmetros que a inversão estimará (a profundidade de cada ponto do polígono).
Lembrando que a inversão busca achar o mínimo da função do ajuste:
$$\phi(\bar{p}) = \sum\limits_{i=1}^N (d^o_i - d_i)^2$$em que $\bar{p}$ é um vetor com os parâmetros (profundidade de cada ponto do polígono), $d^o_i$ é o i-ésimo dado observado e $d_i$ é o dado predito pelos parâmetros.
Essa inversão é do tipo não-linear. Isso quer dizer que precisamos espeficicar um chute inicial para os parâmetros. A inversão é feita em passos a partir desse chute inicial. Vamos usar "todos os parâmetros igual a 1 km".
pontos = 40
topo = bacia.y.min()
p0 = 1000*np.ones(pontos)
xlim = [bacia.x.min(), bacia.x.max()]
misfit = PolygonalBasinGravity(x, z, dado, pontos, bacia.props, top=topo, xlim=xlim).config('levmarq', initial=p0)
A célula abaixo cria algumas funções que usaremos para fazer um gráfico da nossa solução.
def plota_resultado(solver):
fig = mpl.figure(figsize=(12, 4), facecolor='white')
ax = mpl.subplot(2, 1, 1)
mpl.plot(x, solver.data, 'ok')
mpl.plot(x, solver.predicted(), '-r', linewidth=2)
ax.set_ylim(1.2*dado.min(), 1.1*dado.max())
ax.grid(True)
ax.set_xticklabels([])
ax.set_ylabel('Anomalia (mGal)')
ax = mpl.subplot(2, 1, 2)
mpl.polygon(bacia, style='o-k', fill='gray', alpha=0.5, linewidth=2)
mpl.polygon(solver.estimate_, style='o-r', linewidth=2)
ax.set_xlabel('x (km)')
ax.set_ylabel('z (km)')
mpl.set_area(area)
mpl.m2km()
ax.grid(True)
ax.invert_yaxis()
mpl.tight_layout(pad=0)
return fig
Rode a célula abaixo para criar um gráfico interativo. A parte de cima do gráfico mostra os dados observados (pontos) e preditos pela inversão (linha vermelha). A parte de baixo mostra o seu modelo de bacia (preto) e o resultado da inversão (vermelho).
O botão "erro" pode ser movido para mudar o valor do erro que será aplicado aos dados. Cada vez que você mover o botão, o dado é alterado e inversão rodada novamente.
O botão "seed" altera ligeiramente o erro que é aplicado ao dado (mas não o tamanho do erro).
def interativo_sem_vinculo(erro, seed):
misfit.data = utils.contaminate(dado, erro, seed=seed)
misfit.fit()
fig = plota_resultado(misfit)
mpl.close(fig)
display(Image(print_figure(fig, dpi=80)))
misfit.data = dado
widgets.interact(interativo_sem_vinculo,
erro=widgets.FloatSliderWidget(min=0, max=1, step=0.1, value=0),
seed=widgets.IntSliderWidget(min=0, max=10, step=1, value=0))
Vamos utilizar o vínculo de "norma mínima" ou "damping".
Lembre que adicionar vínculos é equivalente a somar uma outra função na nossa função do ajuste $\phi(\bar{p})$:
$$\Gamma(\bar{p}) = \sum\limits_{i=1}^N (d^o_i - d_i)^2 + \mu \sum\limits_{j=1}^M p_j^2$$A função do vínculo é $\sum\limits_{j=1}^M p_j^2$ que impõe que os parâmetros estimados sejam o menor possível. Em outras palavras, eles devem ter norma mínima.
A constante positiva $\mu$ determina qual é o peso do vínculo. $\mu$ cuida do balanço entre ajustar os dados e obedecer aos vínculos.
mi
controla o valor de $\mu$ (um valor de -10
significa que $\mu = 10^{-15}$). Use o botão para aumentar a importância do vínculo de norma mínima.def interativo_damping(erro, seed, mi):
misfit.data = utils.contaminate(dado, erro, seed=seed)
solver = misfit + 10**mi*Damping(misfit.nparams)
solver.fit()
fig = plota_resultado(solver)
mpl.close(fig)
display(Image(print_figure(fig, dpi=80)))
misfit.data = dado
widgets.interact(interativo_damping,
erro=widgets.FloatSliderWidget(min=0, max=1, step=0.1, value=0),
seed=widgets.IntSliderWidget(min=0, max=10, step=1, value=0),
mi=widgets.FloatSliderWidget(min=-15, max=0, step=0.5, value=-15))
Agora vamos utilizar o vínculo de suavidade.
Nesse caso, o vínculo é:
$$\sum\limits_{j=1}^{M - 1} (p_{j+1} - p_j)^2$$O vínculo de suavidade impõe que os parâmetros sejam o mais próximos o possível de seus vizinhos imediatos.
mi
para aumentar a importância do vínculo. Se der uma mensagem de erro, aumente mi
.mi
)def interativo_suave(erro, seed, mu):
misfit.data = utils.contaminate(dado, erro, seed=seed)
solver = misfit + 10**mu*Smoothness1D(misfit.nparams)
solver.fit()
fig = plota_resultado(solver)
mpl.close(fig)
display(Image(print_figure(fig, dpi=80)))
misfit.data = dado
widgets.interact(interativo_suave,
erro=widgets.FloatSliderWidget(min=0, max=5, step=0.1, value=0),
seed=widgets.IntSliderWidget(min=0, max=10, step=1, value=0),
mu=widgets.FloatSliderWidget(min=-15, max=0, step=0.5, value=-15))
Em casos reais não saberemos exatamente qual é a densidade dos sedimentos. Podemos ter uma ideia mas valores de propriedades físicas costumam variar bastante.
Vamos assumir que não sabemos a densidade correta e colocar isso como uma variável que podemos controlar.
def interativo_densidade(seed, mu, densidade):
misfit.data = utils.contaminate(dado, 0.5, seed=seed)
misfit.model['props']['density'] = densidade
solver = misfit + 10**mu*Smoothness1D(misfit.nparams)
solver.fit()
fig = plota_resultado(solver)
mpl.close(fig)
display(Image(print_figure(fig, dpi=80)))
misfit.data = dado
misfit.model['props']['density'] = bacia.props['density']
widgets.interact(interativo_densidade,
seed=widgets.IntSliderWidget(min=0, max=10, step=1, value=0),
mu=widgets.FloatSliderWidget(min=-10, max=0, step=0.5, value=-10),
densidade=widgets.FloatSliderWidget(min=-1000, max=-50, step=20, value=-50))
Para verificar se conseguimos determinar a densidade, vamos imprimir o valor correto abaixo:
print('Densidade correta = {} kg/m3'.format(densidade['density']))