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.
Esse documento que você está usando é 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.
%matplotlib inline
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
from IPython.display import display
import seaborn
import fatiando
from fatiando import mesher, utils, gridder
from fatiando.gravmag import prism, transform
from fatiando.vis import myv, mpl
seaborn.set_context('talk')
print('Versão do Fatiando a Terra: {}'.format(fatiando.__version__))
A Transformada de Fourier nos permite calcular derivadas dos nossos dados. Uma transformação muito utilizada na magnetometria é a Amplitude da Derivada Total, em inglês Total Gradient Amplitude (TGA):
$$ TGA = \sqrt{\left(\frac{\partial \Delta T}{\partial x}\right)^2 + \left(\frac{\partial\Delta T}{\partial y}\right)^2 + \left(\frac{\partial\Delta T}{\partial z}\right)^2} $$O TGA é famoso por concentrar a anomalia magnética de campo total sobre o corpo causador da anomalia.
As derivadas parciais na equação acima podem ser calculadas da anomalia de campo total ($\Delta T$) usando a Transformada de Fourier, como vimos na aula passada. Para lembrar, podemos calcular a derivada de uma função $h(t)$ através de sua transformada $H(f)$
$$ W(f) = H(f) i 2 \pi f $$em que $W(f)$ é a Transformada de Fourier de $\partial h/ \partial t$. Uma conhecendo a transformada, podemos calcular a derivada através da transformada inversa de Fourier
$$ \frac{\partial h}{\partial t} = \int\limits_{-\infty}^{\infty} W(f) e^{i 2 \pi f t} df $$A Amplitude da Derivada Total é erroneamente chamada de Sinal Analítico 3D (3D Analytic Signal Amplitude - ASA). Esse nome é errado (Reid, 2012). Na verdade o que temos é o tamanho (amplitude) do gradiente da anomalia (vetor das derivadas em cada direção).
Um mito persistente é o TGA não depende da direção do campo da Terra ou da magnetização do corpo (inclinação e declinação). Diversos autores já mostraram que esse não é o caso (ver Li, 2006).
Vamos explorar a cara do TGA causado por um modelo simples de 1 prisma alongado. Vocês poderão controlar a inclinação e declinação do campo magnético da Terra e o erro aleatório aplicado ao dado. Assumiremos somente magnetização induzida.
A célula abaixo cria nosso modelo de 1 prisma e gera uma figura 3D. Feche a figura antes de continuar.
shape = (100, 100)
x, y, z = gridder.regular((-5000, 5000, -5000, 5000), shape, z=0)
dx, dy = 500, 5000
cubo = mesher.Prism(-dy/2, dy/2, -dx/2, dx/2, 400, 4000)
cubo_area = cubo.get_bounds()[:4][::-1]
bounds = (-5000, 5000, -5000, 5000, 0, 5000)
scene = myv.figure()
myv.prisms([cubo])
myv.outline(bounds)
myv.wall_bottom(bounds)
myv.wall_west(bounds)
oa = myv.mlab.orientation_axes()
oa.axes.x_axis_label_text = 'N'
oa.axes.y_axis_label_text = 'E'
oa.text_property.color = (0.0, 0.0, 0.0)
myv.show()
def derivada_total(inc, dec, erro):
tf = prism.tf(x, y, z, [cubo], inc, dec, pmag=utils.ang2vec(1, inc, dec))
if erro > 0:
tf = utils.contaminate(tf, erro, seed=0)
total = transform.tga(x, y, tf, shape)
fig, axes = plt.subplots(1, 2, figsize=(14, 9))
for ax, data, title in zip(axes.ravel(), [tf, total], ['Anomalia', 'Derivada Total']):
ax.set_aspect('equal')
plt.sca(ax)
plt.title(title)
mpl.square(cubo_area, style='--w')
scale = np.abs([data.min(), data.max()]).max()
plt.contourf(y.reshape(shape), x.reshape(shape), data.reshape(shape),
30, cmap='RdBu_r', vmin=-scale, vmax=scale)
cb = plt.colorbar(orientation='horizontal', pad=0.01, aspect=50, shrink=0.9)
cb.set_label('nT/m' if title != "Anomalia" else "nT")
ax.set_xticks([])
ax.set_yticks([])
plt.tight_layout()
widgets.interactive(derivada_total,
inc=widgets.IntSlider(min=-90, max=90, step=5, value=45),
dec=widgets.IntSlider(min=-90, max=90, step=5, value=0),
erro=widgets.FloatSlider(min=0, max=20, step=1, value=0))
A anomalia magnética é muito complicada de interpretar. Por isso existem tantos filtros e transformações. A melhor situação é quando a inclinação do campo magnético é $\pm90^\circ$ (nos polos). Nessas ocasiões, a anomalia se concentra sobre o corpo causador.
Felizmente, há um jeito de calcular como a anomalia que medimos ficaria se estivesse nos polos. A técnica chama-se (numa explosão de criatividade) "Redução ao Polo". Um dos jeitos de se calcular a redução ao polo é usando a transformada de Fourier.
É necessário conhecer a direção de magnetização do corpo para aplicar a redução ao polo. Isso é fácil se o corpo tiver somente magnetização induzida pelo campo geomagnético. A magnetização é paralela ao campo da Terra. A situação complica quando há magnetização remanente, aquela que os minerais ferromagnéticos guardam quando se resfriam abaixo da temperatura de Curie. Por isso a redução ao polo não é tão facilmente utilizada quanto o TGA.
Vamos utilizar a mesma simulação que fizemos acima para testar a redução ao polo. Você controlará a inclinação e declinação do campo da Terra (inc
e dec
). Vamos assumir somente magnetização induzida, logo a inclinação e declinação que será utilizada na redução ao polo é a mesma que a da Terra.
shape = (100, 100)
x, y, z = gridder.regular((-5000, 5000, -5000, 5000), shape, z=0)
dx, dy = 500, 5000
cubo = mesher.Prism(-dy/2, dy/2, -dx/2, dx/2, 400, 4000)
cubo_area = cubo.get_bounds()[:4][::-1]
def reducao_polo(erro, inc, dec):
tf = prism.tf(x, y, z, [cubo], inc, dec, pmag=utils.ang2vec(5, inc, dec))
if erro > 0:
tf = utils.contaminate(tf, erro, seed=0)
rtp = transform.reduce_to_pole(x, y, tf, shape, inc, dec, inc, dec)
fig, axes = mpl.subplots(1, 2, figsize=(14, 9))
for ax, data, title in zip(axes.ravel(), [tf, rtp], ['Anomalia', 'Reduzida ao polo']):
ax.set_aspect('equal')
plt.sca(ax)
plt.title(title)
mpl.square(cubo_area, style='--w')
scale = np.abs([data.min(), data.max()]).max()
plt.contourf(y.reshape(shape), x.reshape(shape), data.reshape(shape),
30, cmap='RdBu_r', vmin=-scale, vmax=scale)
cb = plt.colorbar(orientation='horizontal', pad=0.01, aspect=50, shrink=0.9)
cb.set_label('nT')
ax.set_xticks([])
ax.set_yticks([])
plt.tight_layout()
widgets.interactive(reducao_polo,
erro=widgets.FloatSlider(min=0, max=20, step=1, value=0),
inc=widgets.IntSlider(min=-90, max=90, step=5, value=45),
dec=widgets.IntSlider(min=-90, max=90, step=5, value=0))
Em gravimetria e magnetometria, uma transformação que pode ser feita com os dados é a "continuação para cima". Essa transformação nos permite calcular como a anomalia medida seria se estivesse a uma altitude maior. Um dos jeitos de calcular a continuação é através da transformada de Fourier:
$$ H_{up} = H e^{-\Delta z |k|} $$em que $H$ é a transformada dos dados a uma altitude $z$, $H_{up}$ é a transformada dos dados a uma altitude $z + \Delta z$ e $|k|$ são as frequências.
Vamos testar a continuação para cima em uma anomalia causada por dois corpos: um grande e profundo e outro razo e pequeno.
shape = (100, 100)
x, y, z = gridder.regular((-6000, 6000, -6000, 6000), shape, z=0)
modelo = [mesher.Prism(-200, 200, -200, 200, 100, 500),
mesher.Prism(-4000, 4000, -4000, 4000, 1000, 5000)]
inc, dec = 45, -10
mag = utils.ang2vec(1, inc, dec)
modelo[0].props['magnetization'] = 1*mag
modelo[1].props['magnetization'] = 1*mag
bounds = (-6000, 6000, -6000, 6000, 0, 6000)
myv.figure()
myv.prisms(modelo)
myv.outline(bounds)
myv.wall_bottom(bounds)
myv.wall_west(bounds)
oa = myv.mlab.orientation_axes()
oa.axes.x_axis_label_text = 'N'
oa.axes.y_axis_label_text = 'E'
oa.text_property.color = (0.0, 0.0, 0.0)
myv.show()
def continuacao(altitude, erro):
cubo_area = modelo[0].get_bounds()[:4][::-1]
tf = prism.tf(x, y, z, modelo, inc, dec)
if erro > 0:
tf = utils.contaminate(tf, erro, seed=0)
cont = transform.upcontinue(x, y, tf, shape, altitude + 1e-5)
fig, axes = mpl.subplots(1, 2, figsize=(14, 9))
titles = ['Anomalia em 0 m', 'Continuada para {:.0f} m'.format(altitude)]
for ax, data, title in zip(axes.ravel(), [tf, cont], titles):
ax.set_aspect('equal')
plt.sca(ax)
plt.title(title)
for prisma in modelo:
mpl.square(prisma.get_bounds()[:4][::-1], style='--w')
scale = np.abs([data.min(), data.max()]).max()
plt.contourf(y.reshape(shape), x.reshape(shape), data.reshape(shape),
30, cmap='RdBu_r', vmin=-scale, vmax=scale)
cb = plt.colorbar(orientation='horizontal', pad=0.01, aspect=50, shrink=0.9)
cb.set_label('nT')
ax.set_xticks([])
ax.set_yticks([])
plt.tight_layout()
widgets.interactive(continuacao,
altitude=widgets.FloatSlider(min=0, max=4000, step=50, value=0),
erro=widgets.FloatSlider(min=0, max=20, step=1, value=0))