Memoirs on Diﬀerential Equations and Mathematical Physics Volume 53, 2010, 63–97

(1)

Volume 53, 2010, 63–97

J. Elschner, G. C. Hsiao, and A. Rathsfeld

RECONSTRUCTION OF ELASTIC OBSTACLES FROM THE FAR-FIELD DATA OF SCATTERED ACOUSTIC WAVES

(2)

Abstract. We consider the inverse problem for an elastic body emerged in a fluid due to an acoustic wave. The shape of this obstacle is to be reconstructed from the far-field pattern of the scattered wave. For the numerical solution in the two-dimensional case, we compare a simple Newton type iteration method with the Kirsch–Kress algorithm. Our computational tests reveal that the Kirsch–Kress method converges faster for obstacles with very smooth boundaries. The simple Newton method, however, is more stable in the case of not so smooth domains and more robust with respect to measurement errors.

2010 Mathematics Subject Classification. 35R30 76Q05 35J05.

Key words and phrases. Acoustic and elastic waves, inverse scattering, simple Newton iteration, Kirsch–Kress method.

æØ . łª Œ Œª Ł ª æ æ Ł Øº Ø Æ

ł œ æŁ Æ Æ æŁ ª Æ ØæŁ æŒ æŁ Øºø Œ . Ø Æ -

ºŁ º Ø Æ Œ Ł æŒÆ Œ Œ æŁ Ł º æŁ ª Ł

ØºŒ ø Ø . º Œ ºØ Ł Œ Ø ª ª Ø ª Œ æ ºŒ -

ø Ø æŁ ø ª ØºŒ Œ ª Æ - Ł º Ø Ø -

æŁ ØºŒ Œ . æ Øº ªŁ ª łª Œ Œ, ºØ Ø Ł Łæª

Ø ºŒ ª Œ Æ ºŁ ª - Ø ºÆ æ æŒª Łıº

æ º ß Æº . Ø º Ø ª, Ø ª Œ æ ºŒ Ø º-

Æ ø Ł æ º Łæ Ø Ø ª ª , ºÆ ø Æ ºŁ

ª łŒ Ø Ł Łæª , ª æ º ØÆ Æ ºØª

øÆºØ Ł Ø Ø .

(3)

1. Introduction

If an elastic body is subject to an acoustic wave propagating through the surrounding fluid, then an elastic wave is generated inside the body, and the acoustic wave is perturbed (cf. Figure 1). The wave perturbation is characterized by the asymptotics of the scattered field, namely, the far- field pattern. Suppose the material properties of body and surrounding fluid are known. Then the usual inverse problem of obstacle scattering is to determine the shape of the body from measured far-field data generated by plane waves incident from one or from a finite number of directions.

This problem is extremely ill-posed such that regularization techniques are needed for the solution.

Incident Acoustic Wave Elastic Body

Scattered Field Γ

Ω

Compressible Fluid

c= IR^d

Ω (ΩUΓ) Figure 1. Acoustic wave and obstacle.

Clearly, the same numerical methods used for the inverse problems for obstacles with sound-hard and sound-soft boundaries or for penetrable obstacles can be adapted to the scattering by elastic bodies. Among the available numerical methods, in recent years sampling and factorization methods are very popular (cf. e.g. [12]). Without any a priori information about geometrical details like connectivity components or holes, these methods provide good approximations for the shape of the obstacle. The case of acoustic scattering by elastic bodies in [17] is treated by the linear sampling method. Classical methods such as in [3], [15] (cf. [4] for the case of scattering by elastic obstacles) generally require more information on the geometry of the obstacle. For instance, the boundary of the obstacle is required to be homeomorphic to a circle for 2-D and to a sphere for 3-D problems, respectively. Starting from a reasonable initial guess, the parametrization of the obstacle boundary is approximated in a Newton type iteration. Though the accuracy of the reconstructed solution is always

(4)

66 J. Elschner, G. C. Hsiao, and A. Rathsfeld

limited by the ill-posedness, we expect the classical Newton approach to be more accurate than the factorization methods. To avoid the solution of direct problems in each step of iteration, besides the boundaries also the wave field can be included into the components of the iterative solutions. For instance, a method proposed by Kirsch and Kress (cf. e.g. [13], [3], [25] and cf. [5] for the case of scattering by elastic obstacles) represents the waves by potentials with generating layer functions defined over artificial curves.

Note that, for inverse problems in acoustic scattering by elastic obstacles, difficulties with unpleasant eigensolutions of the direct problem, referred to as Jones modes, can be avoided if the Kirsch–Kress method is applied.

In this paper we consider the two-dimensional case and compare the simple Newton method of [4] with the Kirsch–Kress method of [5] for which we present numerical results for the first time. We implement the same parametrization for the approximate boundary curves iterated by both numerical methods. For a simple egg shaped domain and for a nonconvex domain, we apply the Newton method and the Kirsch–Kress algorithm.

The numerical tests show that the Kirsch–Kress method is more accurate due to the better approximation of the fields by potentials in the case of analytic boundaries. Unfortunately, this method is related to an integral equation approach for the direct problem. If the latter integral equation is severely ill-posed, then the Kirsch–Kress algorithm is divergent. Conse- quently, this method diverges if the curves for the potential representations are too far from the boundary curve of the true obstacle or if the latter curve has large Fourier coefficients. In particular, for the reconstruction of the nonconvex obstacle, the Kirsch–Kress method is divergent. To obtain a convergent version of this method, we use a variant with updated curves for the potential representations during the iteration. For transmission and boundary value problems in acoustic scattering, a comparable update of curves has been proposed in [24] (cf. also the curve updates in [21, Chapter 5] and [15]). Furthermore, our numerical examples reveal that the Kirsch–

Kress method is more sensitive with respect to noise in the far-field data, which is also typical for a higher degree of ill-posedness. Finally, we present an example for the reconstruction of an obstacle with Jones modes. Both methods converge for this case.

We start discussing the solution of the direct problem in Section 2. Using the direct solution, we introduce the two numerical schemes for the inverse problem in Section 3. Then we recall the convergence results from [4], [5]. In Section 4 we discuss some details of the implementation. For the least squares problem of the Kirsch–Kress method, we give the formulas for the functional and its gradients in the appendix. Finally, we present the numerical results in Section 5.

2. Direct Problem: Elastic Obstacle in Fluid

Suppose a bounded elastic body is emerged in a homogeneous compressible inviscid fluid. We denote the domain of the body by Ω, its boundary

(5)

ΩR

R Ω

Γ0

Γ

ν ν

Ω^c

Figure 2. Domains.

curve by Γ (cf. Figure 2), and assume that an incoming plane wave is mov- ing in the exterior Ω^c:=R²\Ω toward the body. This wave is scattered by the body and generates an elastic wave inside the body. Mathematically, the acoustic wave is described by the pressure perturbationpover Ω^c and by the displacement functionuon Ω. The displacement fulfills the Navier (time-harmonic Lam´e) equation

∆^∗u(x) +%ω²u(x) = 0, x∈Ω, (1)

∆^∗u(x) :=µ∆u(x) + (λ+µ)∇[∇ ·u(x)].

Here ω is the frequency, % the density of body, and λ, µ are the Lam´e constants. The total pressurepis the sum of the incoming wavep^incand the scattered wavep^swhich satisfies the Helmholtz equation and the radiation condition at infinity

∆p^s(x) +kw2p^s(x) = 0, x∈Ω^c, (2) x

|x|· ∇p^s(x)−ikwp^s(x) =o(|x|^−1/2), |x| → ∞, (3) where k_w² = ω²/c² is the wave number and c the speed of sound. The pressure and the displacement field are coupled through the transmission conditions

u(x)·ν(x) = 1

%fω²

n∂p^s(x)

∂ν +∂p^inc(x)

∂ν o

, x∈Γ, (4)

t[u](x) = −©

p^s(x) +p^inc(x)ª

ν(x), x∈Γ, (5)

t[u](x) := 2µ∂u

∂ν

¯¯

¯Γ+λ[∇ ·u]ν¯

¯Γ+µν×[∇ ×u]¯

¯Γ,

(6)

ν×[∇ ×u]¯

¯Γ:=

µν2(∂x1u2−∂x2u1) ν1(∂x2u1−∂x1u2)

¶ ¯¯¯

¯¯

Γ

.

Here %f is the density of the fluid and ν denotes the unit normal at the points of Γ exterior with respect to Ω.

For numerical computations, we truncate the exterior domain Ω^c to the annular domain ΩR with the outer boundary Γ0 (cf. Figure 2). The Helmholtz equation (2) is solved over ΩR, and a non-local boundary condition is imposed on Γ0(cf. the boundary integral equation techniques in [11]).

Using standard techniques, the boundary value problem can be reformulated in a variational form and solved by the finite element method (cf. [10], [16], [4]). Suppose the boundary Γ of the obstacle is piecewise smooth and choose the auxiliary curve Γ0 such that the corresponding interior domain has no Dirichlet eigenvalue equal tok_ω² for the negative Laplacian. Then existence and uniqueness of the variational solutions as well as the convergence of the finite element method (cf. [4]) can be shown whenever there is no nontrivial solutionu0 of

∆^∗u0(x) +ρω²u0(x) = 0, x∈Ω,

t[u0](x) = 0, x∈Γ, (6) u0(x)·ν = 0, x∈Γ.

Note that nontrivial solutions of (6) are called Jones modes, and a frequency ω, for which the given domain Ω has a nontrivial solution of (6), is called Jones frequency. It is known that domains with Jones frequencies exist but are exceptional. More precisely, Harg´e [8] has shown that the set of domains with Jones frequencies is nowhere dense in a certain metric, and Natroshvili et al. [19] have proved that domains with two non-parallel flat faces have no Jones frequencies. An example of a two-dimensional domain with Jones frequencyω is the disk ΩJ :={x∈R² : |x| < rJ} with rJ = _ω¹p

µ/% r_J⁰. Herer_J⁰ is any of the positive roots of the equationrJ₁⁰(r) =J1(r), andJ1

is the Bessel function of order one. One Jones mode over ΩJ is defined by u0(x) =J1

³ ω

r%

µ|x|´ µ−x2/|x|

x1/|x|

¶

, x∈ΩJ. (7) Note that the smallest positive root ofrJ₁⁰(r) =J1(r) isr⁰_J = 5.135622. . .. Three-dimensional Jones modes are described, e.g., in [19].

Alternatively to the finite element solution, the complete pressure function and the displacement field can be approximated by potentials with sources over auxiliary curves (cf. Figure 3). We introduce the curve Γi

“close” to Γ, but inside Ω, and the curve Γ_e in Ω_R surrounding Γ. We represent the pressure and the displacement by

p^s(x) =£ V_Γ^ac_iϕi

¤(x), x∈Ω^c, u(x) =£ V_Γ^el_eϕ~e

¤(x), x∈Ω (8)

(7)

00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000

11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111 11111111111111111111111111111111111111111111

Γ0

Γ Γ

Γ

e

i

Ω ΩR

Figure 3. Domains and auxiliary curves

with a scalar layer function ϕi and a vector layer ϕ~e. The potentials are defined by

£V_Λ^acp¤ (x) :=

Z

Λ

p(y)G(x, y;kω)dΛy, x∈R², (9)

G(x, y;kw) := i

4H₀⁽¹⁾(kw|x−y|), (10)

£V_Λ^elu¤ (x) :=

Z

Λ

G^el(y, x)u(y)dΛy, x∈R², (11)

G^el(y, x) := 1 µ

µ

G(x, y;ks)δij+ 1 k_s²

∂²¡

G(x, y;k_s)−G(x, y;k_p)¢

∂xi∂xj

¶₂

i,j=1

, where the wave numberskpand ks are defined by%ω²= (λ+ 2µ)k²_p=µk²_s and H₀⁽¹⁾ is the Hankel function of the first kind and of order 0. The layer functions in (8) are chosen such that the corresponding pressure and displacements fields satisfy the transmission conditions (4) and (5). In other words, to get a good approximate solution we have to solve the integral equations

t[V_Γêl_eϕ~e](x) + [V_Γâc_iϕi](x)ν(x) =−pîncν(x), x∈Γ, (12)

%fω²ν(x)·[V_Γêl_eϕ~e](x)−∂ν[V_Γâc_iϕi](x) =∂νpînc(x), x∈Γ. (13) Numerical methods based on the discretization of (9), (11), (12) and (13) are well-known to exhibit high rates of convergence (cf. e.g. Sect. 9.8 in [6]

and [2], [7], [9]). However, for not so simple geometries, an appropriate choice of Γi and Γe and an appropriate quadrature of the integrals is not

(8)

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

x initial_solution

curve in.circle ex.circle

Figure 4. Geometry of scatterer.

trivial. A bad choice may lead to extremely ill-posed equations (12), (13) and to false results.

For a point x tending to infinity, the scattered pressure field p^s(x) is known to have the following asymptotics

p^s(x) = e^ik^w^|x|

|x|^1/2 p^∞³x

|x|

´

+O³ 1

|x|^3/2

´

, |x| → ∞, (14) p^∞(e^it) =− e^iπ/4

√8πkω

Z

Γi

e^−ik^w^y·e^itϕi(y) dΓiy.

The functionF[p^s](t) :=p^∞(e^it) is called the far-field pattern of the scattered field. This is the entity which can be measured.

In order to prepare the numerical results for the inverse problem, we con- clude this section by the computation of the corresponding direct problem.

If we choose the nonconvex domain with boundary curve Γ according to Figure 4, the constants

ω= π

2 kHz, %= 6.75·10⁻⁸ kg/m³, λ= 1.287373095 Pa, µ= 0.66315 Pa, c= 1500 m/s, %f = 2.5·10⁻⁸ kg/m³,

(15)

and the direction of the incoming plane wave equal tov= (1,0)^>, then we get by the finite element method [4] the far-field pattern plotted in Figure 5.

3. Inverse Problem and Iterative Approximation

Now we suppose that the boundary curve Γ of the obstacle is star-shaped and included between the inner curve Γi :={x∈ R² : |x| =ri} and the

(9)

−1 1 3

30

210

60

240

90

270 120

300 150

330

180 0

far−field pattern

Imag.Part Real Part

Figure 5. Far-field pattern.

outer curve Γe:={x∈R²: |x|=re}, i.e., Γ :=©

r(t)e^it: 0≤t≤2πª

, r(t) =ba0+ X∞

j=1

©bajcos(jt) +bbjsin(jt)ª (16) with the constraintri<r(t)< re, 0≤t≤2π. To avoid this constraint, we can use the parametrization Γ = Γ^r

Γ^r:=©

er(t)e^it: 0≤t≤2πª

, er(t) := re+ri

2 +re−ri

π arctan(r(t)). (17) Having in mind this representation, the star-shaped curve is uniquely determined by the real valued function ror, equivalently, by the Fourier coefficients {baj,bbj}. The direct problem of the previous section defines a continuous mapping (cf. [4])

F:H_per^1+ε[0,2π]−→L²_per[0,2π], r7→p^∞,

where p^∞ = F[p^s] is the far-field of the scattered field p^s, and p^s is the pressure part of the solution (p^s, u) to the direct problem (1), (2), (4), (5), and (3) including the interface Γ = Γ^r and a fixed incoming plane wave p^inc. The spaceH_per^1+ε[0,2π] is the periodic Sobolev space of order 1 +ε >1 over the interval [0,2π]. L²[0,2π] is the corresponding Lebesgue space.

Now the inverse problem is the following: For a given far-field patternp^∞, find the shape of the obstacle with boundary rsol such that the scattered field corresponding to the fixed incoming plane wavep^inc has the far-field patternp^∞, i.e., such thatF(rsol) =p^∞. To our knowledge, results on the uniqueness of the solutionrsol are not known yet. For the case of far-field

(10)

data given in all directions of incidence, we refer to the theoretical results in [18], [17].

We now define three different optimization problems equivalent to the inverse problem. Numerical algorithms for the inverse problem can be de- rived simply by applying numerical minimization schemes to the optimization problems. More precisely, the minimization schemes are applied to regularized modifications of the optimization problems.

The first optimization problem is to find a least-squares solution rmin, i.e., a minimizer of the following problem

inf

r∈Hper^1+ε[0,2π]

J_γ(r), J_γ(r) :=°

°F(r)−p^∞°

°²

L²[0,2π].

Since the inverse problem is ill-posed and since the measured far-field data is given with noise, we replace the last optimization problem by

inf

r∈Hper^1+ε[0,2π]J_γ¹(r), J_γ¹(r) :=°

°F(r)−p^∞_noisy°

°²

L²[0,2π]+γkrk²_H1+ε

per[0,2π], (18) where γ is a small positive regularization parameter. As usual this parameter is to be chosen in dependence on the noise level. To guarantee convergence for noise level tending to zero and forγ→0, we suppose

°°p^∞−p^∞_noisy°

°²

L²[0,2π]≤cγ (19)

for a constant c independent of γ. The first numerical algorithm (cf. [4]) consists now in discretizing the mappingF by finite elements and applying a Gauss–Newton method to determine a minimizer of (18). This is a modified Newton method for the operator equationF(rsol) =p^∞_noisy which we shall call the simple Newton iteration.

Theorem 3.1([4]). SupposeΓ0 is chosen such that the corresponding inte- rior domain has no Dirichlet eigenvalue equal tok²_ωfor the negative Lapla- cian. Then we have:

(i) For any γ >0, there is a minimizer r^γ of (18).

(ii) Suppose the far-field patternp^∞ is the exact pattern for a fixed so- lution r^∗ of the inverse problem, i.e.,F(r^∗) =p^∞ and J₀¹(r^∗) = 0.

Then, for ε > 0 and for any set of minimizers r^γ, there exists a subsequence r^γⁿ converging weakly in H_per^1+ε[0,2π] and strongly in H_per^1+ε⁰[0,2π], 0< ε⁰ < ε, to a solution r^∗∗ of (18) withγ = 0 and, therewith, to a solution of the inverse problem.

(iii) If, additionally to the assumptions of (iii), the solution r^∗ of the inverse problem is unique, then we even get that r^γ tends to r^∗ weakly inH_per^1+ε[0,2π] and strongly inH_per^1+ε⁰[0,2π],0< ε⁰ < ε.

Unfortunately, for the first method the computation ofF requires a solution of a direct problem. In particular, if the curve Γ is the boundary of a domain with Jones frequency or close to such a boundary, the direct solution by finite elements is not easy. One way would be to compute with slightly modified frequencies. However, it might be difficult to check whether the

(11)

curve is “close” to the boundary of a domain with Jones frequency and to choose a modified frequency appropriately.

In order to motivate the second numerical method, the Kirsch–Kress algorithm, which corresponds to a third optimization problem, we introduce a second intermediate optimization method first. The plan is to define a method, where a solution of the direct method is not needed. Therefore, besides the unknown curve Γ the pressure p^s and the displacement field uare included into the set of optimization “parameters”. Additionally to the term of the least squares deviation of F[p^s] from p^∞_noisy, new terms are needed which enforce the fulfillment of the equations (1), (2), (4), (5), and (3) at least approximately. Hence, the regularized second optimization problem is to find a minimizer (rmin, umin, pmin) of

inf

r∈Hper^2+ε[0,2π], u∈[H¹(Ω)]², p^s∈H¹(ΩR)J_γ²(r, u, p^s), (20) J_γ²(r, u, p^s) :=°

°F[p^s]−p^∞_noisy°

°²

L²[0,2π])+°

°∆^∗u+%ω²u°

°²

[H⁻¹(Ω)]²+ +°

°∆p^s+kw2p^s°

°²

H⁻¹(ΩR)+ +°

°t[u] +{p^s+p^inc}ν°

°²

[H^−1/2(Γ)]²+ +

°°

°u·ν− 1

%fω² n∂p^s

∂ν +∂p^inc

∂ν o°°

°²

H^−1/2(Γ)+ +

°°

°V_Γ^ac₀[∂νp^s] +1

2[p^s]−K_Γ^ac₀[p^s]

°°

°²

H^1/2(Γ0)+ +c1γkrk²_H2+ε

per[0,2π]+c2γkuk²_H1(Ω)+c3γkp^sk²_H1(ΩR), K_Γ^ac₀[p^s](x) :=

Z

Γ0

∂G(x, y;k_w)

∂ν(y) p^s(y) dΓ0y

where ci > 0, i = 1,2,3, are calibration constants and γ is a small positive regularization parameter. Of course, this is a theoretical optimization problem only. For a numerical realization, the operators should be replaced by those of the variational formulation. However, it is clearly seen that the price for avoiding a solution of the direct problem is an increase in the number of the optimization “parameters”. The numerical solution of the discretized optimization problem (20) is higher dimensional and might be more involved than that for the case of (18).

The third optimization problem is a modification of (20). The optimization “parameters”u andp^s are replaced by the layer functions ϕi and ϕ~e

of the potential representations (8). In other words, in the numerical discretization the finite elements over the domains Ω and ΩR are replaced by lower dimensional boundary elements over the curves Γiand Γe. Instead of the terms inJ_γ² enforcing the conditions (1), (2), (4), (5), and (3), we only

(12)

need terms enforcing (12) and (13). Hence, the regularized third optimization problem is to find a minimizer (rmin, ϕi,min, ~ϕe,min) of

inf

r∈Hper^2+ε[0,2π], ϕi∈H⁻¹(Γi), ~ϕe∈[H⁻¹(Γe)]²

J_γ³(r, ϕ_i, ~ϕ_e), (21) J_γ³(r, ϕi, ~ϕe) :=c

°°

°F£ V_Γ^ac_iϕi

¤−p^∞_noisy

°°

°²

L²[0,2π]+ +γkϕik²_H−1(Γi)+γkϕ~ek²_[H−1(Γe)]²+ +

°°

°t£ V_Γ^el_eϕ~e

¤+£ V_Γ^ac_iϕi

¤ν+p^incν

°°

°²

L²(Γ^r)+ +

°°

°%fω²ν·£ V_Γ^el_iϕ~e

¤−∂ν

£V_Γ^ac_iϕi

¤−∂νp^inc

°°

°²

L²(Γ^r), (22) where γ is a small positive regularization parameter and c a positive calibration constant. We choose the layers ϕi,min and ϕ~e,min in an unusual Sobolev space of negative order to enable approximations by Dirac-delta functionals, i.e., by the method of fundamental solutions. Though the number of optimization parameters in a discretization of (21) is larger than that in a discretization of (18), the objective functionalJ_γ³ is simpler thanJ_γ¹. Applying an optimization scheme like the conjugate gradient method or the Levenberg–Marquardt algorithm to (21), we arrive at the Kirsch–Kress method. Note that the accuracy of the solution of this method is limited by the accuracy of solving the integral equations (12) and (13) with a Tikhonov regularization. To improve this, the curves Γi and Γe can be updated during the iterative solution of the optimization problem (compare the iterative schemes in [21], [24], [15]).

Theorem 3.2([5]). Supposek²_ω is not a Dirichlet eigenvalue for the nega- tive Laplacian in the interior ofΓiand thatp^∞is the exact far-field pattern of a scattered fieldp^s corresponding to some Γ^r^∗. Then we have:

(i) For any γ >0, there is a minimizer (r^γ, ϕ^γ_i, ~ϕ^γ_e)of(21).

(ii) For any set of minimizers (r^γ, ϕ^γ_i, ~ϕ^γ_e), there exists a subsequence (r^γⁿ, ϕ^γ_iⁿ, ~ϕ^γ_eⁿ) such that r^γⁿ converges weakly in H_per^1+ε[0,2π] and strongly inH_per^1+ε⁰[0,2π],0< ε⁰< ε, to a solution r^∗∗ of the inverse problem.

(iii) If, additionally, the solution r^∗ of the inverse problem is unique, then we even get that r^γ tends to r^∗ weakly in H_per^1+ε[0,2π] and strongly inH_per^1+ε⁰[0,2π], 0< ε⁰< ε.

Formulas for the discretization of the optimization problem (21) and for the derivatives of the objective functional are presented in Section 6.

4. Some Details of the Implementation

For the solution of the optimization problems, a lot ofnumerical optimization schemesare available (cf. [20]). Unfortunately, global methods

(13)

which yield the global minimum are often very slow. We recommend gradient based local optimization schemes. They provide local minimizers, i.e. solutions with minimal value of the objective functional in a neighbour- hood of the minimizer. In general, it cannot be guaranteed that the local minimizer is the global minimizer. However, using a good initial guess, the local minimizer will coincide with the global. In particular, we have tested the Gauss–Newton method, the Levenberg–Marquardt algorithm (cf. [14]), and the conjugate gradient method. The last method has been tested for the Kirsch–Kress method to avoid the solution of linear systems in the size of the direct problem.

In order to compute derivatives of the objective functionals in case of the simple Newton iteration the calculus of shape derivatives can be applied.

The derivatives result from solving the finite element system of the direct problem with new right-hand side vectors. This is fast if the finite element system is solved by an LU factorization for sparse systems (cf. [22], [4]).

The derivatives for the Kirsch–Kress method can be obtained by a simple differentiation of the kernel functions in the potential representations. Since the elasticity kernel contains second-order derivatives of the acoustic kernel and since the terms enforcing the transmission conditions contain first-order derivatives of the elastic potential, we need fourth order derivatives of the acoustic kernel. We present the needed formulas in Section 6.

Normally, quadrature rules are needed if the layer functions ϕi and

~

ϕe in the potential representation (8) are approximated by functions of a finite dimensional space. The potential integrals of these functions must be approximated by appropriate quadratures. However, in the case of the Kirsch–Kress method we can approximate the layer functions by linear com- binations of Dirac delta functions

ϕi ∼ϕi,M :=

XM

κ=1

bκδx_i,κ, bκ∈C, xi,κ:=rie^it^κ, tκ:= 2πκ

M , (23)

~

ϕe∼ϕ~e,M :=

XM

κ=1

cκδxe,κ, cκ∈C², xe,κ:=ree^it^κ. (24) This works since the potential operators are smoothing operators from the curve Γe, Γito Γ. Only in the case that Γeor Γiis close to Γ, a trigonometric or spline approximation ofϕiandϕ~etogether with an accurate quadrature must be employed.

Another important issue is the scalingof the optimization scheme. In- deed, the number of necessary iterations depends on the conditioning of the optimization problem. Using an appropriate scaling, the conditioning can be essentially improved. The first choice is, of course, the natural scaling.

The far-field values should be scaled such that the measurement uncertain- ties of the scaled far-field values coincide, and the parameters should be scaled in accordance with the accuracy requirements. A scaling different from the natural one is chosen not to improve the reconstruction operator,

(14)

but to speed up the optimization algorithm. This calibration may include different constants in front of the individual terms in the objective functional (cf. the factorsc andγ in the definition ofJ_γ³) and the replacement of the optimization parameters by the products of these parameters with convenient constants. The constants can be chosen, e.g., to minimize the conditioning of the Jacobian of the mapping that maps the parameters to the far-field values. Alternatively, the constants can be chosen by checking typical test examples with known solution. To improve the conditioning of the optimization in the Kirsch–Kress method, we have replaced the “optimization parameters”r, ϕi, andϕ~e by the parameters

r⁰=r/cr, ϕi0=ϕi/ci, ϕ~⁰_e=ϕ~e/ce. (25) 5. Numerical Results

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

Figure 6. Initial solution and egg shaped domain.

5.1. The curves for the numerical examples and some technical details. We have employed (i) the simple Newton iteration and (ii) the

(15)

Kirsch–Kress method, both with a circle as initial solution, to reconstruct two different obstacles. The first is an easy egg shaped domain (cf. Figure 6) with a boundary given by (17), byr_i= 2,r_e= 6, and by the fast decaying Fourier coefficients

ba0= 0,

ba1=−1, ba2 = 0.1, ba3 = 0.01, ba4 = −0.001, ba5= 0.0001, bb1= 1, bb2 = 0.1, bb3 = 0.01, bb4 = 0.001, bb5= 0.0001.

(26) The second body is the nonconvex obstacle from the end of Section 2 (cf. Fig- ure 7), and its boundary is given by ri = 2, re = 6 and by the Fourier coefficients

b a0= 0, b

a1= 1, ba2 = 0.10, ba3 = 0.040, ba4 = 0.016, ba5=0.008, bb1=−1, bb2 = 0.02, bb3 = −1.500, bb4 = −0.010, bb5=0.008.

(27) Clearly, both obstacles are defined by a truncated Fourier series and are analytic. However, the egg shaped domain is smoother since the nonzero Fourier coefficients have the strong decay property |baj| ≤10^−j and |bbj| ≤ 10^−j. More precisely, the norms

krkr:=

vu

ut|ba0|²+1 2

X∞

j=1

r^−2j|baj|²+1 2

X∞

j=1

r^−2j|bbj|², r >1,

of analytic functions are smaller for the egg shaped domain than for the nonconvex example. Note thatkrkr is the norm

vu ut|ba0|²+

X∞

j=1

r^−2j

¯¯

¯baj−ibbj

2

¯¯

¯²+ X∞

j=1

r^−2j

¯¯

¯baj+ibbj

2

¯¯

¯², of the analytic extension

z=%e^it7→ ba0+ X∞

j=1

h baj−ibbj

2 i

%^je^ijt+ X∞

j=1

h baj+ibbj

2 i

%^−je^−ijt of the function e^it 7→ r(t) = P

j

b

ajcos(jt) +P

j

bbjsin(jt) onto the annular domain{z∈C: 1/r <|z|< r}.

In all computations, we have chosen the physical constants in accordance with (15). The incoming plane wave has been fixed top^inc(x) :=e^i(1,0)^>^·x. Moreover, for all initial curves and all iterative solutions, we have fixed the zeroth Fourier coefficient ba0 to zero. The “measured” far-field data {p^∞(k/M⁰⁰), k = 1, . . . , M⁰⁰}, M⁰⁰ = 80 (cf. Figure 5) has been simulated by the piecewise linear finite element method (FEM) described in Section 2.

To avoid what is called an inverse crime, we have chosen the meshsize of the FEM grid (determined by NETGEN [23]) for the far-field computation by a factor of at least 0.25 smaller than that of the FEM grids involved in the

(16)

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

Figure 7. Initial solution and nonconvex domain.

inverse algorithms. Our tests have revealed that the far-field of the FEM method is more reliable than that computed by the regularized system (12)–

(13). The scaling parametersc, cr, ci, andcefor the Kirsch–Kress method (cf. (25) and the definition of J_γ³) have been determined experimentally such that the reconstruction by Gauss–Newton iteration converges with the smallest number of iteration steps. For example, for the egg shaped domain andM =M⁰ = 44 points of discretization for the approximate integration over Γ, Γi, Γe(cf. the discretized objective functional in (A.6)), these values arec= 4000,cr= 1, ci= 0.1, andce= 0.005.

5.2. Convergence of the simple Newton iteration. The results for the egg shaped domain and for the simple Newton iteration have been similar to those presented in [4], where the constants where slightly different and the obstacle was similar to our nonconvex body. After a small number (≤20) of iterations, the algorithm reconstructs the obstacle. The regularization parameter γ can even be set to zero, which is not surprising since only

(17)

10 unknown real parameters are reconstructed from 160 real measurement values. The left Table 1 exhibits the meshsize h, the number of Gauss–

Newton iterations it, and the accuracy err := ker−er_{F EM}k_L∞[0,2π] of the reconstruction rF EM. The first row contains the accuracy of the initial guess. For the nonconvex obstacle, the results are similar (cf. right Table 1). Most of the computing time is spent on the evaluation of the objective functional including the solution of a direct problem. Therefore, it is not necessary to replace the expensive Gauss–Newton iteration by a different optimization scheme.

h err it

1.2596 0

0.5 0.0759 6

0.25 0.0247 8

0.125 0.00876 8 0.0625 0.00329 10 0.03125 0.00156 10

h err it

1.5733 0 0.25 1.1435 20 0.125 0.00924 17 0.0625 0.00401 15 0.03125 0.00157 18 Table 1. Reconstruction by simple Newton iteration for egg shaped domain (left) and for nonconvex domain (right).

5.3. Convergence of the Kirsch–Kress algorithm. We have started the tests of the Kirsch–Kress method with the nonconvex domain. However, the optimization algorithms did not converge. To fix the problem, we have checked the solution of the corresponding direct problem. We have observed that the far-field of the solution computed by (8), (12), and (13) did not match that of the FEM. Even a Tikhonov regularization in accordance with the last four terms of the functionalJ_γ³did not help. Only a regularization with a truncated singular value decomposition and a well-chosen truncation parameter led to the correct far field. In other words, the reason for the divergence of the Kirsch–Kress method is the high degree of ill-posedness of the system (12), (13). On the other hand, if we commit the inverse crime and take the incorrect far-field data computed by solving (12), (13), then the Kirsch–Kress algorithm does converge.

To show the convergence of the Kirsch–Kress method with FEM generated far-field data, we consider the egg shaped domain. This time the solution curve has a higher degree of smoothness, and the direct solution of (12), (13) together with a Tikhonov regularization yields a far-field solution close to that of the FEM. Table 2 shows that the Kirsch–Kress method converges for the egg shaped domain. Indeed, the table shows the regularization parameter γ, the error ker−erKKkL^∞_per[0,2π] of the Kirsch–Kress reconstruction rKK, and the number of necessary iteration steps. These depend on the number of discretization points M = M⁰ for the approximate integration over Γ, Γi, Γe (cf. the discretized objective functional in

(18)

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

x non-convex domain

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

x non-convex domain

Figure 8. Initial solution with Fourier coefficients ba¹_i,bb¹_i and nonconvex domain with modified curves Γi and Γe.

(A.6)) and on the choice of the optimization method. In particular, we have checked the Gauss–Newton method with experimentally chosen regularization parameterγ (GNw), the Levenberg–Marquardt method with the same regularization parameter (LMw), and the Levenberg–Marquardt method without regularization (LMo). The results show much better approximations than for the simple Newton iteration. Unfortunately, the conjugate gradient method did not converge.

To get convergence of the Kirsch–Kress method also for the nonconvex domain of Figure 7, we have changed the curves Γi and Γe (cf. Figure 8).

If these are closer to the curve Γ^r, then the degree of ill-posedness of the operators in (12), (13) is reduced. We have chosen the initial guess of the Fourier coefficients as

b

a⁰₀= 0.0, b

a⁰₁= 1.3, ba⁰₂=−0.10, ba⁰₃= 0.1, ba⁰₄=−0.05, ba⁰₅= 0.018, bb⁰₁=−0.8, bb⁰₂= 0.05, bb⁰₃=−1.7, bb⁰₄= 0.03, bb⁰₅=−0.020.

(28)

(19)

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

x domain

curve in.curve ex.curve

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

x domain

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

x domain

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

x domain

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

x domain

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

x domain

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

x domain

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

x domain

Figure 9. Kirsch–Kress steps 1-4 to reconstruct nonconvex domain.

(20)

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

x domain

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

x domain

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

x domain

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

x domain

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

x domain

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

x domain

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

x domain

-6 -4 -2 0 2 4 6

-8 -6 -4 -2 0 2 4 6 8

y

x domain

Figure 10. Kirsch–Kress steps 5-8 to reconstruct nonconvex domain.