{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# 20. Image Domain Least-squares migration\nSeismic migration is the process by which seismic data are manipulated to create\nan image of the subsurface reflectivity.\n\nIn one of the previous tutorials, we have seen how the process can be formulated\nas an inverse problem, which requires access to a demigration-migration engine.\nAs performing repeated migrations and demigrations can be very expensive, an\nalternative approach to obtain accurate and high-resolution estimate of the\nsubsurface reflectivity has emerged under the name of image-domain least-squares\nmigration.\n\nIn image-domain least-squares migration, we identify a direct, linear link between\nthe migrated image $\\mathbf{m}$ and the sought after\nreflectivity $\\mathbf{r}$, namely:\n\n\\begin{align}\\mathbf{m} = \\mathbf{H} \\mathbf{r}\\end{align}\n\nHere $\\mathbf{H}$ is the Hessian, which can be written as:\n\n\\begin{align}\\mathbf{H} = \\mathbf{L}^H \\mathbf{L}\\end{align}\n\nwhere $\\mathbf{L}$ is the demigration operator, whilst its adjoint\n$\\mathbf{L}^H$ is the migration operator. In other words, we say that the\nmigrated image can be seen as the result of a pair of demigration/migration of\nthe reflectivity.\n\nWhilst there exists different ways to estimate $\\mathbf{H}$, the approach\nthat we will be using here entails applying demigration and migration to a special\nreflectivity model composed of regularly space scatterers. What we obtain is the\nspatially-varying impulse response of the migration operator, where each filter is\nalso usually referred to as local point spread function (PSF).\n\nOnce these PSFs are computed (an operation that requires one migration and one\ndemigration, much cheaper than what we do in LSM), the migrated image can be deconvolved\nusing the :py:class:`pylops.signalprocessing.NonStationaryConvolve2D` operator.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\nimport pylops\n\nplt.close(\"all\")\nnp.random.seed(0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To start we create a simple model with 2 interfaces (the same we used in\nthe LSM tutorial) and our PSF model with regularly spaced scatteres\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Velocity Model\nnx, nz = 81, 60\ndx, dz = 4, 4\nx, z = np.arange(nx) * dx, np.arange(nz) * dz\nv0 = 1000  # initial velocity\nkv = 0.0  # gradient\nvel = np.outer(np.ones(nx), v0 + kv * z)\n\n# Reflectivity Model\nrefl = np.zeros((nx, nz))\nrefl[:, 30] = -1\nrefl[:, 50] = 0.5\n\n# PSF Reflectivity Model\npsfrefl = np.zeros((nx, nz))\npsfin = (10, 15)\npsfend = (-10, -5)\npsfj = (30, 30)\n\npsfx = np.arange(psfin[0], nx + psfend[0], psfj[0])\npsfz = np.arange(psfin[1], nz + psfend[1], psfj[1])\nPsfx, Psfz = np.meshgrid(psfx, psfz, indexing=\"ij\")\npsfrefl[psfin[0] : psfend[0] : psfj[0], psfin[1] : psfend[-1] : psfj[-1]] = 1\n\n# Receivers\nnr = 51\nrx = np.linspace(10 * dx, (nx - 10) * dx, nr)\nrz = 20 * np.ones(nr)\nrecs = np.vstack((rx, rz))\ndr = recs[0, 1] - recs[0, 0]\n\n# Sources\nns = 51\nsx = np.linspace(dx * 10, (nx - 10) * dx, ns)\nsz = 10 * np.ones(ns)\nsources = np.vstack((sx, sz))\nds = sources[0, 1] - sources[0, 0]\n\nfig, axs = plt.subplots(1, 3, sharey=True, figsize=(10, 5))\naxs[0].imshow(vel.T, cmap=\"summer\", extent=(x[0], x[-1], z[-1], z[0]))\naxs[0].scatter(recs[0], recs[1], marker=\"v\", s=150, c=\"b\", edgecolors=\"k\")\naxs[0].scatter(sources[0], sources[1], marker=\"*\", s=150, c=\"r\", edgecolors=\"k\")\naxs[0].axis(\"tight\")\naxs[0].set_xlabel(\"x [m]\"), axs[0].set_ylabel(\"z [m]\")\naxs[0].set_title(\"Velocity\")\naxs[0].set_xlim(x[0], x[-1])\n\naxs[1].imshow(refl.T, cmap=\"gray\", extent=(x[0], x[-1], z[-1], z[0]))\naxs[1].scatter(recs[0], recs[1], marker=\"v\", s=150, c=\"b\", edgecolors=\"k\")\naxs[1].scatter(sources[0], sources[1], marker=\"*\", s=150, c=\"r\", edgecolors=\"k\")\naxs[1].axis(\"tight\")\naxs[1].set_xlabel(\"x [m]\")\naxs[1].set_title(\"Reflectivity\")\naxs[1].set_xlim(x[0], x[-1])\n\naxs[2].imshow(psfrefl.T, cmap=\"gray_r\", extent=(x[0], x[-1], z[-1], z[0]))\naxs[2].scatter(recs[0], recs[1], marker=\"v\", s=150, c=\"b\", edgecolors=\"k\")\naxs[2].scatter(sources[0], sources[1], marker=\"*\", s=150, c=\"r\", edgecolors=\"k\")\naxs[2].axis(\"tight\")\naxs[2].set_xlabel(\"x [m]\")\naxs[2].set_title(\"PSF Reflectivity\")\naxs[2].set_xlim(x[0], x[-1])\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now create our Kirchhoff modelling object which we will use to model\nand migrate the data, as well as to model and migrate the PSF model.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nt = 151\ndt = 0.004\nt = np.arange(nt) * dt\nwav, wavt, wavc = pylops.utils.wavelets.ricker(t[:41], f0=20)\n\nkop = pylops.waveeqprocessing.Kirchhoff(\n    z,\n    x,\n    t,\n    sources,\n    recs,\n    v0,\n    wav,\n    wavc,\n    mode=\"analytic\",\n    dynamic=False,\n    wavfilter=True,\n    engine=\"numba\",\n)\nkopdyn = pylops.waveeqprocessing.Kirchhoff(\n    z,\n    x,\n    t,\n    sources,\n    recs,\n    v0,\n    wav,\n    wavc,\n    mode=\"analytic\",\n    dynamic=True,\n    wavfilter=True,\n    aperture=2,\n    angleaperture=50,\n    engine=\"numba\",\n)\n\nd = kop @ refl\nmmig = kopdyn.H @ d\n\ndpsf = kop @ psfrefl\nmmigpsf = kopdyn.H @ dpsf\n\nfig, axs = plt.subplots(1, 2, figsize=(10, 6))\naxs[0].imshow(\n    dpsf[ns // 2, :, :].T,\n    extent=(rx[0], rx[-1], t[-1], t[0]),\n    cmap=\"gray\",\n    vmin=-200,\n    vmax=200,\n)\naxs[0].axis(\"tight\")\naxs[0].set_xlabel(\"x [m]\"), axs[0].set_ylabel(\"t [m]\")\naxs[0].set_title(r\"$d_{psf}$\")\naxs[1].imshow(\n    mmigpsf.T, cmap=\"gray\", extent=(x[0], x[-1], z[-1], z[0]), vmin=-200, vmax=200\n)\naxs[1].scatter(Psfx.ravel() * dx, Psfz.ravel() * dz, c=\"r\")\naxs[1].set_xlabel(\"x [m]\"), axs[1].set_ylabel(\"z [m]\")\naxs[1].set_title(r\"$m_{psf}$\")\naxs[1].axis(\"tight\")\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now extract the local PSFs and create the 2-dimensional\nnon-stationary filtering operator\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "psfsize = (21, 21)\npsfs = np.zeros((len(psfx), len(psfz), *psfsize))\n\nfor ipx, px in enumerate(psfx):\n    for ipz, pz in enumerate(psfz):\n        psfs[ipx, ipz] = mmigpsf[\n            int(px - psfsize[0] // 2) : int(px + psfsize[0] // 2 + 1),\n            int(pz - psfsize[1] // 2) : int(pz + psfsize[1] // 2 + 1),\n        ]\n\nfig, axs = plt.subplots(2, 1, figsize=(10, 5))\naxs[0].imshow(\n    psfs[:, 0].reshape(len(psfx) * psfsize[0], psfsize[1]).T,\n    cmap=\"gray\",\n    vmin=-200,\n    vmax=200,\n)\naxs[0].set_title(r\"$m_{psf}$ iz=0\")\naxs[0].axis(\"tight\")\naxs[1].imshow(\n    psfs[:, 1].reshape(len(psfx) * psfsize[0], psfsize[1]).T,\n    cmap=\"gray\",\n    vmin=-200,\n    vmax=200,\n)\naxs[1].set_title(r\"$m_{psf}$ iz=1\")\naxs[1].axis(\"tight\")\nplt.tight_layout()\n\nCop = pylops.signalprocessing.NonStationaryConvolve2D(\n    hs=psfs, ihx=psfx, ihz=psfz, dims=(nx, nz), engine=\"numba\"\n)\n\nmmigpsf = Cop @ refl\n\nfig, axs = plt.subplots(1, 2, figsize=(10, 5))\naxs[0].imshow(\n    mmig.T, cmap=\"gray\", extent=(x[0], x[-1], z[-1], z[0]), vmin=-1e3, vmax=1e3\n)\naxs[0].set_title(r\"$m_{mig}$\")\naxs[0].axis(\"tight\")\naxs[1].imshow(\n    mmigpsf.T, cmap=\"gray\", extent=(x[0], x[-1], z[-1], z[0]), vmin=-1e3, vmax=1e3\n)\naxs[1].set_title(r\"$m_{mig, psf}$\")\naxs[1].axis(\"tight\")\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, we are ready to invert our seismic image for its corresponding\nreflectivity using the :py:func:`pylops.optimization.sparsity.fista` solver.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "minv, _, resnorm = pylops.optimization.sparsity.fista(\n    Cop, mmig.ravel(), eps=1e5, niter=100, eigsdict=dict(niter=5, tol=1e-2), show=True\n)\nminv = minv.reshape(nx, nz)\n\nfig, axs = plt.subplots(1, 2, figsize=(10, 5))\naxs[0].imshow(\n    mmig.T, cmap=\"gray\", extent=(x[0], x[-1], z[-1], z[0]), vmin=-500, vmax=500\n)\naxs[0].set_title(r\"$m_{mig}$\")\naxs[0].axis(\"tight\")\naxs[1].imshow(minv.T, cmap=\"gray\", extent=(x[0], x[-1], z[-1], z[0]), vmin=-1, vmax=1)\naxs[1].set_title(r\"$m_{inv}$\")\naxs[1].axis(\"tight\")\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For a more advanced set of examples of both reflectivity and impedance\nimage-domain LSM head over to this\n[notebook](https://github.com/mrava87/pylops_notebooks/blob/master/developement/LeastSquaresMigration_imagedomainmarmousi.ipynb).\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
}