{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Wavelet estimation\nThis example shows how to use the :py:class:`pylops.avo.prestack.PrestackWaveletModelling` to\nestimate a wavelet from pre-stack seismic data. This problem can be written in mathematical\nform as:\n\n\\begin{align}\\mathbf{d}= \\mathbf{G} \\mathbf{w}\\end{align}\n\nwhere $\\mathbf{G}$ is an operator that convolves an angle-variant reflectivity series\nwith the wavelet $\\mathbf{w}$ that we aim to retrieve.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nimport matplotlib.pyplot as plt\nfrom scipy.signal import filtfilt\n\nimport pylops\nfrom pylops.utils.wavelets import ricker\n\nplt.close('all')\nnp.random.seed(0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's start by creating the input elastic property profiles and wavelet\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "nt0 = 501\ndt0 = 0.004\nntheta = 21\n\nt0 = np.arange(nt0)*dt0\nthetamin, thetamax = 0, 40\ntheta = np.linspace(thetamin, thetamax, ntheta)\n\n# Elastic property profiles\nvp = 1200 + np.arange(nt0) + filtfilt(np.ones(5)/5., 1, np.random.normal(0, 80, nt0))\nvs = 600 + vp/2 + filtfilt(np.ones(5)/5., 1, np.random.normal(0, 20, nt0))\nrho = 1000 + vp + filtfilt(np.ones(5)/5., 1, np.random.normal(0, 30, nt0))\nvp[201:] += 500\nvs[201:] += 200\nrho[201:] += 100\n\n# Wavelet\nntwav = 41\nwavoff = 10\nwav, twav, wavc = ricker(t0[:ntwav//2+1], 20)\nwav_phase = np.hstack((wav[wavoff:], np.zeros(wavoff)))\n\n# vs/vp profile\nvsvp = 0.5\nvsvp_z = np.linspace(0.4, 0.6, nt0)\n\n# Model\nm = np.stack((np.log(vp), np.log(vs), np.log(rho)), axis=1)\n\nfig, axs = plt.subplots(1, 3, figsize=(13, 7), sharey=True)\naxs[0].plot(vp, t0, 'k')\naxs[0].set_title('Vp')\naxs[0].set_ylabel(r'$t(s)$')\naxs[0].invert_yaxis()\naxs[0].grid()\naxs[1].plot(vs, t0, 'k')\naxs[1].set_title('Vs')\naxs[1].invert_yaxis()\naxs[1].grid()\naxs[2].plot(rho, t0, 'k')\naxs[2].set_title('Rho')\naxs[2].invert_yaxis()\naxs[2].grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We create now the operators to model a synthetic pre-stack seismic gather\nwith a zero-phase as well as a mixed phase wavelet.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Create operators\nWavesop = \\\n pylops.avo.prestack.PrestackWaveletModelling(m, theta, nwav=ntwav, wavc=wavc,\n vsvp=vsvp, linearization='akirich')\nWavesop_phase = \\\n pylops.avo.prestack.PrestackWaveletModelling(m, theta, nwav=ntwav, wavc=wavc,\n vsvp=vsvp, linearization='akirich')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's apply those operators to the elastic model and create some synthetic data\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "d = (Wavesop*wav).reshape(ntheta, nt0).T\nd_phase = (Wavesop_phase*wav_phase).reshape(ntheta, nt0).T\n\n#add noise\ndn = d + np.random.normal(0, 3e-2, d.shape)\n\nfig, axs = plt.subplots(1, 3, figsize=(13, 7), sharey=True)\naxs[0].imshow(d, cmap='gray', extent=(theta[0], theta[-1], t0[-1], t0[0]),\n vmin=-0.1, vmax=0.1)\naxs[0].axis('tight')\naxs[0].set_xlabel(r'$\\Theta$')\naxs[0].set_ylabel(r'$t(s)$')\naxs[0].set_title('Data with zero-phase wavelet', fontsize=10)\naxs[1].imshow(d_phase, cmap='gray', extent=(theta[0], theta[-1], t0[-1], t0[0]),\n vmin=-0.1, vmax=0.1)\naxs[1].axis('tight')\naxs[1].set_title('Data with non-zero-phase wavelet', fontsize=10)\naxs[1].set_xlabel(r'$\\Theta$')\naxs[2].imshow(dn, cmap='gray', extent=(theta[0], theta[-1], t0[-1], t0[0]),\n vmin=-0.1, vmax=0.1)\naxs[2].axis('tight')\naxs[2].set_title('Noisy Data with zero-phase wavelet', fontsize=10)\naxs[2].set_xlabel(r'$\\Theta$')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can invert the data. First we will consider noise-free data, subsequently\nwe will add some noise and add a regularization terms in the inversion\nprocess to obtain a well-behaved wavelet also under noise conditions.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "wav_est = Wavesop / d.T.flatten()\nwav_phase_est = Wavesop_phase / d_phase.T.flatten()\nwavn_est = Wavesop / dn.T.flatten()\n\n# Create regularization operator\nD2op = pylops.SecondDerivative(ntwav, dtype='float64')\n\n# Invert for wavelet\nwavn_reg_est, istop, itn, r1norm, r2norm = \\\n pylops.optimization.leastsquares.RegularizedInversion(Wavesop, [D2op], dn.T.flatten(),\n epsRs=[np.sqrt(0.1)], returninfo=True,\n **dict(damp=np.sqrt(1e-4),\n iter_lim=200, show=0))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As expected, the regularization helps to retrieve a smooth wavelet\neven under noisy conditions.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# sphinx_gallery_thumbnail_number = 3\nfig, axs = plt.subplots(2, 1, sharex=True, figsize=(8, 6))\naxs[0].plot(wav, 'k', lw=6, label='True')\naxs[0].plot(wav_est, '--r', lw=4, label='Estimated (noise-free)')\naxs[0].plot(wavn_est, '--g', lw=4, label='Estimated (noisy)')\naxs[0].plot(wavn_reg_est, '--m', lw=4, label='Estimated (noisy regularized)')\naxs[0].set_title('Zero-phase wavelet')\naxs[0].grid()\naxs[0].legend(loc='upper right')\naxs[0].axis('tight')\naxs[1].plot(wav_phase, 'k', lw=6, label='True')\naxs[1].plot(wav_phase_est, '--r', lw=4, label='Estimated')\naxs[1].set_title('Wavelet with phase')\naxs[1].grid()\naxs[1].legend(loc='upper right')\naxs[1].axis('tight')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally we repeat the same exercise, but this time we use a *preconditioner*.\nInitially, our preconditioner is a :py:class:`pylops.Symmetrize` operator\nto ensure that our estimated wavelet is zero-phase. After we chain\nthe :py:class:`pylops.Symmetrize` and the :py:class:`pylops.Smoothing1D`\noperators to also guarantee a smooth wavelet.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Create symmetrize operator\nSop = pylops.Symmetrize((ntwav+1)//2)\n\n# Create smoothing operator\nSmop = pylops.Smoothing1D(5, dims=((ntwav+1)//2,), dtype='float64')\n\n# Invert for wavelet\nwavn_prec_est = \\\n pylops.optimization.leastsquares.PreconditionedInversion(Wavesop, Sop,\n dn.T.flatten(),\n returninfo=False,\n **dict(damp=np.sqrt(1e-4),\n iter_lim=200,\n show=0))\n\nwavn_smooth_est = \\\n pylops.optimization.leastsquares.PreconditionedInversion(Wavesop, Sop*Smop,\n dn.T.flatten(),\n returninfo=False,\n **dict(damp=np.sqrt(1e-4),\n iter_lim=200,\n show=0))\n\nfig, ax = plt.subplots(1, 1, sharex=True, figsize=(8, 3))\nax.plot(wav, 'k', lw=6, label='True')\nax.plot(wav_est, '--r', lw=4, label='Estimated (noise-free)')\nax.plot(wavn_prec_est, '--g', lw=4, label='Estimated (noisy symmetric)')\nax.plot(wavn_smooth_est, '--m', lw=4, label='Estimated (noisy smoothed)')\nax.set_title('Zero-phase wavelet')\nax.grid()\nax.legend(loc='upper right')" ] } ], "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 }