| jupytext | kernelspec | ||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
|
In many applications, acoustic or electromagnetic waves are harnessed to see things; under-water acoustics, radar, sonar, medical ultrasound, ground penetrating radar, seismic imaging, global seismology. In some of these applications the measurements are passive; we record waves that are emitted by the object of investigation and we try to infer its location, size, etc. In active measurements, waves are generated and the response of the object is recorded. An example of such an application is shown in {numref}seismic.
---
height: 150px
name: seismic
---
Acquisition setup of a bat and an active marine seismic survey.
In order to treat the inverse problem, we must first understand the forward problem. In its most rudimentary form, the wave-equation is given by
:label: waveequation
L[c]v(t,x) = q(t,x),
with
and where
To define a unique solution of {eq}waveequation we need to supply boundary and initial conditions. In the applications discussed above one typically considers an unbounded domain (
In scattering experiments it is common to rewrite the wave-equation in terms of an incoming and a scattered wavefield,
:label: scattering
L[c_0]v_i(t,x) = q(t,x), \quad L[c_0]v_s(t,x) = -u(x)\frac{\partial^2 v}{\partial t^2}(t,x).
Under a weak scattering assumption, we may ignore the interaction of scattering by
The measurements are typically taken to be the (scattered) wavefield restricted to
Based on this basic setup, we will discuss three possible inverse problems:
-
Inverse source problem: Recover the source
$q$ from the measurements for known (and constant)$c(x)\equiv c_0$ . -
Inverse scattering: Recover the scattering potential
$u(x)$ from measurements of the scattered field for multiple (known) sources, assuming that$c_0$ is known. -
Waveform tomography: Recover
$c(x)$ from measurements of the total wavefield for multiple (known) sources.
Below you will find some typical examples of inverse problems that come up in practice.
-
Earth-quake localization: An earthquake described by a source term
$w(t)q(x)$ is recorded by multiple seimosgraphs at locations$\Delta = {r_k}_{k=1}^n$ . The goal is to recover$q$ in order to determine the location of the earthquake. -
Passive sonar: Soundwaves emitted from an unidentified target are recorded using an array
$\Delta = {x_0 + rp , | , r\in[-h,h] }$ , where$p\in\mathbb{S}^2$ denotes the orientation of the array and$h$ its width. The goal is to recover the source term$w(t)q(x)$ to determine the origin and nature of the source. -
Radar imaging: Incident plane waves, parametrized by direction
$s \in \Sigma \subset \mathbb{S}^{2}$ , are send into the medium and their reflected response is recorded by an array. The goal is retreive the scattering potential. -
Full waveform inversion: In exploration seismology, the goal is to recover
$c$ from measurements of the total wavefield on the surface:$\Sigma = \Delta = {x ,|, n\cdot x = 0}$ . -
Ultrasound tomography: The goal is to recover
$c$ from the total wavefield for sources and receivers surrounding the object.
A few prominent examples of acquisition setups we will consider here are the following:
-
Point-sources:
$q(t,x) = w(t)\delta(x - s)$ . \item Incoming plane waves with frequency$\omega$ and direction$\xi$ :$v_i(t,x) = \sin(\omega(t - \xi\cdot x/c))$ -
Point-measurements along a hyperplane:
$\Delta = {x \in \mathbb{R}^n, |, n\cdot x = x_0}$ . -
Point-measurements on the sphere:
$\Delta = {x \in \mathbb{R}^n , | , |x| = \rho}$ .
We will further assume that the quantities of interest are compactly supported in both space and time and that the measurement time assumption.
---
height: 150px
name: assumption
---
All information on the compactly supported source term is captured in the light cone defined by $t = c\cdot x$. We assume that $T$ is large enough to capture the complete intersection of the light cone with $\Delta$.
For a general source term we may express the solution as a convolution
where
For
for
for
for
The solution for non-constant
$$ v = g_0q - g_0\left(u\cdot \frac{\partial^2 v}{\partial t^2}\right), $$
where
For non-constant
We need to truncate the spatial domain to
which can be discretized using finite-differences as well.
Ignoring the higher order terms leads to
and
$$ \widetilde{v}{i+1,-J} = \widetilde{v}{i,-J} + \frac{c_{i,-J}\Delta t}{\Delta x}\left(\widetilde v_{i,-J+1} - \widetilde v_{i,-J}\right), $$
$$ \widetilde{v}{i+1,J} = \widetilde{v}{i,J} + \frac{c_{i,J}\Delta t}{\Delta x}\left(\widetilde v_{i,J} - \widetilde v_{i,J-1}\right), $$
with
The accuracy of the approximation follows directly from the higher order terms we left out. Another important aspect is stability, which tells us how errors propagate. Since the equations are all linear, errors will propagate according to the same recursion relation. One way of studying this is von Neumann stability analysis, where we study the behaviour of individual components of the error:
with
Define forward map
with
Data are measured at locations
-
Uniqueness. Can we find sources
$r_0$ for which$|Kr_0| = 0$ ? Yes! First define$w_0(t,x)$ compactly supported in$[0,T] \times \Omega$ so that$w_0(t,x) = 0$ for$x\in P$ and set
then
-
Stability. We can also construct sources that radiate an arbitraliy small amount of energy by picking
$w_{\epsilon}$ such that$|w_{\epsilon}| = \mathcal{O}(\epsilon)$ and$|Lw_{\epsilon}| = \mathcal{O}(1)$ as$\epsilon\downarrow 0$ . Then$K(q + r_{\epsilon}) = d + w_{\epsilon}$ and small perturbation in data leads to large perturbation in the solution.
This will be explored in more detail in the assignments.
Under the weak scatterin assumption, the scattered field is given by
which we measure at
We study a variant of the inverse source problem in which
with measurements on the sphere with radius
We see that
using the time-reversed Green function and evaluating at
To see why this works, we study the normal operator
and
so
For
The kernel
is sometimes refered to as the point-spread function.
Under ideal conditions, i.e, measurements are available for all frequencies
We only get contributions to the integral for directions orthogonal to
Thus with a perfect acquisition we have
This filtering can also be implemented in the Fourier domain by multiplying by
:tags: [hide-input]
import numpy as np
import matplotlib.pyplot as plt
import scipy.sparse as sp
from scipy.ndimage import laplace
def solve(q,c,dt,dx,T=1.0,L=1.0,n=1):
'''Solve n-dim wave equation using Leap-Frog scheme: u_{n+1} = Au_n + Bu_{n-1} + Cq_n'''
# define some quantities
gamma = dt*c/dx
nt = int(T/dt + 1)
nx = int(2*L/dx + 2)
#
q = np.resize(q,(nx**n,nt))
# define matrices
A,B,C,L = getMatrices(gamma,nx,n)
# main loop
u = np.zeros((nx**n,nt))
for k in range(1,nt-1):
u[:,k+1] = A@u[:,k] + L@u[:,k] + B@u[:,k-1] + (dx**2)*C@q[:,k]
return u
def getMatrices(gamma,nx,n):
# setup matrices
l = (gamma**2)*np.ones((3,nx))
l[1,:] = -2*(gamma**2)
l[1,0] = -gamma
l[2,0] = gamma
l[0,nx-2] = gamma
l[1,nx-1] = -gamma
if n == 1:
a = 2*np.ones(nx)
a[0] = 1
a[nx-1] = 1
b = -np.ones(nx)
b[0] = 0
b[nx-1] = 0
c = (gamma)**2*np.ones(nx)
c[0] = 0
c[nx-1] = 0
L = sp.diags(l,[-1, 0, 1],shape=(nx,nx))
else:
a = 2*np.ones((nx,nx))
a[0,:] = 1
a[nx-1,:] = 1
a[:,0] = 1
a[:,nx-1] = 1
a.resize(nx**2)
b = -np.ones((nx,nx))
b[0,:] = 0
b[nx-1,:] = 0
b[:,0] = 0
b[:,nx-1] = 0
b.resize(nx**2)
c = (gamma)**2*np.ones((nx,nx))
c[0,:] = 0
c[nx-1,:] = 0
c[:,0] = 0
c[:,nx-1] = 0
c.resize(nx**2)
L = sp.kron(sp.diags(l,[-1, 0, 1],shape=(nx,nx)),sp.eye(nx)) + sp.kron(sp.eye(nx),sp.diags(l,[-1, 0, 1],shape=(nx,nx)))
A = sp.diags(a)
B = sp.diags(b)
C = sp.diags(c)
return A,B,C,L
:tags: [hide-input]
# parameters
L = 1.0
T = 1.0
dx = 1e-2
dt = .5e-2
nt = int(T/dt + 1)
nx = int(2*L/dx + 2)
c = 1
# define source term
u = np.zeros((nx,nx))
u[nx//2 - 10:nx//2+10,nx//2 - 10:nx//2+10] = 1
q = np.zeros((nx*nx,nt))
q[:,1] = u.flatten()
# forward solve
w_forward = solve(q,c,dt,dx,T=T,L=L,n=2)
# sample wavefield
theta = np.linspace(0,2*np.pi,20,endpoint=False)
xs = 0.8*np.cos(theta)
ys = 0.8*np.sin(theta)
I = np.ravel_multi_index(np.array([[xs/dx + nx//2],[ys//dx + nx//2]],dtype=np.int), (nx,nx))
d = w_forward[I,:]
# define adjoint source
r = np.zeros((nx*nx,nt))
r[I,:] = d
r = np.flip(r,axis=1)
# adjoint solve
w_adjoint = solve(r,c,dt,dx,T=T,L=L,n=2)
# plot
fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(8, 5), sharey=True)
plt.gray()
ax[0,0].plot(xs,ys,'r*')
ax[0,0].imshow(w_forward[:,2].reshape((nx,nx)), extent=(-L,L,-L,L))
ax[0,0].set_title('t = 0')
ax[0,1].plot(xs,ys,'r*')
ax[0,1].imshow(w_forward[:,101].reshape((nx,nx)), extent=(-L,L,-L,L))
ax[0,1].set_title('t = 0.5')
ax[0,2].plot(xs,ys,'r*')
ax[0,2].imshow(w_forward[:,200].reshape((nx,nx)), extent=(-L,L,-L,L))
ax[0,2].set_title('t = 1')
ax[1,0].plot(xs,ys,'r*')
ax[1,0].imshow(w_adjoint[:,200].reshape((nx,nx)), extent=(-L,L,-L,L))
ax[1,1].plot(xs,ys,'r*')
ax[1,1].imshow(w_adjoint[:,101].reshape((nx,nx)), extent=(-L,L,-L,L))
ax[1,2].plot(xs,ys,'r*')
ax[1,2].imshow(w_adjoint[:,2].reshape((nx,nx)), extent=(-L,L,-L,L))
plt.show()
Assuming point sources for the incident field,
so the data ellips.
The adjoint operator here is given by
As with the inverse source problem, applying the adjoint already yields an image. This is called migration in practice.
The foward operator has a nice expression in the Fourier domain when using far field measurements. The derivation is left as an excersise.
---
height: 150px
name: ellips
---
A single data point $f(t,s,r)$ is the integral of the scattering potential along an ellips: $|x - s| + |x - r| = t$. Likewise, a single point of $f(x)$ contributes to $d$ along a hyperola $t = |x - s| + |x - r|$.
Here, we aim to recover
where
for
We can solve such a non-linear equation using Newton's method:
where
We find
where the incident field,
Approximating the inverse using backprojection we obtain
where
+++
+++
- Take
$n=1$ ,$c=1$ and source term$\delta(t)u(x)$ with$u$ square integrable and supported on$[-\frac{1}{2},\frac{1}{2}]$ and measurements at$x = 1$ for$t\in [0,T]$ .
Show that the forward operator can be expressed as
and that the operator is bounded.
:class: tip, dropdown
For $n = 1$ the solution to the wave equation with the given source term is given by
$$v(t,x) = \frac{1}{2}\int_{-1/2}^{1/2} \int \delta(t')u(x') H(t - t' - |x - x'|)\mathrm{d}t'\mathrm{d}x' = \frac{1}{2}\int u(x) H(t - |x - x'|) \mathrm{d}x'.$$
For measurements at $x = 1$ this simplifies to
$$f(t) = v(t,1) = \frac{1}{2}\int_{-1/2}^{1/2} u(x') H(t - |1 - x'|) \mathrm{d}x'.$$
This leads to the desired expression becuase we need $t \geq (1-x')$.
To show that this a bounded operator consider
$$\|f\|^2 = \int_0^T \left(\frac{1}{2}\int_{\max(\min(1-t,1/2),0)}^{1/2} u(x')\mathrm{d}x'\right)^2 \mathrm{d}t \leq \frac{1}{4}\int_0^T (1/2-\max(\min(1-t,1/2),0))\int_{-1/2}^{1/2} u(x')^2 \mathrm{d}x\mathrm{d}t = C \|u\|_2^2,$$
where we have used the Cauchy-Schwarz inequality to bound $\left(\int_a^b u(x) \mathrm{d}x\right)^2 \leq \left(\int_a^b 1 \mathrm{d}x \right)\left(\int_a^b u(x)^2 \mathrm{d}x \right).$
- Show that
$u$ can in principle be reconstructed from$f(t)$ with$T = \frac{3}{2}$ with the following reconstruction formula:
:class: tip, dropdown
For $1/2 \leq t \leq 3/2$ we have $f(t) = (1/2)\int_{1-t}^{1/2} u(x) \mathrm{d}x$, so $f'(t) = u(1-t)/2$.
- Show that
$v_{\epsilon}(x) = \sin(2\pi\epsilon x)$ is an almost non-radiating source in the sense that$|K v_{\epsilon}|/|v_{\epsilon}| = \mathcal{O}(\epsilon^{-1})$ as$\epsilon \rightarrow \infty$ .
:class: tip, dropdown
Integration brings out a factor $\epsilon^{-1}$, whereas the norm of $v_{\epsilon}$ is $\mathcal{O}(1)$ as $\epsilon \rightarrow \infty$.
- Now consider noisy measurements
$f^{\delta}(t) = K u(t) + \delta \sin(t/\delta)$ and show that the error in the reconstruction is of order 1, i.e.,
as
:class: tip, dropdown
Using the reconstruction formula we derived before we see that the particular noise term leads to a constant factor.
- In conclusion, is this inverse source problem well-posed? Why (not)?
:class: tip, dropdown
We conclude that the problem is not well-posed as we cannot uniquely and stably reconstruct the solution.
+++
Consider the inverse scattering problem for
- Show that for
$\rho \gg 1$ , the measurements are given by
where
:class: tip, dropdown
We start from the forward operator
$$f(t,r,s) = Ku(t,r,s) = \int_{\Omega} u(x) \frac{\delta''(t - |x - r| - |x - s|)}{|x-r||x-s|}\mathrm{d}x.$$
After Fourier transform in $t$ we get
$$\widehat{f}(\omega,r,s) = \omega^2\int_{\Omega}u(x)\frac{\exp(\imath\omega(|x-r|+|x-s|))}{|x-r||x-s|}.$$
With the far-field approximation ($\rho = |r| = |s| \gg |x|$) $|x-s| \approx |s| + x\cdot s / |s|$ and $|x-r| \approx |r| + x\cdot r / |r|$ and introducing $\xi = s/|s|$, $\eta = r/|r|$ we get
$$\widehat{f}(\omega,\xi,\eta) \approx \frac{\omega^2 e^{2\imath\omega \rho}}{\rho^2}\int_{\Omega}
\frac{\exp(\imath\omega x\cdot(\xi + \eta))}{\rho^2} u(x)\mathrm{d}x,$$
which can be interpreted as a Fourier transform of $u$ evaluated at wavenumber $\omega (\xi + \eta)$.
- Assuming that measurements are available for
$\omega \in [\omega_0,\omega_1]$ , sketch which wavenumbers of$u$ can be stably reconstructed. In what sense is the problem ill-posed?
:class: tip, dropdown
We can retrieve $u$ from $\widehat{u}$ if we have samples of its Fourier transform everywhere. However, the data only gives us samples $\omega (\xi + \eta)$, where $\xi,\eta \in \mathbb{S}^2$. We can thus recover all samples in the unit sphere with radius $2\omega_1$.
For
-
Determine suitable
$\Delta t$ ,$\Delta x$ and$T$ and generate data for scattering potential$u(x) = \exp(-200|x|^2)$ and an incident plane wave with$\xi = (1,0)$ and$\omega = 10\pi$ . -
Construct a non-scattering potential for this incident field and measurements by using the result you obtained in the previous exercise.
-
Can you construct a non-scattering potential that is invisible from multiple directions?
:tags: [hide-cell]
import numpy as np
import scipy.sparse as sp
import matplotlib.pyplot as plt
def solve(q,c,dt,dx,T=1.0,L=1.0,n=1):
'''Solve n-dim wave equation using Leap-Frog scheme: u_{n+1} = Au_n + Bu_{n-1} + Cq_n'''
# define some quantities
gamma = dt*c/dx
nt = int(T/dt + 1)
nx = int(2*L/dx + 2)
#
q.resize((nx**n,nt))
# define matrices
A,B,C,L = getMatrices(gamma,nx,n)
# main loop
u = np.zeros((nx**n,nt))
for k in range(1,nt-1):
u[:,k+1] = A@u[:,k] + L@u[:,k] + B@u[:,k-1] + (dx**2)*C@q[:,k]
return u
def multiply(u,c,dt,dx,T=1.0,L=1.0,n=1):
# define some quantities
gamma = dt*c/dx
nt = int(T/dt + 1)
nx = int(2*L/dx + 2)
#
u.resize((nx**2,nt))
# define matrices
A,B,C,L = getMatrices(gamma,nx,n)
# main loop
q = np.zeros((nx**n,nt))
for k in range(1,nt-1):
q[k] = (u[:,k+1] - 2*u[:,k] + u[:,k-1] - L@u[:,k])/(dt*c)**2
return u
def sample(xin,xout1,xout2=[]):
'''Spatial sampling by simple interpolation'''
if len(xout2):
n = 2
else:
n = 1
m = len(xout1)
nx = len(xin)
rw = []
cl = []
nz = []
if n == 1:
for k in range(m):
i = 0
while xin[i] < xout1[k]:
i = i + 1
if i < nx - 1:
a = (xout1[k] - xin[i+1])/(xin[i] - xin[i+1])
b = (xout1[k] - xin[i])/(xin[i+1] - xin[i])
rw.append(k)
cl.append(i)
nz.append(a)
rw.append(k)
cl.append(i+1)
nz.append(b)
P = sp.coo_matrix((nz,(rw,cl)),shape=(m,nx))
else:
for k in range(m):
i = 0
j = 0
while xin[i] < xout1[k]:
i = i + 1
while xin[j] < xout2[k]:
j = j + 1
if i < nx - 1 and j < nx - 1:
a = (xout1[k] - xin[i+1])*(xout2[k] - xin[j+1])/(xin[i] - xin[i+1])/(xin[j] - xin[j+1])
b = (xout1[k] - xin[i])*(xout2[k] - xin[j+1])/(xin[i+1] - xin[i])/(xin[j] - xin[j+1])
c = (xout1[k] - xin[i+1])*(xout2[k] - xin[j])/(xin[i] - xin[i+1])/(xin[j+1] - xin[j])
d = (xout1[k] - xin[i])*(xout2[k] - xin[j])/(xin[i+1] - xin[i])/(xin[j+1] - xin[j])
rw.append(k)
cl.append(i+nx*j)
nz.append(a)
rw.append(k)
cl.append(i+1+nx*j)
nz.append(b)
rw.append(k)
cl.append(i+nx*(j+1))
nz.append(c)
rw.append(k)
cl.append(i+1+nx*(j+1))
nz.append(d)
P = sp.coo_matrix((nz,(rw,cl)),shape=(m,nx*nx))
return P
def getMatrices(gamma,nx,n):
# setup matrices
l = (gamma**2)*np.ones((3,nx))
l[1,:] = -2*(gamma**2)
l[1,0] = -gamma
l[2,0] = gamma
l[0,nx-2] = gamma
l[1,nx-1] = -gamma
if n == 1:
a = 2*np.ones(nx)
a[0] = 1
a[nx-1] = 1
b = -np.ones(nx)
b[0] = 0
b[nx-1] = 0
c = (gamma)**2*np.ones(nx)
c[0] = 0
c[nx-1] = 0
L = sp.diags(l,[-1, 0, 1],shape=(nx,nx))
else:
a = 2*np.ones((nx,nx))
a[0,:] = 1
a[nx-1,:] = 1
a[:,0] = 1
a[:,nx-1] = 1
a.resize(nx**2)
b = -np.ones((nx,nx))
b[0,:] = 0
b[nx-1,:] = 0
b[:,0] = 0
b[:,nx-1] = 0
b.resize(nx**2)
c = (gamma)**2*np.ones((nx,nx))
c[0,:] = 0
c[nx-1,:] = 0
c[:,0] = 0
c[:,nx-1] = 0
c.resize(nx**2)
L = sp.kron(sp.diags(l,[-1, 0, 1],shape=(nx,nx)),sp.eye(nx)) + sp.kron(sp.eye(nx),sp.diags(l,[-1, 0, 1],shape=(nx,nx)))
A = sp.diags(a)
B = sp.diags(b)
C = sp.diags(c)
return A,B,C,L
:tags: [hide-cell]
# grid
nt = 201
nx = 101
x = np.linspace(-1,1,nx)
y = np.linspace(-1,1,nx)
t = np.linspace(0,1,nt)
xx,yy,tt = np.meshgrid(x,y,t)
# velocity
c = 1.0
# scattering potential
a = 200;
r = np.exp(-a*xx**2)*np.exp(-a*yy**2);
# incident field, plane wave at 10 Hz.
omega = 2*np.pi*5
ui = np.sin(omega*(tt - (xx + 1)/c))
# solve
us = solve(r*ui,c,t[1]-t[0],x[1]-x[0],n=2)
# sample
xr = 0.8*np.ones(101)
yr = np.linspace(-1,1,101)
P = sample(x,xr,yr)
d = P@us
# plot
ui.resize((nx,nx,nt))
us.resize((nx,nx,nt))
plt.subplot(241)
plt.imshow(ui[:,:,1])
plt.clim((-1,1))
plt.subplot(242)
plt.imshow(ui[:,:,51])
plt.clim((-1,1))
plt.subplot(243)
plt.imshow(ui[:,:,101])
plt.clim((-1,1))
plt.subplot(244)
plt.imshow(ui[:,:,151])
plt.clim((-1,1))
plt.subplot(245)
plt.imshow(us[:,:,1])
plt.clim((-.001,.001))
plt.subplot(246)
plt.imshow(us[:,:,51])
plt.clim((-.001,.001))
plt.subplot(247)
plt.imshow(us[:,:,101])
plt.clim((-.001,.001))
plt.subplot(248)
plt.imshow(us[:,:,151])
plt.clim((-.001,.001))
plt.show()
For
- Express the solution
$v(t,x)$ by using the Green function
The measurements are given by
-
Plot the objective as a function of
$c$ for various values of$\sigma$ . -
Do you think Newton's method will be successful in retrieving the correct
$c$ ?