Inverse Obstacle Scattering for Elastic Waves in Three Dimensions

Consider an exterior problem of the three-dimensional elastic wave equation, which models the scattering of a time-harmonic plane wave by a rigid obstacle. The scattering problem is reformulated into a boundary value problem by introducing a transparent boundary condition. Given the incident field, the direct problem is to determine the displacement of the wave field from the known obstacle; the inverse problem is to determine the obstacle's surface from the measurement of the displacement on an artificial boundary enclosing the obstacle. In this paper, we consider both the direct and inverse problems. The direct problem is shown to have a unique weak solution by examining its variational formulation. The domain derivative is studied and a frequency continuation method is developed for the inverse problem. Numerical experiments are presented to demonstrate the effectiveness of the proposed method.

associated spherical harmonics when studying the transparent boundary condition. Computationally, it is also more intensive.
The rigid obstacle is assumed to be embedded in an open space filled with a homogeneous and isotropic elastic medium. The scattering problem is reduced into a boundary value problem by introducing a transparent boundary condition on a sphere. We show that the direct problem has a unique weak solution by examining its variational formulation. The proofs are based on asymptotic analysis of the boundary operators, the Helmholtz decomposition, and the Fredholm alternative theorem.
The calculation of domain derivatives, which characterize the variation of the wave field with respect to the perturbation of the boundary of an medium, is an essential step for inverse scattering problems. The domain derivatives have been discussed by many authors for the inverse acoustic and electromagnetic obstacle scattering problems [10,16,32]. Recently, the domain derivative is studied in [20] for the elastic wave by using boundary integral equations. Here we present a variational approach to show that it is the unique weak solution of some boundary value problem. We propose a frequency continuation method to solve the inverse problem. The method requires multi-frequency data and proceed with respect to the frequency. At each frequency, we apply the descent method with the starting point given by the output from the previous step, and create an approximation to the surface filtered at a higher frequency. Numerical experiments are presented to demonstrate the effectiveness of the proposed method. A topic review can be found in [4] for solving inverse scattering problems with multi-frequencies to increase the resolution and stability of reconstructions.
The paper is organized as follows. Section 2 introduces the formulation of the obstacle scattering problem for elastic waves. The direct problem is discussed in section 3 where well-posedness of the solution is established. Section 4 is devoted to the inverse problem. The domain derivative is studied and a frequency continuation method is introduced for the inverse problem. Numerical experiments are presented in section 5. The paper is concluded in section 6. To avoid distraction from the main results, we collect in the appendices some necessary notation and useful results on the spherical harmonics, functional spaces, and transparent boundary conditions.
2. Problem formulation. Consider a bounded and rigid obstacle D ⊂ R 3 with a Lipschitz boundary ∂D. The exterior domain R 3 \D is assumed to be filled with a homogeneous and isotropic elastic medium, which has a unit mass density and constant Lamé parameters λ, µ satisfying µ > 0, λ + µ > 0. Let B R = {x ∈ R 3 : |x| < R}, where the radius R is large enough such thatD ⊂ B R . Define Γ R = {x ∈ R 3 : |x| = R} and Ω = B R \D.
Since the obstacle is elastically rigid, we have (2.4) u = 0 on ∂D.
The total field u consists of the incident field u inc and the scattered field v: Subtracting (2.2) from (2.3) yields that v satisfies For any solution v of (2.5), we introduce the Helmholtz decomposition by using a scalar function φ and a divergence free vector function ψ: Substituting (2.6) into (2.5), we may verify that φ and ψ satisfy (2.7) ∆φ + κ 2 p φ = 0, ∆ψ + κ 2 s ψ = 0.
In addition, we require that φ and ψ satisfy the Sommerfeld radiation condition: Using the identity ∇ × (∇ × ψ) = −∆ψ + ∇(∇ · ψ), we have from (2.7) that ψ satisfies the Maxwell equation: It can be shown (cf. [8,Theorem 6.8]) that the Sommerfeld radiation for ψ in (2.8) is equivalent to the Silver-Müller radiation condition: Given u inc , the direct problem is to determine u for the known obstacle D; the inverse problem is to determine the obstacle's surface ∂D from the boundary measurement of u on Γ R . Hereafter, we take the notation of a b or a b to stand for a ≤ Cb or a ≥ Cb, where C is a positive constant whose specific value is not required but should be clear from the context.

Direct scattering problem.
In this section, we study the variational formulation for the direct problem and show that it admits a unique weak solution.
3.1. Transparent boundary condition. We derive a transparent boundary condition on Γ R . Given v ∈ L 2 (Γ R ), it has the Fourier expansion: where {(T m n , V m n , W m n ) : n = 0, 1, . . . , m = −n, . . . , n} is an orthonormal system in L 2 (Γ R ) and v m jn are the Fourier coefficients of v on Γ R . Define a boundary operator which is assumed to have the Fourier expansion: Taking ∂ r of v in (C.6), evaluating it at r = R, and using the spherical Bessel differential equations [34], we get where z n (t) = th n is the spherical Hankel function of the first kind with order n, φ m n and ψ m jn are the Fourier coefficients for φ and ψ on Γ R , respectively. Noting (C.6) and using ∇ · v = ∆φ = 2 (3.5) Comparing (3.2) with (3.5), we have where the matrix 32 = µ n(n + 1)(z n (κ s R) − 1). Here where Λ n = z n (κ p R)(1 + z n (κ s R)) − n(n + 1). Using the above notation and combining (3.6) and (C.10), we derive the transparent boundary condition: Proof. Using the asymptotic expansions of the spherical Bessel functions [34], we may verify that z n (t) = −(n + 1) + 1 16n It follows from straightforward calculations that = n(n + 1) A simple calculation yields which completes the proof by applying Sylvester's criterion. Lemma 3.2. The boundary operator T : Proof. For any given u ∈ H 1/2 (Γ R ), it has the Fourier expansion Let u m n = (u m 1n , u m 2n , u m 3n ) . It follows from (3.7) and the asymptotic expansions of M (n) ij that which completes the proof.

Well-posedness.
Using the transparent boundary condition (3.7), we obtain a boundary value problem for u: where the sesquilinear form b : Here A : B = tr(AB ) is the Frobenius inner product of square matrices A and B.
The following result follows from the standard trace theorem of the Sobolev spaces. The proof is omitted for brevity.
Lemma 3.4. It holds the estimate Lemma 3.5. For any ε > 0, there exists a positive constant C(ε) such that Proof. Let B be the ball with radius R > 0 such thatB ⊂ D. DenoteΩ = B\B . Given u ∈ H 1 ∂D (Ω), letũ be the zero extension of u from Ω toΩ, i.e., The extension ofũ has the Fourier expansioñ A simple calculation yields Sinceũ(R , θ, ϕ) = 0, we haveũ m jn (R ) = 0. For any given ε > 0, it follows from Young's inequality that The proof is completed by noting that Lemma 3.6. It holds the estimate Proof. As is defined in the proof of Lemma 3.5, letũ be the zero extension of u from Ω toΩ. It follows from the Cauchy-Schwarz inequality that Hence we have The proof is completed by noting that Theorem 3.7. The variational problem (3.14) admits a unique weak solution u ∈ H 1 ∂D (Ω). Proof. Using the Cauchy-Schwarz inequality, Lemma 3.2, and Lemma 3.4, we have which shows that the sesquilinear form b(·, ·) is bounded.
It follows from Lemma 3.1 that there exists an N 0 ∈ N such thatM n is positive definite for n > N 0 . The sesquilinear form b can be written as Taking the real part of b, and using Lemma 3.1, Lemma 3.6, Lemma 3.5, we obtain Letting ε > 0 to be sufficiently small, we have C 1 − C 2 ε > 0 and thus Gårding's inequality. Since the injection of H 1 ∂D (Ω) into L 2 (Ω) is compact, the proof is completed by using the Fredholm alternative (cf. [31,Theorem 5.4.5]) and the uniqueness result in Theorem 3.3.
4. Inverse scattering. In this section, we study a domain derivative of the scattering problem and present a continuation method to reconstruct the surface.
4.1. Domain derivative. We assume that the obstacle has a C 2 boundary, i.e., ∂D ∈ C 2 . Given a sufficiently small number h > 0, define a perturbed domain Ω h which is surrounded by ∂D h and Γ R , where Here the function p ∈ C 2 (∂D).
Consider the variational formulation for the direct problem in the perturbed do- Similarly, we may follow the proof of Theorem 3.7 to show that the variational problem (4.1) has a unique weak solution u h ∈ H 1 ∂D h (Ω h ) for any h > 0. Since the variational problem (3.7) is well-posed, we introduce a nonlinear scattering operator: which maps the obstacle's surface to the displacement of the wave field on Γ R . Let u h and u be the solution of the direct problem in the domain Ω h and Ω, respectively. Define the domain derivative of the scattering operator S on ∂D along the direction p as For a given p ∈ C 2 (∂D), we extend its domain toΩ by requiring that p ∈ , J η h and J ξ h are the Jacobian matrices of the transforms η h and ξ h , respectively.
For a test function v h in the domain Ω h , it follows from the transform thatv is a test function in the domain Ω. Therefore, the sesquilinear form b h in (4.2) becomes which gives an equivalent variational formulation of (4.1): Here I is the identity matrix. Following the definitions of the Jacobian matrices, we may easily verify that where the matrix J p = ∇p.
Substituting the above estimates into (4.3)-(4.5), we obtain Hence we have Theorem 4.1. Given p ∈ C 2 (∂D), the domain derivative of the scattering operator S is S (∂D; p) = u | Γ R , where u is the unique weak solution of the boundary value problem:

7)
and u is the solution of the variational problem (3.14) corresponding to the domain Ω.
Proof. Given p ∈ C 2 (∂D), we extend its definition to the domainΩ as before. It follows from the well-posedness of the variational problem (3.14) thatȗ → u in H 1 ∂D (Ω) as h → 0. Taking the limit h → 0 in (4.6) gives which shows that (ȗ − u)/h is convergent in H 1 ∂D (Ω) as h → 0. Denote the limit bẏ u and rewrite (4.8) as First we compute g 1 (p)(u, v). Noting p = 0 on ∂B and using the identity we obtain from the divergence theorem that Noting µ∆u + (λ + µ)∇∇ · u + ω 2 u = 0 in Ω, we have from the integration by parts that Using the integration by parts again yields Let τ 1 (x), τ 2 (x) be any two linearly independent unit tangent vectors on ∂D. Since u = v = 0 on ∂D, we have Using the identities which gives ∂D (p · ∇v) · (ν · ∇u) − (p · ν)(∇u : ∇v) dγ = 0.

Hence we get
Combining the above identities gives Define u =u − p · ∇u. It is clear to note that p · ∇u = 0 on Γ R since p = 0 on Γ R . Hence, we have which shows that u is the weak solution of the boundary value problem (4.7). To verify the boundary condition of u on ∂D, we recall the definition of u and have fromȗ = u = 0 on ∂D that Noting u = 0 on ∂D, we have (4.13) p · ∇u = (p · ν)∂ ν u, which completes the proof by combining (4.12) and (4.13).
Consider the objective function The inverse problem can be formulated as the minimization problem: In order to apply the descend method, we have to compute the gradient of the objective function: We have from Theorem 4.2 that We assume that the scattering data U is available over a range of frequencies ω ∈ [ω min , ω max ], which may be divided into ω min = ω 0 < ω 1 < · · · < ω J = ω max . We now propose an algorithm to reconstruct the Fourier coefficients c i , i = 1, . . . , 6(N + 1) 2 .
Algorithm: Frequency continuation algorithm for surface reconstruction. Denote C (1) k0 = (c 1 , c 2 , . . . , c 6(k0+1) 2 ) and consider the iteration: where τ > 0 and L > 0 are the step size and the number of iterations for every fixed frequency, respectively. 3. Continuation: increase to ω 1 , let k 1 = [ω 1 ], repeat Step 2 with the previous approximation to r j,N as the starting point. More precisely, approximate r j,N by and determine the coefficientsc i , i = 1, . . . , 6(k 1 + 1) 2 by using the descent method starting from the previous result. 4. Iteration: repeat Step 3 until a prescribed highest frequency ω J is reached.

Numerical experiments.
In this section, we present two examples to show the effectiveness of the proposed method. The scattering data is obtained from solving the direct problem by using the finite element method with the perfectly matched layer technique, which is implemented via FreeFem++ [12]. The finite element solution is interpolated uniformly on Γ R . To test the stability, we add noise to the data: where rand are uniformly distributed random numbers in [−1, 1] and δ is the relative noise level, x k are data points. In our experiments, we pick 100 uniformly distributed points x k on Γ R , i.e., K = 100.
In the following two examples, we take λ = 2, µ = 1, R = 1. The radius of the initial R 0 = 0.5. The noise level δ = 5%. The step size in (4.15) is τ = 0.005/k i where k i = [ω i ]. The incident field is taken as a plane compressional wave.
Example 2. Consider a cushion-shaped obstacle: where r(θ, ϕ) = (0.75 + 0.45(cos(2ϕ) − 1)(cos(4θ) − 1)) 1/2 . Figure 5.2(a) shows the exact surface. This example is much more complex than the bean-shaped obstacle due to its multiple concave parts. Multiple incident directions are needed in order to obtain a good result. In this example, the obstacle is illuminated by the compressional wave from 6 directions, which are the unit vectors pointing to the origin from the face centers of the cube. The multiple frequencies are the same as the first example, i.e., the frequency ranges from ω min = 1 to ω max = 5 with ω i = i + 1, i = 0, . . . , 4. For each fixed frequency and incident direction, repeat L = 50 times with previous result as starting points. The step size for the decent method is 0.005/ω i and number of recovered coefficients is 6(ω i +2) 2 for corresponding frequency. 6. Conclusion. In this paper, we have studied the direct and inverse obstacle scattering problems for elastic waves in three dimensions. We develop an exact transparent boundary condition and show that the direct problem has a unique weak solution. We examine the domain derivative of the total displacement with respect to the surface of the obstacle. We propose a frequency continuation method for solving the inverse scattering problem. Numerical examples are presented to demonstrate the effectiveness of the proposed method. The results show that the method is stable and accurate to reconstruct surfaces with noise. Future work includes the surfaces of different boundary conditions and multiple obstacles where each obstacle's surface has a parametric equation. We hope to be able to address these issues and report the progress elsewhere in the future.
For a smooth scalar function u(R, θ, ϕ) defined on Γ R , let be the tangential gradient on Γ R . Define a sequence of vector spherical harmonics: Let L 2 (Ω) = L 2 (Ω) 3 be equipped with the inner product and norm: Denote by H 1 (Ω) the standard Sobolev space with the norm given by It can be verified that H −s (Γ R ) is the dual space of H s (Γ R ) with respect to the inner product Introduce three tangential trace spaces: For any tangential field u ∈ H s t (Γ R ), it can be represented in the series expansion Using the series coefficients, the norm of the space H s t (Γ R ) can be characterized by the norm of the space H −1/2 (curl, Γ R ) can be characterized by n m=−n 1 1 + n(n + 1) |u m 1n | 2 + 1 + n(n + 1)|u m 2n | 2 ; the norm of the space H −1/2 (div, Γ R ) can be characterized by Given a vector field u on Γ R , denote by u Γ R = −e r × (e r × u) the tangential component of u on Γ R . Define the inner product in where v * is the conjugate transpose of v.
Appendix B. Transparent boundary conditions. Recall the Helmholtz decomposition (2.6): where the scalar potential function φ satisfies (2.7) and (2.8): as r → ∞, the vector potential function ψ satisfies (2.9) and (2.10): where r = |x| andx = x/r. In the exterior domain R 3 \B R , the solution φ of (B.1) satisfies where h (1) n is the spherical Hankel function of the first kind with order n and We define the boundary operator T 1 such that where z n (t) = th  Evaluating the derivative of (B.3) with respect to r at r = R and using (B.4), we get the transparent boundary condition for the scalar potential function φ: The following result can be easily shown from (B.4)-(B.5).