01_fwi_theory.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 1: The Full-waveform inversion (FWI) algorithm\n",
    "\n",
    "## Format:\n",
    "\n",
    "The presentation will be structured as follows:\n",
    "- Outline the general steps of the FWI algorithm\n",
    "- Introduction to the domain specific language Devito\n",
    "- Building FWI frameworks with Devito\n",
    "\n",
    "For a more interactive experience, the Devito based material can be run live using [Colab](https://colab.research.google.com/):\n",
    "\n",
    "- After navigating to the website, click on the GitHub tab.\n",
    "- Under Enter a GitHub URL or search by organisation or user enter the following: devitocodes/GW4SHM-Workshop.\n",
    "- You will then see the various notebooks from today's session.\n",
    "- After opening a notebook, uncomment the !pip install devito located at the top of each notebook and execute it (Note: You may see an error the first time you execute the cell, in this case simply execute the cell again)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "Many of you will be familiar with ultrasound imaging techniques. Standard ultrasound scans are seismic based images. Medical FWI images are also seismic based, again utilizing ultrasound. So what are the key differences?\n",
    "\n",
    "<img src=\"figures/us.jpg\" alt=\"Drawing\" style=\"width: 450px;\"/>\n",
    "\n",
    "They key differences lie in the data we collect and how we process it:\n",
    "- Standard ultrasound scans have a single wave-source and the emitting device also (often) contains the receiver.\n",
    "- This means that only direct reflections are utilised to form an image.\n",
    "- This vastly reduces the amount of detail you can resolve.\n",
    "- However, this also means the imaging algorithms involved are numerically cheap and images can be rendered almost instantaneously.\n",
    "\n",
    "<img src=\"figures/us_fwi1.png\" alt=\"Drawing\" style=\"width: 800px;\"/>\n",
    "<img src=\"figures/us_fwi2.png\" alt=\"Drawing\" style=\"width: 800px;\"/>\n",
    "\n",
    "Images from [*Guasch et. al. (2020)*](https://www.nature.com/articles/s41746-020-0240-8)\n",
    "\n",
    "- In FWI algorithms on the other hand, we're trying to resolve the full wave-field, we hence need to collect and process **MUCH** more data.\n",
    "- This leads to more accurate images ! And much higher computational complexity and cost..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FWI\n",
    "\n",
    "- FWI is a computational scheme for generating high-resolution, high-fidelity models of physical properties using finite-frequency waves. (The waves could be electromagnetic, acoustic, elastic or of various other kinds).\n",
    "- The method is used in medical imaging of soft tissues, in non-destructive testing, in petroleum exploration, in earthquake seismology, and to image the interior of the Sun.\n",
    "- FWI is a form of tomography, but conventional tomography assumes that energy travels along infinitely thin geometric ray paths, that there are no finite-frequency wave effects. (That is to say, in FWI we're trying to resolve more physics).\n",
    "- Like other simpler forms of tomography, FWI is a local, iterated inversion scheme that successively improves a starting model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some other applications\n",
    "\n",
    "- Ultrasound waves have frequencies of $20KHz$ or more. Seismic based FWI with much lower frequency waves (10-20Hz or so) is also used for geophysical imaging applications such as in volcanology and geophysical exploration."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<tr>\n",
    "    <td> <img src=\"figures/survey-ship-diagram.png\" alt=\"Drawing\" style=\"width: 450px;\"/> </td>\n",
    "    <td> <img src=\"figures/Marmousi3D.png\" alt=\"Drawing\" style=\"width: 450px;\"/> </td>\n",
    "</tr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Left:** Sketch of offshore seismic survey. **Right:** Example model result for $v_p$.\n",
    "\n",
    "- FWI is also used for structural health monitoring e.g. imaging/(locating) damage in damns, mine shafts and other 'difficult to access' structures. In these applications seismics or ground penetrating radar can be used.\n",
    "\n",
    "As larger and larger simulations become more affordable in the coming years, I predict that more and more fields of science and engineering will make use of FWI."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The FWI algorithm outline\n",
    "\n",
    "The aim of FWI is find a model that minimises some measure of the misfit between a dataset predicted by a model and an observed dataset - this measure is called the *objective function*.\n",
    "\n",
    "A simple geometric analogy, in which the model has just two parameters, is to regard the misfit as being represented by the local height of a two-dimensional error surface, and the two model parameters as representing the $x$ and $y$-coordinates of a point on this surface. FWI then involves starting at some point on this surface, and trying to find the bottom of the deepest valley by heading downhill in a sequence of finite steps. To do this, we have to discover which way is downhill, and how far to step. In real FWI, the model has not just two parameters, but many millions, but the analogy is still appropriate. The algorithm proceeds as follows:\n",
    "1. Calculate the direction of the local gradient $\\nabla_\\mathbf{m}$ of the objective function f with respect to the model parameters - this points uphill\n",
    "    - Using the *starting model* $\\mathbf{m}$ and a known *source* $\\mathbf{s}$, calculate the forward *wavefield* $\\mathbf{u}$ everywhere in the model including the *predicted data* $\\mathbf{p}$ at the receivers.\n",
    "    - At the receivers, subtract the observed data d from the predicted data to obtain the *residual data* $\\delta\\mathbf{d}$.\n",
    "    - Treating the receivers as virtual sources, back-propagate the residual data into the model, to generate the residual wavefield $\\delta\\mathbf{u}$.\n",
    "    - Scale the residual wavefield by the local slowness $1/c$, or squared slowness $1/c^2$, and differentiate it twice in time. (But don't worry too much about this, we can abstract it into the maths).\n",
    "    - At every point in the model, cross-correlate the forward and scaled residual wavefields, and take the zero lag in time to generate the *gradient* for one source.\n",
    "    - Do this for every source, and stack together the results to make the global gradient.\n",
    "2. Find the step length - how far is the bottom of the hill?\n",
    "    - Take a small step and a larger step directly downhill, and calculate the objective function at the current model and in these two new models.\n",
    "    - Assume a linear relationship between changes in the model and changes in the residual data so that there will be a parabolic relationship between changes in the model and changes in the objective function, then fit a parabola through these three points.\n",
    "    - The lowest point on this parabola represents the optimal step length (assuming a locally linear relationship).\n",
    "    - Step downhill by the required amount, and update the model.\n",
    "3. Do it all over again\n",
    "    - Use the new model as the starting model, and repeat steps 1. and 2.\n",
    "    - Repeat this process until the model is 'good enough', that is the model is no longer changing (to some numerical tolerance), or we run out of time, money or patience.\n",
    "\n",
    "This is the basic algorithm. There are several ways to enhance and improve it, but nearly all of these involve a greater computational cost (which is already high)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Wave Equation\n",
    "\n",
    "The wave equation is a simplified model for, i.e. the displacement of a vibrating\n",
    "string (approx. 1D), a membrane (such as a drum skin, approx. 2D) or an elastic solid in 3D (the situation relevant to FWI). That is, the main physics the wave equation is attempting to capture is, broadly speaking, the transfer through space of oscillatory energy (vibrations in time).\n",
    "\n",
    "The simplest wave equation that is commonly used in FWI is:\n",
    "\n",
    "\\begin{equation}\n",
    "  \\frac{1}{c^2}\\frac{\\partial^2 u}{\\partial t^2}-\\nabla^2 u = s.\n",
    "\\label{eq:we0} \\tag{1}\n",
    "\\end{equation}\n",
    "\n",
    "- $u$ is the propagating wavefield measured using some appropriate material property (for example electric field in an EM wave or acoustic pressure in an acoustic wave)\n",
    "- $s$ is the driving source that produces the wavefield\n",
    "- $c$ is the wave speed\n",
    "- Note that $u$ and $s$ vary in space and time, and c varies in space.\n",
    "\n",
    "### Other forms of the wave equation\n",
    "\n",
    "More general forms of the wave equation can be written as\n",
    "\n",
    "\\begin{equation}\n",
    "  \\rho(\\mathbf{x})\\frac{\\partial^2 \\mathbf{u}}{\\partial t^2}(\\mathbf{x},t)-\\nabla\\cdot\\mathbf{\\sigma}(\\mathbf{x},t)=\\mathbf{f}(\\mathbf{x},t),\n",
    "  \\label{eq:awe0} \\tag{7}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\mathbf{u}(\\mathbf{x},t)$ is the displacement field and $\\rho(\\mathbf{x})$, $\\mathbf{\\sigma}(\\mathbf{x},t)$, and $\\mathbf{f}(\\mathbf{x},t)$ represent the material density, stress tensor and an external force density respectively. Depending on the fidelity of model we wish to implement, $\\mathbf{\\sigma}(\\mathbf{x},t)$ can take on many different forms.\n",
    "\n",
    "### Matrix form\n",
    "\n",
    "The wave equation represents a linear relationship between a wavefield $u$ and the source $s$ that generates the wavefield. After discretisation (with for example finite differences) we can therefore write $\\eqref{eq:we0}$ as a matrix equation\n",
    "\n",
    "\\begin{equation}\n",
    "  \\mathbf{A}\\mathbf{u}=\\mathbf{s},\n",
    "\\label{eq:we1} \\tag{15}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\mathbf{u}$ and $\\mathbf{s}$ are column vectors that represent the source and wavefield at discrete points\n",
    "in space and time, and $\\mathbf{A}$ is a matrix that represents the discrete numerical implementation of the operator\n",
    "\n",
    "\\begin{equation}\n",
    "  \\frac{1}{c^2}\\frac{\\partial^2}{\\partial t^2}-\\nabla^2.\n",
    "\\label{eq:we3} \\tag{16}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Objective Function\n",
    "\n",
    "- The central purpose of FWI is to find a physical model of the wave-transmitting medium that minimises the difference between an observed dataset and the same dataset as predicted by the model. \n",
    "- A common way of doing this is to measure the $L_2-norm$ of the difference between our predicted data set `p` and the observed data set `d`.\n",
    "\\begin{equation}\n",
    "  \\phi(\\mathbf{m})=\\frac{1}{2}||\\mathbf{p}-\\mathbf{d}||^2=\\frac{1}{2}||\\delta\\mathbf{d}||^2=\\frac{1}{2}\\delta\\mathbf{d}^{T}\\delta\\mathbf{d}.\n",
    "\\label{eq:oe0} \\tag{18}\n",
    "\\end{equation}\n",
    "- We now want to minimize this functional with respect to the model $\\mathbf{m}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Local Inversion\n",
    "\n",
    "Going through the maths, the model update we need, $\\delta\\mathbf{m}$, to minimize the objective function is:\n",
    "\n",
    "\\begin{equation}\n",
    "  \\delta\\mathbf{m} \\approx - \\left(\\frac{\\partial^2 \\phi}{\\partial\\mathbf{m}^2}\\right)^{-1}\\frac{\\partial \\phi}{\\partial\\mathbf{m}} \\equiv -\\mathbf{H}^{-1}\\nabla_{\\mathbf{m}}\\phi.\n",
    "\\label{li4} \\tag{23}\n",
    "\\end{equation}\n",
    "\n",
    "- $\\nabla_{\\mathbf{m}}\\phi$ is the *gradient* of the objective function $\\phi$ with respect to the model parameters\n",
    "- $\\mathbf{H}$ is the *Hessian* matrix of second differentials, both evaluated at $\\mathbf{m}_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Steepest Descent\n",
    "\n",
    "- The Hessian is in general not practical to compute. We could, e.g., simply replace it a scalar $\\alpha$, the step length, and find a suitable value for this:\n",
    "\n",
    "\\begin{equation}\n",
    "  \\delta\\mathbf{m} = -\\alpha\\frac{\\partial \\phi}{\\partial \\mathbf{m}} = -\\alpha\\nabla_{\\mathbf{m}}\\phi .\n",
    "  \\tag{24}\n",
    "\\end{equation}\n",
    "\n",
    "The method that uses this approach is called the method of *steepest descent*, and in its\n",
    "simplest form it consists of the following steps:\n",
    "1. start from a model $\\mathbf{m}_0$,\n",
    "2. evaluate the gradient of the objective function, $\\nabla_{\\mathbf{m}}\\phi$, for the current model,\n",
    "3. find the step length $\\alpha$,\n",
    "4. subtract $\\alpha$ times the gradient from the current model to obtain a new model,\n",
    "5. iterate from step 2 using the new model until the objective function is sufficiently small (or we run out of patience).\n",
    "\n",
    "To implement this, we need a method of calculating the local gradient.\n",
    "\n",
    "**Note:**\n",
    "\n",
    "Simply ignoring the Hessian is a gross simplification and while it\n",
    "is not normally possible to incorporate its effects fully, there are several possibilities for\n",
    "approximating its effects - L-BFGS is widely used, as are conjugate gradients.\n",
    "\n",
    "Useful links:\n",
    "- **L-BFGS**: https://en.wikipedia.org/wiki/Limited-memory_BFGS\n",
    "- **Conjugate gradient**: https://en.wikipedia.org/wiki/Conjugate_gradient_method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating the gradient\n",
    "\n",
    "In principle, we could find the gradient by perturbing each of the model parameters in turn,\n",
    "and calculating what happens to the objective function each time. For $n$ model parameters,\n",
    "that would require $n+1$ modelling runs, and this is not computationally feasible. Fortunately\n",
    "there is a faster way using a solution to the *adjoint* problem.\n",
    "\n",
    "First, write the gradient in terms of the residual data $\\delta\\mathbf{d}=\\mathbf{p}-\\mathbf{d}$:\n",
    "\n",
    "\\begin{equation}\n",
    "  \\nabla_{\\mathbf{m}}\\phi=\\frac{\\partial \\phi}{\\partial \\mathbf{m}}=\\frac{\\partial}{\\partial \\mathbf{m}}\\left(\\frac{1}{2}\\delta\\mathbf{d}^T\\delta\\mathbf{d}\\right)=\\frac{\\partial (\\mathbf{p}-\\mathbf{d})^T}{\\partial \\mathbf{m}}\\delta\\mathbf{d}=\\left(\\frac{\\partial \\mathbf{p}}{\\partial \\mathbf{m}}\\right)^T\\delta\\mathbf{d}.\n",
    "  \\label{eq:rw0} \\tag{25}\n",
    "\\end{equation}\n",
    "\n",
    "After playing around with some maths (don't worry about the details!) we end up with\n",
    "\n",
    "\\begin{equation}\n",
    "  \\nabla_{\\mathbf{m}}\\phi=-\\mathbf{u}^T\\left(\\frac{\\partial \\mathbf{A}}{\\partial \\mathbf{m}}\\right)^T(\\mathbf{A^{-1}})^T\\mathbf{R}^T\\delta\\mathbf{d}.\n",
    "  \\label{eq:rwe6} \\tag{31}\n",
    "\\end{equation}\n",
    "\n",
    "That is, to find the gradient, we must\n",
    "- Calculate the forward wavefield $\\mathbf{u}$\n",
    "- Differentiate the numerical operator $\\mathbf{A}$ with respect to the model parameters (this is an operation that we can do analytically).\n",
    "- We must compute $(\\mathbf{A}^{-1})^T\\mathbf{R}^T\\delta\\mathbf{d}$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpreting the expression for the gradient\n",
    "\n",
    "- The final term our definition of the gradient, $(\\mathbf{A}^{-1})^T\\mathbf{R}^T\\delta\\mathbf{d}$, turns out to correspond to what we call the residual wave field $\\delta\\mathbf{u}$.\n",
    "- $\\delta\\mathbf{u}$ that is generated by a (virtual) source $\\delta\\mathbf{d}$ located at the receivers, and that is propagated by the operator $\\mathbf{A}^T$ which is the *adjoint* of the operator in the original wave equation.\n",
    "\n",
    "- For one source, the gradient calculated this way requires only two modelling runs rather\n",
    "than the $n+1$ modelling runs that direct methods require.\n",
    "- For multi-source datasets, the full gradient is a sum over all sources.\n",
    "- In practical applications with real datasets, the wave equation will nearly always be modified in various ways to include additional physics, but this does not change the underlying approach."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FWI algorithm summary\n",
    "\n",
    "1. Run our experiment in the lab/field/ambulance and collect the observed data $\\mathbf{d}$. (Data will be collected for many shots).\n",
    "2. Choose a wave-equation and hence come up with an appropriate starting model (initial guess) $\\mathbf{m}_0$.\n",
    "3. For each shot, run the forward numerical simulation and compute $\\mathbf{p}$, the modeled data, and then $\\delta\\mathbf{d}=\\mathbf{p}-\\mathbf{d}$, the residual.\n",
    "4. Using our receivers as virtual sources, back-propagate the residual data from each shot to get the 'residual wavefield' $\\delta\\mathbf{u}$. \n",
    "5. For each shot, multiply the residual wavefield by some stuff we know (the modeled wavefield $u$ and an analytic expression) and sum up the contribution from each shot to get the gradient $\\nabla_{\\mathbf{m}}\\phi$.\n",
    "6. Update out model according to $\\mathbf{m}_1=\\mathbf{m}_0 -\\alpha\\nabla_{\\mathbf{m}}\\phi$.\n",
    "\n",
    "Then, repeat steps 2-6 until we're happy with out result!"
   ]
  }
 ],
 "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}