Stable determination of polyhedral interfaces from boundary data for the Helmholtz equation

# Stable determination of polyhedral interfaces from boundary data for the Helmholtz equation

Elena Beretta  Dipartimento di Matematica "Brioschi", Politecnico di Milano, Milano, Italy (elena.beretta@polimi.it)    Maarten V. de Hoop  Department of Mathematics, Purdue University, West Lafayette, USA (mdehoop@purdue.edu)    Elisa Francini  Dipartimento di Matematica e Informatica “U. Dini”, Università di Firenze, Italy (elisa.francini@unifi.it)    Sergio Vessella  Dipartimento di Matematica e Informatica “U. Dini”, Università di Firenze, Italy (sergio.vessella@unifi.it)
July 30, 2019
###### Abstract

We study an inverse boundary value problem for the Helmholtz equation using the Dirichlet-to-Neumann map. We consider piecewise constant wave speeds on an unknown tetrahedral partition and prove a Lipschitz stability estimate in terms of the Hausdorff distance between partitions.

Keywords. Inverse boundary value problem, Helmholtz equation, Lipschitz stability

MSC: 35R30, 35J08, 35J25

## 1 Introduction

We consider an inverse boundary value problem for the Helmholtz equation

 Δu+ω2q(x)u=0in Ω⊂R3,

where and is the wavespeed. The data are the Dirichlet-to-Neumann map and the objective is to recover the wavespeed. The uniqueness of this inverse problem was established by Sylvester and Uhlmann [20] for . Concerning stability, conditional logarithmic continuous dependence of the wavespeed on the Dirichlet-to-Neumann map has been proven in [2] in the case of wavespeeds in with . We refer to Novikov [13] for a refinement of this stability estimate. The logarithmic rate of stability is optimal [12]. For the inverse conductivity problem the authors of [3] proposed restricting the class of unknown coefficients to a finite dimensional set to obtain Lipschitz stability estimates. The result was extended to complex-valued conductivities in [6]. In this finite dimensional setting, in [4, 5], a Lipschitz stability estimate for the recovery of piecewise constant wavespeeds for a given domain partition from boundary data for the Helmholtz equation, and an estimate for the stability constant in terms of the number of domains in the partition, were obtained.

Here, we study the problem of determining the finite partition from boundary data given a (possibly large) finite set of attainable values for the wavespeed. Due to the severe nonlinearity of the problem the derivation of Lipschitz stability estimates is more subtle. For this reason, we consider a partitoning of the domain with a (regular) unstructured tetrahedral mesh. In fact, an unstructured tetrahedral mesh admits a local refinement and, with piecewise constant wavespeeds, can accurately approximate realistic models in applications. In geophysics, we mention as an example the work of Rüger and Hale [16]. Here, knowledge of a set of attainable values for the wavespeed can be motivated by the general knowledge of relevant rock types. The deformation allows one to adjust the mesh and recover structures in the models. In geodynamics, these structures can be an imprint of the local geology and tectonics [18]. Moreover, one can parametrize major discontinuities at (polyhedral) surfaces by connecting boundaries of subdomains in the partition via a segmentation for example.

In this paper, we establish a Lipschitz stability estimate expressed in terms of the Hausdorff distance between partitions using tetrahedra from the Dirichlet-to-Neumann map. Lipschitz stability estimates provide a framework for optimization, specifically, iterative reconstruction of the wavespeed with a convergence radius determined by the stability constant [7, 8]. The recovery of polyhedral interfaces then becomes a shape optimization. The analysis in [7] makes explicit use of a Landweber iteration. Via successive approximations, and making use of estimates for the corresponding growth of the stability constant, the reconstruction can be cast into a multi-level scheme [8] effectively enlarging the radius of convergence. As an important application, we mention so-called time-harmonic full waveform inversion (FWI) developed in reflection seismology [14, 15, 19, 21] with the goal to image wavespeed variations in Earth’s interior. The data, here, are essentially the single-layer potential operator. However, stability estimates for the Dirichlet-to-Neumann map directly carry over to stability estimates for this operator.

We give an outline of the paper. We first state the main result and the main assumptions (Section LABEL:sec:2). Then we establish a rough stability estimate for the potentials using complex geometrical optics (CGO) solutions following the outline of an estimate in Beretta et al. [5] (Section LABEL:sec:3). The CGO solutions were introduced by Sylvester and Uhlmann [20] in their proof of uniqueness of this inverse boundary value problem. The CGO solutions in our analysis differ slightly from theirs to obtain better constants in the stability estimates as proposed in [17]. We proceed with establishing the recovery of the number of tetrahedra in the mesh from the potential, and with expressing the Hausdorff distance between meshes in terms of the difference of piecewise constant potentials defined on these meshes. Naturally, the information on the Hausdorff distance between meshes can be transformed to information on the vertices of the tetrahedra forming the meshes (Section LABEL:sec:4). The main part of the proof of our result pertains to obtaining a lower bound for the Gateaux derivative of the Dirichlet-to-Neumann map under mesh deformation (Section LABEL:sec:5).

### Notation

We use the Fourier transform convention,

 ^f(ξ)=∫R3f(x)eix⋅ξdx.

If the function is defined on a subset of , it is extended to attaining the value zero. We denote by the inverse Fourier transform of ,

 ˇf(x)=1(2π)3∫R3f(ξ)e−ix⋅ξdξ. \hb@xt@.01(1.1)

We introduce coordinates, , in , where and . We denote the open ball in centered at of radius by , and the open ball in centered at of radius by .

## 2 Assumptions and main result

We let be a bounded domain in such that is connected,

 Ω⊂BR(0) for some R>0, \hb@xt@.01(2.1)

and

 Ω has a Lipschitz boundary with constants r0 and K0, \hb@xt@.01(2.2)

that is, for any point , there exists a rigid transformation of coordinates under which and

 Ω∩{(x′,x3)∈R3:|x′|ψ(x′)},

where is a Lipschitz continuous (level set) function in such that

 ψ(0)=0 and ∥∇ψ∥L∞(B′r0)≤K0.

We consider the boundary value problem for the Helmholtz equation,

 {Δu+ω2qu=0 in Ω,u=ϕ on ∂Ω \hb@xt@.01(2.3)

for , and introduce the Dirichlet-to-Neumann map

 Λq:H1/2(∂Ω)→H−1/2(∂Ω) \hb@xt@.01(2.4)

according to

 ϕ→Λq(ϕ):=∂u∂ν∣∣∣∂Ω. \hb@xt@.01(2.5)

The normal derivative is defined in the weak sense as

for every . In the above, is identified with where denotes the wavespeed. The solution of (LABEL:4.1) exists in and is unique if is not in the Dirichlet spectrum of on .

We introduce , such that and

 ω1≤√λ1(BR)2Q0, \hb@xt@.01(2.6)

where is the first eigenvalue of on . We recall that . (If we detect the spectrum, we substitute the true first eigenfrequency for .) We then assume that

 ω0≤ω≤ω1. \hb@xt@.01(2.7)

### Unstructured tetrahedral mesh

We let be a regular partition of into tetrahedra, namely a collection of closed tetrahedra such that

 ¯¯¯¯Ω=∪Nj=1Tj; \hb@xt@.01(2.8)
 for j≠k either Tj∩Tk=∅ % or it consists of a common vertex, \hb@xt@.01(2.9) a common edge or a common facet;
 the radius of the insphere of each tetrahedron is larger than r1>0. \hb@xt@.01(2.10)

We say that two different tetrahedra of such regular partition are adjacent if they share a common facet.

###### Remark 1

Assumption (LABEL:1.5.5), together with (LABEL:1.1) implies that the tetrahedra of the partition are not degenerate. In particular, there are two positive numbers and (depending on and only) such that

 for each Tj the distance between vertices is % greater than d1 \hb@xt@.01(2.11) and internal angles of triangular facets are greater than α1.

Indeed, we point out that assumptions (LABEL:1.5.5) and (LABEL:1.1) are equivalent to the following

###### Assumption 1

There exists a positive constant such that

 ∣∣Br(P)∩Tj∣∣≥C1r3, \hb@xt@.01(2.12)

for every , every , and .

We show an illustration of a typical model and the assumptions pertaining to the mesh in Figure LABEL:fig:1.

We introduce a finite set of numbers,

 Q={~q1,…,~qL}

representing the possible values which the wavespeed can attain in the domain ,

 Q0=max{|~qj|:j=1,…,L}, \hb@xt@.01(2.13)

and

 c0=min{|~qj−~qk|:j,k=1,…,L,j≠k}. \hb@xt@.01(2.14)
###### Assumption 2

The potentials are piecewise constant and of the form

 q(x)=N∑j=1qjχTj(x) \hb@xt@.01(2.15)

such that is a regular partition of with

 N≤N0 \hb@xt@.01(2.16)

for some ,

 qj∈Q for every j=1,…,N, \hb@xt@.01(2.17)

and

 qj≠qk if Tj is adjacent to Tk. \hb@xt@.01(2.18)

We denote by the norm in defined by

 ∥T∥⋆=sup{⟨Tϕ,ψ⟩:ϕ,ψ∈H1/2(∂Ω), ∥ψ∥H1/2(∂Ω)=∥ϕ∥H1/2(∂Ω)=1}.

We refer to the values of , , , , , , , and as to the a priori data. In the sequel we will introduce a number of constants that we will always denote by and, unless otherwise stated, will depend on a priori data only. The values of these constants might differ from one line to the other.

We state the main result

###### Theorem 2.1

Given a domain satisfying (LABEL:1.1) and (LABEL:1.1.5), a set of values , and , there exist two positive constants and depending on the a priori data and on only such that, for every pair of potentials

 q(0)=N∑j=1q(0)jχT(0)j and q(1)=M∑k=1q(1)kχT(1)k \hb@xt@.01(2.19)

satisfying Assumptions LABEL:Ass:1 and LABEL:Ass:2, if

 ∥Λq(0)−Λq(1)∥⋆≤ε0, \hb@xt@.01(2.20)

then

 N=M \hb@xt@.01(2.21)

and the order of the tetrahedra can be rearranged so that for every we have

 q(0)j=q(1)j, \hb@xt@.01(2.22)

and

 dH(T(0)j,T(1)j)≤C0∥Λq(0)−Λq(1)∥⋆, \hb@xt@.01(2.23)

where denotes the Hausdorff distance.

## 3 A rough stability estimate

We begin with developing a rough stability estimate for the recovery of the potential or wavespeed.

###### Theorem 3.1

Given , , and as in Theorem LABEL:MainTheorem, there exist two positive constants and depending on , , , , , such that, for ,

 ∥q(0)−q(1)∥L2(Ω)≤C2√N0∣∣log(∥Λq(0)−Λq(1)∥⋆)∣∣−1/7. \hb@xt@.01(3.1)

Proof. We proceed as in [5]. Alessandrini’s identity states that

 ω2∫Ω(q(0)−q(1))u0u1dx=⟨(Λ0−Λ1)(u0|∂Ω),u1|∂Ω⟩ \hb@xt@.01(3.2)

for every pair of functions and such that

 Δuk+ω2q(k)uk=0 in Ω for k=0,1,

where we use the shorthand notation, .

We fix and let and be unit vectors in such that is an orthogonal set of vectors. We let be a parameter to be chosen later, and set, for ,

 ζk=⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩(−1)k+1μ√2(√1−|ξ|22μ2η1+(−1)k√2μξ+iη2) if |ξ|μ√2<1,(−1)k+1μ√2((−1)k√2μξ+i√|ξ|22μ2−1η1+η2) if |ξ|μ√2≥1. \hb@xt@.01(3.3)

As can be easily checked,

 ζ0+ζ1=ξ,
 ζk⋅ζk=0 for k=0,1

and

 |ζk|=max{μ,|ξ|√2}. \hb@xt@.01(3.4)

We use here complex geometrical optics (CGO) solutions of the Helmholtz equation and, in particular, the estimates in [17, Theorem 3.8] which are due to [9]. For , there is a solution of

 Δuk+ω2q(k)uk=0 in Ω

of the form

 uk(x)=eix⋅ζk(1+φk(x)), \hb@xt@.01(3.5)

with

 ∥φk∥L2(Ω) ≤ Cω12Q0|ζk|≤Cω12Q0μ, \hb@xt@.01(3.6) ∥∇φk∥L2(Ω) ≤ Cω12Q0,

where .

Inserting (LABEL:e2.3) into (LABEL:e1.2), we get

 ω2∣∣(ˆq(0)−ˆq(1))(ξ)∣∣≤∣∣⟨(Λ0−Λ1)(u0|∂Ω),u1|∂Ω⟩∣∣ +ω2∣∣∣∫Ω(q(0)(x)−q(1)(x))eiξ⋅x(φ0(x)+φ1(x)+φ0(x)φ1(x))dx∣∣∣ ≤∥Λ0−Λ1∥⋆∥u0∥H1(Ω)∥u1∥H1(Ω)+2ω2Q0∣∣∣∫Ω(φ0+φ1+φ0φ1)dx∣∣∣.

Hence,

 ∣∣(ˆq(0)−ˆq(1))(ξ)∣∣2 ≤2ω04∥Λ0−Λ1∥2⋆∥u0∥2H1(Ω)∥u1∥2H1(Ω)+8Q20∣∣∣∫Ω(φ0+φ1+φ0φ1)dx∣∣∣2 ≤2ω04∥Λ0−Λ1∥2⋆∥u0∥2H1(Ω)∥u1∥2H1(Ω)+8Q20|Ω|(∥φ0∥L2(Ω)+∥φ1∥L2(Ω)) +8B20∥φ0∥L2(Ω)∥φ1∥L2(Ω).

With (LABEL:e2.3) and (LABEL:e3.1) we find that there exists a constant depending only on such that, for ,

 ∥uk∥H1(Ω)≤Ce2R(μ+|ξ|), \hb@xt@.01(3.7)

, where . Hence,

 ∣∣(ˆq(0)−ˆq(1))(ξ)∣∣2≤C(e8R(μ+|ξ|)∥Λ0−Λ1∥2⋆+1μ2), \hb@xt@.01(3.8)

where . But then, for ,

 ∥q(0)−q(1)∥2L2(Ω) = ∫|ξ|≤ρ∣∣(ˆq(0)−ˆq(1))(ξ)∣∣2dξ+∫|ξ|>ρ∣∣(ˆq(0)−ˆq(1))(ξ)∣∣2dξ ≤ Cρ3(e8R(μ+ρ)∥Λ0−Λ1∥2⋆+1μ2) +∫|ξ|>ρ∣∣(ˆq(0)−ˆq(1))(ξ)∣∣2dξ.

To estimate the integral in (LABEL:e4.2) we show that for every

 ∥q(0)−q(1)∥2Hs(Ω)≤C√N0, \hb@xt@.01(3.10)

where . Indeed, by [11] we have

 ∥q(0)−q(1)∥2Hs(Ω)≤2(∥q(0)∥2Hs(Ω)+∥q(1)∥2Hs(Ω)) ≤2(N∑j=1|q(0)j|2|T(0)j|1−2s|∂T(0)j|2s+M∑k=1|q(1)k|2|T(1)k|1−2s|∂T(1)k|2s) ≤CN0,

where .

Using (LABEL:e5.1),

 ∫|ξ|>ρ∣∣(ˆq(0)−ˆq(1))(ξ)∣∣2dξ ≤ 1ρ2s∫|ξ|>ρ(1+|ξ|s)2∣∣(ˆq(0)−ˆq(1))(ξ)∣∣2dξ \hb@xt@.01(3.11) ≤ 1ρ2s∥q(0)−q(1)∥2Hs(Ω)≤CN0ρ2s.

Finally, by inserting (LABEL:e5.2) into (LABEL:e4.2), we get that

 ∥q(0)−q(1)∥2L2(Ω)≤CN0{ρ3(e8R(μ+ρ)∥Λ0−Λ1∥2⋆+1μ2)+1ρ2s},

where . We then choose

 ρ=μ23+2s,

and observe that there is a constant depending only on such that, for ,

 ρ3e8R(μ+ρ)≤e18Rμ

so that

 ∥q(0)−q(1)∥2L2(Ω)≤CN0⎛⎜⎝e18Rμ∥Λ0−Λ1∥2⋆+1μ4s3+2s⎞⎟⎠,

where .

We now take

 μ=118R|log∥Λ0−Λ1∥⋆|

and assume that

 ∥Λ0−Λ1∥⋆≤e−18Rc3=:ε1

so that . Then

 ∥q(0)−q(1)∥2L2(Ω)≤CN0(∥Λ0−Λ1∥2⋆+∣∣∣∫log∥Λ0−Λ1∥⋆∣∣∣−α),

where . The claim follows upon choosing .

Next, we establish an estimate for the Haussdorff distance between two domain partitions in terms of the difference of potentials defined on these partitions.

###### Proposition 3.2

Given , and as in Theorem LABEL:MainTheorem, there exists a positive constant depending on , , and such that, if

 ∥q(0)−q(1)∥L2(Ω)≤σ1 \hb@xt@.01(3.12)

then

 N=M \hb@xt@.01(3.13)

and the order of the tetrahedra can be rearranged so that for every

 q(0)j=q(1)j \hb@xt@.01(3.14)

and

 dH(T(0)j,T(1)j)≤∥q(0)−q(1)∥2/3L2(Ω)(c20C1)1/3, \hb@xt@.01(3.15)

where is given by (LABEL:1.3) and by (LABEL:2.2).

Proof. We write

 σ=∥q(0)−q(1)∥L2(Ω). \hb@xt@.01(3.16)

For every we let

 B(0)l={j∈{1,…,N}:q(0)j=~ql} \hb@xt@.01(3.17)

and

 B(1)l={k∈{1,…,M}:q(1)k=~ql}. \hb@xt@.01(3.18)

We note that

 ∥q(0)−q(1)∥2L2(Ω) = L∑l=1⎛⎜ ⎜⎝∑j∈B(0)l∑k∉B(1)l∣∣q(0)j−q(1)k∣∣2∣∣T(0)j∩T(1)k∣∣⎞⎟ ⎟⎠. \hb@xt@.01(3.19)

If and then, by (LABEL:1.3),

 ∣∣q(0)j−q(1)k∣∣≥c0;

hence, by (LABEL:P2.3) and (LABEL:P2.0), we have

 σ2 ≥ c20L∑l=1∑j∈B(0)l∑k∉B(1)l∣∣T(0)j∩T(1)k∣∣ \hb@xt@.01(3.20)

so that

 ∣∣T(0)j∩T(1)k∣∣≤σ2c20 for every j,k such that q(0)j≠q(1)k. \hb@xt@.01(3.21)

By assumption (LABEL:2.1), estimate (LABEL:P3.2) implies that is close to . To make this precise, we introduce

 T(0)j,δ={x∈T(0)j:d(x,∂T(0)j)>δ}

and prove that

 T(1)k∩T(0)j,δσ=∅ \hb@xt@.01(3.22)

with

 δσ=(σ2c20C1)1/3. \hb@xt@.01(3.23)

Indeed, assume that for some , and that there is a point such that

 d(P,∂T(0)j)≥δ, \hb@xt@.01(3.24)

that is, . Using assumption (LABEL:2.1) and (LABEL:2.2) in Remark LABEL:rem1, it then follows that

 ∣∣T(0)j∩T(1)k∣∣≥∣∣Bδ(P)∩T(1)k∣∣≥C1δ3 \hb@xt@.01(3.25)

if . By (LABEL:P3.2)

 C1δ3≤σ2c20. \hb@xt@.01(3.26)

Thus (LABEL:P3.3) holds provided that

 δσ=(σ2c20C1)1/3≤r1,

that is,

 σ≤σ1=√r31c20C1. \hb@xt@.01(3.27)

Now we consider for and for some . Since is a partition of , we can write

 T(0)j,δσ = T(0)j,δσ∩(M⋃k=1T(1)k) = M⋃k=1(T(0)j,δσ∩T(1)k).

Using (LABEL:P3.3),

 T(0)j,δσ∩T(1)k=∅ for k∉B(1)l,

and we then obtain

 T(0)j,δσ=⋃k∈B(1)l(T(0)j,δσ∩T(1)k). \hb@xt@.01(3.28)

If and , then and cannot be adjacent by assumption (LABEL:3.4). This means that there is a unique such that

 T(0)j,δσ∩T(1)k≠∅ \hb@xt@.01(3.29)

and, with (LABEL:P5.1),

 T(0)j,δσ=T(0)j,δσ∩T(1)k⊂T(1)k.

Thus we proved that for every there is a unique index such that

 q(0)j=q(1)¯¯¯k(j) \hb@xt@.01(3.30)

and

 T(0)j,δσ⊂T(1)¯¯¯k(j). \hb@xt@.01(3.31)

In particular, this implies that .

By interchanging the roles of and it follows that , is a permutation on and

 T(0)j,δσ⊂T(1)¯¯¯k(j) and T(1)¯¯¯k(j),δσ⊂T(0)j

that, by (LABEL:P3.4), gives (LABEL:P1.4).

Combining Theorem LABEL:rough and Proposition LABEL:proptetr, we obtain the following logarithmic stability estimate

###### Corollary 3.3

Under the assumptions of Theorem LABEL:rough, there is a constant depending only on the a priori data such that, if

 ∥Λq(0)−Λq(1)∥⋆≤ε2

then

 N=M

and the order of tetrahedra can be rearranged so that

 q(0)j=q(1)j

and

 dH(T(0)j,T(1)j)≤(C22N0c20C1)1/3∣∣log(∥Λq(0)−Λq(1)∥⋆)∣∣−2/21. \hb@xt@.01(3.32)

## 4 Geometric estimates, construction of an intermediate partition and augmenting the domain

Here, we map the information on the Haussdorff distance of tetrahedra in information on the distance between vertices of these tetrahedra. It is straightforward to see that if , , are tetrahedra generated by vertices , , that then

 dH(T(0),T(1))≤min℘max1≤i≤4∣∣P(0)i−P(1)℘(i)∣∣, \hb@xt@.01(4.1)

where denotes a permutation on the set . Moreover, if and satisfies assumption (LABEL:1.5.5) for , then there exists a positive constant , depending on and only, such that

 min℘max1≤i≤4∣∣P(0)i−P(1)℘(i)∣∣≤A1dH(T(0),T(1)). \hb@xt@.01(4.2)

Using Corollary LABEL:cor1 we then obtain

###### Proposition 4.1

Under the assumptions of Theorem LABEL:rough, there is a positive constant such that if

 ∥Λq(0)−Λq(1)∥⋆≤ε3

then for every vertex of (with ) there is a unique vertex of such that

 d(P(0)j,i,P(1)j,i)≤d14 \hb@xt@.01(4.3)

for as in LABEL:2.1.

Proof. It is sufficient to consider , such that

 A1(C22N0c20C1)1/3|log(ε3)|−2/21

and the statement follows.

We introduce a deformation of the tetrahedra forming the partition of . To this end, for each , we define tetrahedra by its vertices,

 P(t)j,i=P(0)j,i+tvj,i for t∈[0,1], \hb@xt@.01(4.4)

where

 vj,i=P(1)j,i−P(0)j,i. \hb@xt@.01(4.5)

The resulting partition is a regular partition of satisfying condition (LABEL:1.5.5). We point out that, by (LABEL:geo1) and (LABEL:geo2), there is a positive constant such that

 A−12(