# An optimal control problem in photoacoustic tomography

Maïtine Bergounioux, Xavier Bonnefond, Thomas Haberkorn, Yannick Privat

► To cite this version:

Maïtine Bergounioux, Xavier Bonnefond, Thomas Haberkorn, Yannick Privat. An optimal control problem in photoacoustic tomography. Mathematical Models and Methods in Applied Sciences, 2014, 24 (12), pp.2525-2548. 10.1142/S0218202514500286. hal-00833867v2

HAL Id: hal-00833867

https://hal.science/hal-00833867v2

Submitted on 1 Mar 2014

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
An optimal control problem in photoacoustic tomography*†

M. Bergounioux‡, X. Bonnefond‡, T. Haberkorn‡, and Y. Privat§

‡Univ. d'Orléans, Labo. MAPMO, CNRS, UMR 6628, Fédération Denis Poisson,
FR 2964, Bat. Math., BP 6759, 45067 Orléans cedex 2, France.

email: maitine.bergounioux,xavier.bonnefond,thomas.haberkorn@univ-orleans.fr

§CNRS, Université Pierre et Marie Curie (Univ. Paris 6), UMR 7598, Laboratoire
Jacques-Louis Lions, F-75005, Paris, France.

email: yannick.privat@upmc.fr

Abstract

This article is devoted to the introduction and study of a photoacoustic tomography model, an imaging technique based on the reconstruction of an internal photoacoustic source distribution from measurements acquired by scanning ultrasound detectors over a surface that encloses the body containing the source under study. In a nutshell, the inverse problem consists in determining absorption and diffusion coefficients in a system coupling a hyperbolic equation (acoustic pressure wave) with a parabolic equation (diffusion of the fluence rate), from boundary measurements of the photoacoustic pressure. Since such kinds of inverse problems are known to be generically ill-posed, we propose here an optimal control approach, introducing a penalized functional with a regularizing term in order to deal with such difficulties. The coefficients we want to recover stand for the control variable. We provide a mathematical analysis of this problem, showing that this approach makes sense. We finally write necessary first order optimality conditions and give preliminary numerical results.

**Keywords:** Photoacoustic tomography; inverse problem; optimal control.

**AMS classification:** 49J20, 35M33, 80A23, 93C20

# 1 Introduction

Photoacoustic imaging constitutes a cutting-edge technology that has drawn considerable attention in the medical imaging area. It uniquely combines the absorption contrast between two media with ultrasound high resolution. Moreover, it is non-ionizing and non-invasive, and is the fastest growing new biomedical method, with clinical applications on the way.

The main idea of the photoacoustic effect is simple. The tissue to be imaged is usually irradiated by a nanosecond-pulsed laser at a given optical wavelength. This energy is converted into heat. Absorption of light by molecules beneath the surface creates a

<footnote>*This work is supported by ANR (AVENTURES - ANR-12-BLAN-BS01-0001-01)</footnote>

<footnote>†The fourth author was partially supported by the ANR project OPTIFORM</footnote>
thermally induced pressure jump that propagates as a sound wave, which can be detected. By detecting the pressure waves, we can localize their heterogeneities (i.e., where light was absorbed) and obtain important information about the studied sample (See Figure 1 below).

<figure><img src="image_2.png" /><figcaption>Figure 1: Realization of a tomograph with integrating transducers, from Patch and Scherzer[35].</figcaption></figure>

This hybrid system uses an electromagnetic input and records ultrasound waves as an output. The electromagnetic energy is distributed at a given time through the object. The induced increase of temperature depends on the local absorption properties. For example, vascularized tissues, as early-stage tumors, absorb more energy. This opens the way to the detection of heterogeneities *via* measurements of the pressure field. Heterogeneities behave like internal acoustic sources, and the signals recorded by pressure detectors outside the medium under study provide information on their distribution.

Note that one speaks of thermo-acoustic tomography (TAT) when the heating is realized by means of microwaves (with wavelengths comparable to 1 m), and of photoacoustic tomography (PAT) when optical heating is used (high-frequency radiation near infrared with sub-µm wavelength). While in TAT waves of radio frequency range are used to trigger the ultrasound signal, in the PAT the frequency lies in the visual or near infrared ranges. In brief, TAT and PAT are two hybrid techniques using electromagnetic waves as an excitation (input) and acoustic waves as an observation (output).

Both techniques lead to an ill-posed inverse problem of the same form which, under simplifying assumptions, entails inversion (in the wide sense) of the spherical Radon transform.

In all generality, the study of such problems leads to a coupled system, constituted by one equation driving the behavior of the acoustic pressure and another one depending on the nature of the problem (PAT or TAT). In many works on the TAT problem, some physical approximations permit to rewrite the direct problem as a single partial differential equation (see the references mentioned below). More precisely, it writes: given the sound speed $c(x)$ and measured data $y_{\text{obs}}$ on $S \subset \mathbb{R}^n$ ($n = 2,3$), find the initial value $u_0(x)$ of
the pressure $y(t, x)$ where $y$ is the solution to the problem

$$ egin{cases} 
\frac{\partial^2 y}{\partial t^2}(t, x) - c^2(x) \Delta y(t, x) = 0, \quad (t, x) \in [0, T] \times \mathbb{R}^n, \\ 
y(0, x) = u_o(x), \quad x \in \mathbb{R}^n, \\ 
\frac{\partial y}{\partial t}(0, x) = 0, \quad x \in \mathbb{R}^n, \\ 
y(t, x) = y_{\text{obs}}(t, x), \quad x \in \mathcal{S}, t \in [0, T]. 
\end{cases} 
\quad (1.1) $$

We also mention [15] where a coupled system is introduced to model the TAT problem. The initial value $u_o$ is the TAT image. This problem is known to be highly ill-posed. In the sequel, we will propose a relevant model for the PAT problem.

In most reconstruction methods in PAT (or TAT), additional assumptions are performed such as conditions on the support of the function to be recovered and/or the observation surface, or a constant sound speed. Notice that a nice overview of the state of the art for the thermo-acoustic inverse problem has been done in [26].

One currently has a choice between three main types of reconstruction procedures for closed observation surfaces, namely the filtered backprojection formulae, eigenfunction expansion methods and time reversal methods.

* The filtered backprojection approach is the most popular [17, 22, 23, 24, 29, 41]. However, it is not clear that backprojection-type formulae could be written for any closed observation surface $S$. In [23], inversion formulae are provided assuming odd dimensions and constant sound speed. Indeed, in this case the Huygens' principle holds. Roughly speaking, it asserts that for any initial source with a compact support, the wave leaves any bounded domain in a finite time. This is no longer true if the spatial dimension is even and/or the sound speed is not constant. All known formulae of filtered backprojection type assume constant sound speed and thus are not available for acoustically inhomogeneous media. In addition, the only closed bounded surface $S$ for which such formulae are known is a sphere. Let us also mention [34] where a reconstruction algorithm in this vein (using the Radon transform) is proposed.

* Expansion series are useful in the case where the Huygens principle is valid. This approach was extended to the constant speed and arbitrary closed observation surface and modified by the use of the eigenfunctions of the Laplacian with Dirichlet conditions on $S$ [10]. It theoretically works for any closed surface and for variable sound speeds [38]. One can also refer to [28, 36].

* The time reversal method (see for example [25, 26]) can be used to approximate the initial pressure when the sound speed inside the object is variable. It works for arbitrary geometries of the closed observation surface $S$. Ammari et al. [3, 5] have performed sharp analysis of these problems both from the modeling and numerical point of view.

We also mention as possible additional techniques those based on finite elements discretization [42].
In this paper we propose to investigate the PAT model and the related inverse problem with an alternative formulation. We use an optimal control approach. Indeed, in our model the function to be recovered is the control function while the pressure is the state function which satisfies a wave equation.

This article is organized as follows: Section 2 is devoted to the description of the mathematical model driving the behaviors of the light transport and the wave pressure. It leads to the introduction of a coupled system of two partial differential equations, respectively of hyperbolic and parabolic type, in Section 2.3. The inverse problem mentioned above is interpreted as an optimal control problem in Section 3.1 and an existence result is provided. Section 3.2 is devoted to the differentiability analysis of the cost function and the computation of its gradient. Necessary first order optimality conditions are then derived in Section 3.3. We end the paper by giving leads for numerical experimentation.

In the sequel, if $E \subset \mathbb{R}^n L^s(E)$ is the space of measurable functions $f$ on $E$ such that $\int_E |f|^s < +\infty$, $s \in [1, +\infty)$, $L^\infty(E)$ is the space of essentially bounded functions on $E$ and $L^\infty(E, [a, b])$ is the set of essentially bounded functions with values in $[a, b]$ (where $a, b \in L^\infty(E)$). For $s, q \in \mathbb{N} \cup \{+\infty\}$, $W^{s,q}(E)$ is the Sobolev-space of $L^q$ function $s$ times differentiable in the distributional sense whose derivatives are in $L^q(E)$ (see [1, 13]). Note that $W^{0,2}(E) = L^2(E)$ and $W^{1,2}(E) = H^1(E)$. The space of functions in $H^1(E)$ whose trace vanishes on $\partial E$ is denoted $H_0^1(E)$. For $s, q \in \mathbb{N} \cup \{+\infty\}$, the set of $C^s$ functionals from $[0, T]$ to $E$ is denoted $C^s(0, T; E)$. The space of functions of bounded variation is denoted $BV(E)$ (see [13]).

# 2 The Photo Acoustic model

## 2.1 Light transport

In the PAT set-up, the tissues to investigate are illuminated with a laser source in the near infrared frequency range. As they propagate in the body, the particles are subject to absorption and diffusion, and are governed by the Boltzmann equation[11]. However, this equation requires the knowledge (in the direct problem) or the reconstruction (in the inverse problem) of a phase function representing the probability of scattering from one direction to another. Since this function depends on numerous factors, it is usually not known, nor is it easy to reconstruct.

Fortunately, the modeling becomes easier when dealing with soft, deep tissues. As a matter of fact, these latter are known to be highly diffusive media, so that the scattering tends to lose anisotropy as the particles go deeper into the tissue.

In this situation, the fluence rate $I$, that is the average of the luminous intensity in all the directions, satisfies the diffusion equation (see for example [11])

$$ \left\{ \begin{array}{ll} \displaystyle \frac{1}{\nu} \frac{\partial I}{\partial t}(t, x) + \mu_a(x) I(t, x) - \operatorname{div}(D \nabla I)(t, x) = S(t, x), & (t, x) \in [0, T] \times \Omega \\ \\ I(0, x) = 0, & x \in \Omega, \end{array} \right. \qquad (2.2) $$

where $\nu$ is the speed of light, $S$ is the incident light source, $\mu_a$ is the absorption coefficient, $D$ is the diffusion coefficient, and $T > 0$ is the duration of the acquisition process.
Here, $\Omega$ stands for the part of the body where the diffusion approximation is relevant. Usually, this domain does not include the tissues next to the surface of the body, since the photons first have a quasi ballistic behavior, which is not consistent with the diffusion approximation. Yet, the scatter sites of the early propagation act as isotropic sources for the diffusion equation (see Figure 2). As a consequence, when the incident light comes from a source point (optical fiber) located at the surface of the body, we can assume that the diffusion approximation holds in the entire body, provided that the source term $S$ in equation (2.2) is correctly chosen. Indeed, in the set up of Figure 2, if the real source $S_{\text{real}}$ has an amplitude $S_{\text{real}}^0$ and is located at $x = 0$, we can define $S$ with an amplitude $S^0 = aS_{\text{real}}^0$ and located at a depth of $\frac{1}{\mu_t}$, where $a$ and $\mu_t$ are, respectively, the transport albedo and the total interaction coefficient [21]. Although they are not precisely known, a model based on reasonable values of these two quantities (for the first layers of the skin) should be precise enough for our purpose.

<figure><img src="image_3.png" /><figcaption>Figure 2: The incident light is scattered at different sites along the $x$ axis.</figcaption></figure>

From now on, we make the assumption that the diffusion approximation holds in the entire body, which means that $\Omega$ stands now for the body. In order to complete this light transport model, we need to find appropriate boundary conditions for Equation (2.2).

Following [27], we start from the Robin condition

$$I(t, x) = AD \frac{\partial I}{\partial \nu}(t, x), \quad \text{a.e.}(t, x) \in [0, T] \times \partial\Omega,$$
where A is related to the internal reflection and can be deduced from the Fresnel reflection coefficients and $\nu$ denotes the outward pointing normal vector. This condition expresses that no photon current goes back into the body from the external medium, and can be satisfied with null Dirichlet conditions on an extrapolated boundary at $2AD$ from $\partial\Omega$ (as in Figure 2)[21]. We will still denote by $\Omega$ this enlarged set so that the boundary condition writes

$$I(t, x) = 0, \text{ a.e. } (t, x) \in [0, T] \times \partial\Omega,$$

The fluence $I$ is now assumed to vanish outside $\Omega$.

## 2.2 Pressure wave

Photoacoustic tomography is a mixed medical imaging technique, meaning that the out-coming signal is not of the same kind as the incoming energy. This is due to the thermo-acoustic effect [24, 37]: the incident light is absorbed by the tissue, and the resulting thermal expansion generates a pressure wave $p^0$ governed by the system

$$\begin{cases} \frac{\partial^2 p^0}{\partial t^2}(t, x) - \operatorname{div}(v_s^2 \nabla p^0)(t, x) = 1_{\Omega}(x)\Gamma(x)\mu_a(x)\frac{\partial I}{\partial t}(t, x), & (t, x) \in [0, T] \times \mathcal{B}, \\ p^0(0, x) = \frac{\partial p^0}{\partial t}(0, x) = 0, & x \in \mathcal{B} \end{cases}$$

where the notation $1_{\Omega}$ stands for the characteristic function of the domain $\Omega$, defined for almost every $x \in \mathcal{B}$ by

$$1_{\Omega}(x) = \begin{cases} 1 & \text{if } x \in \Omega \\ 0 & \text{otherwise.} \end{cases}$$

Here, the Grueneisen coefficient $\Gamma$, coupling the energy absorption to the thermal expansion, is assumed to be known. So is the speed of sound $v_s$, satisfying $v_s \in [v_s^{\min}, v_s^{\max}]$, with $v_s^{\min} > 0$.

The domain $\mathcal{B}$ is the place where the wave propagates. Obviously, it includes $\Omega$ and it has to be bounded in view of numerical simulations. The ball $\mathcal{B}$ is chosen large enough in such a way that $p^0$ vanishes on $\partial\mathcal{B}$ during the recording process. The size of $\mathcal{B}$ depends consequently on the location of the recording equipment and the duration $T$ of the acquisition. In other words, the reflected wave coming from $\partial\mathcal{B}$ doesn't have time to reach the acquisition equipment before time $T$.

**Remark 2.1.** Actually, the source term of the wave equation should be $\Gamma(x)\frac{\partial^2 T}{\partial t^2}(t, x)$, where $T$ is solution of some heat equation

$$\frac{\partial T}{\partial t}(t, x) - K\Delta T(t, x) = \mu_a(x)I(t, x),$$

(see [33] for example). Nevertheless, in most applications, the illumination is so short that no thermal conduction can occur during the thermoacoustic coupling, so that the constant $K$ can be considered to vanish, leading to Equation (2.2). Most of well known works on photoacoustic tomography make the further assumption that the illumination approaches a Dirac pulse, so that the thermoacoustic coupling becomes instantaneous and that the fluence rate $I$ does not depend on time anymore. Yet, it is possible to achieve
photoacoustic imaging with a continuous-wave laser illumination (see e.g. [31]), though there is a loss of resolution. This is the reason why we didn't make the assumption of a laser pulse illumination. The effect of the thermal conduction should be discussed in a future work.

Without any loss of generality, and without loss of information, it is more convenient to work with the new state p defined by

$$p(t, x) = \int_{0}^{t} p^{0}(s, x)ds.$$

According to Leibniz rule, this latter satisfies

$$\left\{
\begin{array}{ll}
\displaystyle \frac{\partial^2 p}{\partial t^2}(t, x) - \operatorname{div}(v_s^2 \nabla p)(t, x) = 1_{\Omega}(x)\Gamma(x)\mu_a(x)I(t, x), & (t, x) \in [0, T] \times \mathcal{B}, \\
\\
p(t, x) = 0, & (t, x) \in [0, T] \times \partial\mathcal{B}, \\
\\
p(0, x) = \frac{\partial p}{\partial t}(0, x) = 0, & x \in \mathcal{B}.
\end{array}
\right.$$

## 2.3 The direct problem

The effectiveness of photoacoustic tomography relies on the relation between inhomogeneities of the biological tissues and variations of the coefficients $\mu_a$ and $D$. Depending on the frequency range of the illumination (usually in the red or near infrared region), the gray level mapping of the absorptivity can achieve useful functional and structural imaging through, for instance, quantification of oxygen saturation or hemoglobin content [32, 40].

These considerations suggest to define $\mu := (\mu_a, D)$ as the **control variable** that we want to identify. Let $\mu_a^{\min} < \mu_a^{\max}$ and $D^{\min} < D^{\max}$ denote positive real numbers. The minimal (natural) assumptions on $\mu_a$ and $D$ are

$$\mu_a \in [\mu_a^{\min}, \mu_a^{\max}] \text{ and } D \in [D^{\min}, D^{\max}] \text{ a.e. in } \mathcal{B}, \qquad (2.3)$$

so that these maps lie in $L^{\infty}(\mathcal{B})$. We recall that $\Omega$ and $\mathcal{B}$ are two bounded open sets of $\mathbb{R}^d$ ($d \ge 2$), with $\mathcal{C}^1$-boundaries, satisfying $\Omega \subset \mathcal{B}$. The set $\Omega$ being the (extrapolated) body, we may assume that $\mu_a$ and $D$ are known on $\mathcal{B}\setminus\Omega$.

Introduce the set $\mathcal{Q}$ and its boundary $\Sigma$ defined by

$$\mathcal{Q} = (0, T) \times \Omega \text{ and } \Sigma = (0, T) \times \partial\Omega.$$

Since there are two variables to reconstruct, we might need at least two sets of data. This idea has been explored in a slightly different context in [14]. Following this work, we assume that the experiment is repeated with different light sources, denoted by $(S_k)_{1 \le k \le s}$ with $s \ge 2$ and each $S_k$ in $L^{\infty}(\mathcal{Q})$.

Provided that the frequency of the sources $S_k$ doesn't change, the coefficients $\mu_a$ and $D$ remain the same. However, the fluence rate $I$ and the acoustic signal $p$ may change with $k$. Then, we may define $I_k$ and $p_k$, for $k \in \{1, \dots, s\}$ as the solutions of the two **state equations**
$$ egin{cases} 
\frac{\partial^2 p_k}{\partial t^2}(t, x) - \operatorname{div}(v_s^2 \nabla p_k)(t, x) = \mathbb{1}_{\Omega}(x)\Gamma(x)\mu_a(x)I_k(t, x), & (t, x) \in (0, T) \times \mathcal{B}, \\ 
p_k(t, x) = 0, & (t, x) \in (0, T) \times \partial\mathcal{B}, \\ 
p_k(0, x) = \frac{\partial p_k}{\partial t}(0, x) = 0, & x \in \mathcal{B}, 
\end{cases} 
\quad (2.4) $$

and

$$ egin{cases} 
\frac{1}{\nu} \frac{\partial I_k}{\partial t}(t, x) + \mu_a(x) I_k(t, x) - \operatorname{div}(D \nabla I_k)(t, x) = S_k(t, x), & (t, x) \in (0, T) \times \Omega, \\ 
I_k(0, x) = 0, & x \in \Omega, \\ 
I_k(t, x) = 0 & x \in \mathcal{B} \setminus \Omega \\ 
I_k(t, x) = 0, & (t, x) \in \Sigma. 
\end{cases} 
\quad (2.5) $$

The photoacoustic tomography model is completely described by the coupling of equations (2.5) and (2.4), in which $I_k$ is extended to 0 on $\mathcal{B} \setminus \Omega$. We first mention that this system is well-posed, in other terms that (2.4)-(2.5) has a unique solution under standard assumptions. The following theorem is standard and its proof can be found for example in [20].

**Theorem 2.1.** Let $\Omega$ be a bounded connected open set of $\mathbb{R}^d$ with $C^1$ boundary, $\Gamma \in L^\infty(\mathcal{B})$, $v_s \in L^\infty(\mathcal{B}, [v_s^{\min}, v_s^{\max}])$. Assume that the assumptions (2.3) hold. Then,

i) Equation (2.5) has a unique solution $I_k$ such that

$$ 
\begin{aligned}
I_k \in C^0(0, T; L^2(\Omega)) \cap L^2(0, T; H_0^1(\Omega)), \\ 
\frac{\partial I_k}{\partial t} \in L^2(0, T; H^{-1}(\Omega)).
\end{aligned} 
$$

ii) Equation (2.4) has a unique solution $p_k$ such that

$$ p_k \in C(0, T; H_0^1(\mathcal{B})) \cap C^1(0, T; L^2(\mathcal{B})). $$

**Remark 2.2.** Even if they are reasonable in this setting, the assumptions made earlier on the variables $\mu_a$, $D$ and $v_s$ are not sharp, neither are the regularity results stated here. Nevertheless, our purpose does not require stronger statements.

The last step to complete the description of the direct model is the formalization of data acquisition. In PAT, ultrasonic transducers are placed in a neighborhood of the body and record the resulting pressure wave $p^0$ for all times in $[0, T]$. Let us denote by $\omega$ the set of the locations of these transducers, which can be either finite, discrete or (ideally) some hypersurface of $\mathbb{R}^d$. Assume for example that

$$ \omega = \bigcup_{i=1}^{N} \{x_i\}, $$

where each point $x_i$ belongs to $\mathcal{B} \setminus \Omega$. Unfortunately, this choice of acquisition set do not allow to apply some classic techniques of optimal control, such as Stokes' theorem. To
overcome this difficulty, we propose to thicken the set $\omega$ into a union of non empty open sets of $\mathbb{R}^d$. Namely, we replace in the sequel the set $\omega$ by the set $\omega_\varepsilon$ defined for $\varepsilon > 0$ by

$$ \omega_{\varepsilon} = \bigcup_{x \in \omega} B(x, \varepsilon), \qquad (2.6) $$

where $B(x, \varepsilon)$ denotes the open ball with radius $\varepsilon$ centered at $x$. It is illustrated on Figure 3.

<figure><img src="image_4.png" /><figcaption>Figure 3: The dots are $\omega$, the blue (light gray) balls are $\omega_\varepsilon$.</figcaption></figure>

We thus make the assumption that the pressure $p_k^0$ is known on $[0, T] \times \omega_\varepsilon$. Still defining the state variables $p_k$ as $\int_0^t p_k^0(t, x) dx$, the PAT data are given by

$$ \{p_k(t, x) | 1 \le k \le s, t \in [0, T], x \in \omega_\varepsilon\}. $$

Actually, we don't have access to such an information (we only record $p_k^0$ on $\omega$). Nevertheless, once the space discretization step is set to $\delta x$, $\varepsilon$ can be set to $\delta x/2$, so that the thickened data have the same discrete counterpart as the actual data.

Next Section is devoted to the sensitivity analysis of these state equations.

## 2.4 Sensitivity analysis

In this section, we will omit the indices $k$ when we refer to $p_k$ and $I_k$, the solutions of (2.4)-(2.5), for the sake of clarity. Define $U_{ad}$, the set of admissible controls $\mu = (\mu_a, D)$ as

$$ U_{ad} = \{(\mu_a, D) \in [L^\infty(\mathcal{B})]^2 | \mu_a \in [\mu_a^{\min}, \mu_a^{\max}] \text{ and } D \in [D^{\min}, D^{\max}] \text{ a.e. in } \mathcal{B}\}. \quad (2.7) $$

Using Theorem 2.1, we define the maps
$$ I: \begin{array}{rcl} \mathcal{U}_{ad} & \to & C^0(0, T; L^2(\Omega)) \cap L^2(0, T; H_0^1(\Omega)) \\ (\mu_a, D) & \mapsto & I(\mu_a, D) \end{array} \qquad (2.8) $$

where $I(\mu_a, D)$ satisfies (2.5) and

$$ p: \begin{array}{rcl} \mathcal{U}_{ad} & \to & C^0(0, T; H_0^1(\mathcal{B})) \\ (\mu_a, D) & \mapsto & p(\mu_a, D), \end{array} \qquad (2.9) $$

where $p(\mu_a, D)$ is the solution to (2.4).

The following theorem constitutes the main result of this section.

**Theorem 2.2.** The operator $p$ defined by (2.9) is weakly continuous with respect to $\mu$ and strongly continuous with respect to $D$ from $\mathcal{U}_{ad}$ to $L^2(\mathcal{B})$. More precisely, let $(\mu_a^n, D^n)_{n \in \mathbb{N}} \in \mathcal{U}_{ad}^N$ be such that $\mu_a^n$ weakly converges to $\mu_a^*$ in $L^2(\Omega)$ and $D^n$ strongly converges to $D^*$ in $L^2(\Omega)$ for some $(\mu_a^*, D^*)$ in $\mathcal{U}_{ad}$; then the sequence $(p^n)_{n \in \mathbb{N}}$ defined by $p^n = p(\mu_a^n, D^n)$ strongly converges up to a subsequence to $p^* = p(\mu_a^*, D^*)$ in $L^2(\mathcal{Q})$.

*Proof.* Let $(\mu_a^n, D^n)_{n \in \mathbb{N}} \in [\mathcal{U}_{ad}]^N$ and $(\mu_a^*, D^*)$ in $\mathcal{U}_{ad}$ be such that $\mu_a^n \to \mu_a^*$ (weakly) in $L^2(\Omega)$ and $D^n \to D^*$ (strongly) in $L^2(\Omega)$, as $n \to +\infty$.

**Step 1:** let us first prove that the sequence $(I^n)_{n \in \mathbb{N}}$ defined by $I^n = I(\mu_a^n, D^n)$ weakly converges in $L^2(0, T; H_0^1(\Omega))$ to some $I^* \in L^2(0, T; H_0^1(\Omega))$.

Indeed, note first that the weak formulation of System (2.5) writes: for every $\Psi \in L^2(0, T; H_0^1(\Omega))$,

$$ \int_{\mathcal{Q}} \frac{1}{\nu} \frac{\partial I^n}{\partial t} \Psi + D^n \nabla I^n \nabla \Psi + \mu_a^n I^n \Psi - S \Psi = 0. \qquad (2.10) $$

Taking $\Psi = I^n$ yields:

$$ \int_{\mathcal{Q}} (D^n |\nabla I^n|^2 + \mu_a^n |I^n|^2) \le \int_{\mathcal{Q}} S I^n, $$

since $I^n(0) = 0$. Using both Cauchy-Schwarz inequality and the fact that $D^n$ and $\mu_a^n$ lie in $\mathcal{U}_{ad}$ for every $n \in \mathbb{N}$ yields

$$ \min(\mu_a^{\min}, D^{\min}) \int_{0}^{T} \|I^n\|_{H^1(\Omega)}^2 \le \|S\|_{L^2(\mathcal{Q})} \sqrt{\int_{0}^{T} \|I^n\|_{H^1(\Omega)}^2}, \qquad (2.11) $$

It thus follows that the sequence $(I^n)_{n \in \mathbb{N}}$ is uniformly bounded in $L^2(0, T; H_0^1(\Omega))$ and weakly converges (up to a subsequence) to some $I^*$ in $L^2(0, T; H_0^1(\Omega))$. This convergence is actually strong in $L^2(\mathcal{Q})$ by virtue of the Aubin-Lions lemma (see[30] Theorem 5.1 p. 58).

**Step 2:** Let us now prove that $I^* = I(\mu_a^*, D^*)$. Notice first that the combination of inequalities (2.10) and (2.11) ensures that $\frac{\partial I^n}{\partial t}$ is uniformly bounded in $L^2(0, T; H^{-1}(\Omega))$. As a consequence, $\frac{\partial I^*}{\partial t}$ belongs to $L^2(0, T; H^{-1}(\Omega))$ and $I^*$ is thus the limit of a subsequence of $(I^n)_{n \in \mathbb{N}}$ in $C^0(0, T; L^2(\Omega))$. It ensures that $I^*(0)$ is well-defined and vanishes.
With a slight notational abuse, we still denote by $(I^n)_{n \in \mathbb{N}}$ the subsequence introduced above. Using a pointwise version of equation (2.10), we claim that for almost every time $t \in [0, T]$, one has

$$ \forall \Psi \in H_0^1(\Omega), \quad \int_{\Omega} \frac{1}{\nu} \frac{\partial I^n}{\partial t}(t, \cdot) \Psi + D^n \nabla I^n(t, \cdot) \nabla \Psi + \mu_a^n I^n(t, \cdot) \Psi = \int_{\Omega} S(t, \cdot) \Psi. $$

The weak convergence of $\frac{\partial I^n}{\partial t}$ to $\frac{\partial I^*}{\partial t}$ yields

$$ \int_{\Omega} \frac{\partial I^n}{\partial t}(t, \cdot) \Psi \to \int_{\Omega} \frac{\partial I^*}{\partial t}(t, \cdot) \Psi \quad \text{as } n \to +\infty. $$

Moreover, using the strong convergence of $(I^n(t, \cdot))_{n \in \mathbb{N}}$ to $I^*(t, \cdot)$ in $L^2(\Omega)$ for almost every $t \in [0, T]$, one gets

$$ \int_{\Omega} \mu_a^n I^n(t, \cdot) \Psi \to \int_{\Omega} \mu_a^* I^*(t, \cdot) \Psi \quad \text{as } n \to +\infty. $$

Note that

$$ \langle (D^* - D^n) \nabla I^*(t, \cdot), \nabla \Psi \rangle_{L^2(\Omega)} \to 0 \quad \text{as } n \to +\infty. $$

and

$$ \langle D^n \nabla (I^* - I^n)(t, \cdot), \nabla \Psi \rangle_{L^2(\Omega)} \to 0 \quad \text{as } n \to +\infty, $$

since the sequence $(D^n)_{n \in \mathbb{N}}$ strongly converges to $D^*$ as $n \to +\infty$.

Combining these results with the following decomposition

$$ \langle D^* \nabla I^*(t, \cdot) - D^n \nabla I^n(t, \cdot), \nabla \Psi \rangle_{L^2(\Omega)} = \langle D^n \nabla (I^* - I^n)(t, \cdot), \nabla \Psi \rangle_{L^2(\Omega)} \\ + \langle (D^* - D^n) \nabla I^*(t, \cdot), \nabla \Psi \rangle_{L^2(\Omega)} . $$

shows that $I^* = I(\mu_a^*, D^*)$.

**Step 3:** let us finally show that the sequence $(p^n)_{n \in \mathbb{N}}$ defined by $p^n = p(\mu_a^n, D^n)$ strongly converges, up to a subsequence, to $p^* = p(\mu_a^*, D^*)$ in $L^2(\mathcal{B})$ as $n \to +\infty$. The source term in (2.4) is $h_n = 1_{\Omega} \Gamma \mu_a^n I^n$. For every $\Psi \in L^2([0, T] \times \mathcal{B})$ we get

$$ \langle h_n - h^*, \Psi \rangle_{L^2([0,T] \times \mathcal{B})} = \langle h_n - h^*, \Psi \rangle_{L^2(\mathcal{Q})} \\ = \langle I^n - I^*, \Gamma \mu_a^n \Psi \rangle_{L^2(\mathcal{Q})} + \langle \mu_a^n - \mu_a^*, \Gamma I^* \Psi \rangle_{L^2(\mathcal{Q})}, $$

with $h^* = 1_{\Omega} \Gamma \mu_a^n I^*$. Since $(I^n)_{n \in \mathbb{N}}$ strongly converges to $I^*$ in $L^2(\mathcal{Q})$, the sequence $(h_n)_{n \in \mathbb{N}}$ weakly converges to $h^*$ in $L^2([0, T] \times \mathcal{B})$ (recall that $I^n$ and $I^*$ vanish on $\mathcal{B} \setminus \Omega$). Since $p^n$ is solution to the wave equation (2.4), it implies that the sequence $(p^n)_{n \in \mathbb{N}}$ strongly converges up to a subsequence to $p^*$ in $L^2([0, T] \times \mathcal{B})$. $\square$

**Remark 2.3.** It would be interesting to investigate the robustness of the optimal control approach with respect to some uncertainty of the coefficients of the diffusion and acoustic equations. Theorem 2.2 provides yet some hints, and the issue of modeling the uncertainties and investigate their influence on the optimal control solution will be the object of a future study.
# 3 The inverse problem

## 3.1 Formulation as an optimal control problem

Let us now consider the open set $\omega_{\epsilon}$ defined by (2.6), where we will measure the outgoing pressures $p_k^{\text{obs}}$ at every time. In the sequel, $\epsilon$ is fixed and we rather write $\omega$ instead of $\omega_{\epsilon}$, for the sake of simplicity. The asymptotic behavior of the solutions as $\epsilon \to 0$ will be investigated in a forthcoming work.

Physically, it seems relevant to minimize a least square functional with respect to $\mu_a$ and $D$. We choose to add a penalization term in order to ensure the existence of an optimal control.

Let us define the functional $J$ by

$$J(\boldsymbol{\mu}) = \mathcal{F}(\boldsymbol{\mu}) + f(\boldsymbol{\mu}), \qquad (3.12)$$

for every $\boldsymbol{\mu} = (\mu_a, D) \in \mathcal{U}_{ad}$, where $f(\boldsymbol{\mu})$ stands for a regularizing term and $\mathcal{F}$ is a least square functional with respect to the measured pressure data. We set

$$F_k(\boldsymbol{\mu}) = \frac{1}{2} \int_{[0,T] \times \omega} (p_k(t,x) - p_k^{\text{obs}}(t,x))^2 dx dt$$

where $p_k^{\text{obs}}$ is the measured pressure (observed state) on $\omega$ when the source signal is $S_k$. Fix $\alpha \ge 0$ and $\beta \ge 0$. Assuming that we perform $s$ experiments, we define

$$\mathcal{F}(\boldsymbol{\mu}) = \sum_{k=1}^{s} F_k(\boldsymbol{\mu}) = \sum_{k=1}^{s} \frac{1}{2} \int_{[0,T] \times \omega} (p_k(t,x) - p_k^{\text{obs}}(t,x))^2 dx dt$$

and

$$f(\boldsymbol{\mu}) = \begin{cases} \alpha \int_{\Omega} (B\mu_a)^2(x)dx + \beta TV(D) & \text{if } D \in \text{BV}(\Omega) \\ +\infty & \text{otherwise.} \end{cases}$$

Here $BV(\Omega)$ denotes the space of functions of bounded variation[2, 13], $TV(D)$ stands for the total variation of $D$, and $B: L^2(\Omega) \to L^2(\Omega)$ is an invertible linear operator.

**Remark 3.1.** The operator $B$ is usually the $L^2(\Omega)$ identity operator. However one can decide to focus on specific frequencies of $\mu_a$ and $B$ can be chosen as a pass-band filter. Following[16], $B$ can be chosen as a mollifier operator for example.

The choice of the total variation as a regularization term is a technical choice that fits the physical meaning. Indeed, strong $L^2$ convergence of the $D$ part of minimizing sequences is needed to use Theorem 2.2 and obtain an existence result. The TV term seems to be the weakest one that provides such a convergence while respecting the physical requirements since discontinuities (and contours) are preserved.

The original inverse problem to perform parameter identification can be viewed as the following optimal control problem

$$(\mathcal{P}) \quad \min_{\boldsymbol{\mu} \in \mathcal{U}_{ad}} J(\boldsymbol{\mu}),$$

where the admissible set $\mathcal{U}_{ad}$ is defined by (2.7).

Notice that the values of the coefficients $\mu_a$ and $D$ on $B\setminus\Omega$ are already known and that $\mathcal{U}_{ad}$ is a closed convex subset of $L^2(B) \times L^2(B)$.
**Theorem 3.1** (Existence of an optimal control $\boldsymbol{\mu}$). Assume that $\alpha \ge 0$ and $\beta > 0$. Then, Problem (P) has at least a solution $\bar{\boldsymbol{\mu}} = (\bar{\mu}_a, \bar{D})$.

*Proof.* Let $(\mu_a^n, D^n)_{n \in \mathbb{N}}$ be a minimizing sequence. Since $(\mu_a^n)_{n \in \mathbb{N}}$ is bounded in $L^\infty(\Omega)$ (and in $L^2(\Omega)$), it weakly converges (up to a subsequence) to some $\bar{\mu}_a$ in $L^2(\Omega)$ as $n \to +\infty$. The sequence $(D^n)_{n \in \mathbb{N}}$ is bounded in $L^\infty$ as well, so it is bounded in $L^1$. From the boundedness of $(TV(D^n))_{n \in \mathbb{N}}$, we deduce that $(D^n)$ is bounded in $BV$ and weakly converges up to a subsequence to some $\bar{D}$ in $BV(\Omega)$ as $n \to +\infty$. The space $BV(\Omega)$ being compactly embedded in $L^1(\Omega)$, the sequence $(D^n)_{n \in \mathbb{N}}$ strongly converges to $\bar{D}$ for the $L^1$-topology. Since the sequence $(D^n)_{n \in \mathbb{N}}$ is uniformly bounded in $L^\infty(\Omega)$, we get the strong convergence of $(D^n)_{n \in \mathbb{N}}$ to $\bar{D}$ for the $L^2$-topology.

Using Theorem 2.2, we conclude that the sequence $(p_k^n)_{n \in \mathbb{N}}$ defined by $p_k^n = p_k(\mu_a^n, D^n)$ strongly converges (up to a subsequence) to $\bar{p}_k = p_k(\bar{\mu}_a, \bar{D})$ in $L^2([0, T] \times \mathcal{B})$, for every $k \in \{1, \dots, s\}$ as $n \to +\infty$. The lower semicontinuity of every $F_k$ with respect to the $L^2$-convergence and the lower semicontinuity of $f$ with respect to the $L^1$-convergence implies that the pair $(\bar{\mu}_a, \bar{D})$ is a solution of Problem (P). □

**Remark 3.2.** We are not able to prove uniqueness by now. As already mentioned, it seems necessary to get more that one data set, that is $s \ge 2$ [14]. Moreover, we will have to assume $\alpha > 0$. Notice that, in [4, 7], the authors get a unique solution by reducing the model.

## 3.2 Computation of the cost functional derivative

In order to write the necessary first order optimality conditions for Problem (P), we first compute the derivative of $\mathcal{F}$ with respect to the control variable $\boldsymbol{\mu} = (\mu_a, D)$. Since $\mathcal{F} = \sum_{k=1}^s F_k$, it suffices to compute the derivative of $F_k$. For the sake of clarity and readability, we will omit the index $k$ in the sequel.

In order to write the optimality conditions in the most simple way, let us notice that $L^\infty(\mathcal{B}) \subset L^2(\mathcal{B})$ so that we can endow $\mathcal{U}_{ad}$ with the usual hilbertian structure of $L^2(\mathcal{B})$.

Let $\boldsymbol{\mu} \in \mathcal{U}_{ad}$ and $\xi = (\xi_{\mu_a}, \xi_D) \in L^2(\Omega) \times L^2(\Omega)$ be an admissible perturbation of $\boldsymbol{\mu}$. In the sequel, if $\boldsymbol{\mu} \in \mathcal{U}_{ad} \mapsto g(\boldsymbol{\mu})$ is a Gâteaux-differentiable functional at $\boldsymbol{\mu}$ in direction $\xi$, we will indifferently denote by $\langle dg(\boldsymbol{\mu}), \xi \rangle$ or $\dot{g}(\boldsymbol{\mu})$ the Gâteaux derivative of $g$ at $\boldsymbol{\mu}$ in direction $\xi$, that is

$$ \dot{g}(\boldsymbol{\mu}) = \langle dg(\boldsymbol{\mu}), \xi \rangle = \lim_{t \searrow 0} \frac{g(\boldsymbol{\mu} + t\xi) - g(\boldsymbol{\mu})}{t}. $$

A calculus of variation standard analysis permits to show, applying shrewdly the implicit function theorem, that the functional $F$ is differentiable at $\boldsymbol{\mu}$ in direction $\xi$. Its derivative writes

$$ \langle dF(\boldsymbol{\mu}), \xi \rangle = \int_{[0,T] \times \omega} (p(t,x) - p^{\text{obs}}(t,x)) \dot{p}(t,x) dx dt, \quad (3.13) $$
where $\dot{p}$ is the solution of the system

$$ egin{cases} 
\frac{\partial^2 \dot{p}}{\partial t^2} - \operatorname{div}(v_s^2 \nabla \dot{p}) = 1_{\Omega} \Gamma \xi_{\mu_a} I + 1_{\Omega} \Gamma \mu_a \dot{I} & \text{in } [0, T] \times \mathcal{B} \\ 
\dot{p}(0, \cdot) = \frac{\partial \dot{p}}{\partial t}(0, \cdot) = 0 & \text{in } \mathcal{B} \\ 
\dot{p} = 0 & \text{on } [0, T] \times \partial \mathcal{B} 
\end{cases} 
\quad (3.14) $$

and $I$ is solution of the following system

$$ egin{cases} 
\frac{1}{\nu} \frac{\partial \dot{I}}{\partial t} + \mu_a \dot{I} + \xi_{\mu_a} I - \operatorname{div}(D \nabla \dot{I}) - \operatorname{div}(\xi_D \nabla I) = 0 & \text{in } \mathcal{Q} \\ 
\dot{I}(0, \cdot) = 0 & \text{in } \Omega \\ 
\dot{I} = 0 & \text{on } \Sigma 
\end{cases} 
\quad (3.15) $$

Since the expression (3.13) does not permit to express the first order optimality conditions easily, it is convenient to introduce some adjoint states to rewrite this derivative into a more workable expression. For that purpose, let us define $q_1$ and $q_2$ as the respective solutions of the systems

$$ egin{cases} 
\frac{\partial^2 q_1}{\partial t^2} - \operatorname{div}(v_s^2 \nabla q_1) = (p - p^{\text{obs}})1_{\omega_\epsilon} & \text{in } [0, T] \times \mathcal{B} \\ 
q_1(T, \cdot) = \frac{\partial q_1}{\partial t}(T, \cdot) = 0 & \text{in } \mathcal{B} \\ 
q_1 = 0 & \text{on } [0, T] \times \partial \mathcal{B} 
\end{cases} 
\quad (3.16) $$

and

$$ egin{cases} 
-\frac{1}{\nu} \frac{\partial q_2}{\partial t} + \mu_a q_2 - \operatorname{div}(D \nabla q_2) = \Gamma \mu_a q_1 & \text{in } \mathcal{Q} \\ 
q_2(T, \cdot) = 0 & \text{on } \Omega \\ 
q_2 = 0 & \text{on } \Sigma. 
\end{cases} 
\quad (3.17) $$

It is standard that under the assumptions of Theorem 2.1, System (3.17) has a unique solution

$$ q_2 \in C^0(0, T; L^2(\Omega)) \cap L^2(0, T; H_0^1(\Omega)) $$

and System (3.16) has a unique solution

$$ q_1 \in C^0(0, T; H_0^1(\mathcal{B})) \cap C^1(0, T; L^2(\mathcal{B})). $$

Let us now compute the derivative of $F$ at $\mu$ in the direction $\xi$.

**Proposition 3.1.** For every $\xi = (\xi_a, \xi_D) \in L^2(\Omega) \times L^2(\Omega)$, the functional $F$ is Gâteaux-differentiable at $\mu = (\mu_a, D)$ in the direction $\xi$ and

$$ egin{align} 
\langle dF(\boldsymbol{\mu}), \xi \rangle_{L^2(\Omega)} &= \int_{\Omega} \nabla F(\boldsymbol{\mu})(x) \xi(x) dx \nonumber \\ 
&= \int_{\Omega} \left( \frac{\partial F}{\partial \mu_a}(\mu_a, D)(x) \xi_{\mu_a} + \frac{\partial F}{\partial D}(\mu_a, D)(x) \xi_D(x) \right) dx \tag{3.18} 
\end{align} $$
where

$$ \nabla F(\boldsymbol{\mu}) = \left( \frac{\partial F}{\partial \mu_a}(\boldsymbol{\mu}), \frac{\partial F}{\partial D}(\boldsymbol{\mu}) \right) = \left( \int_0^T (1_{\Omega} \Gamma q_1 - q_2) I, - \int_0^T \nabla q_2 \cdot \nabla I \right). $$

*Proof.* Using integration by parts and Green's formula, one gets

$$ \begin{aligned} \langle dF(\mu_a), \xi \rangle_{L^2(\Omega)} &= \int_{[0,T] \times \omega} (p - p^{\text{obs}}) \dot{p} \\ &= \int_{[0,T] \times \mathcal{B}} \dot{p} \left( \frac{\partial^2 q_1}{\partial t^2} - \text{div}(v_s^2 \nabla q_1) \right) \\ &= \int_{[0,T] \times \mathcal{B}} q_1 \left( \frac{\partial^2 \dot{p}}{\partial t^2} - \text{div}(v_s^2 \nabla \dot{p}) \right) \\ &= \int_{[0,T] \times \mathcal{B}} q_1 (1_{\Omega} \Gamma \xi_{\mu_a} I + 1_{\Omega} \Gamma \mu_a \dot{I}) \\ &= \int_{\mathcal{Q}} \Gamma q_1 \xi_{\mu_a} I + \int_{\mathcal{Q}} \Gamma q_1 \mu_a \dot{I}. \end{aligned} $$

Let us now rewrite the term $\int_{\mathcal{Q}} \Gamma q_1 \mu_a \dot{I}$. One has

$$ \begin{aligned} \int_{\mathcal{Q}} \Gamma q_1 \mu_a \dot{I} &= \int_{\mathcal{Q}} \left( -\frac{1}{\nu} \frac{\partial q_2}{\partial t} + \mu_a q_2 - \text{div}(D \nabla q_2) \right) \dot{I} \\ &= \int_{\mathcal{Q}} \left( \frac{1}{\nu} \frac{\partial \dot{I}}{\partial t} + \mu_a \dot{I} - \text{div}(D \nabla \dot{I}) \right) q_2 \\ &= \int_{\mathcal{Q}} (-\xi_{\mu_a} I + \text{div}(\xi_D \nabla I)) q_2 \end{aligned} $$

We finally get

$$ \int_{\mathcal{Q}} \Gamma \mu_a q_1 \dot{I} = \int_{\mathcal{Q}} (-\xi_{\mu_a} q_2 I - \xi_D \nabla q_2 \cdot \nabla I) \quad (3.19) $$

so that

$$ \langle dF(\mu_a), \xi \rangle_{L^2(\Omega)} = \int_{\Omega} ((\Gamma q_1 - q_2)I\xi_{\mu_a} - \nabla q_2 \cdot \nabla I\xi_D) . $$

□

We deduce from the previous result the following expression of the derivative of $\mathcal{F}$.

**Theorem 3.2.** For every $\xi = (\xi_a, \xi_D) \in L^2(\Omega) \times L^2(\Omega)$, the functional $\mathcal{F}$ is Gâteaux-differentiable at $\boldsymbol{\mu} = (\mu_a, D)$ in the direction $\xi$ and

$$ \langle d\mathcal{F}(\boldsymbol{\mu}), \xi \rangle_{L^2(\Omega) \times L^2(\Omega)} = \int_{\Omega} \nabla \mathcal{F}(\boldsymbol{\mu}) \cdot \xi $$

where

$$ \nabla F(\boldsymbol{\mu}) = \sum_{k=1}^{s} \left( \int_{0}^{T} (1_{\Omega} \Gamma q_{1}^{k} - q_{2}^{k}) I_{k}, - \int_{0}^{T} \nabla q_{2}^{k} \cdot \nabla I_{k} \right), $$
and, for every $k \in \{1, \dots, s\}$,

$$ egin{cases} 
\frac{\partial^2 q_1^k}{\partial t^2} - \operatorname{div}(v_s^2 \nabla q_1^k) = (p_k - p_k^{\text{obs}})1_{\omega_e} & \text{in } [0, T] \times \mathcal{B} \\ 
q_1^k(T, \cdot) = \frac{\partial q_1^k}{\partial t}(T, \cdot) = 0 & \text{in } \mathcal{B} \\ 
q_1^k = 0 & \text{on } [0, T] \times \partial \mathcal{B} 
\end{cases} 
\quad (3.20) $$

$$ egin{cases} 
-\frac{1}{\nu} \frac{\partial q_2^k}{\partial t} + \mu_a q_2^k - \operatorname{div}(D \nabla q_2^k) = \Gamma \mu_a q_1^k & \text{in } \mathcal{Q} \\ 
q_2^k(T, \cdot) = 0 & \text{on } \Omega \\ 
q_2^k = 0 & \text{on } \Sigma 
\end{cases} 
\quad (3.21) $$

### 3.3 First order optimality conditions for Problem (P)

Assume that $\bar{\mu} = (\bar{\mu}_a, \bar{D})$ is an optimal solution to problem (P). Introduce the so called indicator function of the set $\mathcal{U}_{ad}$, denoted $\iota_{\mathcal{U}_{ad}}$ and defined by

$$ \iota_{\mathcal{U}_{ad}}(x) = \begin{cases} 0 & \text{if } x \in \mathcal{U}_{ad} \\ +\infty & \text{otherwise.} \end{cases} $$

The regularization function $f$ is not Gâteaux differentiable because of the Total Variation term. However the subdifferential $\partial TV$ is well known[18] and we get

$$ 0 \in \partial TV(\boldsymbol{\mu}) \iff \boldsymbol{\mu} \in \partial TV^*(0), $$

where the total variation conjugate functional $TV^*$ is the indicator function $\iota_{\bar{K}}$ of $\bar{K}$ with

$$ K = \{ \operatorname{div} \varphi \mid \varphi \in C_c^1(\Omega, \mathbb{R}^2), \|\varphi\|_{\infty} \le 1 \}. $$

This gives useful algorithms to compute the total variation subgradients (see [18, 39] for example).

Writing (P) as

$$ \min_{\boldsymbol{\mu} \in [L^\infty(\Omega)]^2} \mathcal{F}(\boldsymbol{\mu}) + f(\boldsymbol{\mu}) + \iota_{\mathcal{U}_{ad}}(\boldsymbol{\mu}), $$

the classical optimality condition reads

$$ 0 \in \partial(\mathcal{F}(\bar{\boldsymbol{\mu}}) + f(\bar{\boldsymbol{\mu}}) + \iota_{\mathcal{U}_{ad}}(\bar{\boldsymbol{\mu}})). $$

Using standard computational rules[19] and decoupling the first order optimality conditions on $\mu_a$ and $D$ yields:

i) Equation on $\mu_a$. For every $\mu_a \in L^\infty(\Omega)$ such that $\mu_a \in [\mu_a^{\min}, \mu_a^{\max}]$,

$$ \left\langle \frac{\partial \mathcal{F}}{\partial \mu_a}(\bar{\mu}_a, \bar{D}) + 2\alpha B^* B \bar{\mu}_a, \mu_a - \bar{\mu}_a \right\rangle_{L^2(\Omega)} \ge 0, \quad (3.22) $$
ii) Equation on D.

$$ -\frac{\partial \mathcal{F}}{\partial D}(\bar{\mu}_a, \bar{D}) \in \partial \mathcal{TV}(\bar{D}) + \partial \iota_{[D^{\min}, D^{\max}]}, \qquad (3.23) $$

where $B^*$ is the $L^2$-adjoint operator of $B$.

With the previous computations, equation (3.22) writes

$$ \forall \mu_a \in L^\infty(\Omega), \text{ s.t. } \mu_a \in [\mu_a^{\min}, \mu_a^{\max}], \quad \left\langle \sum_{k=1}^{s} (1_\Omega \Gamma q_1^k - q_2^k) I^k + 2\alpha B^* B \bar{\mu}_a, \mu_a - \bar{\mu}_a \right\rangle_{L^2(\Omega)} \ge 0, \qquad (3.24) $$

while equation (3.23) becomes

$$ \exists \delta^* \in \partial \mathcal{TV}(\bar{D}), \quad \forall D \in L^\infty(\Omega) \text{ s.t. } D \in [D^{\min}, D^{\max}], \quad \left\langle \sum_{k=1}^{s} \nabla I^k \cdot \nabla q_2^k - \delta^*, D - \bar{D} \right\rangle_{L^2(\Omega)} \ge 0. \qquad (3.25) $$

The following theorem summarizes these optimality conditions.

**Theorem 3.3.** Assume $\bar{\mu} = (\bar{\mu}_a, \bar{D})$ is an optimal solution to Problem (P). Then, there exists $q_1^k, q_2^k, k = 1, \dots, s$ and $\delta^* \in \partial \mathcal{TV}(\bar{D})$ such that

* The 2s state equations (2.4) for the pression and (2.5) for the fluence are satisfied (with s sources $S_k, k \in \{1, \dots, s\}$)

* The 2s adjoint state equations (3.20)-(3.21) are satisfied by $q_1^k$ and $q_2^k$ respectively, for $k \in \{1, \dots, s\}$.

* Inequations (3.24) and (3.25) hold.

**Remark 3.3.** In the very case where $D$ is constant and/or known, we are only interested in $\mu_a$. The (reduced) optimality system writes then: 2s state equations (2.4) and (2.5), 2s adjoint state equations (3.20)-(3.21) and (3.24).

# 4 Numerical experiments

The approach we use leads to an optimality system that can be solved numerically. However, the solving of this coupled optimal control problem raises some issues like the non differentiability of the BV regularization and the difference of speed scale (sound versus light) of the two equations. There are many numerical difficulties that we do not detail by now.

* One of them is related to the non differentiability of the BV regularization term. One can use dedicated algorithms as Chambolle's projection algorithm[18] for example. However, we need to study precisely how such a projection algorithm can be combined with the solving of the adjoint system. Another way could be to simply use a differentiable approximation of the total variation as $\mathcal{TV}(D) \approx \sqrt{\|\nabla D\|^2 + \varepsilon}$, where $\varepsilon > 0$ is small.
* The second numerical delicate issue is the coupled state (and adjoint state) system solving. The full coupled system has different time-scales. Indeed, the wave equation has a characteristic time related to the speed of sound while the one of the fluence is related to the speed of light. As discussed before, if the light source is a short pulse (with respect to the speed of sound), then we can modify the fluence equation by replacing it with the stationary one (dropping the $\partial_t$ term). This would then replace the source term in the wave equation by a Dirac. Such a modification would only simplify the derivation of the necessary condition as well as the numerical implementation. On the other hand, when considering a long lasting light source, we need to solve the 2 concurrent unstationary equations. Since this doesn't solve the time-scale difference problem, an idea is to use an asymptotic development to capture the main features of the RTE solution and plug it into the wave equation. This implies that we have to consider appropriate methods (multiscale, asymptotic methods or adaptive meshes) to perform competitive numerical computations. In addition, this has to be compared with existing methods.

Those interesting issues will be addressed in a forthcoming paper.

To illustrate the control approach and show that it is a relevant alternative method to the classical ones (that we mentioned in the introduction), we briefly present numerical experiments to compute a simple TAT model. Shortly speaking, we assume that the fluence equation is not useful any longer and consider equation (1.1) as a good model for TAT (as usual in TAT papers). More precisely we want to recover the source $u$ which drives the following equation

$$ 
\left\{ 
\begin{array}{l} 
\left( \frac{\partial^2 p}{\partial t^2} - \operatorname{div}(v^2 \nabla p) \right) (t, x) = 0, \quad (t, x) \in [0, T] \times \mathcal{B}, \\ 
p(t, x) = 0, \quad (t, x) \in [0, T] \times \partial \mathcal{B}, \\ 
p(0, x) = u, \quad \frac{\partial p}{\partial t}(0, x) = 0, \quad x \in \mathcal{B} 
\end{array} 
\right. 
\qquad (4.26) 
$$

from measurements $p_{obs}$ on the boundary of $\Omega$ where $\Omega$ is the 0-centered ball of radius $1/\sqrt{2}$ (so as to be outside the square containing the phantom). We consider the case where the observation surface is not closed (half a sphere) and the sound speed $v$ is not constant. We consider a 2D problem. Though we should deal with the 3D problem, the 2D one is still interesting, since it covers the case where detectors are lineic[34]. With the previous notations, we assume that $k = 1$ (only one source), $B$ is a mollifier[16] and $\beta = 0$ (total variation not included). The control problem writes

$$ (P_{\varepsilon}) \quad \left\{ 
\begin{array}{l} 
\min \frac{1}{2} \|p(u) - p_{obs}\|_{L^2([0,T]\times\omega)}^2 + \frac{\alpha}{2} \|Bu\|_{L^2([0,T]\times\Omega)}^2 \\ 
u \in L^2([0,T] \times \mathcal{B}), 
\end{array} 
\right. 
\qquad (4.27) $$

where $p(u)$ is the solution to (4.26) and $u$ is supported in $\Omega \subset \mathcal{B}$. This uncoupled system gives rise to the same kind of optimality system as the coupled system, except that there is only a slight modification of adjoint equation (3.20) and necessary condition (3.22) to consider. Moreover, we do not introduce bounds on $u$ so no projection has to be performed.

The tests have been done using Shepp-Logan phantom. As our purpose is to illustrate the relevance of our approach we do not focus on code and/or optimization methods so
that we do not report CPU time for example. The known speed of sound is supposed to be 1 outside the Shepp-Logan phantom and in [0.95, 1.05] inside. This choice of variation of speed represents the real variations between soft tissues and water (where the body to be reconstructed would be submerged).

Other methods than the optimal control approach could be used to solve this reconstruction problem. Most of those methods have originally been devised so as to deal with constant sound speed and a closed domain of observation and later on have been adapted to less stringent assumptions. For example, the time-reversal method[36, 38] has been adapted to variable speed of sound and open observation domain. A method based on the eigenvalues of the Laplacian has been extended to deal with open observation domain[29].

We chose to solve the optimality system by means of the conjugate gradient algorithm. The forward and backward problems are solved by means of a leapfrog discretization scheme on a staggered grid. In order to avoid handling large grids (due to the size of $B$), we use an appropriate PML (Perfectly Matched Layer) technique[12].

All the computations are performed on a standard computer using the *Scilab* software. We use the algorithm on the 512 by 512 pixels Shepp-Logan phantom, given on Figure 4.

<figure><img src="image_5.png" /><figcaption>Figure 4: The Shepp-Logan phantom. Left: only the phantom (512×512 pixels). Right: phantom and detectors distribution (number of detectors will vary).</figcaption></figure>

We investigate 5 different cases

* data are not corrupted or a white gaussian noise (SNR=0.15) is added to the simulated data;

* the number of detectors is 500, 50 or 10 (only for noiseless data) with a uniform angular sampling on the right half circle (as shown on Figure 4).

All the results correspond to a zero initialization of $u$ and 10 iterations (20 for 10 detectors) of the conjugate gradient. When tackling noiseless data, we use $\alpha = 0.1$ while we use $\alpha = 0.4$ for noisy data. The data ($p_{obs}$) are simulated with the known speed of sound
while the forward and backward problems use a noisy speed with an added white gaussian noise (SNR=0.02 so as to be less than the amplitude of the speed variation).

Figure 5 shows the results of the solving for uncorrupted data with 500, 50 and 10 detectors. For the noiseless data, we see that the reconstruction enables the recovery of

<figure><img src="image_6.png" /><figcaption>Figure 5: Noiseless data. Left: 500 detectors. Center: 50 detectors. Right: 10 detectors</figcaption></figure>

all the features of the Shepp-Logan phantom when using 500 or 50 detectors. A close inspection of the heavily dense (500) detectors population shows that the right side of the phantom is better reconstructed than the left side, which is to be expected since the detectors are in the right half plane. On the 50 detectors solving, we see some wave like artifacts inside the phantom but also in the outside, coming from the sparsely distributed detectors. The over pessimistic result with 10 detectors is only displayed because it clearly shows that the properly reconstructed boundaries are only the ones tangential to circles originating from the detectors.

Figure 6 shows the results for noisy data with 500 and 50 detectors. For the noisy

<figure><img src="image_7.png" /><figcaption>Figure 6: Noisy data. Left: 500 detectors. Right: 50 detectors.</figcaption></figure>

data, both reconstructions are slightly blurrier than for the noiseless data, as expected
since the regularization term is stronger. Even with sparsely distributed detectors, we can
still recover the different elements of the phantom.

5 Conclusion

We have presented a new model for PAT phenomena involving a coupled system. The optimal control approach we use seems promising in this context. Techniques we use are classical but they provide flexibility to study the inverse problem, especially in the case where the direct problem is not easy to handle ( involving many nonlinear equations for example). This leads to optimality systems that can be solved numerically, though there remain challenging numerical issues.

Many open questions remain, as the uniqueness of the solution of the optimal control problem for example. Identifiability issues have to be addressed as well. On the other hand, we have to investigate precisely the behavior of the solutions when the observation set $\omega_e$ measure reduces to 0: this case corresponds to pointwise sensors. At last, we go on developing precise models taking into account physical phenomena that are usually neglected, especially in the TAT context where the direct problem can be described both by a pressure equation and a Maxwell equation.

References

[1] R.A. Adams and J. Fournier, *Sobolev spaces*. Pure and Applied Mathematics Series, 140, Academic Press (2003).

[2] L. Ambrosio, N. Fusco and D. Pallara, *Functions of bounded variation and free discontinuity problems*. Oxford mathematical monographs, Oxford University Press (2000).

[3] H. Ammari, E. Bossy, V. Jugnon and H. Kang, *Mathematical modeling in photoacoustic imaging of small absorbers*. SIAM Rev. **52** (2010), no. 4, 677-695

[4] H. Ammari, E. Bossy, V. Jugnon and H. Kang, *Reconstruction of the optical absorption coefficient of a small absorber from the absorbed energy density*. SIAM J. Appl. Math. **71** (2011), no. 3, 676-693.

[5] H. Ammari, E. Bretin, V. Jugnon and A. Wahab, *Photoacoustic imaging for attenuating acoustic media*. Chapter in a Lecture Notes in Mathematics Volume, Springer-Verlag, 2011

[6] H. Ammari, E. Bretin, J. Garnier and V. Jugnon, *Coherent interferometry alogorithms for photoacousting imaging*. SIAM J. Numeri Anal. **50**, (2012), no. 5, 2259-2280.

[7] H. Ammari, H. Kang and S. Kim, *Sharp estimates for the Neumann functions and applications to quantitative photo-acoustic imaging in inhomogeneous media*. J. Differential Equations **253** (2012), no. 1, 41-72.
[8] H. Ammari, J. Garnier, Josselin, W. Jing and L. Hoang Nguyen, *Quantitative thermo-acoustic imaging: an exact reconstruction formula*. J. Differential Equations **254** (2013), no. 3, 1375-1395.

[9] H. Ammari, E. Bretin, J. Garnier and A. Wahab, *Time reversal in attenuating acoustic media*. Contemporary Mathematics, Vol. 548, (2011), 151-164

[10] M.A Anastasio, J.Zhang, D. Modgil and P.J La Rivière, *Application of inverse source concepts to photoacoustic tomography*, Inverse Problems **23** (2007), no. 6, 21-35

[11] Arridge, S. R., *Optical tomography in medical imaging*, Inverse Problems **15** (1999), no. 2, R41

[12] J. P. Berenger, *A perfectly matched layer for the absorption of electromagnetic waves*, J. Comput. Phys., pp. 185-200, Oct. 1994.

[13] H. Attouch, G. Buttazzo and G. Michaille, *Variational analysis in Sobolev and BV spaces : applications to PDEs and optimization*. MPS-SIAM series on optimization, (2006)

[14] G. Bal and K. Ren, *Multi-source quantitative photoacoustic tomography in a diffusive regime*, Inverse Problems, **27** (2011), 075003, 20 p.

[15] G. Bal, K. Ren, G. Uhlmann, and T. Zhou, *Quantitative thermo-acoustics and related problems*, Inverse Problems, **27** (2011), 055007, 15 p.

[16] X. Bonnefond and P. Maréchal, *A variational approach to the inversion of some compact operators*. Pacific Journal of Optimization, 5(1) (2009), 97-110

[17] P. Burgholzer, J. Bauer-Marschallinger, H. Grün, M. Haltmeier and G. Paltauf, *Temporal back-projection algorithms for photoacoustic tomography with integrating line detectors*, Inverse Problems **23** (2007), no. 6, 65-80

[18] A. Chambolle, *An algorithm for total variation minimization and applications*, Journal of Mathematical Imaging and Vision **20** (2004), 89-97

[19] I. Ekeland and R. Temam, *Convex Analysis and Variational problems*, SIAM Classic in Applied Mathematics, **28** (1999)

[20] L.C. Evans, *Partial Differential Equations*, Graduate Studies in Mathematics, **19**, AMS (1998).

[21] T.J. Farrell, M.S. Patterson and B. Wilson, *A diffusion theory model of spatially resolved, steady-state, diffuse reflectance for the noninvasive determination of tissue optical properties in vivo*, Med. Phys., **19** (1992), 879-888.

[22] D. Finch, M. Haltmeier and R. Rakesh, *Inversion of spherical means and the wave equation in even dimensions*, SIAM J. Appl. Math. **68**, (2007), no. 2, 392-412

[23] D. Finch, S. Patch and R. Rakesh, *Determining a function from its mean values over a family of spheres*, SIAM J. Math. Anal. **35**, (2004), 1213-40
[24] M. Haltmeier, T. Schuster and O. Scherzer, *Filtered backprojection for thermoacoustic computed tomography in spherical geometry*, Math. Meth. Appl. Sci. **28**, (2005), no. 16, 1919-1937

[25] Y. Hristova, P. Kuchment and L.Nguyen, *Reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and inhomogeneous media*, Inverse Problems, **24**, (2008), no.5, 055006, 25 p.

[26] Y. Hristova, *Time reversal in thermoacoustic tomography: an error estimate*, Inverse Problems, **25**, (2009), no. 5, 055008, 14 p.

[27] M. Keizer, W.M. Star and P.R. Storchi, *Optical diffusion in layered media*, App. Opt., **27** (1988), 1820-1824.

[28] K.P. Köstliand and P.C.Beard, *Two-dimensional photoacoustic imaging by use of Fourier-transform image reconstruction and a detector with an anisotropic response*, App. Opt., **42** (2003), pp. 1899-1908.

[29] L. A. Kunyansky, *Thermoacoustic tomography with detectors on an open curve: an efficient reconstruction algorithm*, Inverse Problems, **24**, (2008), no. 5, 055021, 18 p.

[30] J.L. Lions, *Quelques méthodes de résolution des problèmes aux limites non linéaires*, Gauthier-Villars, Paris, 1969

[31] K. Maslov and L.V. Wang, *Photoacoustic imaging of biological tissue with intensity modulated continuous-wave laser*, Journal of Biomedical Optics, **13**, (2008), no. 2, 024006, 5 p.

[32] J. Mobley and T. Vo-Dinh, in *Biomedical Photonic Handbook*, edited by T. Vo-Dinh, CRC Press (2003), Chap. 2.

[33] P.M. Morse and M. Feshbach, *Theoretical acoustics*, Reading, MA, Hardwood Academic (1991).

[34] G. Paltauf, R. Nuster, M. Haltmeier and P. Burgholzer, *Experimental evaluation of reconstruction algorithms for limited view photoacoustic tomography with line detectors*, Inverse Problems **23**, (2007), no. 6, S81

[35] S. K. Patch and O.Scherzer, *Photo- and Thermo-Acoustic Imaging Introduction*, Inverse Problems, **23**, (2007), no. 6, 1-10

[36] J.Qian, P. Stefanov, G. Uhlmann and H. Zhao, *An efficient Neumann-series based algorithm for Thermoacoustic and Photoacoustic Tomography with variable sound speed*, SIAM J. Imaging Sci. 4 (2011), no. 3, 850-883

[37] O. Scherzer, M. Grasmair, H. Grossauer, M. Haltmeier and F. Lenzen, *Variational Methods in Imaging*, Springer, 2008

[38] P. Stefanov and G.Uhlmann, *Thermoacoustic tomography with variable sound speed*, Inverse Problems 25 (2009), no. 7, 075011, 16 p.
[39] P. Weiss, L. Blanc-Féraud and G. Aubert, *Efficient schemes for total variation minimization under constraints in image processing*, SIAM J. Sci. Comput., 31(3), (2009), 2047-2080

[40] M. Xu and L.V. Wang, *Photoacoustic imaging in biomedicine*, Rev. Sci. Instrum., **77** (2006), no 4

[41] M. Xu and L.V. Wang, *Universal back-projection algorithm for photoacoustic computed tomography*, Physical Review E, **71** (2005), 1-7

[42] Z. Yuan and H. Jiang, *Three-dimensional finite-element-based photoacoustic tomography: Re-construction algorithm and simulations*, Medical Physics, **34** (2007), no. 538