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
from fatiando.gravmag import prism
seaborn.set_context('talk')
print('Versão do Fatiando a Terra: {}'.format(fatiando.__version__))
A primeira tarefa servirá para ilustrar o conceito de Transformada de Fourier que vimos na aula teórica.
Como vimos em aula, uma função $h(t)$ pode ser escrita como:
$$ h(t) = \int\limits_{-\infty}^{\infty} H(f)\ e^{i 2 \pi f t}df $$em que $f$ é a frequência. $e^{i 2 \pi f t}$ são "coisas" que oscilam com frequência $f$ pois
$$ e^{i 2 \pi f t} = \cos(2 \pi f t) + i\sin(2 \pi f t) $$A equação da transformada nos diz que a função $h(t)$ pode ser escrita como uma "soma" infinita de oscilações de frequências diferentes. A oscilação de frequência $f$ possui uma amplitude $H(f)$. A função que descreve as amplitudes, $H(f)$, é a Transformada de Fourier. Podemos calcular a transformada $H(f)$ da função utilizando a fórmula abaixo:
$$ H(f) = \int\limits_{-\infty}^{\infty} h(t)\ e^{-i 2 \pi f t}dt $$$H(f)$ é um número complexo. É comum fazer um gráfico do módulo da transformada $|H(f)|$ pela frequência. Esse gráfico é chamado de espectro de amplitudes.
Rode a célula abaixo para gerar uma figura interativa.
def dois_senos(f1, A1, f2, A2):
sample = 1/200
t = np.arange(0, 10, sample)
n = t.size
s = A1*np.sin(2*np.pi*f1*t) + A2*np.sin(2*np.pi*f2*t)
f = np.fft.fftfreq(n, sample)[:n//2]
S = np.fft.fft(s)[:n//2]/n
fig = plt.figure(figsize=(12, 6))
ax = plt.subplot(211)
plt.title(r'$h(t) = A_1 \sin(2\pi f_1 t) + A_2 \sin(2\pi f_2 t)$', fontsize=16)
plt.plot(t, s, '-k')
plt.xlim(0, 5)
plt.ylim(-200, 200)
ax.set_xlabel('tempo (s)', labelpad=-15)
plt.ylabel('Amplitude')
ax = plt.subplot(212)
plt.title(r'$|H(f)|$', fontsize=16)
plt.vlines(f, [0], np.abs(S), colors='k', linewidth=3)
plt.hlines(0, f.min(), f.max(), colors='k', linewidth=5)
plt.xlim(0, 50)
plt.ylim(0, 60)
ax.set_xlabel('frequencia (Hz)', labelpad=-10)
plt.ylabel('Amplitude')
plt.tight_layout()
widgets.interactive(dois_senos,
f1=widgets.IntSlider(min=1, max=40, step=1, value=1),
A1=widgets.IntSlider(min=0, max=100, step=10, value=100),
f2=widgets.IntSlider(min=1, max=40, step=1, value=10),
A2=widgets.IntSlider(min=0, max=100, step=10, value=0))
Na gravimetria e magnetometria fazemos diversas operações que dependem das derivadas das anomalias. Precisamos de um jeito de calcular essas derivadas sem sabermos a função que produz a anomalia. Felizmente, as derivadas parciais podem ser calculadas usando a Transformada de Fourier.
Sabemos que a transformada de Fourier da derivada $W(f)$ é
$$ W(f) = H(f)\ i 2 \pi f $$Ou seja, podemos calcular a transformada da derivada através da transformada dos dados. Uma vez tendo a transformada da derivada, podemos calcular a derivada fazendo a transformada inversa:
$$ \dfrac{\partial h}{\partial t} = \int\limits_{-\infty}^{\infty} W(f)\ e^{i 2 \pi f t}df $$Rode a célula abaixo para gerar uma figura interativa.
def dois_senos_derivada(f2, A2):
f1, A1 = 1, 100
sample = 1/200
t = np.arange(0, 10, sample)
n = t.size
h = A1*np.sin(2*np.pi*f1*t) + A2*np.sin(2*np.pi*f2*t)
f = np.fft.fftfreq(n, sample)
H = np.fft.fft(h)
# A linha abaixo calcula a transformada da derivada
W = H*1j*2*np.pi*f
# e essa calcula a derivada através da transformada inversa
w = np.fft.ifft(W)
f = f[:n//2]
Habs = np.abs(H[:n//2]/n)
Wabs = np.abs(W[:n//2]/n)
fig, axes = plt.subplots(2, 2, figsize=(14, 7))
# Plot the function and spectrum
ax = axes[0][0]
ax.set_title(r'$h(t) = A_1 \sin(2\pi f_1 t) + A_2 \sin(2\pi f_2 t)$', fontsize=16)
ax.plot(t, h, '-k')
ax.set_xlim(0, 2.5)
ax.set_ylim(-200, 200)
ax.set_xlabel('t (s)', labelpad=-15)
ax.set_ylabel('Amplitude')
ax = axes[1][0]
ax.set_title(r'$|H(f)|$', fontsize=16)
ax.vlines(f, [0], Habs, colors='k', linewidth=3)
ax.hlines(0, f.min(), f.max(), colors='k', linewidth=5)
ax.set_xlim(0, 50)
ax.set_ylim(0, 60)
ax.set_xlabel('f (Hz)', labelpad=-10)
ax.set_ylabel('Amplitude')
# Plot the derivative and spectrum
ax = axes[0][1]
ax.set_title(r'$\partial h / \partial t$', fontsize=16)
ax.plot(t, w, '-k')
ax.set_xlim(0, 2.5)
ax.set_xlabel('t (s)', labelpad=-15)
ax = axes[1][1]
ax.set_title(r'$|W(f)|$', fontsize=16)
ax.vlines(f, [0], Wabs, colors='k', linewidth=3)
ax.hlines(0, f.min(), f.max(), colors='k', linewidth=5)
ax.set_xlim(0, 50)
ax.set_xlabel('f (Hz)', labelpad=-10)
plt.tight_layout()
widgets.interactive(dois_senos_derivada,
f2=widgets.IntSlider(min=1, max=40, step=1, value=10),
A2=widgets.IntSlider(min=0, max=50, step=10, value=0))
Vamos utilizar a tranformada de Fourier para calcular a derivada da anomalia de campo total ($\Delta T$) de um dique magnetizado pelo campo da Terra. Nosso dique está localizado em x = 0 e terá 100 m de espessura e 1 km de profundidade.
Rode a célula abaixo para gerar uma figura interativa.
areacubo = [-50, 50, 50, 1050]
cubo = mesher.Prism(areacubo[0], areacubo[1], -50000, 50000, areacubo[2], areacubo[3])
x1, x2 = -2500, 2500
xp = np.linspace(x1, x2, 200)
yp = np.zeros_like(xp)
zp = np.zeros_like(xp)
n = xp.size
sample = xp[1] - xp[0]
f = np.fft.fftfreq(n, sample)
def anomalia_derivada_cubo(inc, erro):
cubo.addprop('magnetization', utils.ang2vec(1, inc, 0))
tf = prism.tf(xp, yp, zp, [cubo], inc, 0)
if erro > 0:
tf = utils.contaminate(tf, erro)
TF = np.fft.fft(tf)
# A linha abaixo calcula a transformada da derivada
W = TF*1j*2*np.pi*f
# e essa calcula a derivada através da transformada inversa
w = np.fft.ifft(W)
espectro_tf = np.abs(TF[:n//2]/n)
espectro_deriv = np.abs(W[:n//2]/n)
fig, axes = plt.subplots(2, 2, figsize=(14, 7))
# Plot the function and spectrum
ax = axes[0][0]
ax.set_title(r'$\Delta T$', fontsize=16)
ax.plot(xp, tf, '-k')
ax.set_xlim(xp.min(), xp.max())
ax.set_xlabel('x (m)')
ax.set_ylabel(r'nT')
ax = axes[1][0]
ax.set_title(r'$|H(f)|$', fontsize=16)
ax.vlines(f[:n//2], [0], espectro_tf, colors='k', linewidth=3)
ax.hlines(0, 0, f.max(), colors='k', linewidth=5)
ax.set_xlim(0, f.max())
ax.set_xlabel('f (1/m)')
ax.set_ylabel('Amplitude')
# Plot the derivative and spectrum
ax = axes[0][1]
ax.set_title(r'$\partial \Delta T / \partial x$', fontsize=16)
ax.plot(xp, w, '-k')
ax.set_xlim(xp.min(), xp.max())
ax.set_xlabel('x (m)')
ax.set_ylabel(r'nT/m')
ax = axes[1][1]
ax.set_title(r'$|W(f)|$', fontsize=16)
ax.vlines(f[:n//2], [0], espectro_deriv, colors='k', linewidth=3)
ax.hlines(0, 0, f.max(), colors='k', linewidth=5)
ax.set_xlim(0, f.max())
ax.set_xlabel('f (1/m)')
plt.tight_layout()
widgets.interactive(anomalia_derivada_cubo,
inc=widgets.IntSlider(min=-90, max=90, step=15, value=-45),
erro=widgets.IntSlider(min=0, max=20, step=1, value=0))