Reconstruction Algorithms for Photoacoustic Tomography in Heterogeneous Damping Media

Abstract

In this article, we study several reconstruction methods for the inverse source problem of photoacoustic tomography with spatially variable sound speed and damping. The backbone of these methods is the adjoint operators, which we thoroughly analyze in both the $L^2$- and $H^1$-settings. They are casted in the form of a nonstandard wave equation. We derive the well posedness of the aforementioned wave equation in a natural functional space and also prove the finite speed of propagation. Under the uniqueness and visibility condition, our formulations of the standard iterative reconstruction methods, such as Landweber’s and conjugate gradients (CG), achieve a linear rate of convergence in either $L^2$- or $H^1$-norm. When the visibility condition is not satisfied, the problem is severely ill posed and one must apply a regularization technique to stabilize the solutions. To that end, we study two classes of regularization methods: (i) iterative and (ii) variational regularization. In the case of full data, our simulations show that the CG method works best; it is very fast and robust. In the ill-posed case, the CG method behaves unstably. Total variation regularization method (TV), in this case, significantly improves the reconstruction quality.

Appendix

1.1 Proof of Theorem 4

Let $B_R$ denote the ball of radius R centered at the origin and $R:=R_0 + c_+ T$, where $R_0$ satisfies $\varOmega \subset B_{R_0}$. Let $H_0^1(B_R)$ be the closure of $C_0^\infty (B_R)$ with respect to the norm

$$\begin{aligned} \Vert f\Vert _{H_0^1(B_R)} = \left[ \int _{B_R} |\nabla f|^2 \mathrm{d}x \right] ^{1/2}. \end{aligned}$$

Our proof is divided into two steps:

Step 1 There exists a weak solution q of (5) on $B_R$. That is,

• i’) $q \in L^2([0,T];H_0^1(B_R))$, $q' \in L^2([0,T];L^2(B_R))$, $q'' \in L^2([0,T];H^{-1}(B_R))$,

• ii’) $q(0)=0$ and $q'(0) =0$, and

• iii’) for any function $\phi \in H_0^1(B_R)$

  $$\begin{aligned}&\int _{B_R} c^{-2}(x) \, q_{tt}(x,t) \, \phi (x) \,\mathrm{d}x \\&\qquad + \int _{B_R} a(x) \, q_{t}(x,t) \, \phi (x) \,\mathrm{d}x \\&\qquad + \int _{B_R} \nabla q(x,t) \, \nabla \phi (x) \,\mathrm{d}x \\&\quad = - \int _{\partial \varOmega } \, g(y,t) \, \phi (y) \, \mathrm{d}y \,\quad \text{ a.e } t \in [0,T]. \end{aligned}$$

Step 2 The solution q in Step 1 satisfies: $q(x,t) =0$ for all $(x,t) \in \varOmega ^c \times [0,T]$ such that $dist(x,\partial \varOmega ) \ge c_+ t$.

Once both steps are proved, the solution q of Eq. (5) is just the trivial extension of q into $[0,T] \times \mathbb R^d$. Let us now proceed to prove those steps.

Proof of Step 1 Let $\{\phi _{k}\}_k$ be an orthogonal basis of $H_0^1(B_R)$.² For any integer N, we define

$$\begin{aligned} q_N(x,t) = \sum _{i=1}^N d_i(t) \phi _i(x) \end{aligned}$$

to be a solution of the system

$$\begin{aligned}&\int _{B_R} c^{-2}(x) \, q_{N,tt}(x,t) \, \phi _i(x) \,\mathrm{d}x \nonumber \\&\qquad + \int _{B_R} a(x) \, q_{N,t}(x,t) \, \phi _i(x) \,\mathrm{d}x \nonumber \\&\qquad +\int _{B_R} \nabla q_N(x,t) \, \nabla \phi _i(x) \,\mathrm{d}x \nonumber \\&\quad = - \int _{\partial \varOmega } g(y,t) \, \phi _i(y) \, \mathrm{d}y, \quad i =1,\dots ,N, \end{aligned}$$

(13)

together with the initial condition $q_N(x,0) = q_{N,t}(x,0) =0$. Since the above system is a standard linear ODE system for $(d_1,\dots , d_N)$, $q_N$ uniquely exists. Multiplying each equation by $d_i'(t)$ and summing them up, we obtain:

$$\begin{aligned}&\int _{B_R} c^{-2}(x) q_{N,tt} (x,t) q_{N,t}(x,t) \,\mathrm{d}x \\&\qquad + \int _{B_R} a(x) \, [q_{N,t}(x,t)]^2 \,\mathrm{d}x \\&\qquad +\int _{B_R} \nabla q_N(x,t) \, \nabla q_{N,t} \,\mathrm{d}x\\&\quad = - \int _{\partial \varOmega } g(y,t) \, q_{N,t} (y,t) \, \mathrm{d}y. \end{aligned}$$

This implies

$$\begin{aligned}&\frac{1}{2} \frac{\mathrm{d}}{\mathrm{d}t} \left[ \int _{B_R} c^{-2}(x) |q_{N,t}(x,t)|^2 \,\mathrm{d}x + \int _{B_R} |\nabla q_N(x,t)|^2 \mathrm{d}x \right] \\&\quad \le - \int _{\partial \varOmega } g(y,t) \, q_{N,t}(y,t) \, \mathrm{d}y. \end{aligned}$$

Taking the integration of both sides with respect to t and using the initial conditions for $q_N$:

$$\begin{aligned}&\frac{1}{2} \left[ \Vert q_{N,t}(\,\cdot \,, t)\Vert ^2 + \Vert q_N(\,\cdot \,,t)\Vert ^2_{H_0^1(B_R)} \right] \le \\&\quad - \int _{\partial \varOmega } g(y,t) \, q_N(y,t) \, \mathrm{d}y + \int _0^t \int _{\partial \varOmega } g_t(y,t) \, q_N(y,t) \, \mathrm{d}y. \end{aligned}$$

Bounding the first term of the right-hand side, we obtain

$$\begin{aligned}&\frac{1}{2} \left[ \Vert q_{N,t}(\,\cdot \,,t)\Vert ^2 + \Vert q_N(\,\cdot \,,t)\Vert ^2_{H_0^1(B_R)} \right] \\&\quad \le \Vert g(\,\cdot \,,t)\Vert _{H^{-1/2}(\partial \varOmega )} \Vert q_N(\,\cdot \,,t)\Vert _{H^{1/2}(\partial \varOmega )} \\&\qquad + \int _0^t \Vert g_t(\,\cdot \,,t)\Vert ^2_{H^{-1/2}(\partial \varOmega )} + \int _0^t \Vert q_N(\,\cdot \,,t)\Vert ^2_{H^{1/2}(\partial \varOmega )}. \end{aligned}$$

Now, Young’s inequality gives

$$\begin{aligned}&\frac{1}{2} \left[ \Vert q_{N,t}(\,\cdot \,,t)\Vert ^2 + \Vert q_N(\,\cdot \,,t)\Vert ^2_{H_0^1(B_R)} \right] \\&\quad \le A \Vert g(\,\cdot \,,t)\Vert ^2_{H^{-1/2}(\partial \varOmega )} \\&\quad + \frac{1}{2A}\Vert q_N(\,\cdot \,,t)\Vert ^2_{H^{1/2}(\partial \varOmega )} + \int _0^t \Vert g_t(\,\cdot \,,t)\Vert ^2_{H^{-1/2}(\partial \varOmega )} \\&\quad + \int _0^t \Vert q_N(\,\cdot \,,t)\Vert _{H^{1/2}(\partial \varOmega )}, \end{aligned}$$

where $A>0$ can be any constant, whose value will be specified later. Noting that $\Vert q_N(\,\cdot \,,t)\Vert _{H^{1/2}(\partial \varOmega )} \le C \Vert q_N(\,\cdot \,,t)\Vert _{H_0^{1}(B_R)}$ we obtain by choosing A big enough

$$\begin{aligned}&\frac{1}{2} \left[ \Vert q_{N,t}(\,\cdot \,,t)\Vert ^2 + \Vert q_N(\,\cdot \,,t)\Vert ^2_{H_0^1(B_R)} \right] \\&\quad \le A \Vert g(\,\cdot \,,t)\Vert ^2_{H^{-1/2}(\partial \varOmega )} \\&\qquad + \frac{1}{4} \Vert q_N(\,\cdot \,,t)\Vert ^2_{H^{1}_0(B_R)} + \int _0^t \Vert g_t(\,\cdot \,,t)\Vert ^2_{H^{-1/2}(\partial \varOmega )} \\&\qquad +\,C \int _0^t \Vert q_N(\,\cdot \,,t)\Vert ^2_{H^{1}_0(B_R)}. \end{aligned}$$

Here and in the sequel, C is a generic constant whose value may vary from one place to another. Therefore,

$$\begin{aligned}&\Vert q_{N,t}(\,\cdot \,,t)\Vert ^2 + \Vert q_N(\,\cdot \,,t)\Vert ^2_{H_0^1(B_R)} \le C \big ( \Vert g(\,\cdot \,,t)\Vert ^2_{H^{-1/2}(\partial \varOmega )} \\&\quad + \int _0^T \Vert g_t(\,\cdot \,,t)\Vert ^2_{H^{-1/2}(\partial \varOmega )}\\&\quad + \int _0^t \Vert q_N(\,\cdot \,,t)\Vert ^2_{H_0^{1}(B_R)}\big ), \quad t \in [0,T]. \end{aligned}$$

Let $E_N(t) := \int _0^t \Vert q_{N,t}(\,\cdot \,,t)\Vert ^2 + \Vert q_N(\,\cdot \,,t)\Vert ^2_{H_0^{1}(B_R)}$. We arrive at

$$\begin{aligned}&E_N'(t) - C E_N(t) \\&\quad \le C \big ( \Vert g(\,\cdot \,,t)\Vert ^2_{H^{-1/2}(\partial \varOmega )} \\&\qquad + \Vert g_t\Vert ^2_{L^2([0,T],H^{-1/2}(\partial \varOmega ))} \big ),~ t \in [0,T]. \end{aligned}$$

From the Grownwall’s inequality, we obtain

$$\begin{aligned}&E_N(T) \le C( \Vert g\Vert ^2_{L^2([0,T],H^{-1/2}(\partial \varOmega ))} \nonumber \\&\quad +\,\Vert g_t\Vert ^2_{L^2([0,T],H^{-1/2}(\partial \varOmega ))} ). \end{aligned}$$

(14)

Since C is a constant independent of N, $\{q_{N}\}$ and $\{q_{N,t}\}$ are bounded sequences in $L^2([0,T],H^1_0(B_R))$ and $L^2([0,T];L^2(B_R))$, respectively. After possibly passing over to subsequences, we obtain $q_N \rightharpoonup q$ in $L^2([0,T];H^1_0(B_R))$ and $q_{N,t} \rightharpoonup q_1$ in $L^2([0,T];L^2(B_R))$. It is easy to show that $q_1=q'$. Since $\{\phi _k\}$ is a basis of $H_0^1(B_R)$, from (13), we obtain for any $v \in L^2([0,T];H_0^1(\varOmega ))$:

$$\begin{aligned}&\lim _{N \rightarrow \infty } \int _0^T \int _{\mathbb R^d} c^{-2}(x) \, q_{N,tt}(x,t) \, v(x,t) \,\mathrm{d}x \mathrm{d}t \\&\qquad + \int _0^T \int _{\mathbb R^d} a(x) \, q_{t}(x,t) \, v(x,t) \,\mathrm{d}x \mathrm{d}t \\&\qquad + \int _0^T \int _{\mathbb R^d} \nabla q(x,t) \, \nabla v(x,t) \,\mathrm{d}x \\&\quad = - \int _{\partial \varOmega } g(y,t) \, v(y,t) \, \mathrm{d}y. \end{aligned}$$

That is, $q_{N,tt}$ converges to an element in $L^2([0,T], H^{-1}(B_R))$. That is, $q_{tt} \in L^2([0,T], H^{-1}(B_R))$ and

$$\begin{aligned}&\int _0^T \int _{\mathbb R^d} c^{-2}(x) \, q_{tt}(x,t) \, v(x,t) \,\mathrm{d}x \,\mathrm{d}t\\&\qquad + \int _0^T \int _{\mathbb R^d} a(x) \, q_{t}(x,t) \, v(x,t) \,\mathrm{d}x \mathrm{d}t \\&\qquad + \int _0^T \int _{\mathbb R^d} \nabla q(x,t) \, \nabla v(x,t) \,\mathrm{d}x \mathrm{d}t \\&\quad = - \int _0^T \int _{\partial \varOmega } g(y,t) \, v(y,t) \, \mathrm{d}y \, \mathrm{d}t. \end{aligned}$$

Let $\phi \in H_0^1(B_R)$. For any $t_0 \in (0,T)$, choosing³$v(x,t) = \phi (x) \chi _{[t_0-\epsilon , t_0+ \epsilon ]}(t)$, we obtain

$$\begin{aligned}&\int _{t_0-\epsilon }^{t_0+\epsilon } \int _{\mathbb R^d} c^{-2}(x) \, q_{tt}(x,t) \, \phi (x) \,\mathrm{d}x \,\mathrm{d}t \\&\qquad + \int _{t_0-\epsilon }^{t_0+ \epsilon } \int _{\mathbb R^d} a(x) \, q_{t}(x,t) \, \phi (x) \,\mathrm{d}x \mathrm{d}t \\&\qquad + \int _{t_0-\epsilon }^{t_0+\epsilon } \int _{\mathbb R^d} \nabla q(x,t) \, \nabla \phi (x) \,\mathrm{d}x \mathrm{d}t \\&\quad = - \int _{t_0 - \epsilon }^{t_0+\epsilon } \int _{\partial \varOmega } g(y,t) \, \phi (y) \, \mathrm{d}y \, \mathrm{d}t. \end{aligned}$$

Dividing both sides by $2 \epsilon$ and send $\epsilon \rightarrow 0$, we obtain

$$\begin{aligned}&\int _{B_R} c^{-2}(x) \, q_{tt}(x,t_0) \, \phi (x) \,\mathrm{d}x \\&\qquad +\int _{B_R} a(x) \, q_{t}(x,t_0) \, \phi (x) \,\mathrm{d}x\\&\qquad + \int _{B_R} \nabla q(x,t_0) \, \nabla \phi (x) \,\mathrm{d}x \\&\quad = - \int _{\partial \varOmega } \, g(y,t_0) \, \phi (y) \, \mathrm{d}y \,\quad \text{ a.e } t_0 \in [0,T] \end{aligned}$$

This finishes the proof of Step 1, since ii’) easily follows from the fact that $q_N (\,\cdot \,,0) =0$ and $q_{N,t} (\,\cdot \,,0) =0$.

Proof of Step 2 We first prove the result in the case $q' \in L^2([0,T],H^1(\varOmega ))$ and $q'' \in L^2([0,T],L^2(\varOmega ))$. Let $(x_0,t_0) \in (B_R \setminus \varOmega ) \times [0,T]$ such that $dist(x_0,\partial \varOmega )> c_+ t_0$. There is $\epsilon _0>0$ such that for each $t \in [0,t_0]$, we have $B(x_0, (c_++ \epsilon _0) (t_0-t)) \cap \partial \varOmega =\emptyset$. We also denote $\mathcal {O}_t = B(x_0, c (t_0-t)) \cap B_R$ and

$$\begin{aligned} E(t)= & {} \frac{1}{2} \int _{\mathcal {O}_t} c^{-2}(x) |q_t(x,t)|^2\\&+\,|\nabla q(x,t)|^2 \mathrm{d}x, \quad 0 \le t \le t_0. \end{aligned}$$

Then,

$$\begin{aligned} \frac{\mathrm{d}}{\mathrm{d}t}E(t)= & {} -\frac{c_+}{2} \int _{\partial \mathcal {O}_t \setminus \partial B_R} c^{-2}(x) |q_t(x,t)|^2\\&+\,|\nabla q(x,t) |^2 d\sigma (x)\\&+\int _{\mathcal {O}_t} c^{-2}(x) q_t(x,t) \, q_{tt}(x,t) \\&+\,\nabla q(x,t) \nabla q_t(x,t) \, \mathrm{d}x. \end{aligned}$$

Taking integration by parts for the second integral gives the following formula of $\frac{\mathrm{d}}{\mathrm{d}t}E(t)$:

$$\begin{aligned}&-\frac{c_+}{2} \int _{\partial \mathcal {O}_t \setminus \partial \varOmega _R} \big [ c^{-2}(x) |q_t(x,t)|^2 + |\nabla q(x,t)|^2 - 2 \partial _\nu q(x,t)\\&\quad \times \frac{q_t(x,t)}{c_+} \big ] d\sigma (x) + \int _{\mathcal {O}_t} \big [c^{-2}(x) q_{tt}(x,t)\\&\quad - \varDelta q(x,t) \big ] q_t(x,t) \, \mathrm{d}x. \end{aligned}$$

Noting that the integrand of the first term on the right-hand side is nonnegative, we arrive to

$$\begin{aligned} \frac{\mathrm{d}}{\mathrm{d}t}E(t) \le \int _{\mathcal {O}_t} \big [c^{-2}(x) q_{tt}(x,t) - \varDelta q(x,t) \big ] q_t(x,t) \, \mathrm{d}x. \end{aligned}$$

Let us recall that for any function $\phi \in H_0^1(B_R)$

$$\begin{aligned}&\int _{B_R} c^{-2}(x) \, q_{tt}(x,t) \, \phi (x) \,\mathrm{d}x + \int _{B_R} a(x) \, q_{t}(x,t) \, \phi (x) \,\mathrm{d}x \\&\quad + \int _{B_R} \nabla q(x,t) \, \nabla \phi (x) \,\mathrm{d}x = -\int _{\partial \varOmega } \, g(y,t) \, \phi (x) \, \mathrm{d}y. \end{aligned}$$

For $0<\epsilon <\epsilon _0$, we choose $\varphi _\epsilon \in C^\infty (\mathbb R^d)$ be a nonnegative function such that $\varphi _\epsilon \equiv 1$ on $B_{(x_0, c_+(t_0 -t))}$ and $\varphi _\epsilon \equiv 0$ outside of $B_{(x_0,(c_+ + \epsilon ) (t_0-t))}$ and $\lim _{\epsilon \rightarrow 0} \varphi _\epsilon = \chi _{B_{x_0,c_+(t_0-t)}}$ on $L^2(\mathbb R^d)$. Choosing $\phi (x)= q_t(x,t) \varphi _\epsilon (x)$,⁴ we obtain

$$\begin{aligned}&\int _{B_R} c^{-2}(x) \, q_{tt}(x,t) \, q_t(x,t) \varphi _\epsilon (x) \,\mathrm{d}x\\&\quad + \int _{B_R} a(x) \, q_{t}(x,t) \, q_t(x,t) \varphi _\epsilon (x) \,\mathrm{d}x \\&\quad + \int _{B_R} \nabla q(x,t) \, \nabla [v_t(x,t) \varphi _\epsilon (x)] \,\mathrm{d}x =0. \end{aligned}$$

Taking integration by parts for the last integral and combine it with the first integral, we obtain

$$\begin{aligned}&\int _{B_R} \left[ c^{-2}(x) \, q_{tt}(x,t) - \varDelta q(x,t)\right] \, q_t(x,t) \varphi _\epsilon (x) \,\mathrm{d}x \\&\quad + \int _{B_R} a(x) \, q^2_{t}(x,t) \varphi _\epsilon (x) \,\mathrm{d}x =0. \end{aligned}$$

Therefore,

$$\begin{aligned} \int _{B_R} \left[ c^{-2}(x) \, q_{tt}(x,t) - \varDelta q(x,t)\right] \, q_t(x,t) \varphi _\epsilon (x) \,\mathrm{d}x \le 0. \end{aligned}$$

Taking the limit as $\epsilon \rightarrow 0$, we obtain

$$\begin{aligned} \int _{\mathcal {O}_t} \left[ c^{-2}(x) \, q_{tt}(x,t) - \varDelta q(x,t)\right] \, q_t(x,t) \,\mathrm{d}x \le 0. \end{aligned}$$

We obtain $\frac{E(t)}{\mathrm{d}t} \le 0.$ Noting that $E(0) =0$, we arrive at $E(t) =0$ for all $t \in [0,t_0]$. Therefore, $q(x,t) = 0$ on $\mathcal {O}_t$ for all $t \in [0,t_0]$. Since this is correct for all $(x_0,t_0) \in \varOmega ^c \times [0,T]$ such that $dist(x_0, \partial \varOmega ) > c_+ t_0$, It is now easy to see $q(x,t) = 0$ for all $(x,t) \in \varOmega ^c$ such that $dist(x,\partial \varOmega ) \ge c_+ t$.

In general, we do not have the required regularity for the above proof. However, consider $Q(x,t) = \int _0^t q(x,\tau ) d \tau$. Then, Q satisfies the same equation (with a different jump function) and the required regularity. The above proof then shows that $Q(x,t) =0$ for all $(x,t) \in \varOmega ^c \times [0,T]$ such that $dist(x,\partial \varOmega ) \ge c_+ t$. It implies the same result for q(x, t). This finishes proof of Step 2.

Finishing the Proof Now extending q into $\mathbb R^d \times [0,T]$ by zero on $(\mathbb R^d \setminus B_R) \times [0,T]$, we can easily prove that q is a weak solution on $\mathbb R^d \times [0,T]$. Moreover, q satisfies the finite speed of propagation (i). Finally, estimate (6) follows from (14). The uniqueness of q is simple (see, e.g., proof of Theorem A.2 in [8]), we leave the details to the reader.

1.2 A k-Space Method for the Damped Wave Equation

In this subsection, we briefly describe the k-space method as we use it to numerically compute the solution of the wave equation, which is required for evaluating the forward operator $\mathbf W$ and its adjoint $\mathbf W ^*$. For the case $a=0$, several methods for numerically solving the underlying acoustic wave equation have been used in PAT. This includes finite difference methods [10, 47, 57], finite element methods [8] as well as Fourier spectral and k-space methods [14, 29, 60]. We now extend the k-space method to the case $a \ne 0$ because this method does not suffer from numerical dispersion [13].

Consider the solution $p :\mathbb R^d \times (0, T) \rightarrow \mathbb R$ of the damped wave equation

$$\begin{aligned}&[ c^{-2} \, \partial _{tt} + a \, \partial _{t} - \varDelta ] p = s&\text{ on } \mathbb R^d \times (0,T), \end{aligned}$$

(15)

$$\begin{aligned}&p(0) = f&\text{ on } \mathbb R^d, \end{aligned}$$

(16)

$$\begin{aligned}&p_t(0) = - c^2 \, a \, f&\text{ on } \mathbb R^d. \end{aligned}$$

(17)

Here, $s :\mathbb R^d \times (0,T) \rightarrow \mathbb R$ is a given source term and $f :\mathbb R^d \rightarrow \mathbb R$ the given initial pressure. To derive the k-space method, one first rewrites (15) in the form

$$\begin{aligned} {[}\partial _{tt} - c_0^{2} \varDelta ] p = (1 - c_0^{2} / c^{2}) p_{tt}- c_0^{2} a \, p_{t} + c_0^{2} s \end{aligned}$$

(18)

where $c_0>0$ is a suitable constant; we take $c_0 =c_+ := \max \left\{ c(x) :x \in \mathbb R^2\right\}$.

The k-space method is derived from (18) by introducing the auxiliary functions v(x, t) and r(x, t) such that $v_{tt}(x,t) = (1 - c_0^{2} / c^{2}(x) ) p_{tt}(x,t)$ and $r_{tt}(x,t) = c_0^{2} a(x) p_{t}(x,t)$. Such an approach shows that (18) is equivalent to the following system of equations,

$$\begin{aligned} {[}\partial _{tt} - c_0^{2} \varDelta ] w&= c_0^2 \, s + c_0^2 \varDelta v - c_0^2 \, \varDelta r, \end{aligned}$$

(19)

$$\begin{aligned} v&= \left( c^2 / c_0^2 - 1 \right) \, (w - r) \end{aligned}$$

(20)

$$\begin{aligned} p&= v+ w - r \end{aligned}$$

(21)

$$\begin{aligned} r(t)&= c_0^2 a \int _{0}^t p(s) \mathrm{d}s. \end{aligned}$$

(22)

Interpreting $c_0^2 \varDelta v(x,t) - c_0^2 \, \varDelta r(x,t)$ as an additional source term, (19) is a standard wave equation with constant sound speed $c_0$. This suggests the time stepping formula

$$\begin{aligned}&w(x,t + h_t) =2 w(x,t) - w(x,t - h_t)\nonumber \\&\qquad - 4 \mathcal F_\xi ^{-1} \Bigl [ \sin (c_0 \vert \xi \vert h_t/2)^2 \nonumber \\&\quad \times \mathcal F_x [w(x,t) + v(x,t) - r(x,t) ] - (c_0h_t/2)^2 \nonumber \\&\quad \times {{\,\mathrm{sinc}\,}}(c_0 \vert \xi \vert h_t/2)^2 \mathcal F_x [s(x,t) ] \Bigr ], \end{aligned}$$

(23)

where $\mathcal F_x$ and $\mathcal F_\xi ^{-1}$ denote the Fourier and inverse Fourier transforms in the spatial variable x and the spatial frequency variable $\xi$, respectively, and $h_t > 0$ is a time stepping size.

The resulting k-space method for solving (15) is summarized in Algorithm 1.

Algorithm 1

(The k-space method) For given initial pressure f(x) and source term s(x, t) approximate the solution p(x, t) of (15) as follows:

• (1) Set $t = 0$ and define initial conditions

  • $r(x,0) = 0$;

  • $v(x,0) = (1-c_0^2 / c^2 (x)) f(x)$;

  • $w(x,0) = c_0^2 / c^2 (x) f(x)$;

  • $w(x,-h_t) = (1 + h_t c_0^2 a(x) ) w(x,0)$.

• (2) Compute $w(x,t + h_t)$ by evaluating (23);

• (3) Make the updates

  • $v(x,t + h_t) :=\left( c^2(x)/c_0^2- 1\right) \, ( w(x,t+h_t) - r(x,t) )$;

  • $p(x,t + h_t) :=v(x,t + h_t) + w(x,t + h_t) - r(x,t)$;

  • $r(x,t + h_t) :=r(x,t) + c_0^2 a(x) p(x,t + h_t) h_t$;

• (4) Set $t \leftarrow t+h_t$ and go back to (3).

Algorithm 1 can directly be used to evaluate the forward operator $\mathbf W f$ by taking $s(x,t) =0$ and restricting the solution to the measurement surface $S_R$, that is, $\mathbf W f = p|_{S_R \times (0, T)}$. Recall that the adjoint operator is given by $\mathbf W ^*g = q_t(0)$, where $q:\mathbb R^2 \times (0,T) \rightarrow \mathbb R$ satisfies the adjoint wave equation

$$\begin{aligned}&[c^{-2} \, \partial _{tt} - \varDelta ] q = - \delta _{S_R} \, g&\text{ on } \mathbb R^2 \times (0,T) \end{aligned}$$

(24)

$$\begin{aligned}&q_t(T) = q(T) = 0&\text{ on } \mathbb R^d. \end{aligned}$$

(25)

By substituting $t \leftarrow T- t$ and taking $s(x,t) = g(x,T-t) \, \delta _{S}(x)$ as source term in 15, Algorithm 1 can also be used to evaluate the $\mathbf W ^*$. In the partial data case where measurements are made on a subset $S \subsetneq S_R$ only, the adjoint can be implemented by taking the source $s(x,t) = \chi (x,t) \, g(x,T-t)\, \delta _{S_R}(x)$ with an appropriate window function $\chi (x,t)$. In order to use all available data, in our implementations we take the window function to be equal to one on the observation part S and zero outside. This choice of the window function is known to create streak artifacts into the picture [7, 21, 46]. However, as we see in our simulations, the artifacts fade away quickly after several iterations when the problem is well posed.