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. Continue reading.Collapse content.

Inverse problem of seismic imaging

The linearized inverse problem of seismic imaging entails solving \[\begin{equation} \label{eqn_for_n} Fn(t;\mathbf{x_r},\mathbf{x_s})= y(t;\mathbf{x_r},\mathbf{x_s}) \end{equation}\] 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 \[\begin{equation} \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}\end{equation}\] 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 \[\begin{equation}\label{imaging_op} \Lambda:=K F^\dagger \Psi F \end{equation} \] 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. Continue reading.Collapse content.