{ "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 6: Inversão para estimar o embasamento de uma bacia" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Preparação\n", "\n", "Como na [prática 5](http://www.leouieda.com/geofisica1/lessons/gravimetria/pratica5.html), esse documento que você está usando é um [IPython notebook](http://ipython.org/notebook.html). Esta prática usará a parte de inversão 2D da biblioteca [Fatiando a Terra](http://fatiando.org) de modelagem geofísica. Para isso, precisamos carregar (`import`) os módulos que vamos usar 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": [ "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" ] }, { "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: Gere um modelo de uma bacia\n", "\n", "Faça um modelo de uma bacia sedimentar e passe para o colega. Algumas instruções:\n", "* O topo da bacia deve ser reto e estar próximo da superfície (z = 0).\n", "* Desenhe a bacia no **sentido horário**.\n", "* **Não coloque erro no dado**. Faremos isso depois.\n", "* Sua bacia deve conter tanto variações abruptas (falhas) quanto variações suaves.\n", "* Faça uma bacia **realista**. Use a geologia que foi entulhada em sua cabeça até agora." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "area = [0, 100e3, 0, 10e3]\n", "x = np.linspace(area[0], area[1], 80)\n", "z = -np.ones_like(x)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "meu_modelo = Moulder(area, x, z, density_range=[-1000, 0])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "meu_modelo.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rode a célula abaixo para salvar o seu modelo. **Mude o nome do arquivo** e troque com um colega. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "meu_modelo.save('modelo_bacia')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tarefa 2: Faça a inversão sem utilizar vínculos" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Primeiro precisamos pegar o modelo do colega e extrair os dados (sem ruído), o polígono 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", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "modelo = Moulder.load('modelo_bacia')\n", "modelo.plot(figsize=(12, 6))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dado = modelo.predicted\n", "bacia = modelo.model[0]\n", "densidade = bacia.props.copy()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A célula abaixo prepara a inversão. `pontos` define quantos pontos terá o polígono que descreverá a bacia. Esse será o número de **parâmetros** que a inversão estimará (a profundidade de cada ponto do polígono).\n", "\n", "Lembrando que a inversão busca achar o **mínimo da função do ajuste**:\n", "\n", "$$\\phi(\\bar{p}) = \\sum\\limits_{i=1}^N [d^o_i - d_i(\\bar{p})]^2$$\n", "\n", "em que $\\bar{p}$ é um vetor com os parâmetros (profundidade de cada ponto do polígono), $d^o_i$ é o i-ésimo dado observado e $d_i$ é o dado predito pelos parâmetros.\n", "\n", "Essa inversão é do tipo **não-linear**. Isso quer dizer que precisamos espeficicar um chute inicial para os parâmetros. A inversão é feita em passos a partir desse chute inicial. Vamos usar \"todos os parâmetros igual a 1 km\"." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "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)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A célula abaixo cria algumas funções que usaremos para fazer um gráfico da nossa solução." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "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" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rode a célula abaixo para criar um **gráfico interativo**. A parte de cima do gráfico mostra os dados observados (pontos) e preditos pela inversão (linha vermelha). A parte de baixo mostra o seu modelo de bacia (preto) e o resultado da inversão (vermelho).\n", "\n", "O botão \"erro\" pode ser movido para mudar o valor do erro que será aplicado aos dados. Cada vez que você mover o botão, o dado é alterado e inversão rodada novamente.\n", "\n", "O botão \"seed\" altera ligeiramente o erro que é aplicado ao dado (mas não o tamanho do erro). \n", "\n", "### Tarefas e perguntas\n", "\n", "1. Execute a célula abaixo para rodar a inversão e criar o gráfico interativo.\n", "2. Aplique erro no dado. Se der uma mensagem de erro, não se preocupe. Quando isso acontece significa que o problema **não possui solução única**. Diminua o erro para voltar a obter uma solução.\n", " * O que acontece com a estimativa da inversão?\n", " * A inversão é capaz de recuperar o modelo original que gerou os dados?\n", " * A inversão recupera o modelo original se o erro for zero? Por que?\n", "3. Use o botão \"seed\" para mudar ligeiramente o erro aplicado.\n", " * O que acontece com a estimativa?\n", " * Essa inversão é estável?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "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))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tarefa 3: Utilizando o vínculo \"damping\" (ou norma mínima)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Vamos utilizar o vínculo de \"norma mínima\" ou [\"damping\"](http://en.wikipedia.org/wiki/Tikhonov_regularization). \n", "\n", "Lembre que adicionar vínculos é equivalente a somar uma outra função na nossa função do ajuste $\\phi(\\bar{p})$:\n", "\n", "$$\\Gamma(\\bar{p}) = \\sum\\limits_{i=1}^N [d^o_i - d_i(\\bar{p})]^2 + \\mu \\sum\\limits_{j=1}^M p_j^2$$\n", "\n", "A função do vínculo é $\\sum\\limits_{j=1}^M p_j^2$ que impõe que os parâmetros estimados sejam o menor possível. Em outras palavras, eles devem ter **norma mínima**.\n", "\n", "A constante positiva $\\mu$ determina qual é o peso do vínculo. $\\mu$ cuida do balanço entre ajustar os dados e obedecer aos vínculos." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Tarefas e perguntas\n", "\n", "1. Execute a célula abaixo para rodar a inversão e criar o gráfico interativo.\n", "2. Aplique erro ao dado.\n", "3. O botão `mi` controla o valor de $\\mu$ (um valor de `-10` significa que $\\mu = 10^{-15}$). Use o botão para aumentar a importância do vínculo de norma mínima.\n", " * O que acontece se $\\mu$ for muito pequeno? O que acontece com o ajuste? Com $\\mu$ pequeno o problema é estável?\n", " * O que acontece se $\\mu$ for muito grande? O que acontece com o ajuste? Com $\\mu$ grande o problema é estável?\n", " * Por que a bacia estimada tende para esse valor quando $\\mu$ é grande?\n", "4. Determine um valor \"ideal\" para $\\mu$. Por ideal quero dizer: é o menor possível para tornar o problema **estável e ajustar os dados**.\n", " * A inversão é capaz de recuperar a falha? Por que?\n", " * E a região suave?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "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))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tarefa 4: Estabilizando com o vínculo de suavidade" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Agora vamos utilizar o vínculo de suavidade.\n", "\n", "Nesse caso, o vínculo é:\n", "\n", "$$\\sum\\limits_{j=1}^{M - 1} (p_{j+1} - p_j)^2$$\n", "\n", "O vínculo de suavidade impõe que os parâmetros sejam o mais próximos o possível de seus vizinhos imediatos." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Tarefas e perguntas\n", "\n", "1. Execute a célula abaixo para rodar a inversão e criar o gráfico interativo.\n", "2. Aplique erro ao dado.\n", "3. Use o botão `mi` para aumentar a importância do vínculo. 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 é estável?\n", " * O que acontece se $\\mu$ for muito grande? O que acontece com o ajuste? Com $\\mu$ grande o problema é estável?\n", " * Por que a bacia estimada tende para esse valor quando $\\mu$ é grande?\n", "4. Determine um valor \"ideal\" para $\\mu$ (menor possível para tornar o problema **estável e ajustar os dados**).\n", " * A inversão é capaz de recuperar a falha? Por que?\n", " * E a região suave?\n", "5. Aumente o erro aplicado ao dado verifique se o problema continua estável. (Se der erro, não se preocupe e aumente `mi`)\n", " * Se o problema deixar de ser estável, o que deve ser feito?\n", " * Como o erro influencia a qualidade da solução?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "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))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tarefa 5: E se você não souber a densidade?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Em casos reais não saberemos exatamente qual é a densidade dos sedimentos. Podemos ter uma ideia mas valores de propriedades físicas costumam variar bastante.\n", "\n", "Vamos assumir que não sabemos a densidade correta e colocar isso como uma variável que podemos controlar." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Tarefas e perguntas\n", "\n", "1. Ache uma estimativa estável com uma densidade baixa. Lembre que o requisito mínimo é ajustar os dados.\n", "2. Ache uma estimativa estável com uma densidade alta.\n", " * O que acontece com o resultado da inversão (estável!) quando você aumenta a densidade?\n", " * Se não soubermos a densidade, há algum jeito de estimarmos a profundidade do embasamento corretamente? Que outra informação poderia nos ajudar a fazer isso?\n", "3. Tente determinar a densidade correta." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "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))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Para verificar se conseguimos determinar a densidade, vamos imprimir o valor correto abaixo:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print('Densidade correta = {} kg/m3'.format(densidade['density']))" ] } ], "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 }