{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 6 - Sísmica de reflexão: Correção de Normal Moveout (NMO) e empilhamento\n", "\n", "Na aula passada, vimos que uma seção CMP é crucial para conseguirmos uma estimativa da velocidade. Nessa prática, veremos como utilizar a correção de NMO para fazer a **análise de velocidades**. Também veremos a técnica do empilhamento e como ela melhora a razão sinal-ruído do dado.\n", "\n", "\n", "Utilizaremos as simulações de ondas da biblioteca [Fatiando a Terra](http://www.fatiando.org). Essas simulações utilizam o [método de diferenças finitas](http://en.wikipedia.org/wiki/Finite_difference_method) para calcular soluções da equação da onda." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Objetivos\n", "\n", "* Utilizar a correção de NMO para estimar velocidades.\n", "* Ver como o empilhamento melhora a qualidade dos dados." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Questão para entregar\n", "\n", "
\n", "\n", "Explique como é feita a análise de velocidades. As velocidades determinadas são as velocidades verdadeiras das camadas? Por que? Qual é vantagem de realizar o empilhamento? \n", "\n", "
\n", "\n", "### Regras para a resposta\n", "\n", "* Coloque **nome, data e o número da prática** em sua resposta. \n", "* A resposta pode ter no **máximo 1 página** (não uma folha).\n", "* **Execute o notebook** antes de responder. As simulações abaixo foram feitas para te ajudar.\n", "* **Pense e organize** sua resposta andtes de começar a escrever." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Instruções\n", "\n", "Esse documento é um [Jupyter notebook](http://jupyter.org/), um documento interativo que mistura texto (como esse), código (como abaixo), e o resultado de executar o código (números, texto, figuras, videos, etc).\n", "\n", "O notebook te fornecerá exemplos interativos que trabalham os temas abordados no questionário. Utilize esses exemplos para responder as perguntas.\n", "\n", "As células com números ao lado, como `In [1]:`, são código [Python](http://python.org/). Algumas dessas células não produzem resultado e servem de preparação para os exemplos interativos. Outras, produzem gráficos interativos. **Você deve executar todas as células, uma de cada vez**, mesmo as que não produzem gráficos.\n", "\n", "**Para executar uma célula**, clique em cima dela e aperte `Shift + Enter`. O foco (contorno verde ou cinza em torno da célula) deverá passar para a célula abaixo. Para rodá-la, aperte `Shift + Enter` novamente e assim por diante. Você pode executar células de texto que não acontecerá nada." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setup\n", "\n", "Rode as células abaixo para carregar os módulos necessários para essa prática." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline\n", "from __future__ import division, print_function\n", "import math\n", "from multiprocessing import Pool, cpu_count\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import ipywidgets as ipw\n", "from fatiando.seismic import RickerWavelet, FDAcoustic2D\n", "from fatiando import utils\n", "from fatiando.vis import anim_to_html\n", "from fatiando.vis.mpl import seismic_image, seismic_wiggle" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulação de um CMP para um modelo de duas camadas\n", "\n", "Rode as células abaixo para simular uma seção CMP para um modelo de duas camadas. A célula abaixo cria nosso modelo." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "shape = (150, 200)\n", "spacing = 10\n", "extent = [0, shape[1]*spacing, shape[0]*spacing, 0]\n", "\n", "# Velocidades\n", "v1, v2 = 4000, 5000\n", "\n", "densidade = np.ones(shape)*1600\n", "velocidade = np.ones(shape)*v1\n", "l1 = 40\n", "densidade[l1:,:] = 1800\n", "velocidade[l1:,:] = v2\n", "l2 = 100\n", "densidade[l2:,:] = 2000\n", "velocidade[l2:,:] = 8000" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Em seguida, precisamos definir onde serão localizadas as fontes e os receptores em nossa simulação. Vamos aproveitar e calcular também os espassamentos (offsets) dos receptores. Lembre-se: offset é a distância da conte ao receptor." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fontes = np.array(list(reversed(range(55, shape[1]//2 - 3, 3))))\n", "recep = np.array([shape[1] - s for s in fontes])\n", "offsets = (recep - fontes)*spacing\n", "print(\"Utilizando {} fontes e {} receptores.\".format(len(fontes), len(recep)))\n", "print('Fontes: {}'.format(fontes*spacing))\n", "print('Receptores: {}'.format(recep*spacing))\n", "print('Offsets: {}'.format(offsets))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Vamos rodar as simulações que precisamos (uma por fonte). **A barra de progresso não irá aparecer** pois vamos rodar as simulações em paralelo para agilizar o processo." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def simul_shot(fonte, its=1200, verbose=False):\n", " shot = FDAcoustic2D(velocidade, densidade, spacing=spacing, taper=0.005, padding=50, verbose=verbose)\n", " shot.add_point_source((0, fonte), RickerWavelet(1, 120))\n", " shot.run(its)\n", " return shot" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%time\n", "print('Simulando...', end=' ')\n", "pool = Pool(processes=cpu_count())\n", "shots = pool.map(simul_shot, fontes)\n", "pool.close()\n", "print('Pronto.')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Vamos dar uma olhada na simulação de 1 tiro para ter uuma ideia do que acontece." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "shots[0].animate(embed=True, dpi=60, fps=8, every=20, cutoff=0.5, cmap='Greys')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Vamos extrair os dados CMP da nossa simulação." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dt = shots[0].dt\n", "times = np.linspace(0, dt*shots[0].simsize, shots[0].simsize)\n", "CMP = np.empty((shots[0].simsize, len(recep)))\n", "for i, s in enumerate(shots):\n", " CMP[:, i] = s[:, 0, recep[i]]\n", " \n", "fig = plt.figure(figsize=(10, 9))\n", "ax = plt.subplot(111)\n", "ax.set_title('CMP')\n", "ax.set_xlabel('offset (m)')\n", "ax.set_ylabel('times (s)')\n", "cutoff = 0.1\n", "ax.imshow(CMP, extent=[offsets.min(), offsets.max(), times.max(), times.min()], \n", " aspect=1000, cmap='Greys', vmin=-cutoff, vmax=cutoff, interpolation='nearest')\n", "fig.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Para pensar\n", "\n", "* Tente identificar a onda direta e as reflexões principais na simulação e no CMP.\n", "* O que é o outro evento que aparece no CMP (acima da segunda reflexão)?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Análise de velocidades" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Agora que temos nosso CMP, podemos aplicar a correção de NMO e fazer a nossa análise de velocidades. Rode a célula abaixo para produzir uma figura interativa." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def nmo_correction(CMP, times, offsets, v):\n", " nmo = np.zeros_like(CMP)\n", " for i, t0 in enumerate(times):\n", " for j, o in enumerate(offsets):\n", " t = np.sqrt(t0**2 + o**2/v[i]**2)\n", " k = int(math.floor(t/dt))\n", " if k < times.size - 1:\n", " # Linear interpolation of the amplitude\n", " y0, y1 = CMP[k, j], CMP[k + 1, j]\n", " x0, x1 = times[k], times[k + 1]\n", " nmo[i, j] = y0 + (y1 - y0)*(t - x0)/(x1 - x0)\n", " return nmo\n", "\n", "def analise_velocidades(CMP):\n", " def aplica_nmo(v1, v2):\n", " v = v1*np.ones_like(times)\n", " v[times > 0.35] = v2\n", " nmo = nmo_correction(CMP, times, offsets, v)\n", "\n", " fig = plt.figure(figsize=(12, 6))\n", " ax = plt.subplot(121)\n", " ax.set_title('CMP')\n", " ax.set_xlabel('offset (m)')\n", " ax.set_ylabel('times (s)')\n", " cutoff = 0.1\n", " ax.imshow(CMP, extent=[offsets.min(), offsets.max(), times.max(), times.min()], \n", " aspect=1000, cmap='Greys', vmin=-cutoff, vmax=cutoff, interpolation='nearest') \n", " ax.grid()\n", "\n", " ax = plt.subplot(122)\n", " ax.set_title('Corrigido de NMO')\n", " ax.set_xlabel('offset (m)')\n", " ax.set_ylabel('times (s)')\n", " cutoff = 0.1\n", " ax.imshow(nmo, extent=[offsets.min(), offsets.max(), times.max(), times.min()], \n", " aspect=1000, cmap='Greys', vmin=-cutoff, vmax=cutoff, interpolation='nearest') \n", " ax.grid()\n", "\n", " fig.tight_layout()\n", " widget = ipw.interactive(aplica_nmo, \n", " v1=ipw.FloatSlider(min=2000, max=6000, step=100, value=2000),\n", " v2=ipw.FloatSlider(min=2000, max=6000, step=100, value=2000))\n", " return widget\n", "\n", "analise_velocidades(CMP)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Figura acima\n", "\n", "* O painel da esquerda mostra nosso CMP original.\n", "* O painel da direita mostra o CMP após aplicação da correção NMO usando as velocidades especificadas.\n", "* Você pode controlar as velocidades utilizadas na correção NMO do nosso CMP." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Para pensar\n", "\n", "* **Determine a velocidade das duas reflexões**.\n", "* Essa velocidade é a velocidade real das camadas? Por que?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Empilhamento\n", "\n", "O CMP que simulamos acima não é muito realista pois não está contaminado com qualquer tipo de ruído. Então vamos sacanear o problema adicionando ruído aleatório nos nossos dados." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ruido = 0.1\n", "CMP_ruido = CMP + np.random.uniform(-ruido, ruido, CMP.shape)\n", "fig = plt.figure(figsize=(10, 9))\n", "ax = plt.subplot(111)\n", "ax.set_title(u'CMP com ruído aleatório')\n", "ax.set_xlabel('offset (m)')\n", "ax.set_ylabel('times (s)')\n", "cutoff = 0.1\n", "ax.imshow(CMP_ruido, extent=[offsets.min(), offsets.max(), times.max(), times.min()], \n", " aspect=1000, cmap='Greys', vmin=-cutoff, vmax=cutoff, interpolation='nearest')\n", "fig.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Uma técnica que nos possibilita diminuir a influência do ruído aleatório é o **empilhamento**. A ideia é somar diversos dados que possuam um sinal coerente no meio de ruído aleatório.\n", "\n", "Rode a célula abaixo para gerar uma figura interativa que ilustra esse conceito." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def empilhamento(ruido):\n", " N = 500\n", " M = 10\n", " if ruido <= 0:\n", " dados_ruido = np.array([np.zeros(N) for i in range(M)])\n", " else:\n", " dados_ruido = np.array([np.random.uniform(-ruido, ruido, size=N) for i in range(M)])\n", " x = np.arange(N)\n", " sinal = 1*utils.gaussian(x, 250, 2)\n", " dados = dados_ruido + sinal\n", " plt.figure(figsize=(10, 6))\n", " plt.subplot(121)\n", " plt.title(u'Dados com ruído e um sinal não-aleatório')\n", " for i, d in enumerate(dados):\n", " plt.plot(d + i + 1, x, '-k')\n", " plt.xlim(0, len(dados) + 1)\n", " plt.xlabel('# do dado')\n", " plt.grid()\n", " ax = plt.subplot(143)\n", " plt.title('Empilhamento')\n", " plt.plot(dados.sum(0), x, '-k')\n", " plt.grid()\n", " plt.xlim(-15, 15)\n", " plt.tight_layout()\n", "ipw.interactive(empilhamento, ruido=ipw.FloatSlider(min=0, max=0.5, step=0.1, value=0))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Figura acima\n", "\n", "* O painel da esquerda mostra 10 conjuntos de dados com 1 sinal coerente no meio.\n", "* O painel da direita mostra o resultado do empilhamento desses dados.\n", "* Você pode controlar a quantidade de ruído aleatório que é inserido nos dados da esquerda." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Para pensar\n", "\n", "* Com o máximo de ruído você consegue enxergar o sinal nos dados originais? E no empilhamento?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Empilhamento após análise de velocidades\n", "\n", "Para poder fazer o empilhamento, precisamos ter um sinal coerente entre todos os traços. Isso significa que a nossa reflexão precisa acontecer no mesmo tempo em todos os traços. Felizmente, é exatamente isso que obtemos após a análise de velocidades e correção NMO.\n", "\n", "Rode a célula abaixo para gerar a figura interativa da análise de velocidades. Dessa vez vamos motrar também o resultado do empilhamento da seção corrigida de NMO." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def nmo_empilhamento(v1, v2):\n", " v = v1*np.ones_like(times)\n", " v[times > 0.35] = v2\n", " nmo = nmo_correction(CMP_ruido, times, offsets, v)\n", " stack = np.atleast_2d(nmo.sum(1)).T\n", " \n", " fig = plt.figure(figsize=(12, 6))\n", " ax = plt.subplot(131)\n", " ax.set_title('CMP')\n", " ax.set_xlabel('offset (m)')\n", " ax.set_ylabel('times (s)')\n", " cutoff = 0.1\n", " ax.imshow(CMP_ruido, extent=[offsets.min(), offsets.max(), times.max(), times.min()], \n", " aspect=2500, cmap='Greys', vmin=-cutoff, vmax=cutoff, interpolation='nearest') \n", " ax.grid()\n", " \n", " ax = plt.subplot(132)\n", " ax.set_title('Corrigido de NMO')\n", " cutoff = 0.1\n", " ax.imshow(nmo, extent=[offsets.min(), offsets.max(), times.max(), times.min()], \n", " aspect=2500, cmap='Greys', vmin=-cutoff, vmax=cutoff, interpolation='nearest') \n", " ax.grid()\n", " \n", " ax = plt.subplot(165)\n", " ax.set_title(u'Primeiro traço do NMO')\n", " seismic_wiggle(np.atleast_2d(nmo[:, 0]).T, dt=dt, scale=1)\n", " ax.set_xlim(-1, 1)\n", " ax.grid()\n", " \n", " ax = plt.subplot(166)\n", " ax.set_title('Empilhamento')\n", " seismic_wiggle(stack, dt=dt, scale=1)\n", " ax.set_xlim(-10, 10)\n", " ax.grid()\n", " \n", " fig.tight_layout(pad=0, w_pad=0)\n", "\n", "ipw.interactive(nmo_empilhamento, \n", " v1=ipw.FloatSlider(min=2000, max=6000, step=100, value=2000),\n", " v2=ipw.FloatSlider(min=2000, max=6000, step=100, value=2000))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Figura acima\n", "\n", "* O painel 1 mostra o CMP contaminado com ruído aleatório.\n", "* O painel 2 mostra o resultado da correção de NMO.\n", "* O painel 3 mostra o primeiro traço da seção corrigida de NMO (sem empilhamento).\n", "* O painel 4 mostra o resultado do empilhamento da seção corrigida de NMO." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Para pensar\n", "\n", "* Como o resultado do empilhamento muda quando você utiliza as velocidades correta para NMO?\n", "* Como o traço empilhado se compara a um traço normal do NMO (sem empilhamento)?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## License and information\n", "\n", "**Course website**: https://github.com/leouieda/geofisica2\n", "\n", "**Note**: This notebook is part of the course \"Geofísica 2\" 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)" ] } ], "metadata": { "anaconda-cloud": {}, "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.12" }, "widgets": { "state": { "53a28390fd914b219c74e4f5a8971f3c": { "views": [ { "cell_index": 20 } ] }, "b7d6ad3786ed4264970e7ab865a251d5": { "views": [ { "cell_index": 30 } ] }, "bbd18e6621ae44f3a2b49d2aa18e28a0": { "views": [ { "cell_index": 26 } ] } }, "version": "1.2.0" } }, "nbformat": 4, "nbformat_minor": 0 }