An inverse obstacle problem for the wave equation in a finite time domain

We consider an inverse obstacle problem for the acoustic transient wave equation. More precisely, we wish to reconstruct an obstacle characterized by a Dirichlet boundary condition from lateral Cauchy data given on a subpart of the boundary of the domain and over a finite interval of time. We first give a proof of uniqueness for that problem and then propose an “exterior approach” based on a mixed formulation of quasi-reversibility and a level set method in order to actually solve the problem. Some 2D numerical experiments are provided to show that our approach is effective.

1. Introduction. We address an inverse obstacle problem for the wave equation, defined as follows. Let G be an open, bounded and connected domain of R d , d ≥ 2, with Lipschitz boundary. Throughout the paper, we define an obstacle as an open domain O G which is formed by a collection of a finite number of disjoint connected Lipschitz domains (O is not necessarily connected), and such that Ω := G \ O is connected. A generic point in Q := Ω × (0, T ), for T > 0, will be denoted (x, t). Let Γ be a non-empty open subset of ∂G. The notations are illustrated on Figure 1. Given a pair of data (g 0 , g 1 ) on Γ × (0, T ), the inverse obstacle problem consists in finding a domain O (independent of time t) and some function u(x, t) in the space (1) H 1,1 (Q) := {u ∈ L 2 (0, T ; H 1 (Ω)), u ∈ H 1 (0, T ; L 2 (Ω))}, following the notation of [19], and such that where ν is the outward unit normal to Ω. The inverse problem that we address is a geometric inverse problem, the obstacle to retrieve being characterized by a Dirichlet boundary condition. We also emphasize that the data are restricted to  (2) arises in the following practical situation. The medium G is at rest for negative time t. At initial time t = 0, an experimenter generates a known pulse g 1 on the subpart Γ of the surface of the medium G and measures the response g 0 of the medium on Γ during the time interval (0, T ), the complementary part of the surface being inaccessible. The goal is to retrieve the unknown obstacle O from those measurements (g 0 , g 1 ) on Γ × (0, T ).
To the best of our knowledge, there are very few papers dealing with the effective reconstruction of the obstacle: the only one we have found is [11], in which an optimization technique is applied. However, the obstacles which are seeked in the numerical experiments of [11] are a priori known to be circles in 2D, which are characterized by only three real parameters. We mention that in [18] the authors also want to retrieve a surface defined by a Dirichlet boundary condition with the help of a single incident wave, however with measurements in an infinite time interval. Needless to say that there are many contributions in the case of several incident waves, in particular in the case of infinitely many. In this vein, let us mention [1], which is based on the topological gradient, [22], which is based on the boundary control method, or [8], which relies on the Linear Sampling Method.
A crucial issue induced by considering measurements in a finite time domain (0, T ) is uniqueness. More precisely, the natural question is: what is the minimal value of T so that the obstacle O is uniquely determined by the Cauchy data (g 0 , g 1 ) on Γ × (0, T )? From our understanding this question is still open. In [15], however, the author provides some minimal values of T which guarantee uniqueness, as soon as a bound on a certain length characterizing the obstacle is assumed. Considering that our geometric setting is slightly different from the one studied in [15] and also for the sake of self-containment, we give a detailed uniqueness proof in 2D. However the main concern of this paper is the effective reconstruction of the obstacle in 2D from the data when uniqueness holds. Our strategy is to apply the so-called "exterior approach", first introduced in [2] in the case of the Laplace equation and then extended in [3,4,5] to the case of the Stokes system, heat equation and Helmholtz equation, respectively. We remark that an inverse problem such as (2) is both ill-posed and non-linear. The idea consists in addressing these two issues separately by coupling a quasi-reversibility method and a level set method. An initial guess of the obstacle being given, the exterior solution is updated from the lateral Cauchy data by solving a quasi-reversibility problem. An initial guess of the solution being given, the obstacle is updated by solving a level set problem. Our approach is hence iterative, with one step of quasi-reversibility and one step of level set at each iteration, until some stopping criterion is reached. Quasi-reversibility goes back to [17]: it is a Tikhonov regularization of some linear ill-posed PDE problem which is directly in the form of a weak formulation. This weak formulation can hence be discretized with the help of a Finite Element Method. The main drawback of the original method presented in [17] and later in [16] is the fact that the weak formulation associated with an ill-posed second-order problem corresponds to a fourth-order well-posed problem. As a consequence, the justification of such quasireversibility method requires some additional regularity for the exact data and the discretization requires some additional regularity for the approximation space. This is why we developed in [6] an alternative mixed formulation of quasi-reversibility which preserves the order of the original problem. Such mixed formulation has been recently generalized in [5] and we propose to apply it once again in the case of the wave equation with lateral Cauchy data. Concerning the level set method, we reuse the method introduced in [2], which is based on a simple Poisson equation and has shown its efficiency in many situations ever since. It should be noted that, due to the specific form of the uniqueness results in the case of the wave equation in a finite time domain, both the quasi-reversibility method and the level set method are more difficult to justify than in the previous situations of elliptic or parabolic equations.
This paper is organized as follows. The next section shows in which sense all the boundary conditions in problem (2) have to be understood. Section 3 is dedicated to some uniqueness results. Firstly it concerns a well-known unique continuation result which is crucial in what follows and the proof of which is postponed in an Appendix, for the sake of self-containment. Secondly it establishes a uniqueness result for the inverse problem (2) in dimension 2. Section 4 describes the mixed quasi-reversibility method while Section 5 describes the level set method. By merging these two methods we derive the "exterior approach" algorithm, which is detailed in Section 6. Lastly, some numerical results which illustrate the feasibility of this approach to solve problem (2) are presented in Section 7.
2. Some comments on boundary conditions. It is important to note that, due to the absence of boundary conditions on the boundary ∂G \ Γ, the standard regularity results for the wave equation are useless. This explains why, in order to justify that the boundary conditions in problem (2) are meaningful, we assumed that the solution u belongs to the unsual space H 1,1 (Q) defined by (1). Let us give a meaning to the traces (u, ∂ ν u) of u on Σ := Γ × (0, T ) and the traces (u, ∂ t u) of u on S 0 = Ω × {0}. It is readily seen that H 1,1 (Q) coincides with H 1 (Q), since Q can be viewed as a bounded Lipschitz domain of R d+1 . As a result, the trace of u ∈ V on ∂Q is well-defined in H 1/2 (∂Q), and in particular the traces u| Σ and u| S0 are well-defined in H 1/2 (Σ) and H 1/2 (S 0 ) as the sets of restrictions of functions in H 1/2 (∂Q) to Γ × (0, T ) and Ω × {0}, respectively. Furthermore, let us denote, for u ∈ H 1 (Q), v := (∇u, −∂ t u) ∈ (L 2 (Q)) d+1 . If u solves the wave equation We conclude that v ∈ H div (Q): this enables one to define v · ν d+1 | ∂Q ∈ H −1/2 (∂Q), which is the dual space of H 1/2 (∂Q). Here ν d+1 denotes the outward unit normal to the (d + 1)-dimensional domain Q. As a result, the normal derivatives ∂ ν u| Σ and ∂ t u| S0 are well defined in H −1/2 (Σ) and H −1/2 (S 0 ) as the restrictions of distributions in H −1/2 (∂Q) to Γ × (0, T ) and Ω × {0}, respectively. Let us recall that H −1/2 (Σ) and H −1/2 (S 0 ) are the dual spaces ofH 1/2 (Σ) andH 1/2 (S 0 ), the latter being defined as the subsets of functions in H 1/2 (Σ) and H 1/2 (S 0 ) which once extended by 0 on the whole boundary ∂Q are still in H 1/2 (∂Q).

3.1.
A basic unique continuation result. In this section, we recall a useful unique continuation result for the wave equation. To this aim, we need to define the geodesic distance in Ω.
If Ω ⊂ R d is a connected open domain, for x, y ∈ Ω, the geodesic distance between x and y is defined by where g is a continuous path in Ω of length (g). Here, the length of g is defined as where | · | stands for the Euclidian distance and the sup is taken over all decompositions of [0, 1] into an arbitrary (finite) number of intervals.
Now, let us introduce the constant D(Ω, Γ) denoting the largest geodesic distance between some point x in Ω and Γ, which will play a crucial role in the sequel. It is not difficult to see that we may have D(Ω, x 0 ) = +∞, even if Ω is a bounded domain. However we have D(Ω, x 0 ) < +∞ as soon as Ω is a bounded Lipschitz domain, which is a consequence of the fact that for x 0 ∈ Ω, the application x → d Ω (x, x 0 ) has a continuous extension in Ω for the Euclidian topology [23]. In this Lipschitz case, the definition of D(Ω, x 0 ) can hence be extended to the case when x 0 ∈ ∂Ω. We clearly have We conclude that for a bounded Lipschitz domain, D(Ω, Γ) < +∞. Now let us state the following classical theorem, which concerns unique continuation for the wave equation in the presence of lateral Cauchy data, and which is proved in Appendix for the sake of self-containment. Assume that u ∈ H 1 (Q) satisfies the system Then u vanishes in the subdomain Q 0 of Q defined by In particular, if T > D(Ω, Γ), then u vanishes in the subdomain Ω×(0, T −D(Ω, Γ)).   Proof. Given x ∈ ∂O, let us consider some y ∈ ∂O such that

Uniqueness in
Such point y exists because the functionỹ → d ∂O (x,ỹ) is continuous on the compact set ∂O. In dimension 2, there are only two paths g 1 and g 2 joining x and y on ∂O. It happens that (g 1 ) = (g 2 ). Indeed, if we had for example (g 1 ) > (g 2 ), then d ∂O (x, y) = (g 2 ). There would exist someỹ ∈ ∂O such that the two pathsg 1 and g 2 joining x toỹ satisfy and this would contradict the fact that y maximizesỹ → d ∂O (x,ỹ) on ∂O. Since (g 1 ) = (g 2 ), we have P (O) = (g 1 ) + (g 2 ) = 2 (g 2 ) = 2 sup y∈∂O d ∂O (x,ỹ).
We see that supỹ ∈∂O d ∂O (x,ỹ) is independent of x and then P (O) = 2D(∂O).
We will also need the following lemma, which is proved in [7] (see theorem IX.17 and Remark 20).
Let us now state our uniqueness result for problem (2) in two dimensions: it establishes a minimal time T which guarantees uniqueness given an imposed bound P on the perimeter of the unknown obstacle.
Theorem 3.6. For d = 2, let us consider two obstacles O 1 , O 2 such as defined in the introduction and corresponding functions u 1 , u 2 which satisfy (2) with data (g 0 , g 1 ) on Γ. We assume in addition that u i ∈ L 2 (0, T ; C 0 (Ω i )), i = 1, 2, and that for all x 0 ∈ Γ and all sufficiently small t 0 > 0, the function g 0 (x 0 , ·) is not identically zero in the interval (0, t 0 ). Let us denote by P an upper bound of P (O 1 ) and P (O 2 ) and D = D(G, Γ). If we assume that Proof. Let us define Ω 12 as the connected component of Let us consider some point x ∈ O 12 . For any ε > 0 there exists some The path g cuts the boundary of O 12 at some point x 12 such that the restrictiong of g joining x 12 to x 0 belongs to Ω 12 . Without loss of generality, the point x 12 belongs to the boundary of O 1 , so that x 12 is renamed Now let us consider the set R = O 12 \ O 2 and assume that R is not empty. The boundary of R is partitioned into a subpart of ∂O 2 and Γ 12 = ∂Ω 12 ∩ ∂O 1 . For any pointx 1 of Γ 12 , we have (5) and (6) imply We first propagate the vanishing property of function u from point x 0 to point x 1 . By using the second part of Proposition 5 in Ω 12 , we obtain from the above estimate of d Ω12 (x 1 , x 0 ) that for T > P + 2D, the function u(x 1 , ·) vanishes in the interval (0, T − D − P − 3ε). Sincex 1 was an arbitrarily chosen point of Γ 12 , we obtain that u vanishes in Γ 12 × (0, T − D − P − 3ε). We remark that any point of the boundary of R is either on ∂O 2 or on ∂O 12 ∩ ∂O 1 , we hence have . Because of the assumption on the continuity of u 2 , we obtain from Lemma 3.5 that u 2 ∈ L 2 (0, T − D − P − 3ε; H 1 0 (R)). Using the initial condition for u 2 we conclude that We now propagate the vanishing property of function u 2 from point x 1 to point x 0 . If we reuse the second part of Proposition 5 in Ω 12 as well as the estimate of d Ω12 (x 1 , x 0 ) given by (5), we obtain that for T > P + 2D, u 2 vanishes at point x 0 in the time interval (0, T − 2D − P − 5ε). Since ε is arbitrarily small we obtain a contradiction to the assumption on g 0 . We conclude that R is empty. Since The same reasoning is then applied with some point x ∈ O 2 to obtain the same contradiction, hence Remark 1. Here we have used the fact that although O 12 is generally not a Lipschitz domain, P (O 12 ) is still defined since its boundary is still rectifiable and  Proof. Let us consider some x ∈ Ω. For arbitrarily small ε > 0, one may find some and since the inequality is valid for all ε and all x ∈ Ω, we obtain D(Ω, Γ) ≤ D(G, Γ). Otherwise, let us consider some path g joining x to x 0 in G such that (g) ≤ d G (x, x 0 ) + ε. Clearly, there exists a curve h in Ω joining x 0 to x consisting of two parts: one part coincides with a sub-path of g, the other part coincides with a subpart of the boundary of O. The length of the first part is less than (g) while the length of the second part is less than D(∂O), so that (h) ≤ D(∂O) + (g). Hence, Problem (7) corresponds to problem (2) in the case when the obstacle O, and hence the complementary domain Ω, is known. However, it is crucial in view of the "exterior approach" introduced hereafter, to regularize such problem without taking into account the boundary condition on ∂O in the problem (7). Let us assume existence of a solution u ∈ H 1 (Q) to problem (7), which means in some sense that (g 0 , g 1 ) are exact data. Then from Theorem 3.3 u is uniquely defined in We now describe our mixed formulation of quasi-reversibility. Let us introduce the setsΓ := ∂Ω \ Γ,Σ =Γ × (0, T ) and S T = Ω × {T }. Let us also introduce The spaces V , V 0 andṼ 0 are endowed with the same norm · given by That · is actually a norm in these three spaces is a consequence of Poincaré's inequality. We define H For what follows we need the following weak characterization of the solutions to problem (7).
is a solution to problem (7) if and only if u ∈ V g and for all µ ∈Ṽ 0 , where the last integral means duality pairing between H Proof. We first consider some u ∈ H 1 (Q) which solves problem (7). Then u ∈ V g . For µ ∈ H 1 (Q), we use the following integration by parts formula, where the bracket ·, · ∂Q has the meaning of duality between H −1/2 (∂Q) and H 1/2 (∂Q). By considering now µ ∈Ṽ 0 , we use the fact that div d+1 v = 0 in Q, that ∂ t u = 0 on S 0 , µ = 0 on S T , ∂ ν u = g 1 on Σ and µ = 0 onΣ to obtain where the bracket ·, · Σ has the meaning of duality between H −1/2 S0 (Σ) and The weak formulation (11) is achieved. Conversely, let us consider some u ∈ V g which satisfies (11). By taking µ ∈ C ∞ 0 (Q), it follows that u solves the wave equation in the distributional sense in Q. Hence formula (12) is valid, and we obtain that for all µ ∈Ṽ 0 v · ν d+1 , µ ∂Q = g 1 , µ Σ , so that by using the fact that µ = 0 onΣ and where Σ 0 = Σ ∪ S 0 and the brackets ·, · Σ0 have the meaning of duality between We conclude that v · ν d+1 =g 1 on Σ 0 , then v · ν d+1 = g 1 on Σ and v · ν d+1 = 0 on S 0 , that is ∂ ν u = g 1 on Σ and ∂ t u = 0 on S 0 . As a conclusion, problem (7) is satisfied by u.
Due to Lemma 4.1, problem (7) is a particular instance of the abstract framework described in [5], which we recall here. We consider three Hilbert spaces V , M and H, endowed with the scalar products (·, ·) V , (·, ·) M and (·, ·) H and corresponding norms · V , · M and · H . We denote A : V → H a continuous onto operator while for some g ∈ H, we denote V g = {u ∈ V, Au = g}, which is an affine space. The corresponding vector space is denoted V 0 , equipped with the norm · V . For a continuous bilinear form b on V × M and a linear form m on M , let us consider the abstract weak formulation: find u ∈ V g such that for all µ ∈ M , (13) b(u, µ) = m(µ).
The bilinear form b is said to satisfy the inf − sup property on V 0 × M if Assumption 4.2. There exists α > 0 such that The bilinear form b is said to satisfy the solvability property on V 0 × M if According to the Brezzi-Nečas-Babuška theorem (see, for example, [12]), the problem (13) ε The following theorem is proved in [5].
where u m is the unique solution to the minimization problem Our problem (11) coincides with the general problem (13) by choosing spaces V as specified by (8), M =Ṽ 0 and H = H 1/2 S0 (Σ) as specified by (9) and (10) (14), which can be directly applied to our problem, consists in the following problem for some real ε > 0: for (g 0 , g  (Σ), problem (16) has a unique solution (u ε , λ ε ) ∈ V g ×Ṽ 0 . If in addition we assume there exists u ∈ H 1 (Q) satisfying problem (7) for where u m is the unique solution of minimal norm · to problem (7). In particular, u m coincides with u in Q 0 defined by (4) and if T > D(Ω, Γ), u m coincides with u in subdomain Ω × (0, T − D(Ω, Γ)).

5.
The level set method. We present a simple level set method and show it enables us to identify the obstacle O assuming that the function u which solves (2) is known in the uniqueness domain Q 0 defined by (4).

5.1.
Preliminaries. This section presents, in a slightly more general case, some results already given in [2]. Let us consider a function U ∈ H 1 (G) such that, for some obstacle O G formed by a collection of a finite number of disjoint connected open Lipschitz domains such that Ω = G \ O is connected, we have where v ω := w ω + U and w ω is the unique solution w in H 1 0 (ω) of problem ∆w = f − ∆U . Here supp denotes the essential support of a function.

Remark 2.
With additional regularity assumptions, we can give a simpler definition of v ωn and ω n+1 . Indeed, in the case when ω n is a Lipschitz domain, the function v ωn is the unique solution to the problem ∆v = f in ω n v = U on ∂ω n , while when v ωn is a continuous function, the domain ω n+1 is defined by Our level set method essentially relies on the weak maximum principle (see, for example, [13]). In what follows, for a functional space V , V + (resp. V − ) denotes the subset of measurable functions v ∈ V such that v ≥ 0 (resp. v ≤ 0) almost everywhere. For example, the function u of Lemma 5.1 satisfies u ∈ H 1 0 (Ω) − . We begin with the following proposition.
Proof. Let us assume that O ⊂ ω n G. We have to prove that O ⊂ ω n+1 . By using Lemma 5.1 and the fact that f − ∆U ≥ 0, we obtain w ωn ≤ 0 in ω n , that is v ωn = w ωn + U ≤ U in ω n . Since U ≤ 0 in O ⊂ ω n , we obtain v ωn ≤ 0 in O. Hence O ⊂ ω n \ supp(sup(v ωn , 0)) = ω n+1 . The proof is complete. Proposition 1 immediately implies the following proposition by using [14].
The definition of the Hausdorff distance for open domains and various relative properties are detailed in [14]. We need the following lemma, which is proved in [14] (Corollary 3.1.12).
We now consider the following assumption, which concerns the continuity of the Dirichlet solution to the Laplace equation with respect to the domain and is extensively analyzed in [14]. Remark 3. Some sufficient conditions on the sequence (ω n ) that enable us to fulfill Assumption 5.3 are given in [14]. In particular, for d = 2, a sufficient condition to have Assumption 5.3 is that all the sets G \ ω n , n ∈ N, are connected, which is a consequence ofŠverák's Theorem (see Theorem 3.4.14 in [14]). Note, however, that we don't know if the sets G \ ω n are connected from the general definition of the ω n .
Let us denote K n and L n the supports of ϕ n and ψ n respectively. We have K n ⊂ G \ O and L n ⊂ ω ∞ , hence K n ∩ L n ⊂ R. For x ∈ R \ (K n ∩ L n ), either ϕ n (x) = 0 and ψ n (x) ≥ 0, or ϕ n ≥ 0 and ψ n = 0, hence θ n (x) = 0. This implies that supp(θ n ) ⊂ K n ∩ L n , that is θ n in compactly supported in R. Since θ n converges to U in H 1 (R), U ∈ H 1 0 (R). in Ω (20)  Proof. First, that u ∈ H 1 (Q) implies
The following theorem shows that under Assumption 5.3 the sequence of domains (ω n ) converges to the true obstacle O.
Theorem 5.5. We consider some obstacle O and some function u which satisfy problem (2) with T > T (O) (O and u are both uniquely defined). We assume in addition that u ∈ L 2 (0, T ; C 0 (Ω)), that for all x 0 ∈ Γ and all sufficiently small t 0 > 0, the function g 0 (x 0 , ·) is not identically zero in the interval (0, t 0 ). Let U ∈ H 1 (G) and f ∈ H −1 (G) satisfy (20) and (18).
Let ω 0 denote an open domain such that O ⊂ ω 0 G and (ω n ) denotes the decreasing sequence of open domains defined by (19).
With additional Assumption 5.3, we have Proof. We already know from Proposition 2 that O ⊂ ω ∞ . Let us denote R = ω ∞ \O and let us assume that R = ∅. We shall find a contradiction. From Proposition 3 we have U ∈ H 1 0 (R). We remark that u = 0 in ∂R × (0, T − T (O)/2), because and u(·, t) ∈ C 0 (Ω) for almost all t. By using the initial condition in (2), we get u = 0 in R × (0, T − T (O)/2). By the same unique continuation argument that was used in the proof of Theorem 3.6 and since T > T (O), we obtain that there exists some x 0 ∈ Γ such that u(x 0 , ·) vanishes in the interval (0, , which contradicts our assumption on g 0 . We conclude that  6. Description of the "exterior approach". In this section, we deduce from the results of Sections 4 and 5 an algorithm to approximately solve the problem (2) presented in the introduction, that is to retrieve the obstacle O from the lateral Cauchy data (g 0 , g 1 ) on Γ × (0, T ). In Theorem 5.5 we have shown convergence of the level set method to the true obstacle O if the field defined by (20) is based on the true solution u. In practice such function is unknown but following Theorem 4.5, it can be approximated in the subdomain Ω × (0, T − D(Ω, Γ)) with the help of the solution u ε of the quasi-reversibility problem (16). This is why we propose the following algorithm in the continuous framework. Step 2: the domain O n being given, solve the quasi-reversibility problem (16) in Ω n × (0, T ) with Ω n := G \ O n for some selected parameter ε > 0. The solution is denoted (u n , λ n ). 3.
Step 3: the function u n being given, solve the non-homogeneous Dirichlet problem

Go back to
Step 2 until the stopping criteria is reached. Alternatively, it can be tempting, in view of Remark 4, to replace T (O n ) by 2D(Ω n , Γ) in the above algorithm. This is what we will do later in our numerical experiments.  7.1. Validation of the quasi-reversibility method. We begin with some preliminary numerical experiments in 2D with the mixed formulation of quasi-reversibility (16) to solve problem (7) when the obstacle O is supposed to be known. The computation domain Ω × (0, T ), with Ω = G \ O, is then fixed. The domain G is here the disc of center (0, 0) and radius R = 5 in R 2 , while Γ is the whole boundary of G. We discretize the space/time domain Q = Ω × (0, T ) with the help of prismatic finite elements, that is tensor products of P 1 triangular finite elements in dimension 2 for the spatial domain Ω and P 1 finite elements in dimension 1 for the time interval (0, T ). We refer to [4] for a description of such finite element and an analysis of the induced discretization error in a slightly different case (heat equation instead of wave equation). Computations are obtained with the Finite Element Library XLife++ [21]. In all our computations related to the inverse problem, the time interval is 4/15 (whatever T is) while the mesh size with respect to space is such that the number of vertices on ∂G is either 164 (in Figures 5, 8 and 9) or 92 (in Figures 3, 4, 6 and 7). In order to build some artificial data (g 0 , g 1 ) on Γ × (0, T ), we compute the solution to the following forward problem for some Neumann data g. This forward computation is based on a classical finite difference (leap-frog) scheme with respect to time and a finite element method with respect to space, in such a way that the meshes and the discretizations of the forward and inverse problems are different. Our artificial data are then given by (g 0 = u| Γ×(0,T ) , g 1 = g). We can then use these lateral Cauchy data (g 0 , g 1 ) as inputs in the mixed formulation of quasi-reversibility (16), in which we set ε = 0.001. We present two cases.
In the first case, the obstacle O is the disc of center (0, 0) and radius r = 2.5, T = 5.5 while g(x, t) = 10 sin 2 (t), which means it is a full radial case. For that geometry, the uniqueness domain defined by (4) reduces to (23) Q In the Figure 3, we have plotted the discrepancy u ε −u between the quasi-reversibility solution and the exact solution for several fixed times t ∈ (0, T ) as a function of the radial distance |x|. From Theorem 4.5, the function u − u ε is supposed to be small in Q 0 for small ε. It can be seen that Figure 3 is consistent with (23) since R − T = −0.5. In the second case, the obstacle is the union of two discs, one centered at (2.5, 1) of radius 1 and one centered at (−2.5, 0) of radius 1.5, T = 15 and g(x, t) = (20 + 10 cos(2θ)) sin 2 (t). Here θ is the polar angle, that is some point x ∈ ∂G has Euclidian coordinates (5 cos θ, 5 sin θ). On Figure 4 we have plotted the quasireversibility solution u ε as well as the discrepancy u ε − u on the meshed space/time domain Q = Ω×(0, T ). The visible slice corresponds to t = T /2. It can be seen that for t = T /2, the error is small near the outer boundary of Ω and is concentrated near the inner boundary. This is due to the fact that the intersection of the uniqueness domain Q 0 and the plane t = T /2 is formed by the set {x ∈ Ω, d Ω (x, Γ) < T /2} and to the fact that the inverse problem in Q 0 is exponentially ill-posed.

7.2.
Validation of the exterior approach. We are now in a position to test the exterior approach, that is the algorithm of Section 6 to solve the inverse obstacle problem (2) in the case when the unknown obstacle is the union of the two discs as before. Note that the cylindrical domain G × (0, T ) is meshed once and for all with the help of the prismatic finite elements described above. At iteration n, the obstacle O n and its complementary domain Ω n are polygonal domains based on the mesh of G. The initial guess O 0 is chosen as the disk centered at (0, 0) and of radius 4.1. In the step 2 of the algorithm at iteration n, the 3D quasi-reversibility problem is solved in domain Ω n × (0, T ) with the help of the tensorized P 1 ⊗ P 1 finite elements, while in step 3, the 2D level set problem is solved in domain O n with simple P 1 finite elements. Let us discuss some critical parameters of the algorithm of Section 6.
1. About the computation of T (O n ): in view of Remark 4, we shall use 2D(Ω n , Γ) instead of T (O n ) in (21), since 2D(Ω n , Γ) may be much smaller than T (O n ). The distance D(Ω n , Γ) can be computed numerically by solving a forward wave propagation problem in Ω n × (0, +∞): Due to the unit speed of propagation, the minimal time of full impact on Ω n (with the pulse shift of 2 taken into account) coincides with D(Ω n , Γ). In practice, we observe that the only values of u which are involved in (21) are taken on ∂Ω n \ ∂G = ∂O n , which leads us to compute instead of D(Ω n , Γ). In our numerical experiments, we hence replace T (O n ) by 2D n in the algorithm. More precisely, time integration stops once all the nodes of the polygonal line ∂O n have been reached by the pulse, let say when |u| > 0.2 × 10 = 2 (20% of the initial pulse amplitude) in order to take the numerical dispersion effects into account. 2. About the choice of f : in practice the distribution f in (21) is chosen as a constant spatial function which must be sufficiently large in view of (18). Following our algorithm, strictly speaking such constant f does not change during the iterations n. However, it may be more efficient to iteratively decrease the value of f in step 3 in order to prevent the algorithm from a preliminary stop due to the fact that the evolution of the level sets may occur on a scale that cannot be resolved numerically because the mesh size is not small enough. This is why we introduce a refinement of step 3 which in some sense finds the optimal value of f . Note that such refinement also includes a stopping criterion. In the following refined step 3, we have set M n+1 = sup x∈∂On V n (x) for n ≥ 0 (M 0 = +∞) and f is a large constant. 3. About the choice of ε: the quasi-reversibility regularization parameter is ε = 0.001 in the following numerical experiments. In practice, this value should be set in accordance with the amplitude of the noise which contaminates the data (g 0 , g 1 ). For instance, the Morozov's discrepancy principle could be applied to the mixed formulation of quasi-reversibility as done in [5] in the simpler case of Helmholtz equation. However, as shown in [5], such procedure requires to take into account the Dirichlet data g 0 in a weak way instead of in a strong way. Since such work seems quite significant, we intend to address it in a future contribution. Before testing the complete algorithm, we wish to test the level set method only. More precisely, in order to illustrate Theorem 5.5, we only apply step 3 of the algorithm of Section 6 by using the exact solution u in Ω n × (0, T − D n ) instead of the quasi-reversibility solution u n . From Theorem 5.5, the sequence of domains O n is supposed to converge to the true domain O in the sense of Hausdorff distance. This convergence can be checked, up to the spatial mesh size, on the Figure 5. In this computation, the final time is T = 25, which is much larger than 2D(Ω, Γ) = 10, that is the minimal time for convergence in Theorem 5.5. The boundaries of O 0 and some selected intermediate obstacles O n until stationarity of the algorithm are represented on the picture. Now let us try experiments with the complete algorithm of Section 6. In Figure  6, we analyze the identification results for exact data and increasing final time. We have tested three values of T , namely T = 2D(Ω, Γ) = 10, T = 15 and T = 25, keeping the same discretization time step. We observe that when T increases, the identification results improve: not only the final obstacle is better retrieved but the number of iterations to achieve the final obstacle is smaller. In Figure 7, we analyze the identification results for fixed final time T = 25 and noisy data. We artificially perturb both Dirichlet data g 0 and Neumann data g 1 by some pointwise random  noise which is then rescaled in such a way that we exactly control the L 2 relative error δ (see [2] for more details). In other words, the noisy data (g δ 0 , g δ 1 ) satisfy g δ 0 − g 0 L 2 (Γ×(0,T )) = δ g 0 L 2 (Γ×(0,T )) , g δ 1 − g 1 L 2 (Γ×(0,T )) = δ g 1 L 2 (Γ×(0,T )) . The three pictures of 7 correspond to δ = 0 (exact data), δ = 2% and δ = 5%. In Figure 7. Two discs and noisy data. Top left: δ = 0 (exact data). Top right: δ = 0.02. Bottom: δ = 0.05 picture 8, for exact data and final time T = 25, we consider the difficult case of partial data. More precisely, data (g 0 , g 1 ) are only known for π/2 ≤ θ ≤ 2π, that is Γ is now formed by 3/4-th of the circle ∂G. We consider obstacles formed by only one disc, either located far away from the subpart of ∂G on which we have no data, or close to that subpart. Clearly, the identification result is better in the first case. Finally, in the last Figure 9, we consider a boomerang-shape obstacle, which is highly non-convex, in the case of full data (Γ = ∂G) and T = 15, either with δ = 0 (no noise) and δ = 2%.
Remark 6. In difficult configurations like in Figure 9, the identification would be better by improving the mixed formulation of quasi-reversibility in different directions, for example local mesh refinement (like in [9]) and iterations of the mixed formulation (like in [10] and [4]).
We are now able to prove Theorem 3.3.