{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# 02. Post-stack inversion\nThis tutorial focuses on extending post-stack seismic inversion to GPU\nprocessing. We refer to the equivalent `PyLops tutorial `_\nfor a more detailed description of the theory.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# sphinx_gallery_thumbnail_number = 2\nimport numpy as np\nimport torch\nimport matplotlib.pyplot as plt\nfrom scipy.signal import filtfilt\nfrom pylops.utils.wavelets import ricker\n\nimport pylops_gpu\nfrom pylops_gpu.utils.backend import device\n\ndev = device()\nplt.close('all')\nnp.random.seed(10)\ntorch.manual_seed(10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We consider the 1d example. A synthetic profile of acoustic impedance\nis created and data is modelled using both the dense and linear operator\nversion of :py:class:`pylops_gpu.avo.poststack.PoststackLinearModelling`\noperator. Both model and wavelet are created as numpy arrays and converted\ninto torch tensors (note that we enforce ``float32`` for optimal performance\non GPU).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# model\nnt0 = 301\ndt0 = 0.004\nt0 = np.arange(nt0)*dt0\nvp = 1200 + np.arange(nt0) + \\\n filtfilt(np.ones(5)/5., 1, np.random.normal(0, 80, nt0))\nrho = 1000 + vp + \\\n filtfilt(np.ones(5)/5., 1, np.random.normal(0, 30, nt0))\nvp[131:] += 500\nrho[131:] += 100\nm = np.log(vp*rho)\n\n# smooth model\nnsmooth = 100\nmback = filtfilt(np.ones(nsmooth)/float(nsmooth), 1, m)\n\n# wavelet\nntwav = 41\nwav, twav, wavc = ricker(t0[:ntwav//2+1], 20)\n\n# convert to torch tensors\nm = torch.from_numpy(m.astype('float32'))\nmback = torch.from_numpy(mback.astype('float32'))\nwav = torch.from_numpy(wav.astype('float32'))\n\n# dense operator\nPPop_dense = \\\n pylops_gpu.avo.poststack.PoststackLinearModelling(wav / 2, nt0=nt0,\n explicit=True)\n\n# lop operator\nPPop = pylops_gpu.avo.poststack.PoststackLinearModelling(wav / 2, nt0=nt0)\n\n# data\nd_dense = PPop_dense * m.flatten()\nd = PPop * m.flatten()\n\n# add noise\ndn_dense = d_dense + \\\n torch.from_numpy(np.random.normal(0, 2e-2, d_dense.shape).astype('float32'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now estimate the acoustic profile from band-limited data using either\nthe dense operator or linear operator.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# solve dense\nminv_dense = \\\n pylops_gpu.avo.poststack.PoststackInversion(d, wav / 2, m0=mback, explicit=True,\n simultaneous=False)[0]\n\n# solve lop\nminv = \\\n pylops_gpu.avo.poststack.PoststackInversion(d_dense, wav / 2, m0=mback,\n explicit=False,\n simultaneous=False,\n **dict(niter=500))[0]\n\n# solve noisy\nmn = \\\n pylops_gpu.avo.poststack.PoststackInversion(dn_dense, wav / 2, m0=mback,\n explicit=True, epsI=1e-4,\n epsR=1e0, **dict(niter=100))[0]\n\n\nfig, axs = plt.subplots(1, 2, figsize=(6, 7), sharey=True)\naxs[0].plot(d_dense, t0, 'k', lw=4, label='Dense')\naxs[0].plot(d, t0, '--r', lw=2, label='Lop')\naxs[0].plot(dn_dense, t0, '-.g', lw=2, label='Noisy')\naxs[0].set_title('Data')\naxs[0].invert_yaxis()\naxs[0].axis('tight')\naxs[0].legend(loc=1)\naxs[1].plot(m, t0, 'k', lw=4, label='True')\naxs[1].plot(mback, t0, '--b', lw=4, label='Back')\naxs[1].plot(minv_dense, t0, '--m', lw=2, label='Inv Dense')\naxs[1].plot(minv, t0, '--r', lw=2, label='Inv Lop')\naxs[1].plot(mn, t0, '--g', lw=2, label='Inv Noisy')\naxs[1].set_title('Model')\naxs[1].axis('tight')\naxs[1].legend(loc=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We move now to a 2d example. First of all the model is loaded and\ndata generated.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# model\ninputfile = '../testdata/avo/poststack_model.npz'\n\nmodel = np.load(inputfile)\nm = np.log(model['model'][:, ::3])\nx, z = model['x'][::3]/1000., model['z']/1000.\nnx, nz = len(x), len(z)\n\n\n# smooth model\nnsmoothz, nsmoothx = 60, 50\nmback = filtfilt(np.ones(nsmoothz)/float(nsmoothz), 1, m, axis=0)\nmback = filtfilt(np.ones(nsmoothx)/float(nsmoothx), 1, mback, axis=1)\n\n# convert to torch tensors\nm = torch.from_numpy(m.astype('float32'))\nmback = torch.from_numpy(mback.astype('float32'))\n\n\n# dense operator\nPPop_dense = \\\n pylops_gpu.avo.poststack.PoststackLinearModelling(wav / 2, nt0=nz,\n spatdims=nx, explicit=True)\n\n# lop operator\nPPop = pylops_gpu.avo.poststack.PoststackLinearModelling(wav / 2, nt0=nz,\n spatdims=nx)\n\n# data\nd = (PPop_dense * m.flatten()).reshape(nz, nx)\nn = torch.from_numpy(np.random.normal(0, 1e-1, d.shape).astype('float32'))\ndn = d + n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally we perform different types of inversion\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# dense inversion with noise-free data\nminv_dense = \\\n pylops_gpu.avo.poststack.PoststackInversion(d, wav / 2, m0=mback,\n explicit=True,\n simultaneous=False)[0]\n\n# dense inversion with noisy data\nminv_dense_noisy = \\\n pylops_gpu.avo.poststack.PoststackInversion(dn, wav / 2, m0=mback,\n explicit=True, epsI=4e-2,\n simultaneous=False)[0]\n\n# spatially regularized lop inversion with noisy data\nminv_lop_reg = \\\n pylops_gpu.avo.poststack.PoststackInversion(dn, wav / 2, m0=minv_dense_noisy,\n explicit=False,\n epsR=5e1, epsI=1e-2,\n **dict(niter=80))[0]\n\nfig, axs = plt.subplots(2, 4, figsize=(15, 9))\naxs[0][0].imshow(d, cmap='gray',\n extent=(x[0], x[-1], z[-1], z[0]),\n vmin=-0.4, vmax=0.4)\naxs[0][0].set_title('Data')\naxs[0][0].axis('tight')\naxs[0][1].imshow(dn, cmap='gray',\n extent=(x[0], x[-1], z[-1], z[0]),\n vmin=-0.4, vmax=0.4)\naxs[0][1].set_title('Noisy Data')\naxs[0][1].axis('tight')\naxs[0][2].imshow(m, cmap='gist_rainbow',\n extent=(x[0], x[-1], z[-1], z[0]),\n vmin=m.min(), vmax=m.max())\naxs[0][2].set_title('Model')\naxs[0][2].axis('tight')\naxs[0][3].imshow(mback, cmap='gist_rainbow',\n extent=(x[0], x[-1], z[-1], z[0]),\n vmin=m.min(), vmax=m.max())\naxs[0][3].set_title('Smooth Model')\naxs[0][3].axis('tight')\naxs[1][0].imshow(minv_dense, cmap='gist_rainbow',\n extent=(x[0], x[-1], z[-1], z[0]),\n vmin=m.min(), vmax=m.max())\naxs[1][0].set_title('Noise-free Inversion')\naxs[1][0].axis('tight')\naxs[1][1].imshow(minv_dense_noisy, cmap='gist_rainbow',\n extent=(x[0], x[-1], z[-1], z[0]),\n vmin=m.min(), vmax=m.max())\naxs[1][1].set_title('Trace-by-trace Noisy Inversion')\naxs[1][1].axis('tight')\naxs[1][2].imshow(minv_lop_reg, cmap='gist_rainbow',\n extent=(x[0], x[-1], z[-1], z[0]),\n vmin=m.min(), vmax=m.max())\naxs[1][2].set_title('Regularized Noisy Inversion - lop ')\naxs[1][2].axis('tight')\n\nfig, ax = plt.subplots(1, 1, figsize=(3, 7))\nax.plot(m[:, nx//2], z, 'k', lw=4, label='True')\nax.plot(mback[:, nx//2], z, '--r', lw=4, label='Back')\nax.plot(minv_dense[:, nx//2], z, '--b', lw=2, label='Inv Dense')\nax.plot(minv_dense_noisy[:, nx//2], z, '--m', lw=2, label='Inv Dense noisy')\nax.plot(minv_lop_reg[:, nx//2], z, '--g', lw=2, label='Inv Lop regularized')\nax.set_title('Model')\nax.invert_yaxis()\nax.axis('tight')\nax.legend()\nplt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, if you want to run this code on GPUs, take a look at the following `notebook\n`_\nand obtain more and more speed-up for problems of increasing size.\n\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.12" } }, "nbformat": 4, "nbformat_minor": 0 }