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, euler, sphere, polyprism
from fatiando import utils, gridder, mesher
from fatiando.vis import mpl, myv
import fatiando
mpl.rc('lines', linewidth=2)
mpl.rc('font', size=12)
print('Versão do Fatiando a Terra: {}'.format(fatiando.__version__))
Nessa tarefa vamos analisar os efeitos do erro aleatório nos resultados da deconvolução. Para isso vamos usar um modelo simples com uma única esfera.
Rode a célula abaixo para produzir a figura interativa.
inc, dec = 45, 15
modelo = [mesher.Sphere(1000, 1000, 400, 300, {'magnetization': utils.ang2vec(5, inc, dec)})]
area = [0, 2000, 0, 2000]
shape = (100, 100)
x, y, z = gridder.regular(area, shape, z=0)
anomalia_pura = sphere.tf(x, y, z, modelo, inc, dec)
def fonte_simples(erro):
n_janelas = 20
tamanho_janela = 600
if erro > 0:
anomalia = utils.contaminate(anomalia_pura, erro, seed=0)
else:
anomalia = anomalia_pura
dx = fourier.derivx(x, y, anomalia, shape)
dy = fourier.derivy(x, y, anomalia, shape)
dz = fourier.derivz(x, y, anomalia, shape)
prob = euler.Classic(x, y, z, anomalia, dx, dy, dz, 3)
solver = euler.MovingWindow(prob, windows=(n_janelas, n_janelas),
size=(tamanho_janela, tamanho_janela))
solver.fit()
filtrado = [e for e in solver.estimate_ if e[2] > 0]
# Graficos
mpl.figure(figsize=(12, 4.5))
mpl.subplot(1, 2, 1)
mpl.title(u'Anomalia com erro de {} nT'.format(erro))
mpl.axis('scaled')
mpl.contourf(y, x, anomalia, shape, 20, cmap=mpl.cm.RdBu_r)
if filtrado:
x0, y0, z0 = np.transpose(filtrado)
mpl.scatter(y0, x0, c=z0, s=80, cmap=mpl.cm.RdYlGn, vmin=0, vmax=1000)
mpl.colorbar(pad=0).set_label('Profundidade (metros)')
mpl.grid()
mpl.m2km()
mpl.subplot(1, 2, 2)
mpl.title('Derivada em z')
mpl.axis('scaled')
mpl.contourf(y, x, dz, shape, 20, cmap=mpl.cm.RdBu_r)
mpl.colorbar(pad=0).set_label('nT/m')
mpl.square((0, tamanho_janela, 0, tamanho_janela), '--k', linewidth=2, label='Janela')
mpl.legend(loc='upper right')
mpl.m2km()
mpl.tight_layout(pad=0)
return filtrado
w = widgets.interactive(fonte_simples,
erro=widgets.FloatSliderWidget(min=0, max=10, step=1, value=0))
w
Rode a célula abaixo para gerar uma figura 3D com o modelo da esfera e as soluções da deconvolução de Euler.
myv.figure()
for e in modelo:
myv.points([[e.x, e.y, e.z]], color=(0, 0, 1), size=e.radius, opacity=0.3)
myv.points(w.result, size=30, color=(1, 0, 0))
bounds = area + [0, 1000]
myv.axes(myv.outline(bounds), ranges=[b*0.001 for b in bounds], fmt='%.1f')
myv.wall_bottom(bounds)
myv.wall_north(bounds)
myv.show()
Agora que vocês pegaram o jeito da deconvolução, vamos ver como ela se comporta para um modelo com duas esferas. Vamos analisar como as soluções são influenciadas pelo procedimento de dividir a área em janelas. Os parâmetros que controlam esse procedimento são o número de janelas e o seu tamanho.
Rode a célula abaixo para produzir a figura interativa.
inc, dec = 45, 15
modelo = [mesher.Sphere(1500, 1500, 400, 300, {'magnetization': utils.ang2vec(5, inc, dec)}),
mesher.Sphere(2500, 2500, 400, 300, {'magnetization': utils.ang2vec(5, -15, -20)})]
area = [0, 4000, 0, 4000]
shape = (100, 100)
x, y, z = gridder.regular(area, shape, z=0)
anomalia_pura = sphere.tf(x, y, z, modelo, inc, dec)
anomalia = utils.contaminate(anomalia_pura, 5, seed=0)
dx = fourier.derivx(x, y, anomalia, shape)
dy = fourier.derivy(x, y, anomalia, shape)
dz = fourier.derivz(x, y, anomalia, shape)
prob = euler.Classic(x, y, z, anomalia, dx, dy, dz, 3)
def duas_fontes(n_janelas, tamanho_janela):
solver = euler.MovingWindow(prob, windows=(n_janelas, n_janelas),
size=(tamanho_janela, tamanho_janela),
keep=0.15)
solver.fit()
filtrado = [e for e in solver.estimate_ if e[2] > 0]
mpl.figure(figsize=(7, 6))
mpl.title(u'Deconvolução com {} x {} janelas de {} m'.format(n_janelas, n_janelas, tamanho_janela))
mpl.axis('scaled')
mpl.contourf(y, x, anomalia, shape, 20, cmap=mpl.cm.RdBu_r)
if filtrado:
x0, y0, z0 = np.transpose(filtrado)
mpl.scatter(y0, x0, c=z0, s=80, cmap=mpl.cm.RdYlGn, vmin=0, vmax=1000)
mpl.colorbar(pad=0).set_label('Profundidade (metros)')
mpl.square((0, tamanho_janela, 0, tamanho_janela), '--k', linewidth=2, label='Janela')
mpl.legend(loc='upper right')
mpl.m2km()
mpl.tight_layout(pad=0)
return filtrado
w = widgets.interactive(duas_fontes,
n_janelas=widgets.IntSliderWidget(min=5, max=50, step=5, value=30),
tamanho_janela=widgets.FloatSliderWidget(min=200, max=3000, step=100, value=1000))
w
n_janelas
controla o número de janelas em ambas as direções. n_janelas = 10
significa que serão utilizadas 10 x 10 = 100 janelas para fazer a deconvolução.tamanho_janela
controla o tamanho das janelas (em metros). n_janelas = 10
e tamanho_janela = 1000
significa que a área será coberta por 10 x 10 = 100 janelas com 1000 x 1000 metros cada. Note que as janelas se sobrepõe.Rode a célula abaixo para gerar uma figura 3D com o modelo da esfera e as soluções da deconvolução de Euler.
myv.figure()
for e in modelo:
myv.points([[e.x, e.y, e.z]], color=(0, 0, 1), size=e.radius, opacity=0.3)
myv.points(w.result, size=30, color=(1, 0, 0))
bounds = area + [0, 1000]
myv.axes(myv.outline(bounds), ranges=[b*0.001 for b in bounds], fmt='%.1f')
myv.wall_bottom(bounds)
myv.wall_north(bounds)
myv.show()
Para complicar ainda mais, vamos fazer um teste com um modelo mais complexo: um sill, um batólito e um dique. Nesse exercício, vamos explorar a influência que o índice estrutural tem sobre as soluções.
Esse modelo é o que foi utilizado em: Uieda, Oliveira Jr., and Barbosa (2014), Geophysical tutorial: Euler deconvolution of potential-field data, The Leading Edge, 33(4), 448-450..
Rode a célula abaixo para produzir a figura interativa.
area = [5000, 25000, 5000, 25000]
bounds = [0, 30000, 0, 30000, 0, 5000]
shape = (100, 100)
x, y, z = gridder.regular(area, shape, z=-300)
inc, dec = -15, 30
dike = mesher.PolygonalPrism(
[[78.12500000000136, 17968.750000000004],
[15156.250000000002, 19843.750000000004],
[29843.75, 20781.250000000004],
[29843.75, 21015.625000000004],
[15078.125000000002, 20156.250000000004],
[156.25000000000136, 18281.250000000004]], 0, 5000,
props={'magnetization': utils.ang2vec(10, inc, dec)})
sill = mesher.PolygonalPrism(
[[10000.000000000002, 8828.125000000004],
[10156.250000000002, 9062.500000000004],
[11328.125000000002, 10703.125000000004],
[11875.000000000002, 11796.875000000004],
[11953.125000000002, 12890.625000000004],
[11406.250000000002, 13593.750000000004],
[10156.250000000002, 13984.375000000004],
[9375.000000000002, 13984.375000000004],
[9218.750000000002, 13984.375000000004],
[8437.500000000002, 12812.500000000004],
[7968.750000000002, 11796.875000000004],
[7968.750000000002, 10078.125000000004],
[7968.750000000002, 9765.625000000004],
[8593.750000000002, 9218.750000000004]], 1000, 1500,
props={'magnetization': utils.ang2vec(10, inc, dec)})
batolito = mesher.PolygonalPrism(
[[18046.875, 9921.875000000004],
[19062.5, 10234.375000000004],
[19765.625, 11093.750000000004],
[20078.125, 12500.000000000004],
[19843.75, 13515.625000000004],
[18281.25, 14140.625000000004],
[16640.625, 14140.625000000004],
[15781.250000000002, 13671.875000000004],
[15703.125000000002, 11718.750000000004],
[15859.375000000002, 10390.625000000004],
[16250.000000000002, 9843.750000000004]], 500, 4000,
props={'magnetization': utils.ang2vec(2, inc, dec)})
modelo = [dike, sill, batolito]
anomalia = utils.contaminate(polyprism.tf(x, y, z, modelo, inc, dec), 5)
dx = fourier.derivx(x, y, anomalia, shape)
dy = fourier.derivy(x, y, anomalia, shape)
dz = fourier.derivz(x, y, anomalia, shape)
def complexo(indice):
tamanho_janela = 3000
prob = euler.Classic(x, y, z, anomalia, dx, dy, dz, indice)
n_janelas = 50
solver = euler.MovingWindow(prob, windows=(n_janelas, n_janelas),
size=(tamanho_janela, tamanho_janela),
keep=0.3)
solver.fit()
filtrado = [e for e in solver.estimate_ if e[2] > 0]
mpl.figure(figsize=(7, 6))
mpl.title(u'Índice estrutural {}'.format(indice))
mpl.axis('scaled')
mpl.contourf(y, x, anomalia, shape, 20, cmap=mpl.cm.RdBu_r)
if filtrado:
x0, y0, z0 = np.transpose(filtrado)
mpl.scatter(y0, x0, c=z0, s=80, cmap=mpl.cm.RdYlGn, vmin=0, vmax=3000)
mpl.colorbar(pad=0).set_label('Profundidade (metros)')
mpl.square((5e3, 5e3 + tamanho_janela, 5e3, 5e3 + tamanho_janela), '--k', linewidth=2, label='Janela')
mpl.legend(loc='upper left')
for b in modelo:
mpl.polygon(b, xy2ne=True)
mpl.m2km()
mpl.tight_layout(pad=0)
return filtrado
w = widgets.interactive(complexo,
indice=widgets.IntSliderWidget(min=1, max=4, step=1, value=3))
w
indice
controla o índice estrutural.O índice estrutural (IS) representa o decaimento do campo magnético (ou gravitacional) com a distância. Esse fator só é definido para fontes ideias, como esferas, diques retos e verticais, etc, e varia para cada uma dessas fontes. Segundo Reid et al. (1990), para o caso magnético:
Fonte | IS |
---|---|
Dique | 1 |
Sill | 1 |
Esfera | 3 |
Cilindro vertical | 2 |
Referências
Reid, A. B., J. M. Allsop, H. Granser, A. J. Millett, and I. W. Somerton (1990), Magnetic interpretation in three dimensions using Euler deconvolution, Geophysics, 55(1), 80–91, doi:10.1190/1.1442774.
Rode a célula abaixo para gerar uma figura 3D com o modelo da esfera e as soluções da deconvolução de Euler.
myv.figure()
myv.polyprisms(modelo, opacity=0.3)
myv.points(w.result, size=150, color=(1, 0, 0))
myv.axes(myv.outline(bounds), ranges=[b*0.001 for b in bounds], fmt='%.1f')
myv.wall_bottom(bounds)
myv.wall_north(bounds)
myv.show()