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.
Objetivos:
Esse documento que você está usando é um IPython notebook. É um documento interativo que mistura texto (como esse), código (como abaixo), e o resultado de executar o código (que pode ser números, texto, figuras, videos, etc). Esta prática usará a biblioteca Fatiando a Terra de modelagem geofísica 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.
%matplotlib inline
from __future__ import division
from IPython.html import widgets
import numpy as np
from fatiando.gravmag import prism, fourier
from fatiando import utils, gridder, mesher
from fatiando.vis import mpl
import fatiando
mpl.rc('lines', linewidth=2)
mpl.rc('font', size=12)
:0: FutureWarning: IPython widgets are experimental and may change in the future.
print('Versão do Fatiando a Terra: {}'.format(fatiando.__version__))
Versão do Fatiando a Terra: 0.3-237-g65602ba
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 contínua 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 $e^{i 2 \pi f t}$ são "coisas" que oscilam com frequência $f$. A equação acima 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$ é ponderada por uma amplitude $H(f)$. A função que descreve as amplitudes, $H(f)$, é a Transformada de Fourier:
$$ 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(freq1, amp1, freq2, amp2):
sample = 1/200
t = np.arange(0, 10, sample)
n = t.size
s = amp1*np.sin(2*np.pi*freq1*t) + amp2*np.sin(2*np.pi*freq2*t)
f = np.fft.fftfreq(n, sample)[:n//2]
S = np.fft.fft(s)[:n//2]/n
fig = mpl.figure(figsize=(12, 5))
ax = mpl.subplot(211)
mpl.title(r'$h(t) = {:.0f}\ \sin(2\pi\ {:.0f}\ t) + {:.0f}\ \sin(2\pi\ {:.0f}\ t)$'.format(amp1, freq1, amp2, freq2),
fontsize=18)
mpl.plot(t, s, '-k')
mpl.xlim(0, 5)
mpl.ylim(-200, 200)
mpl.grid()
ax.set_xlabel('tempo (s)', labelpad=-15)
mpl.ylabel('Amplitude')
ax = mpl.subplot(212)
mpl.plot(f, np.abs(S), '-k')
mpl.xlim(0, 50)
mpl.ylim(0, 50)
mpl.grid()
ax.set_xlabel('frequencia (Hz)', labelpad=-10)
mpl.ylabel(r'$|H(f)|$', fontsize=14)
mpl.tight_layout(h_pad=0)
widgets.interactive(dois_senos,
freq1=widgets.IntSliderWidget(min=1, max=40, step=1, value=1),
amp1=widgets.IntSliderWidget(min=0, max=100, step=10, value=0),
freq2=widgets.IntSliderWidget(min=1, max=40, step=1, value=1),
amp2=widgets.IntSliderWidget(min=0, max=100, step=10, value=100))
amp1 = 0
e fixe freq2
em algum valor. Varie somente amp2
.amp2
varia?amp1 = 0
e amp2=100
. Aumente freq2
.amp2 = 100
e freq2=1
. Aumente amp1
e freq1
.A Transformada de Fourier é muito útil na geofísica para fazer transformações nos dados. Em gravimetria e magnetometria, uma transformação comum é 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.
Rode a célula abaixo para gerar uma figura interativa.
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
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 = fourier.upcontinue(x, y, tf, shape, altitude)
fig, axes = mpl.subplots(1, 2, figsize=(14, 5.5))
for ax, data, title in zip(axes.ravel(), [tf, cont], ['Anomalia em 0 m', 'Continuada para {:.0f} m'.format(altitude)]):
ax.set_aspect('equal')
mpl.sca(ax)
mpl.title(title)
scale = np.abs([data.min(), data.max()]).max()
mpl.contourf(y, x, data, shape, 30, cmap=mpl.cm.RdBu_r, vmin=-scale, vmax=scale)
mpl.colorbar(pad=0).set_label('nT')
mpl.m2km()
mpl.tight_layout(h_pad=0, w_pad=0)
widgets.interactive(continuacao,
altitude=widgets.FloatSliderWidget(min=0, max=5000, step=50, value=0),
erro=widgets.FloatSliderWidget(min=0, max=20, step=1, value=0))
Outra operação que a Transformada de Fourier nos permite fazer é calcular derivadas dos nossos dados. Uma transformação muito utilizada na magnetometria é a Derivada Total (DT):
$$ DT = \sqrt{\left(\frac{\partial T}{\partial x}\right)^2 + \left(\frac{\partial T}{\partial y}\right)^2 + \left(\frac{\partial T}{\partial z}\right)^2} $$A DT é famosa 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 ($T$) usando a Transformada de Fourier. Derivando a equação de $h(t)$ acima:
$$ \frac{\partial h}{\partial t} = \frac{\partial}{\partial t}\left(\int\limits_{-\infty}^{\infty} H(f)\ e^{i 2 \pi f t}df\right) $$Como a única coisa que depende de $t$ é $e^{i 2 \pi f t}$, podemos passar a derivada para dentro da integral:
$$ \frac{\partial h}{\partial t} = \int\limits_{-\infty}^{\infty} H(f)\ \frac{\partial e^{i 2 \pi f t}}{\partial t}df $$e aplicando a regra da cadeia
$$ \frac{\partial h}{\partial t} = \int\limits_{-\infty}^{\infty} H(f)\ i 2 \pi f\ e^{i 2 \pi f t} df $$$\frac{\partial h}{\partial t}$ também é uma função e pode ser escrita como a soma de infinitas oscilações. Se chamarmos a transformada de Fourier de uma função $h$ de $F\{h\}$, podemos escrever a derivada de $h$ como:
$$ \frac{\partial h}{\partial t} = \int\limits_{-\infty}^{\infty} F\left\{\frac{\partial h}{\partial t}\right\}e^{i 2 \pi f t} df $$Se compararmos as duas equações acima, vemos que
$$ F\left\{\frac{\partial h}{\partial t}\right\} = H(f)\ i 2 \pi f $$Ou seja, a transformada de Fourier da derivada é $i 2 \pi f$ vezes a transformada de Fourier da função. Uma vez tendo a transformada da derivada, podemos calcular a derivada.
Rode a célula abaixo para gerar uma figura interativa.
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 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 = fourier.ansig(x, y, tf, shape)
fig, axes = mpl.subplots(1, 2, figsize=(14, 6))
for ax, data, title in zip(axes.ravel(), [tf, total], ['Anomalia', 'Derivada Total']):
ax.set_aspect('equal')
mpl.sca(ax)
mpl.title(title)
mpl.square(cubo_area)
scale = np.abs([data.min(), data.max()]).max()
mpl.contourf(y, x, data, shape, 30, cmap=mpl.cm.RdBu_r, vmin=-scale, vmax=scale)
if title != "Anomalia":
mpl.colorbar(pad=0).set_label('nT/m')
else:
mpl.colorbar(pad=0).set_label('nT')
mpl.m2km()
mpl.tight_layout(h_pad=0, w_pad=0)
widgets.interactive(derivada_total,
inc=widgets.IntSliderWidget(min=-90, max=90, step=5, value=45),
dec=widgets.IntSliderWidget(min=-90, max=90, step=5, value=0),
erro=widgets.FloatSliderWidget(min=0, max=20, step=1, value=0))
inc
e dec
) e o erro colocado nos dados originais (da esquerda).Vimos na prática 1 que a anomalia magnética é muito complicada de interpretar. A melhor situação é quando a inclinação do campo magnético é $\pm90^\circ$ (nos polos).
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 a Derivada Total.
Rode a célula abaixo para gerar uma figura interativa.
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(inc, dec, erro):
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 = fourier.reduce_to_pole(x, y, tf, shape, inc, dec)
fig, axes = mpl.subplots(1, 2, figsize=(14, 6))
for ax, data, title in zip(axes.ravel(), [tf, rtp], ['Anomalia', 'Reduzida ao polo']):
ax.set_aspect('equal')
mpl.sca(ax)
mpl.title(title)
mpl.square(cubo_area)
scale = np.abs([data.min(), data.max()]).max()
mpl.contourf(y, x, data, shape, 30, cmap=mpl.cm.RdBu_r, vmin=-scale, vmax=scale)
mpl.colorbar(pad=0).set_label('nT')
mpl.m2km()
mpl.tight_layout(h_pad=0, w_pad=0)
widgets.interactive(reducao_polo,
inc=widgets.IntSliderWidget(min=-90, max=90, step=5, value=45),
dec=widgets.IntSliderWidget(min=-90, max=90, step=5, value=0),
erro=widgets.FloatSliderWidget(min=0, max=20, step=1, value=0))
inc
e dec
) e o erro colocado nos dados originais (da esquerda).