Keith W Kelly
A scalable algorithm for inverse medium problems with multiple sources

201 E 24th St
Austin
TX 78712
keith@ices.utexas.edu
George Biros

We consider the problem of acoustic scattering as described by the free-space, time-harmonic scalar wave equation given by

$$\Delta v(\mathbf{r}) + k^2(\mathbf{r})v(\mathbf{r}) = -s(\mathbf{r}),$$ (0.1)

along with radiation boundary conditions. Here, $r$ is a point in $\mathbb{R}^3$, $s$ is the source term, and $k$ is the wavenumber. Our formulation is based on potential theory. First we write $k$ as $k_0 + k(r) + \eta(r)$, where $k_0$ is a constant value, $k(r)$ is a known function of space, and $\eta$ is the unknown medium perturbation. Then, given the sources $s(\mathbf{r})$ and detectors outside the perturbation $\eta$ we solve an inverse problem based on (1) for $\eta$. We observe the total field only at some (finite dimensional set of) sensor positions. We use multiple events, where we consider independent activations or different sources at different times. This problem is common in many areas of science and engineering such as seismic imaging, subsurface imaging, impedance tomography, non-destructive evaluation, and diffuse optical tomography.

We describe an algorithm for efficiently solving the inverse scattering problem for the low-frequency time-harmonic scalar wave equation with multi-point illumination. The following approach is based on the FaIMS algorithm described by Chaillat and Biros in [1]. First, we compress the number of incident fields (computed by solving the variable coefficient Helmholtz equation with point sources) using a randomized QR factorization to compute a low rank approximation. By compressing the incident fields, we greatly reduce the problem size and by using a randomized QR factorization we can compute an approximation to within a specified tolerance, ensuring that there is not significant information loss. After compressing the incident fields, we can compute the medium perturbation $\eta$ by solving $\Delta \phi(\mathbf{r}) + k_0^2\phi(\mathbf{r}) = -\eta(\pmb{r})v(\pmb{r})$, where the scattered field $\phi$ is measured as the difference between the incident field and the measured total field. We transform the Helmholtz equation above into a linear integral equation to obtain

$$\phi(\pmb{r}') = \int_{\mathbb{R}^3} G(\pmb{r},\pmb{r}')\eta(\pmb{r})(\phi(\pmb{r}) + u(\pmb{r})) d\pmb{r},$$ (0.2)

where $G$ is the Green's function for the Helmholtz operator with wavenumber $k_0$. Using the Born approximation we can eliminate the dependence of the right hand side of the Lippmann-Schwinger equation on the scattered field. In order to solve for $\eta$ we use the conjugate gradient method on the normal equations with a randomized QR factorization based preconditioner. To improve the accuracy of the computed solution for slightly larger perturbations of the background medium, we use an iterated Born approximation scheme.

We have extended the results of [1] in a few different ways. In the original paper the background medium was assumed to have a constant wavenumber except for the perturbation, i.e. $k(r)$ was zero, but here we allow for variable background medium. We expand the types of scatterers allowed; the original paper allowed only point scatterers while we allow for multiple continuous scattering objects. We also make use use of the iterated Born approximation in our solver. Finally, we achieve superior scalability to the approached referenced in the FaIMS paper. We rapidly evaluate the integral for the forward scatterer using a fast multipole method for volume potentials in parallel. Furthermore, we use PETSc for in our implementation for its fast iterative solvers.