{ "metadata": { "name": "", "signature": "sha256:f893cdc570cc8bb6c666cf61f7896a24c6af1286fad151d276bc17520a93ca1c" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "**Course website**: http://www.leouieda.com/geofisica1\n", "\n", "**Note**: This notebook is part of the course \"Geof\u00edsica 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\u00e1tica 6: Invers\u00e3o para estimar o embasamento de uma bacia" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Prepara\u00e7\u00e3o\n", "\n", "Como na [pr\u00e1tica 5](http://www.leouieda.com/geofisica1/lessons/gravimetria/pratica5.html), esse documento que voc\u00ea est\u00e1 usando \u00e9 um [IPython notebook](http://ipython.org/notebook.html). Esta pr\u00e1tica usar\u00e1 a parte de invers\u00e3o 2D da biblioteca [Fatiando a Terra](http://fatiando.org) de modelagem geof\u00edsica. Para isso, precisamos carregar (`import`) os m\u00f3dulos que vamos usar e tamb\u00e9m o m\u00f3dulo [numpy](http://www.numpy.org/).\n", "\n", "O notebook \u00e9 divido em c\u00e9lulas (como esta). Para editar o conte\u00fado de uma c\u00e9lula, clique nela (clique nesta para editar esse texto). Para executar uma c\u00e9lula, aperte `Shift + Enter`. Execute as duas c\u00e9lulas abaixo." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from IPython.core.pylabtools import print_figure\n", "from IPython.html import widgets\n", "from IPython.display import Image, display\n", "from fatiando.gravmag import talwani\n", "from fatiando.gravmag.basin2d import PolygonalBasinGravity\n", "from fatiando.gravmag.interactive import Moulder\n", "from fatiando.inversion.regularization import Damping, Smoothness1D, Smoothness2D\n", "from fatiando.vis import mpl\n", "from fatiando import utils\n", "import fatiando\n", "import numpy as np" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "print('Vers\u00e3o do Fatiando a Terra: {}'.format(fatiando.__version__))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tarefa 1: Gere um modelo de uma bacia\n", "\n", "Fa\u00e7a um modelo de uma bacia sedimentar e passe para o colega. Algumas instru\u00e7\u00f5es:\n", "* O topo da bacia deve ser reto e estar pr\u00f3ximo da superf\u00edcie (z = 0).\n", "* Desenhe a bacia no **sentido hor\u00e1rio**.\n", "* **N\u00e3o coloque erro no dado**. Faremos isso depois.\n", "* Sua bacia deve conter tanto varia\u00e7\u00f5es abruptas (falhas) quanto varia\u00e7\u00f5es suaves.\n", "* Fa\u00e7a uma bacia **realista**. Use a geologia que foi entulhada em sua cabe\u00e7a at\u00e9 agora." ] }, { "cell_type": "code", "collapsed": false, "input": [ "area = [0, 100e3, 0, 10e3]\n", "x = np.linspace(area[0], area[1], 80)\n", "z = -np.ones_like(x)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "meu_modelo = Moulder(area, x, z, density_range=[-1000, 0])" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "meu_modelo.run()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rode a c\u00e9lula abaixo para salvar o seu modelo. **Mude o nome do arquivo** e troque com um colega. " ] }, { "cell_type": "code", "collapsed": false, "input": [ "meu_modelo.save('modelo_bacia')" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tarefa 2: Fa\u00e7a a invers\u00e3o sem utilizar v\u00ednculos" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Primeiro precisamos pegar o modelo do colega e extrair os dados (sem ru\u00eddo), o pol\u00edgono que descreve a bacia e a densidade correta. \n", "\n", "**Mude o nome do arquivo** `'modelo_bacia'` abaixo para o nome do arquivo do colega." ] }, { "cell_type": "code", "collapsed": false, "input": [ "modelo = Moulder.load('modelo_bacia')\n", "modelo.plot(figsize=(12, 6))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "dado = modelo.predicted\n", "bacia = modelo.model[0]\n", "densidade = bacia.props.copy()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A c\u00e9lula abaixo prepara a invers\u00e3o. `pontos` define quantos pontos ter\u00e1 o pol\u00edgono que descrever\u00e1 a bacia. Esse ser\u00e1 o n\u00famero de **par\u00e2metros** que a invers\u00e3o estimar\u00e1 (a profundidade de cada ponto do pol\u00edgono).\n", "\n", "Lembrando que a invers\u00e3o busca achar o **m\u00ednimo da fun\u00e7\u00e3o do ajuste**:\n", "\n", "$$\\phi(\\bar{p}) = \\sum\\limits_{i=1}^N (d^o_i - d_i)^2$$\n", "\n", "em que $\\bar{p}$ \u00e9 um vetor com os par\u00e2metros (profundidade de cada ponto do pol\u00edgono), $d^o_i$ \u00e9 o i-\u00e9simo dado observado e $d_i$ \u00e9 o dado predito pelos par\u00e2metros.\n", "\n", "Essa invers\u00e3o \u00e9 do tipo **n\u00e3o-linear**. Isso quer dizer que precisamos espeficicar um chute inicial para os par\u00e2metros. A invers\u00e3o \u00e9 feita em passos a partir desse chute inicial. Vamos usar \"todos os par\u00e2metros igual a 1 km\"." ] }, { "cell_type": "code", "collapsed": false, "input": [ "pontos = 40\n", "topo = bacia.y.min()\n", "p0 = 1000*np.ones(pontos)\n", "xlim = [bacia.x.min(), bacia.x.max()]\n", "misfit = PolygonalBasinGravity(x, z, dado, pontos, bacia.props, top=topo, xlim=xlim).config('levmarq', initial=p0)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A c\u00e9lula abaixo cria algumas fun\u00e7\u00f5es que usaremos para fazer um gr\u00e1fico da nossa solu\u00e7\u00e3o." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def plota_resultado(solver):\n", " fig = mpl.figure(figsize=(12, 4), facecolor='white')\n", " ax = mpl.subplot(2, 1, 1)\n", " mpl.plot(x, solver.data, 'ok')\n", " mpl.plot(x, solver.predicted(), '-r', linewidth=2)\n", " ax.set_ylim(1.2*dado.min(), 1.1*dado.max())\n", " ax.grid(True)\n", " ax.set_xticklabels([])\n", " ax.set_ylabel('Anomalia (mGal)')\n", " ax = mpl.subplot(2, 1, 2)\n", " mpl.polygon(bacia, style='o-k', fill='gray', alpha=0.5, linewidth=2)\n", " mpl.polygon(solver.estimate_, style='o-r', linewidth=2)\n", " ax.set_xlabel('x (km)')\n", " ax.set_ylabel('z (km)')\n", " mpl.set_area(area)\n", " mpl.m2km()\n", " ax.grid(True)\n", " ax.invert_yaxis()\n", " mpl.tight_layout(pad=0)\n", " return fig" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rode a c\u00e9lula abaixo para criar um **gr\u00e1fico interativo**. A parte de cima do gr\u00e1fico mostra os dados observados (pontos) e preditos pela invers\u00e3o (linha vermelha). A parte de baixo mostra o seu modelo de bacia (preto) e o resultado da invers\u00e3o (vermelho).\n", "\n", "O bot\u00e3o \"erro\" pode ser movido para mudar o valor do erro que ser\u00e1 aplicado aos dados. Cada vez que voc\u00ea mover o bot\u00e3o, o dado \u00e9 alterado e invers\u00e3o rodada novamente.\n", "\n", "O bot\u00e3o \"seed\" altera ligeiramente o erro que \u00e9 aplicado ao dado (mas n\u00e3o o tamanho do erro). \n", "\n", "### Tarefas e perguntas\n", "\n", "1. Execute a c\u00e9lula abaixo para rodar a invers\u00e3o e criar o gr\u00e1fico interativo.\n", "2. Aplique erro no dado. Se der uma mensagem de erro, n\u00e3o se preocupe. Quando isso acontece significa que o problema **n\u00e3o possui solu\u00e7\u00e3o \u00fanica**. Diminua o erro para voltar a obter uma solu\u00e7\u00e3o.\n", " * O que acontece com a estimativa da invers\u00e3o?\n", " * A invers\u00e3o \u00e9 capaz de recuperar o modelo original que gerou os dados?\n", " * A invers\u00e3o recupera o modelo original se o erro for zero? Por que?\n", "3. Use o bot\u00e3o \"seed\" para mudar ligeiramente o erro aplicado.\n", " * O que acontece com a estimativa?\n", " * Essa invers\u00e3o \u00e9 est\u00e1vel?" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def interativo_sem_vinculo(erro, seed):\n", " misfit.data = utils.contaminate(dado, erro, seed=seed)\n", " misfit.fit()\n", " fig = plota_resultado(misfit)\n", " mpl.close(fig) \n", " display(Image(print_figure(fig, dpi=80)))\n", " misfit.data = dado\n", " \n", "widgets.interact(interativo_sem_vinculo, \n", " erro=widgets.FloatSliderWidget(min=0, max=1, step=0.1, value=0), \n", " seed=widgets.IntSliderWidget(min=0, max=10, step=1, value=0))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tarefa 3: Utilizando o v\u00ednculo \"damping\" (ou norma m\u00ednima)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Vamos utilizar o v\u00ednculo de \"norma m\u00ednima\" ou [\"damping\"](http://en.wikipedia.org/wiki/Tikhonov_regularization). \n", "\n", "Lembre que adicionar v\u00ednculos \u00e9 equivalente a somar uma outra fun\u00e7\u00e3o na nossa fun\u00e7\u00e3o do ajuste $\\phi(\\bar{p})$:\n", "\n", "$$\\Gamma(\\bar{p}) = \\sum\\limits_{i=1}^N (d^o_i - d_i)^2 + \\mu \\sum\\limits_{j=1}^M p_j^2$$\n", "\n", "A fun\u00e7\u00e3o do v\u00ednculo \u00e9 $\\sum\\limits_{j=1}^M p_j^2$ que imp\u00f5e que os par\u00e2metros estimados sejam o menor poss\u00edvel. Em outras palavras, eles devem ter **norma m\u00ednima**.\n", "\n", "A constante positiva $\\mu$ determina qual \u00e9 o peso do v\u00ednculo. $\\mu$ cuida do balan\u00e7o entre ajustar os dados e obedecer aos v\u00ednculos." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Tarefas e perguntas\n", "\n", "1. Execute a c\u00e9lula abaixo para rodar a invers\u00e3o e criar o gr\u00e1fico interativo.\n", "2. Aplique erro ao dado.\n", "3. O bot\u00e3o `mi` controla o valor de $\\mu$ (um valor de `-10` significa que $\\mu = 10^{-15}$). Use o bot\u00e3o para aumentar a import\u00e2ncia do v\u00ednculo de norma m\u00ednima.\n", " * O que acontece se $\\mu$ for muito pequeno? O que acontece com o ajuste? Com $\\mu$ pequeno o problema \u00e9 est\u00e1vel?\n", " * O que acontece se $\\mu$ for muito grande? O que acontece com o ajuste? Com $\\mu$ grande o problema \u00e9 est\u00e1vel?\n", " * Por que a bacia estimada tende para esse valor quando $\\mu$ \u00e9 grande?\n", "4. Determine um valor \"ideal\" para $\\mu$. Por ideal quero dizer: \u00e9 o menor poss\u00edvel para tornar o problema **est\u00e1vel e ajustar os dados**.\n", " * A invers\u00e3o \u00e9 capaz de recuperar a falha? Por que?\n", " * E a regi\u00e3o suave?" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def interativo_damping(erro, seed, mi):\n", " misfit.data = utils.contaminate(dado, erro, seed=seed)\n", " solver = misfit + 10**mi*Damping(misfit.nparams)\n", " solver.fit()\n", " fig = plota_resultado(solver)\n", " mpl.close(fig) \n", " display(Image(print_figure(fig, dpi=80)))\n", " misfit.data = dado\n", " \n", "widgets.interact(interativo_damping, \n", " erro=widgets.FloatSliderWidget(min=0, max=1, step=0.1, value=0), \n", " seed=widgets.IntSliderWidget(min=0, max=10, step=1, value=0),\n", " mi=widgets.FloatSliderWidget(min=-15, max=0, step=0.5, value=-15))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tarefa 4: Estabilizando com o v\u00ednculo de suavidade" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Agora vamos utilizar o v\u00ednculo de suavidade.\n", "\n", "Nesse caso, o v\u00ednculo \u00e9:\n", "\n", "$$\\sum\\limits_{j=1}^{M - 1} (p_{j+1} - p_j)^2$$\n", "\n", "O v\u00ednculo de suavidade imp\u00f5e que os par\u00e2metros sejam o mais pr\u00f3ximos o poss\u00edvel de seus vizinhos imediatos." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Tarefas e perguntas\n", "\n", "1. Execute a c\u00e9lula abaixo para rodar a invers\u00e3o e criar o gr\u00e1fico interativo.\n", "2. Aplique erro ao dado.\n", "3. Use o bot\u00e3o `mi` para aumentar a import\u00e2ncia do v\u00ednculo. Se der uma mensagem de erro, aumente `mi`.\n", " * O que acontece se $\\mu$ for muito pequeno? O que acontece com o ajuste? Com $\\mu$ pequeno o problema \u00e9 est\u00e1vel?\n", " * O que acontece se $\\mu$ for muito grande? O que acontece com o ajuste? Com $\\mu$ grande o problema \u00e9 est\u00e1vel?\n", " * Por que a bacia estimada tende para esse valor quando $\\mu$ \u00e9 grande?\n", "4. Determine um valor \"ideal\" para $\\mu$ (menor poss\u00edvel para tornar o problema **est\u00e1vel e ajustar os dados**).\n", " * A invers\u00e3o \u00e9 capaz de recuperar a falha? Por que?\n", " * E a regi\u00e3o suave?\n", "5. Aumente o erro aplicado ao dado verifique se o problema continua est\u00e1vel. (Se der erro, n\u00e3o se preocupe e aumente `mi`)\n", " * Se o problema deixar de ser est\u00e1vel, o que deve ser feito?\n", " * Como o erro influencia a qualidade da solu\u00e7\u00e3o?" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def interativo_suave(erro, seed, mu):\n", " misfit.data = utils.contaminate(dado, erro, seed=seed)\n", " solver = misfit + 10**mu*Smoothness1D(misfit.nparams)\n", " solver.fit()\n", " fig = plota_resultado(solver)\n", " mpl.close(fig) \n", " display(Image(print_figure(fig, dpi=80)))\n", " misfit.data = dado\n", "\n", "widgets.interact(interativo_suave, \n", " erro=widgets.FloatSliderWidget(min=0, max=5, step=0.1, value=0), \n", " seed=widgets.IntSliderWidget(min=0, max=10, step=1, value=0),\n", " mu=widgets.FloatSliderWidget(min=-15, max=0, step=0.5, value=-15))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tarefa 5: E se voc\u00ea n\u00e3o souber a densidade?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Em casos reais n\u00e3o saberemos exatamente qual \u00e9 a densidade dos sedimentos. Podemos ter uma ideia mas valores de propriedades f\u00edsicas costumam variar bastante.\n", "\n", "Vamos assumir que n\u00e3o sabemos a densidade correta e colocar isso como uma vari\u00e1vel que podemos controlar." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Tarefas e perguntas\n", "\n", "1. Ache uma estimativa est\u00e1vel com uma densidade baixa. Lembre que o requisito m\u00ednimo \u00e9 ajustar os dados.\n", "2. Ache uma estimativa est\u00e1vel com uma densidade alta.\n", " * O que acontece com o resultado da invers\u00e3o (est\u00e1vel!) quando voc\u00ea aumenta a densidade?\n", " * Se n\u00e3o soubermos a densidade, h\u00e1 algum jeito de estimarmos a profundidade do embasamento corretamente? Que outra informa\u00e7\u00e3o poderia nos ajudar a fazer isso?\n", "3. Tente determinar a densidade correta." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def interativo_densidade(seed, mu, densidade):\n", " misfit.data = utils.contaminate(dado, 0.5, seed=seed)\n", " misfit.model['props']['density'] = densidade\n", " solver = misfit + 10**mu*Smoothness1D(misfit.nparams)\n", " solver.fit()\n", " fig = plota_resultado(solver)\n", " mpl.close(fig) \n", " display(Image(print_figure(fig, dpi=80)))\n", " misfit.data = dado\n", " misfit.model['props']['density'] = bacia.props['density']\n", "\n", "widgets.interact(interativo_densidade, \n", " seed=widgets.IntSliderWidget(min=0, max=10, step=1, value=0),\n", " mu=widgets.FloatSliderWidget(min=-10, max=0, step=0.5, value=-10),\n", " densidade=widgets.FloatSliderWidget(min=-1000, max=-50, step=20, value=-50))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Para verificar se conseguimos determinar a densidade, vamos imprimir o valor correto abaixo:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "print('Densidade correta = {} kg/m3'.format(densidade['density']))" ], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }