{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# 07. Post-stack inversion\nEstimating subsurface properties from band-limited seismic data represents an\nimportant task for geophysical subsurface characterization.\n\nIn this tutorial, the :py:class:`pylops.avo.poststack.PoststackLinearModelling`\noperator is used for modelling of both 1d and 2d synthetic post-stack seismic\ndata from a profile or 2d model of the subsurface acoustic impedence.\n\n\\begin{align}d(t, \\theta=0) = \\frac{1}{2} w(t) * \\frac{\\mathrm{d}\\ln \\text{AI}(t)}{\\mathrm{d}t}\\end{align}\n\nwhere $\\text{AI}(t)$ is the acoustic impedance profile and $w(t)$ is\nthe time domain seismic wavelet. In compact form:\n\n\\begin{align}\\mathbf{d}= \\mathbf{W} \\mathbf{D} \\mathbf{ai}\\end{align}\n\nwhere $\\mathbf{W}$ is a convolution operator, $\\mathbf{D}$ is a\nfirst derivative operator, and $\\mathbf{ai}$ is the input model.\nSubsequently the acoustic impedance model is estimated via the\n:py:class:`pylops.avo.poststack.PoststackInversion` module. A two-steps\ninversion strategy is finally presented to deal with the case of noisy data.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n\n# sphinx_gallery_thumbnail_number = 4\nimport numpy as np\nfrom scipy.signal import filtfilt\n\nimport pylops\nfrom pylops.utils.wavelets import ricker\n\nplt.close(\"all\")\nnp.random.seed(10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's start with a 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.avo.poststack.PoststackLinearModelling`\noperator.\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) + filtfilt(np.ones(5) / 5.0, 1, np.random.normal(0, 80, nt0))\nrho = 1000 + vp + filtfilt(np.ones(5) / 5.0, 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# dense operator\nPPop_dense = pylops.avo.poststack.PoststackLinearModelling(\n wav / 2, nt0=nt0, explicit=True\n)\n\n# lop operator\nPPop = pylops.avo.poststack.PoststackLinearModelling(wav / 2, nt0=nt0)\n\n# data\nd_dense = PPop_dense * m.ravel()\nd = PPop * m\n\n# add noise\ndn_dense = d_dense + np.random.normal(0, 2e-2, d_dense.shape)" ] }, { "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 = pylops.avo.poststack.PoststackInversion(\n d, wav / 2, m0=mback, explicit=True, simultaneous=False\n)[0]\n\n# solve lop\nminv = pylops.avo.poststack.PoststackInversion(\n d_dense,\n wav / 2,\n m0=mback,\n explicit=False,\n simultaneous=False,\n **dict(iter_lim=2000)\n)[0]\n\n# solve noisy\nmn = pylops.avo.poststack.PoststackInversion(\n dn_dense, wav / 2, m0=mback, explicit=True, epsR=1e0, **dict(damp=1e-1)\n)[0]\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)\nplt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see how inverting a dense matrix is in this case faster than solving\nfor the linear operator (a good estimate of the model is in fact obtained\nonly after 2000 iterations of lsqr). Nevertheless, having a linear operator\nis useful when we deal with larger dimensions (2d or 3d) and we want to\ncouple our modelling operator with different types of spatial regularizations\nor preconditioning.\n\nBefore we move onto a 2d example, let's consider the case of non-stationary\nwavelet and see how we can easily use the same routines in this case\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# wavelet\nntwav = 41\nf0s = np.flip(np.arange(nt0) * 0.05 + 3)\nwavs = np.array([ricker(t0[:ntwav], f0)[0] for f0 in f0s])\nwavc = np.argmax(wavs[0])\n\nplt.figure(figsize=(5, 4))\nplt.imshow(wavs.T, cmap=\"gray\", extent=(t0[0], t0[-1], t0[ntwav], -t0[ntwav]))\nplt.xlabel(\"t\")\nplt.title(\"Wavelets\")\nplt.axis(\"tight\")\n\n# operator\nPPop = pylops.avo.poststack.PoststackLinearModelling(wavs / 2, nt0=nt0, explicit=True)\n\n# data\nd = PPop * m\n\n# solve\nminv = pylops.avo.poststack.PoststackInversion(\n d, wavs / 2, m0=mback, explicit=True, **dict(cond=1e-10)\n)[0]\n\nfig, axs = plt.subplots(1, 2, figsize=(6, 7), sharey=True)\naxs[0].plot(d, t0, \"k\", lw=4)\naxs[0].set_title(\"Data\")\naxs[0].invert_yaxis()\naxs[0].axis(\"tight\")\naxs[1].plot(m, t0, \"k\", lw=4, label=\"True\")\naxs[1].plot(mback, t0, \"--b\", lw=4, label=\"Back\")\naxs[1].plot(minv, t0, \"--r\", lw=2, label=\"Inv\")\naxs[1].set_title(\"Model\")\naxs[1].axis(\"tight\")\naxs[1].legend(loc=1)\nplt.tight_layout()" ] }, { "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.0, model[\"z\"] / 1000.0\nnx, nz = len(x), len(z)\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# dense operator\nPPop_dense = pylops.avo.poststack.PoststackLinearModelling(\n wav / 2, nt0=nz, spatdims=nx, explicit=True\n)\n\n# lop operator\nPPop = pylops.avo.poststack.PoststackLinearModelling(wav / 2, nt0=nz, spatdims=nx)\n\n# data\nd = (PPop_dense * m.ravel()).reshape(nz, nx)\nn = np.random.normal(0, 1e-1, d.shape)\ndn = d + n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally we perform 4 different inversions:\n\n* trace-by-trace inversion with explicit solver and dense operator with\n noise-free data\n\n* trace-by-trace inversion with explicit solver and dense operator\n with noisy data\n\n* multi-trace regularized inversion with iterative solver and linear operator\n using the result of trace-by-trace inversion as starting guess\n\n .. math::\n J = ||\\Delta \\mathbf{d} - \\mathbf{W} \\Delta \\mathbf{ai}||_2 +\n \\epsilon_\\nabla ^2 ||\\nabla \\mathbf{ai}||_2\n\n where $\\Delta \\mathbf{d}=\\mathbf{d}-\\mathbf{W}\\mathbf{AI_0}$ is\n the residual data\n\n* multi-trace blocky inversion with iterative solver and linear operator\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# dense inversion with noise-free data\nminv_dense = pylops.avo.poststack.PoststackInversion(\n d, wav / 2, m0=mback, explicit=True, simultaneous=False\n)[0]\n\n# dense inversion with noisy data\nminv_dense_noisy = pylops.avo.poststack.PoststackInversion(\n dn, wav / 2, m0=mback, explicit=True, epsI=4e-2, simultaneous=False\n)[0]\n\n# spatially regularized lop inversion with noisy data\nminv_lop_reg = pylops.avo.poststack.PoststackInversion(\n dn,\n wav / 2,\n m0=minv_dense_noisy,\n explicit=False,\n epsR=5e1,\n **dict(damp=np.sqrt(1e-4), iter_lim=80)\n)[0]\n\n# blockiness promoting inversion with noisy data\nminv_lop_blocky = pylops.avo.poststack.PoststackInversion(\n dn,\n wav / 2,\n m0=mback,\n explicit=False,\n epsR=[0.4],\n epsRL1=[0.1],\n **dict(mu=0.1, niter_outer=5, niter_inner=10, iter_lim=5, damp=1e-3)\n)[0]\n\nfig, axs = plt.subplots(2, 4, figsize=(15, 9))\naxs[0][0].imshow(d, cmap=\"gray\", extent=(x[0], x[-1], z[-1], z[0]), vmin=-0.4, vmax=0.4)\naxs[0][0].set_title(\"Data\")\naxs[0][0].axis(\"tight\")\naxs[0][1].imshow(\n dn, cmap=\"gray\", extent=(x[0], x[-1], z[-1], z[0]), vmin=-0.4, vmax=0.4\n)\naxs[0][1].set_title(\"Noisy Data\")\naxs[0][1].axis(\"tight\")\naxs[0][2].imshow(\n m,\n cmap=\"gist_rainbow\",\n extent=(x[0], x[-1], z[-1], z[0]),\n vmin=m.min(),\n vmax=m.max(),\n)\naxs[0][2].set_title(\"Model\")\naxs[0][2].axis(\"tight\")\naxs[0][3].imshow(\n mback,\n cmap=\"gist_rainbow\",\n extent=(x[0], x[-1], z[-1], z[0]),\n vmin=m.min(),\n vmax=m.max(),\n)\naxs[0][3].set_title(\"Smooth Model\")\naxs[0][3].axis(\"tight\")\naxs[1][0].imshow(\n minv_dense,\n cmap=\"gist_rainbow\",\n extent=(x[0], x[-1], z[-1], z[0]),\n vmin=m.min(),\n vmax=m.max(),\n)\naxs[1][0].set_title(\"Noise-free Inversion\")\naxs[1][0].axis(\"tight\")\naxs[1][1].imshow(\n minv_dense_noisy,\n cmap=\"gist_rainbow\",\n extent=(x[0], x[-1], z[-1], z[0]),\n vmin=m.min(),\n vmax=m.max(),\n)\naxs[1][1].set_title(\"Trace-by-trace Noisy Inversion\")\naxs[1][1].axis(\"tight\")\naxs[1][2].imshow(\n minv_lop_reg,\n cmap=\"gist_rainbow\",\n extent=(x[0], x[-1], z[-1], z[0]),\n vmin=m.min(),\n vmax=m.max(),\n)\naxs[1][2].set_title(\"Regularized Noisy Inversion - lop \")\naxs[1][2].axis(\"tight\")\naxs[1][3].imshow(\n minv_lop_blocky,\n cmap=\"gist_rainbow\",\n extent=(x[0], x[-1], z[-1], z[0]),\n vmin=m.min(),\n vmax=m.max(),\n)\naxs[1][3].set_title(\"Blocky Noisy Inversion - lop \")\naxs[1][3].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.plot(minv_lop_blocky[:, nx // 2], z, \"--y\", lw=2, label=\"Inv Lop blocky\")\nax.set_title(\"Model\")\nax.invert_yaxis()\nax.axis(\"tight\")\nax.legend()\nplt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That's almost it. If you wonder how this can be applied to real data,\nhead over to the following [notebook](https://github.com/equinor/segyio-notebooks/blob/master/notebooks/pylops/01_seismic_inversion.ipynb)\nwhere the open-source [segyio](https://github.com/equinor/segyio) library\nis used alongside pylops to create an end-to-end open-source seismic\ninversion workflow with SEG-Y input data.\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.9.15" } }, "nbformat": 4, "nbformat_minor": 0 }