{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "**Course website**: http://www.leouieda.com/geofisica1\n", "\n", "**Note**: This notebook is part of the course \"Geofísica 1\" of Geology program of the \n", "[Universidade do Estado do Rio de Janeiro](http://www.uerj.br/). \n", "All content can be freely used and adapted under the terms of the \n", "[Creative Commons Attribution 4.0 International License](http://creativecommons.org/licenses/by/4.0/).\n", "\n", "![Creative Commons License](https://i.creativecommons.org/l/by/4.0/88x31.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Prática 4 - Deconvolução de Euler\n", "\n", "Objetivos:\n", "\n", "* Explorar as soluções obtidas com a Deconvolução de Euler.\n", "* Analisar os efeitos do ruído aleatório nas soluções obtidas.\n", "* Adquirir a intuição da sensibilidade da Deconvolução aos parâmetros que a controlam." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Preparação\n", "\n", "Esse documento que você está usando é um [IPython notebook](http://ipython.org/notebook.html). É 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](http://fatiando.org) de modelagem geofísica e também o módulo [numpy](http://www.numpy.org/).\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline\n", "from __future__ import division\n", "from IPython.html import widgets\n", "import numpy as np\n", "from fatiando.gravmag import prism, fourier, euler, sphere, polyprism\n", "from fatiando import utils, gridder, mesher\n", "from fatiando.vis import mpl, myv\n", "import fatiando\n", "mpl.rc('lines', linewidth=2)\n", "mpl.rc('font', size=12)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print('Versão do Fatiando a Terra: {}'.format(fatiando.__version__))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tarefa 1: Fonte simples (uma esfera)\n", "\n", "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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rode a célula abaixo para produzir a figura interativa." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "inc, dec = 45, 15\n", "modelo = [mesher.Sphere(1000, 1000, 400, 300, {'magnetization': utils.ang2vec(5, inc, dec)})]\n", "area = [0, 2000, 0, 2000]\n", "shape = (100, 100)\n", "x, y, z = gridder.regular(area, shape, z=0)\n", "anomalia_pura = sphere.tf(x, y, z, modelo, inc, dec)\n", "def fonte_simples(erro):\n", " n_janelas = 20 \n", " tamanho_janela = 600\n", " if erro > 0:\n", " anomalia = utils.contaminate(anomalia_pura, erro, seed=0)\n", " else:\n", " anomalia = anomalia_pura\n", " dx = fourier.derivx(x, y, anomalia, shape)\n", " dy = fourier.derivy(x, y, anomalia, shape)\n", " dz = fourier.derivz(x, y, anomalia, shape)\n", " prob = euler.Classic(x, y, z, anomalia, dx, dy, dz, 3)\n", " solver = euler.MovingWindow(prob, windows=(n_janelas, n_janelas), \n", " size=(tamanho_janela, tamanho_janela))\n", " solver.fit() \n", " filtrado = [e for e in solver.estimate_ if e[2] > 0]\n", " # Graficos\n", " mpl.figure(figsize=(12, 4.5))\n", " mpl.subplot(1, 2, 1)\n", " mpl.title(u'Anomalia com erro de {} nT'.format(erro))\n", " mpl.axis('scaled')\n", " mpl.contourf(y, x, anomalia, shape, 20, cmap=mpl.cm.RdBu_r)\n", " if filtrado:\n", " x0, y0, z0 = np.transpose(filtrado)\n", " mpl.scatter(y0, x0, c=z0, s=80, cmap=mpl.cm.RdYlGn, vmin=0, vmax=1000)\n", " mpl.colorbar(pad=0).set_label('Profundidade (metros)')\n", " mpl.grid()\n", " mpl.m2km()\n", " mpl.subplot(1, 2, 2)\n", " mpl.title('Derivada em z')\n", " mpl.axis('scaled')\n", " mpl.contourf(y, x, dz, shape, 20, cmap=mpl.cm.RdBu_r)\n", " mpl.colorbar(pad=0).set_label('nT/m')\n", " mpl.square((0, tamanho_janela, 0, tamanho_janela), '--k', linewidth=2, label='Janela')\n", " mpl.legend(loc='upper right')\n", " mpl.m2km()\n", " mpl.tight_layout(pad=0)\n", " return filtrado\n", "w = widgets.interactive(fonte_simples, \n", " erro=widgets.FloatSliderWidget(min=0, max=10, step=1, value=0))\n", "w" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sobre a figura\n", "\n", "* A figura da esquerda mostra a anomalia magnética de campo total gerada pela esfera.\n", "* O centro da esfera está a 400 metros de profundidade.\n", "* Os círculos coloridos representam as soluções da deconvolução de Euler aplicada a anomalia.\n", "* As cores representam a profundidade da solução (barra de cor).\n", "* A figura da direita mostra a derivada da anomalia na direção z (vertical) calculada com a Transformada de Fourier.\n", "* As derivadas em x, y e z são utilizadas na deconvolução. Somente a derivada em z é mostrada aqui para facilitar.\n", "* O quadrado preto tracejado representa a área da janela utilizada na deconvolução.\n", "* O índice estrutural utilizado na deconvolução é 3.\n", "* Somente 20% das soluções são mantidas (as melhores). Soluções localizadas acima dos dados também são descartadas." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rode a célula abaixo para **gerar uma figura 3D** com o modelo da esfera e as soluções da deconvolução de Euler.\n", "* Essa figura 3D reflete os resultados mostrados na figura acima. \n", "* Quando terminar, feche a janela da figura. Enquanto não fechar não poderá rodar mais nada aqui.\n", "* Rode a célula de novo gerar a figura novamente.\n", "* Se alterar a figura acima, a figura 3D será atualizada quando a gerar novamente.\n", "* Para salvar uma foto da figura clicando no disquete." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "myv.figure()\n", "for e in modelo:\n", " myv.points([[e.x, e.y, e.z]], color=(0, 0, 1), size=e.radius, opacity=0.3)\n", "myv.points(w.result, size=30, color=(1, 0, 0))\n", "bounds = area + [0, 1000]\n", "myv.axes(myv.outline(bounds), ranges=[b*0.001 for b in bounds], fmt='%.1f')\n", "myv.wall_bottom(bounds)\n", "myv.wall_north(bounds)\n", "myv.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Tarefas e perguntas\n", "\n", "1. Mantenha o erro em zero e gere a figura 3D.\n", " * Onde estão localizadas as soluções da deconvolução de Euler?\n", "2. Aumente o erro. Se lembre de olhar tanto a figura em planta (o mapa) quanto o 3D.\n", " * Como o erro influencia nas soluções da deconvolução?\n", " * Por que a deconvolução é tão sensível ao erro?\n", " * A distribuição das soluções em 3D possui alguma semelhança com a forma do corpo (esfera)? " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tarefa 2: Múltiplas fontes (duas esferas)\n", "\n", "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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rode a célula abaixo para produzir a figura interativa." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "inc, dec = 45, 15\n", "modelo = [mesher.Sphere(1500, 1500, 400, 300, {'magnetization': utils.ang2vec(5, inc, dec)}),\n", " mesher.Sphere(2500, 2500, 400, 300, {'magnetization': utils.ang2vec(5, -15, -20)})]\n", "area = [0, 4000, 0, 4000]\n", "shape = (100, 100)\n", "x, y, z = gridder.regular(area, shape, z=0)\n", "anomalia_pura = sphere.tf(x, y, z, modelo, inc, dec)\n", "anomalia = utils.contaminate(anomalia_pura, 5, seed=0)\n", "dx = fourier.derivx(x, y, anomalia, shape)\n", "dy = fourier.derivy(x, y, anomalia, shape)\n", "dz = fourier.derivz(x, y, anomalia, shape)\n", "prob = euler.Classic(x, y, z, anomalia, dx, dy, dz, 3)\n", "def duas_fontes(n_janelas, tamanho_janela):\n", " solver = euler.MovingWindow(prob, windows=(n_janelas, n_janelas), \n", " size=(tamanho_janela, tamanho_janela),\n", " keep=0.15)\n", " solver.fit() \n", " filtrado = [e for e in solver.estimate_ if e[2] > 0]\n", " mpl.figure(figsize=(7, 6))\n", " mpl.title(u'Deconvolução com {} x {} janelas de {} m'.format(n_janelas, n_janelas, tamanho_janela))\n", " mpl.axis('scaled')\n", " mpl.contourf(y, x, anomalia, shape, 20, cmap=mpl.cm.RdBu_r)\n", " if filtrado:\n", " x0, y0, z0 = np.transpose(filtrado)\n", " mpl.scatter(y0, x0, c=z0, s=80, cmap=mpl.cm.RdYlGn, vmin=0, vmax=1000)\n", " mpl.colorbar(pad=0).set_label('Profundidade (metros)')\n", " mpl.square((0, tamanho_janela, 0, tamanho_janela), '--k', linewidth=2, label='Janela')\n", " mpl.legend(loc='upper right')\n", " mpl.m2km()\n", " mpl.tight_layout(pad=0)\n", " return filtrado\n", "w = widgets.interactive(duas_fontes, \n", " n_janelas=widgets.IntSliderWidget(min=5, max=50, step=5, value=30),\n", " tamanho_janela=widgets.FloatSliderWidget(min=200, max=3000, step=100, value=1000))\n", "w" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sobre a figura\n", "\n", "* A anomalia é gerada por duas esferas com magnetizações distintas.\n", "* Os dados foram contaminados com erro de 5 nT.\n", "* O quadrado preto tracejado representa a área da janela utilizada na deconvolução.\n", "* O índice estrutural utilizado na deconvolução é 3.\n", "* Somente 15% das soluções são mantidas (as melhores). Soluções localizadas acima dos dados também são descartadas.\n", "* O botão `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. \n", "* O botã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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rode a célula abaixo para **gerar uma figura 3D** com o modelo da esfera e as soluções da deconvolução de Euler.\n", "* Essa figura 3D reflete os resultados mostrados na figura acima. \n", "* Quando terminar, feche a janela da figura. Enquanto não fechar não poderá rodar mais nada aqui.\n", "* Rode a célula de novo gerar a figura novamente.\n", "* Se alterar a figura acima, a figura 3D será atualizada quando a gerar novamente.\n", "* Para salvar uma foto da figura clicando no disquete." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "myv.figure()\n", "for e in modelo:\n", " myv.points([[e.x, e.y, e.z]], color=(0, 0, 1), size=e.radius, opacity=0.3)\n", "myv.points(w.result, size=30, color=(1, 0, 0))\n", "bounds = area + [0, 1000]\n", "myv.axes(myv.outline(bounds), ranges=[b*0.001 for b in bounds], fmt='%.1f')\n", "myv.wall_bottom(bounds)\n", "myv.wall_north(bounds)\n", "myv.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Tarefas e perguntas\n", "\n", "1. Mantenha o tamanho das janelas em 1000 metros. Mude o número de janelas. Não se esqueça de olhar a figura 3D.\n", " * O número de soluções muda?\n", " * Como varia a dispersão das soluções obtidas?\n", "2. Mantenha o número de janelas fixo em 30. Mude o tamanho das janelas.\n", " * O que acontece quando o tamanho das janelas é pequeno? Por que?\n", " * E quando o tamanho é grande? Por que?\n", " * Qual é o tamanho que produz as soluções mais concentradas em cada uma das 2 anomalias? \n", " * Qual é a relação tamanho com as anomalias?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tarefa 3: Fontes não-ideais e o índice estrutural\n", "\n", "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.\n", "\n", "Esse modelo é o que foi utilizado em:\n", "[Uieda, Oliveira Jr., and Barbosa (2014), Geophysical tutorial: Euler deconvolution of potential-field data, The Leading Edge, 33(4), 448-450.](http://wiki.seg.org/wiki/Euler_deconvolution_of_potential_field_data_%28tutorial%29).\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rode a célula abaixo para produzir a figura interativa." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "area = [5000, 25000, 5000, 25000]\n", "bounds = [0, 30000, 0, 30000, 0, 5000]\n", "shape = (100, 100)\n", "x, y, z = gridder.regular(area, shape, z=-300)\n", "inc, dec = -15, 30\n", "dike = mesher.PolygonalPrism(\n", " [[78.12500000000136, 17968.750000000004],\n", " [15156.250000000002, 19843.750000000004],\n", " [29843.75, 20781.250000000004],\n", " [29843.75, 21015.625000000004],\n", " [15078.125000000002, 20156.250000000004],\n", " [156.25000000000136, 18281.250000000004]], 0, 5000,\n", " props={'magnetization': utils.ang2vec(10, inc, dec)})\n", "sill = mesher.PolygonalPrism(\n", " [[10000.000000000002, 8828.125000000004],\n", " [10156.250000000002, 9062.500000000004],\n", " [11328.125000000002, 10703.125000000004],\n", " [11875.000000000002, 11796.875000000004],\n", " [11953.125000000002, 12890.625000000004],\n", " [11406.250000000002, 13593.750000000004],\n", " [10156.250000000002, 13984.375000000004],\n", " [9375.000000000002, 13984.375000000004],\n", " [9218.750000000002, 13984.375000000004],\n", " [8437.500000000002, 12812.500000000004],\n", " [7968.750000000002, 11796.875000000004],\n", " [7968.750000000002, 10078.125000000004],\n", " [7968.750000000002, 9765.625000000004],\n", " [8593.750000000002, 9218.750000000004]], 1000, 1500,\n", " props={'magnetization': utils.ang2vec(10, inc, dec)})\n", "batolito = mesher.PolygonalPrism(\n", " [[18046.875, 9921.875000000004],\n", " [19062.5, 10234.375000000004],\n", " [19765.625, 11093.750000000004],\n", " [20078.125, 12500.000000000004],\n", " [19843.75, 13515.625000000004],\n", " [18281.25, 14140.625000000004],\n", " [16640.625, 14140.625000000004],\n", " [15781.250000000002, 13671.875000000004],\n", " [15703.125000000002, 11718.750000000004],\n", " [15859.375000000002, 10390.625000000004],\n", " [16250.000000000002, 9843.750000000004]], 500, 4000,\n", " props={'magnetization': utils.ang2vec(2, inc, dec)})\n", "modelo = [dike, sill, batolito]\n", "anomalia = utils.contaminate(polyprism.tf(x, y, z, modelo, inc, dec), 5)\n", "dx = fourier.derivx(x, y, anomalia, shape)\n", "dy = fourier.derivy(x, y, anomalia, shape)\n", "dz = fourier.derivz(x, y, anomalia, shape)\n", "def complexo(indice):\n", " tamanho_janela = 3000\n", " prob = euler.Classic(x, y, z, anomalia, dx, dy, dz, indice)\n", " n_janelas = 50\n", " solver = euler.MovingWindow(prob, windows=(n_janelas, n_janelas), \n", " size=(tamanho_janela, tamanho_janela),\n", " keep=0.3)\n", " solver.fit() \n", " filtrado = [e for e in solver.estimate_ if e[2] > 0]\n", " mpl.figure(figsize=(7, 6))\n", " mpl.title(u'Índice estrutural {}'.format(indice))\n", " mpl.axis('scaled')\n", " mpl.contourf(y, x, anomalia, shape, 20, cmap=mpl.cm.RdBu_r)\n", " if filtrado:\n", " x0, y0, z0 = np.transpose(filtrado)\n", " mpl.scatter(y0, x0, c=z0, s=80, cmap=mpl.cm.RdYlGn, vmin=0, vmax=3000)\n", " mpl.colorbar(pad=0).set_label('Profundidade (metros)')\n", " mpl.square((5e3, 5e3 + tamanho_janela, 5e3, 5e3 + tamanho_janela), '--k', linewidth=2, label='Janela')\n", " mpl.legend(loc='upper left')\n", " for b in modelo:\n", " mpl.polygon(b, xy2ne=True)\n", " mpl.m2km()\n", " mpl.tight_layout(pad=0)\n", " return filtrado\n", "w = widgets.interactive(complexo, \n", " indice=widgets.IntSliderWidget(min=1, max=4, step=1, value=3))\n", "w" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sobre a figura\n", "\n", "* A anomalia foi gerada pelos 3 corpos e contaminada com erro de 5 nT.\n", "* As linhas pretas sólidas representam a projeção dos 3 corpos na superfície.\n", "* Nenhum dos corpos aflora.\n", "* Somente 30% das soluções são mantidas (as melhores). Soluções localizadas acima dos dados também são descartadas.\n", "* O número de janelas é de 50 x 50 e o tamanho de 3000 x 3000 metros.\n", "* O botão `indice` controla o índice estrutural." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Índice estrutural\n", "\n", "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:\n", "\n", "| Fonte | IS |\n", "|-------------------|-----|\n", "| Dique | 1 |\n", "| Sill | 1 |\n", "| Esfera | 3 |\n", "| Cilindro vertical | 2 |\n", "\n", "**Referências**\n", "\n", "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](http://dx.doi.org/doi:10.1190/1.1442774)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rode a célula abaixo para **gerar uma figura 3D** com o modelo da esfera e as soluções da deconvolução de Euler.\n", "* Essa figura 3D reflete os resultados mostrados na figura acima. \n", "* Quando terminar, feche a janela da figura. Enquanto não fechar não poderá rodar mais nada aqui.\n", "* Rode a célula de novo gerar a figura novamente.\n", "* Se alterar a figura acima, a figura 3D será atualizada quando a gerar novamente.\n", "* Para salvar uma foto da figura clicando no disquete." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "myv.figure()\n", "myv.polyprisms(modelo, opacity=0.3)\n", "myv.points(w.result, size=150, color=(1, 0, 0))\n", "myv.axes(myv.outline(bounds), ranges=[b*0.001 for b in bounds], fmt='%.1f')\n", "myv.wall_bottom(bounds)\n", "myv.wall_north(bounds)\n", "myv.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Tarefas e perguntas\n", "\n", "1. Mantenha o índice em 3 (esfera). Gere a figura 3D.\n", " * Qual dos 3 corpos mais se aproxima de uma esfera?\n", " * As soluções da deconvolução estão localizadas no centro desse corpo, como acontecia com a esfera da Tarefa 1?\n", " * As soluções tem alguma semelhança com a forma do corpo (todo o volume, não só as bordas)?\n", " * Como são as soluções para os outros 2 corpos?\n", " * Elas acertam as profundidades desses 2 corpos?\n", "2. Diminua o índice estrutural. Lembre de olhar a figura 3D.\n", " * Como variam as soluções (localização horizontal, profundidade, dispersão, etc)?\n", " * Qual seria o índice estrutural mais correto para o dique e para o sill?\n", " * Onde estão localizadas as soluções da deconvolução para esse índice?" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.9" } }, "nbformat": 4, "nbformat_minor": 0 }