The quasi-reversibility method to numerically solve an inverse source problem for hyperbolic equations

We consider an inverse source problem arising from biomedical imaging based on the photo-acoustic and thermo-acoustic effects, which are named as photo-acoustic tomography (PAT) and thermo-acoustic tomograpy (TAT) respectively. In PAT, [1, 2], non-ionizing laser pulses are sent to a biological tissue under inspection (for instance, woman’s breast in mamography). A part of the energy will be absorbed and converted into heat, causing a thermal expansion and a subsequence ultrasonic wave propagating in space. The ultrasonic pressures on a surface around the tissue are measured. The experimental set up for TAT, [3], is similar to PAT except the use of microwave other than laser pulses. Finding some initial information of the pressures from these measurements yields structure inside this tissue.

Due to the important real-world applications, the inverse source problem PAT/TAT has been studied intensively. There are several methods to solve them available. In the case when the waves propagate in the free space, one can find explicit reconstruction formulas in [4, 5, 6, 7], the time reversal method [8, 9, 10, 11, 12], the quasi-reversibility method [13] and the iterative methods [14, 15, 16]. The publications above study PAT/TAT for simple models for non-damping and isotropic media. The reader can find publications about PAT/TAT for more complicated model involving a damping term or attenuation term [17, 18, 19, 20, 21, 22, 23, 24, 25]. The time reversibility requires an approximation of the wave at a “final stage” in the whole medium. This “internal data” might be known assuming that the reflection of the wave on the measured surface is negligible. However, there are many circumstances where this assumption is no longer true. For example, when the biological tissue under consideration is located inside glass, the waves reflects as in a resonant cavity [26, 27, 28]. In this case, measuring or approximating final stage of the wave inside the tissue is impossible. We draw the reader’s attention to [29] for a method to solve PAT/TAT in a bounded cavity with wave reflection. Our contribution in this paper is to apply the quasi-reversibility method to solve the inverse source problem of PAT/TAT for damping and nonhomogeneous media. In this case, the model is a full hyperbolic equation in a bounded domain. By this, the reflection of the waves at the measurement sites is allowed.

The uniqueness and the stability for the inverse source problem for general hyperbolic equations in the damping and inhomogeneous media can be proved by using a Carleman estimate. These important results are the extensions of the uniqueness and stability for the inverse problem for a simple hyperbolic equation in [13] in the non-damping case. However, since the arguments are very similar to the ones in that paper using Carleman estimate for hyperbolic operators, for the brevity, we do not write the proof here.

Instead of using the direct optimal control method, to find the initial value of general hyperbolic equations, we derive a system of elliptic partial differential equations, which are considered as an approximation model for our method. Solution of this system is the vector consisting of Fourier coefficients of the wave with respect to a special basis. This system and the given Cauchy boundary data can be solved by the quasi-reversibility method. The quasi-reversibility method was first proposed by Lattès and Lions [30] in 1969. Since then it has been studied intensively [31, 32, 33, 34, 13, 35, 36, 37]. The application of Carleman estimates for proofs of convergence of those minimizers was first proposed in [38] for Laplace’s equation. In particular, [39] is the first publication where it was proposed to use Carleman estimates to obtain Lipschitz stability of solutions of hyperbolic equations with lateral Cauchy data. We draw the reader’s attention to the paper [40] that represents a survey of the quasi-reversibility method. Using a Carleman estimate, we prove Lipschitz-like convergence rate of regularized solutions generated by the quasi-reversibility method to the true solution of that Cauchy problem.

It seems, in theory, that our method of approximation works for any orthonormal basis of $L^{2}$. However, this observation is not true in the numerical sense. That means, the special basis we use to establish the approximation model is crucial. The basis we use in this paper was first introduced by Klibanov in [48], called $\{\Psi_{n}\}_{n\geq 1}$. It has a very important property that $\Psi^{\prime}_{n}$ is not identically $0$ in an open interval while other bases; for e.g., trigonometric functions and orthonormal polynomials, do not. In this paper, we prove numerically that choosing this basis is optimal for our method.

As mentioned, we establish in this paper the Lipschitz convergence of the quasi-reversibility method. Our main contribution to this field is to relax a technical condition on the noise. In our previous works [41, 42, 43, 44, 45] and references therein, we assumed that the noise contained in the boundary data can be “smoothly extended” as a function defined on the domain $\Omega$. This condition implies that the noise must be smooth. Motivated by the fact that this assumption is not always true, we employ a Carleman estimate involving the boundary integrals to obtain the new convergence without imposing this “extension” condition of the noise function.

The paper is organized as follows. In Section 2, we state the inverse problems under consideration and derive an approximation model whose solution directly yields their solutions. In Section 3, we introduce some auxiliary results and prove the Carleman estimate, which plays an important role in our analysis. In Section 4, we implement the quasi-reversibility method to solve the system of elliptic equations and prove the convergence of the solution as the noise level tends to $0$. Section 5 is for the numerical studies. Section 6 is to provide some numerical results for $1D-$ problem and to show the significance of the used orthonormal basis. Finally, section 7 is for the concluding remarks.

Let $\Omega$ be a smooth and bounded domain in $\mathbb{R}^{d}$, where $d\geq 1$ is the spatial dimension, and $T$ be a positive number. Let $1\leq{\bf a}\in C^{1}(\Omega)$, and $a,{\bf c}\in L^{\infty}(\Omega)$ be functions defined on $\Omega$. Let ${\bf b}$ be a $d$-dimensional vector valued function in $L^{\infty}(\Omega,\mathbb{R}^{d}).$ Define the elliptic operator

for all $\phi\in C^{2}(\overline{\Omega})$. Let $p\in H^{2}_{0}(\Omega)$ represent a source, we consider the problems of solving $u({\bf x},t)\in H^{2}(\Omega\times(0,T))$ generated by the source $p({\bf x})\in C_{0}^{2}(\Omega)$ and subjected to either homogeneous Dirichlet or Neumann boundary condition, which are given by

or

respectively. Our interest is to determine the source $p({\bf x})$, ${\bf x}\in\Omega$, from some boundary observation of the wave. More precisely, the inverse source problems are formulated as:

The unique solvability of problems (2.2) and (2.3) can be obtained by Garlerkin approximations and energy estimates as in Chapter 7, Section 7.2 in [46]. We now focus on our approach for solving these two inverse problems.

Let $\{\Psi_{n}\}_{n\geq 1}$ be an orthonormal basis of $L^{2}(0,T).$ The function $u({\bf x},t)$ can be represented as:

where

Consider

for some cut-off number $N$. This number $N$ is chosen numerically such that ${u^{\small N}}$ well-approximates the function $u$, see Section 5.1 for more details. We have,

Plugging (2.8) and (2.9) into the first equation of problem (2.2), we get

for all $({\bf x},t)\in\Omega\times[0,T]$.

Since $\{\Psi_{n}\}_{n\geq 1}$ is an orthonormal basis of $L^{2}$, multiplying $\Psi_{m}(t)$ to both sides of (2.10) and then integrating the resulting equation with respect to $t$, for each $m\in\{1,\dots,N\},$ yields

where

Furthermore, from (2.7) we have

Therefore, the Cauchy data of $u_{n}$, for all $n=1,\ldots,N$ on the boundary $\partial\Omega$ are determined by:

1. 1.

   Regarding Problem 2.1

   $$\left\{\begin{array}[]{rcl}\partial_{\nu}u_{n}({\bf x})&=&\displaystyle\int_{0}^{T}f_{1}({\bf x},t)\Psi_{n}(t)dt,\\ u_{n}({\bf x})&=&0\end{array}\right.$$

   (2.12)

2. 2.

   Regarding Problem 2.2

   $$\left\{\begin{array}[]{rcl}u_{n}({\bf x})&=&\displaystyle\int_{0}^{T}f_{2}({\bf x},t)\Psi_{n}(t)dt,\\ \partial_{\nu}u_{n}({\bf x})&=&0.\end{array}\right.$$

   (2.13)

The system of elliptic partial differential equations (2.11) together with Cauchy data either (2.12) or (2.13) is our approximation model, see Remark 2.1. It allows to determine coefficients $u_{n}$, for all $n=1,\ldots,N$, and then the approximation $u^{N}({\bf x},t)$ of $u({\bf x},t)$. The source term will be given by $u^{N}({\bf x},0).$ In summary, the numerical method for solving Problem 2.1 and Problem 2.1 is described in Algorithm 1 and Algorithm 2 below respectively.

As mentioned in Remark 2.2 that solving (2.11) and (2.12) for Problem 2.1 and solving (2.11) and (2.13) for Problem 2.2 are interesting when the given data in (2.4) and (2.5) contain noise. We employ the quasi-reversibility method to do so. Let $\epsilon$ be a small positive number playing the role of the regularization parameter. To solve (2.11) and (2.12) we minimize the following mismatch functional

where ${\bf u}=(u_{1},\dots,u_{N})$ is in $H^{2}(\Omega)^{N}$ satisfying ${\bf u}({\bf x})=0$ for all ${\bf x}\in\partial\Omega$. To solve (2.11) and (2.13), we minimize the following mismatch functional

where ${\bf u}=(u_{1},\dots,u_{N})$ is in $H^{2}(\Omega)^{N}$ satisfying $\partial_{\nu}{\bf u}({\bf x})=0$ for all ${\bf x}\in\partial\Omega$.

In the following sections, we prove that $J_{1}$ has a unique minimizer and that minimizer converges to the true solution of (2.11) and (2.12) as the noise level tends to $0$. The results for $J_{2}$ can be obtained in the same manner. We introduce some auxiliary results in the next section.

Let $X$ be a number in $(0,1)$. We denote

Then, there exists $\alpha\in(0,1/2)$ such that the function

We have the following lemma.

Lemma 3.3 is a direct consequence of [49, Chapter 4, §1, Lemma 3] in which the function $\phi$ is independent of the time variable. Let $\tilde{R}$ be a positive number such that $\Omega\subset B(0,\tilde{R})$, where $B(0,\tilde{R})$ denotes a ball of the center at $0$ and the radius $\tilde{R}$. Let $p$ and $q$ be two positive numbers such that

where ${\bf e}_{1}$ is the unit direction vector of $x_{1}$ axis. For any ${\bf x}\in\Omega$, we define $\tilde{\bf x}:=p{\bf x}+q{\bf e}_{1}$, then (3.4) yields $\tilde{\bf x}\in\widetilde{\Omega}$. By modifying constant $C$ in Lemma 3.3 (using $Cp^{2}$ instead of $C$) we have that the Lemma 3.3 holds true in domain $\Omega$. From now on, we apply Lemma 3.3 for all function in the space $H^{2}({\Omega})$. The following result plays an important role in our analysis.

The Carleman estimate (3.5) plays a key role for us to estimate the error of the solution to the inverse problem assuming that the given data contains noise.

As mentioned in section 2, we only prove the convergence of the quasi-reversibility method to solve (2.11) with Cauchy boundary data (2.12). In this case, the objective functional $J_{1}$, now named as $J$, see (2.14), is written as

subject to the constraint $u_{1}=\ldots=u_{N}=0$ on $\partial\Omega$. Let

be a closed subspace of $H^{2}(\Omega)^{N}$. It is clear that $J$ is strictly convex in ${\bf H_{0}}$. We now prove that $J$ is has a unique minimizer in ${\bf H_{0}}$.

We now assume that the measured data contain noise, with a noise level $\delta$. We next study the convergence of ${\bf u}^{\min}$ as noise level $\delta$ tends to 0. Let denote $f_{1}^{\delta}({\bf x},t)$ the noisy data and $f_{1}^{*}({\bf x},t)$ the corresponding noiseless data, $({\bf x},t)\in\partial\Omega\times[0,T]$. By noise level $\delta$, we mean

Since the truncation number $N$ is a finite number, we can write

where $C$ is a generic constant depending only on $N,{\bf a},{\bf b}$ and $c$. For each $m\in\{1,\dots,N\}$, define

for all ${\bf x}\in\partial\Omega.$ The following theorem guarantees the Lipschitz stability of the reconstructed method with respect to noise.