# Project C1 • Local inversion for linear seismic imaging

### Principal investigators

 PD Dr. Peer Kunstmann (7/2015 - ) Prof. Dr. Andreas Rieder (7/2015 - )

### Project summary

In this project we reconstruct subsurface material parameters, for example under the earth's or the ocean's surface. To this end, sources on the surface excite pressure waves going through the material below. Then, the geological inner structure scatters these waves and those parts returning to the surface are recorded by receivers. By using these measurements we image these material parameters.
For this purpose, we develop, analyze, and test a migration scheme based on newly introduced imaging operators in two and three dimensions.
From a mathematical point of view we consider a nonlinear parameter identification problem for the elastic wave equation. This problem is solved by a multi-stage process which starts with determining the wave speed from a simpler model: The acoustic wave equation. By linearization (Born approximation) we get Beylkin's model for linear seismic imaging based on the generalized Radon transform where the sound speed is averaged over reflection isochrones connecting sources and receivers (microphones) by points of equal travel time. Thus, the specific instance of the generalized Radon transform depends on the scanning geometry, that is, on the positions of sources and microphones.

## Inverse problem of seismic imaging

The linearized inverse problem of seismic imaging entails solving $\$$\\label{eqn_for_n} Fn(t;\\mathbf{x_r},\\mathbf{x_s})= y(t;\\mathbf{x_r},\\mathbf{x_s}) \$$$ for $n$ where $y$ are pre-processed seismic measurements for several instants of time $t$ and for several source positions $\\mathbf{x_s}$ and receiver positions $\\mathbf{x_r}$. The sought-for quantity $n$ is a dimensionless reflectivity capturing the high frequency variations of the true sound speed $\\nu$, that is, $\\begin{equation*}\\frac{1}{\\nu^2(\\mathbf{x})}= \\frac{1+n(\\mathbf{x})}{c^2(\\mathbf{x})}\\end{equation*}$ for $\\mathbf{x}\\in\\mathbb{R}^d$ where $d\\in\{2,3\}$ is the spatial dimension and $c$ a smooth and known background speed satisfying the geometric optics approximation, that is, no multiple reflections or caustics occur. The involved operator $\$$\\label{eq:genRadon} F w(t,\\mathbf{x_r},\\mathbf{x_s})=\\int \\frac{w(\\mathbf{x})}{c^2(\\mathbf{x})}a(\\mathbf{x},\\mathbf{x_s})a(\\mathbf{x_r},\\mathbf{x})\\delta\\big(t-\\tau(\\mathbf{x},\\mathbf{x_s})-\\tau(\\mathbf{x},\\mathbf{x_r})\\big)\\,\\mathrm{d}\\mathbf{x}\$$$ is a generalized Radon transform (GRT) in which the travel time $\\tau$ and the amplitude $a$ solve the eikonal and transport equation. The processed data $y$ of \\eqref{eqn_for_n} are given by $\\begin{equation*} y(t,\\mathbf{x_r},\\mathbf{x_s})=\\int_0^t \\!\\!(t-r)^{d-2}(u-\\widetilde{u})(r;\\mathbf{x_r},\\mathbf{x_s})\\,\\mathrm{d}r\\end{equation*}$ where $d\\in\{2,3\}$ is the spatial dimension, $u$ are the recorded seismic measurements (reflected wave field) and $\\widetilde{u}$ is the computed reference solution with respect to the known background wave speed $c$.
In the following, we restrict ourselves to two spatial dimensions. However, we obtained similar results for the three dimensional situation [GKQR18b].
For the moment, we work under the following assumptions. We assume the background velocity $c$ to be constant. So, we use $c=1$ for simplicity. Further, with $X=\\mathbb{R}^2_+=\{\\mathbf{x}=(x_1,x_2)\\in\\mathbb{R}^2\\,|\\, x_2>0\}$ we consider $n\\in L^2(X)$ with compact support. At least, for the common offset data acquisition geometry with offset $\\alpha\\ge0$ we work with sources and receivers parameterized by $s\\in\\mathbb{R}$ via $\\begin{equation*} \\mathbf{x_s}(s)=(s-\\alpha,0)^\\top\\quad\\text{ and }\\quad\\mathbf{x_r}(s)=(s+\\alpha,0)^\\top.\\end{equation*}$ Using these assumptions we get $\\tau(\\mathbf{x},\\mathbf{y})=|\\mathbf{x}-\\mathbf{y}|\\quad\\text{and}\\quad a(\\mathbf{x},\\mathbf{y})= \\frac{1}{\\sqrt{|\\mathbf{x}-\\mathbf{y}|}}.$ Thus, the 2D generalized Radon transform \\eqref{eq:genRadon} integrates over ellipses and may be written as $\\begin{equation*} Fn(s,t)=\\int_X A(s,\\mathbf{x})n(\\mathbf{x})\\delta\\big(t-\\varphi(s,\\mathbf{x})\\big)\,\\mathrm{d}\\mathbf{x},\\quad t>2\\alpha, \\end{equation*}$ with $\\varphi(s,\\mathbf{x}):= \\vert \\mathbf{x_s}(s)- \\mathbf{x}\\vert+\\vert \\mathbf{x_r}(s)- \\mathbf{x}\\vert\\quad\\text{and}\\quad A(s,\\mathbf{x})=\\frac{1}{\\sqrt{\\vert \\mathbf{x_s}(s)- \\mathbf{x}\\vert\\,\\vert \\mathbf{x_r}(s)- \\mathbf{x}\\vert}}.$ Further, let $S_0\\subset\\mathbb{R}$ be the bounded open set containing the parameters for the source receiver pairs. For this setting we introduce an imaging operator.

## Our imaging operator

Under the assumptions above we consider the linearized inverse problem \\eqref{eqn_for_n} on the data space $Y=S_0\\times (2\\alpha,\\infty)$. Here, we do not try to reconstruct $n$ directly (there is no inversion formula in general, even in the full data case, $S_0=\\mathbb{R}$). Instead, we define imaging operators of the form $\$$\\label{imaging_op} \\Lambda:=K F^\\dagger \\Psi F \$$$ where $\\Psi\\colon Y\\to [0,\\infty)$ is a smooth compactly supported cut-off function and $K$ is a properly supported pseudo-differential operator on $X$ of non-negative order $m$ with symbol $k$. Further, $F^\\dagger$ is the generalized backprojection operator $F^\\dagger u(\\mathbf{x})= \\int_{S_0}W(s,\\mathbf{x})u(s,\\varphi(s,\\mathbf{x}))\\,\\mathrm{d} s$ for $u\\in\\mathcal{D}(Y)$ where $W$ is a smooth positive weight. Then, the composition $F^\\dagger\\Psi F$ is well-defined for distributions with compact support. Note that the formal $L^2$-adjoint $F^*$ has the weight $W=A$.

## Motivation and Results for explicit choices of $\\Lambda$

Our imaging operator $\\Lambda$ is a pseudo-differential operator of order $m-1$. Thus, the reconstruction $\\Lambda n=KF^\\dagger\\Psi g$ from measured data $g$ (in an ideal setting: $g=Fn$) shows some singularities of $n$. The cut-off $\\Psi$ adds no non-smooth artifacts, so it has no influence on the singularities. For these reasons we chose $\\Lambda$ as in \\eqref{imaging_op}. In directions where $\\Lambda$ is microlocally elliptic the singularities are even emphasized. To study its (microlocal) ellipticity we computed the top order symbol of $\\Lambda$.
Our first choice $F^\\dagger=F^\\ast$ and $K=\\Delta$ with the Laplacian $\\Delta$ resulted in an operator $\\Lambda$ of order $1$. Further studies of the symbol led us to choose $F^\\dagger = F^*$ and $K=\\Delta M +\\alpha\\operatorname{Id}$ with $M$ being the multiplication operator by $x_2$. This operator, also of order $1$, reconstructs jumps in $n$ having the same height but being located at different depths with the same intensities almost independent of $\\alpha$. For the numerical reconstructions from elliptic means $g=Fn$ see the figures above. Here, on the left-hand side $\\Delta F^*\\Psi g$ is used, on the right-hand side $(\\Delta M +\\alpha\\operatorname{Id})F^*\\Psi g$. How we generated those figures, is explained in the next paragraph.

## Numerical experiments: Method of the approximate inverse

The structure of $\\Lambda$ suggest to use the concept of approximate inverse for imaging. This means, we calculate a smoothed version of $\\Lambda n$ by using duality arguments. More details...
Let $\{e_{\\gamma}\}_{\\gamma>0}$ be a family of mollifiers, that is, $e_\\gamma\\to\\delta$ as $\\gamma\\to0$. Further, let $e_{\\gamma}$ be compactly supported in the ball about the origin with radius $\\gamma$. Then, for $\\mathbf{p}\\in X$ and $0<\\gamma < p_2$, we have $\\Lambda n \\ast e_\\gamma ( \\mathbf{p})= \\langle \\Psi g ,\\widetilde{\\psi}_{\\mathbf{p},\\gamma}\\rangle_{L^2(Y)}\\quad\\text{for } g=Fn$ where the reconstruction kernel $\\widetilde\\psi_{\\mathbf{p},\\gamma}$ is given by $\\widetilde{\\psi}_{\\mathbf{p},\\gamma}= F(M\\Delta+\\alpha \\operatorname{Id})e_{\\gamma}(\\,\\cdot\\, - \\mathbf{p}).$

We calculated the reconstruction kernel for a certain family of mollifiers semi-analytically. So, by taking the inner product with the data $g$ we get the smoothed version $\\Lambda n \\ast e_\\gamma ( \\mathbf{p})$ of $\\Lambda n(\\mathbf{p})$. This approach even carries over to other scanning geometries like common source and common midpoint.

## Publications

1. , , , and . Microlocal analysis of imaging operators for effective common offset seismic reconstruction. Inverse Problems, 34(11):114001, 24, November . [preprint] [bibtex]

2. , , , and . Approximate inverse for the common offset acquisition geometry in 2D seismic imaging. Inverse Problems, 34(1):014002, 25, January . [preprint] [bibtex]

## Preprints

1. , , , and . Imaging with the elliptic Radon transform in 3D from an analytical and numerical perspective. CRC 1173 Preprint 2020/13, Karlsruhe Institute of Technology, April . [bibtex]

2. , , , and . Beyond Kirchhoff migration in 2D seismic imaging. Oberwolfach Report, Mathematisches Forschungsinstitut Oberwolfach, May . [bibtex]

## Theses

1. . Seismic imaging with the elliptic Radon transform in 3D: analytical and numerical aspects. PhD thesis, Karlsruhe Institute of Technology (KIT), November . [bibtex]