Extended sampling method for inverse elastic scattering problems using one incident wave

# Extended sampling method for inverse elastic scattering problems using one incident wave

J. Liu Department of Mathematical Sciences, Jinan University, Guangzhou, 130012, China (liujuan@jnu.edu.cn).    X. Liu NCMIS and Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China. (xdliu@amt.ac.cn).    J. Sun Department of Mathematical Sciences, Michigan Technological University, Houghton, MI 49931, U.S.A. (jiguangs@mtu.edu).
###### Abstract

We consider the inverse elastic scattering problems using the far field data due to one incident plane wave. A simple method is proposed to reconstruct the location and size of the obstacle using different components of the far field pattern. The method sets up linear ill-posed integral equations for sampling points in the domain of interrogation and uses the (approximate) solutions to compute indicators. Using the far field patterns of rigid disks as the kernels of the integral equations and moving the measured data to the right hand side, the method has the ability to process limited aperture data. Numerical examples show that the method can effectively determine the location and approximate the support of the obstacle with little a priori information.

## 1 Introduction

Motivated by applications in non-destructive testings, medical imaging and seismic exploration, the inverse elastic scattering problems have received significant attention recently [12, 3, 5, 1, 14, 18, 17, 4]. There exist two main groups of methods to recover the location and shape of an obstacle using the scattering data. The first group are the iterative methods to minimize some cost functions [4, 18]. These methods usually need to solve the forward scattering problems. The second group are the non-iterative or direct methods, e.g., the linear sampling method [8, 3], the factorization method [2, 6, 13], the reciprocity gap method [9, 20, 11], the range test method [21], the reverse time migration [7] and the direct sampling method [15, 17]. These methods do not solve any forward scattering problem and thus are fast in general.

In this paper, we consider the inverse elastic scattering problems using the far field pattern due to one incident plane wave. A simple method, called extended sampling method (ESM) [19], is proposed to reconstruct the location and size of the obstacle. The method is based on an idea similar to that of the linear sampling method, which has been studied extensively in the literature [8, 10].

The linear sampling method uses the full aperture far field pattern, i.e., the far field pattern of all incident and observation directions. It sets up linear ill-posed integral equations for the sampling points in the domain of interrogation and uses the (approximate) solutions to reconstruct the location and shape of the unknown obstacle. The kernel of the integral equations is the (measured) full aperture far field pattern. ESM also sets up integral equations for the sampling points in the domain of interrogation and uses the (approximate) solutions to compute some indicators. However, the kernels of these integral equations for ESM are the full aperture far field patterns of rigid disks, which can be computed easily in advance. The measured far field pattern is moved to the right hand side of the equation. This arrangement gives ESM the ability to process the far field pattern due to a single incident plane wave.

The rest of the paper is organized as follows. In Section 2, the direct elastic scattering problem is presented. In Section 3, we first derive a relation between the scattered field of the elastic scattering problems and the Sommerfeld radiation solutions of Helmholtz equations. This relation is used to transform the inverse problem with the compressional part or shear part of the far field pattern into an inverse acoustic scattering problem. Then a reconstruction algorithm is proposed based on ESM. In Section 4, using the far field patterns of elastic waves for rigid disks and the related translation property, a new far field equation is proposed. The behavior of the solutions is analyzed. An ESM algorithm for the inverse problem using the new equation is presented. In Section 5, numerical examples are provided.

## 2 Preliminaries and the direct scattering problem

We begin with the notations used throughout this paper. Vectors are written in bold to distinguish from scalars. For a vector , let and be obtained by rotating anticlockwise. Denote by for simplicity. In addition to the usual differential operators and , we will make use of and .

The propagation of time-harmonic waves in an isotropic homogeneous medium with Lame constants , (, ) and density is governed by the Navier equation

 Δ∗u+ρω2u=0, (2.1)

where denotes the displacement field and denotes the circular frequency. The differential operator . In this paper, we assume for simplicity.

The solution of (2.1) can be decomposed as

 u=up+us,

where

are known as the compressional part of associated with the wave number , and the shear part of associated with the wave number .

The direct elastic scattering problem for an obstacle is as follows. Given a bounded domain of class and an incident field such that and is a solution of (2.1) in a neighborhood of , find the scattered field such that

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩Δ∗u+ω2u=0,   in R2∖¯¯¯¯¯D,limr→∞√r(∂up/∂r−ikpup)=0,   r=|x|,limr→∞√r(∂us/∂r−iksus)=0,   r=|x|. (2.3)

It is well known that every radiation solution to the Navier equation has an asymptotic behavior of the form

 u(x)=eikp|x|√|x|u∞p(^x)^x+eiks|x|√|x|u∞s(^x)^x⊥+O(|x|−3/2),   |x|→∞, (2.4)

uniformly in all directions , where and are analytic functions on the unit circle

 S:={^x|^x∈R2,|^x|=1}.

Throughout the paper, is defined as the far field pattern of , and and are defined as compressional part and shear part of the far-field pattern, respectively. For the well-posedness of the above direct scattering problem, one needs to impose suitable conditions on , which depend on the physical properties of the scatterer. The scattered field satisfies

• the Dirichlet boundary condition

 u=−uinc   on ∂D,

for a rigid body;

• the Neumann boundary condition

 Tνu=−Tνuinc   on ∂D,

for a cavity, where denotes the surface traction operator and is the unit outward normal to ;

• the impedance boundary condition

 Tνu+iσu=−Tνuinc−iσuinc   on ∂D,

with some real-valued parameter .

The inverse scattering problem of interests is, using the far-field pattern of all observation directions due to one incident wave, to reconstruct the location and approximate the support of the scatterer without knowing the physical properties of the scatterer. More specifically, the following two inverse elastic scattering problems will be considered:

• Determine the location and size of the scatterer from the knowledge of compressional part or shear part of the far field pattern due to one incident wave.

• Determine the location and size of the scatterer from the knowledge of the far field pattern due to one incident wave.

## 3 Extended sampling method for IP-P

By building a relation between the scattered solution of elastic scattering problems and the radiating solutions of Helmholtz equations, we can extend ESM for inverse acoustic scattering problems to solve IP-P.

### 3.1 Radiating solutions of Helmholtz equations

Recall that is the scattered field of (2.3) and is the corresponding far field pattern. The following theorem shows that and are related with the radiating solutions and their far field patterns of some Helmholtz equations, respectively.

###### Theorem 3.1.

Let be a solution of (2.3). Then and are the radiating solutions of

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩Δϕ+k2pϕ=0,   in R2∖¯¯¯¯¯D,Δψ+k2sψ=0,   in R2∖¯¯¯¯¯D,limr→∞√r(∂ϕ/∂r−ikpϕ)=0,    r=|x|,limr→∞√r(∂ψ/∂r−iksψ)=0,   r=|x|. (3.1)

Furthermore, for the far field patterns and of and , respectively,

 ϕ∞(^x)=1ikpu∞p(^x),    ψ∞(^x)=1iksu∞s(^x). (3.2)
###### Proof.

It can be easily verified that and satisfy the first and second Helmholtz equations of (3.1), respectively.

From (2.4), the radiating solution of (2.3) in polar coordinates has an asymptotic behavior

 u(r,θ) = (u1(r,θ)u2(r,θ)) =

Consequently, we have that

 ∂∂ru1(r,θ)=ikpeikpr√ru∞p(θ)cosθ−ikseiksr√ru∞s(θ)sinθ+O(r−3/2),
 ∂∂ru2(r,θ)=ikpeikpr√ru∞p(θ)sinθ+ikseiksr√ru∞s(θ)cosθ+O(r−3/2).

Using these two equations, we obtain

 div u = ∂u1∂x1+∂u2∂x2 = (∂u1∂rcosθ−1r∂u1∂θsinθ)+(∂u2∂rsinθ+1r∂u2∂θcosθ) = (∂u1∂rcosθ+∂u2∂rsinθ)+(1r∂u2∂θcosθ−1r∂u1∂θsinθ) = ikpeikpr√ru∞p(θ)+O(r−3/2)

and

 div⊥ u = ∂u2∂x1−∂u1∂x2 = (∂u2∂rcosθ−1r∂u2∂θsinθ)−(∂u1∂rsinθ+1r∂u1∂θcosθ) = (∂u2∂rcosθ−∂u1∂rsinθ)−(1r∂u1∂θcosθ+1r∂u2∂θsinθ) = ikseiksr√ru∞s(θ)+O(r−3/2).

Thus

 ϕ(x)=−1k2pdiv u=1ikpeikp|x|√|x|u∞p(^x)+O(|x|−3/2), (3.3)
 ψ(x)=−1k2sdiv⊥u=1ikseiks|x|√|x|u∞s(^x)+O(|x|−3/2), (3.4)

which imply that and satisfy the Sommerfeld radiation conditions in (3.1).

Since and are radiating solutions of the Helmholtz equations, the far field patterns and have the following asymptotic expansions (see [10])

 ϕ(x)=eikp|x|√|x|{ϕ∞(^x)+O(1|x|)},   |x|→∞,
 ψ(x)=eiks|x|√|x|{ψ∞(^x)+O(1|x|)},   |x|→∞.

Then from (3.3) and (3.4), (3.2) is proved. ∎

###### Remark 3.2.

If the obstacle is a rigid body, the solution of the elastic scattering problem (2.3) satisfies the Dirichlet boundary condition . From (2.2), and satisfy

Denote by the unit tangent vector and by the unit outward normal vector on . By straightforward calculation, and satisfy the coupled boundary conditions

 ∂ϕ∂ν+∂ψ∂τ=−ν⋅uinc,   ∂ϕ∂τ−∂ψ∂ν=−τ⋅uinc,   on ∂D.

From Theorem 3.1, satisfies the following Helmholtz equations with coupled boundary conditions

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩Δϕ+k2pϕ=0,   in R2∖¯¯¯¯¯D,Δψ+k2sψ=0,   in R2∖¯¯¯¯¯D,∂ϕ∂ν+∂ψ∂τ=−ν⋅uinc,   on ∂D,∂ϕ∂τ−∂ψ∂ν=−τ⋅uinc,   on ∂D,limr→∞√r(∂ϕ/∂r−ikpϕ)=0,    r=|x|,limr→∞√r(∂ψ/∂r−iksψ)=0,   r=|x|, (3.6)

which has a unique solution (see [22]). On the other hand, from (2.2), if solves (3.6), direct computation shows that satisfies the elastic scattering problem (2.3) with Dirichlet boundary condition.

### 3.2 Extended sampling method

Denote by and the scattered field and far field pattern of the unknown scatterer due to an incident wave, respectively. The inverse problem IP-P is to determine the location and approximate support of the scatterer from the knowledge of the compressional part or the shear part . In this subsection, for simplicity we use to represent either or .

Let be a sound soft disk centered at with radius and let , solve

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩Δu+k2tu=0,   in R2∖¯¯¯¯¯¯Bz,u=−eiktx⋅d,   on ∂Bzlimr→∞√r(∂u/∂r−iktu)=0,   r=|x|,

where is the incident direction. Denote the far field pattern of (for its series expansion, see e.g., [10, 19]). Then for IP-P, we introduce the far field equation

 ∫SuBz∞(^x;kt,d)gz(d)ds(d)=1iktU∞t(^x),^x∈S. (3.7)

The following theorem is the main result for (3.7).

###### Theorem 3.3.

Let be a sound soft disk centered at with radius . Let be an obstacle. Assume that does not coincide with any zero of the Bessel functions . Then the following results hold for the far field equation (3.7):

• If , for a given , there exists a function such that

 ∥∥∥∫SuBz∞(^x;kt,d)gεz(d)ds(d)−1iktU∞t(^x)∥∥∥L2(S)<ε (3.8)

and the Herglotz wave function

 vgεz(x):=∫∂Bzeiktx⋅dgεz(d)ds(d),x∈Bz

converges to the solution of the Helmholtz equation with on as .

• If , every that satisfies (3.8) for a given is such that

 limε→0∥vgεz∥H1(Bz)=∞.
###### Proof.

Denote the right hand side of the far field equation (3.7) by . From Theorem 3.1, we know that is the far field pattern of

 ϕt(x):=⎧⎪⎨⎪⎩−1k2p% div U,   t=p,−1k2sdiv⊥U,   t=s,

which is a radiating solution of the Helmholtz equation outside . The rest of proof is exactly the same as the that of Theorem 3.1 in [19]. ∎

Now we are ready to present the extended sampling method (ESM) for IP-P. Let be a domain containing . For a sampling point , we consider the linear ill-posed integral equation (3.7). By Theorem 3.3, one expects that is relatively large when is outside and relatively small when is inside . Consequently, a reconstruction of the location and support of can be obtained by plotting for all sampling points .

• The Extended Sampling Method for IP-P

• Generate a set of sampling points for which contains .

• For each sampling point ,

• compute and set up a discrete version of (3.7);

• use the Tikhonov regularization to compute an approximate solution of (3.7).

• Find the global minimum point for .

• Choose to be the reconstruction for .

As in [19], one can use a multilevel technique to find a suitable radius of the sampling disks .

• The Multilevel ESM for IP-P

• Choose the sampling disks with a large radius . Generate a proper set of sampling points. Using ESM, determine the global minimum point for and an approximation for .

• For

• Let and generate a proper set of sampling points.

• Find the minimum point for and an approximation for . If , go to Step 3.

• Choose and to be the location and approximate support of , respectively.

## 4 Extended sampling method for IP-F

In this section, a novel far field equation is introduced. The series expansion for the far field patterns of rigid disks, which serve as the kernels of the integrals, will be studied in detail. Then an extended sampling method for IP-F is proposed.

### 4.1 Far field pattern for rigid disks

Let be a rigid disk centered at the origin with radius . Let and be the radiating solution and far field pattern of the elastic scattering problem (2.3) for due to an incident wave , respectively. From Remark 3.2, we know that , where is the unique solution of

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩Δϕ+k2pϕ=0,   in R2∖¯¯¯¯B,Δψ+k2sψ=0,   in R2∖¯¯¯¯B,∂ϕ∂ν+∂ψ∂τ=−ν⋅uinc,   on ∂B,∂ϕ∂τ−∂ψ∂ν=−τ⋅uinc,   on ∂B,limr→∞√r(∂ϕ/∂r−ikpϕ)=0,    r=|x|,limr→∞√r(∂ψ/∂r−iksψ)=0,   r=|x|.

Denote by the Hankel function of the first kind of order . In polar coordinates , and can be written as (see [10])

 ϕ(r,θ)=∞∑n=−∞anH(1)|n|(kpr)einθ,   r≥R, (4.1)
 ψ(r,θ)=∞∑n=−∞bnH(1)|n|(ksr)einθ,   r≥R. (4.2)

Since and satisfy boundary conditions and on , using the fact that and , we have

 ∞∑n=−∞kpanH(1)′|n|(kpR)einθ+iR∞∑n=−∞nbnH(1)|n|(ksR)einθ=−ν(θ)⋅uinc(R,θ),
 iR∞∑n=−∞nanH(1)|n|(kpR)einθ−∞∑n=−∞ksbnH(1)′|n|(ksR)einθ=−τ(θ)⋅uinc(R,θ).

Multiplying these two equations by , and integrating with respect to , we obtain a linear system for the coefficients and ,

 kpH(1)′|n|(kpR)an+inRH(1)|n|(ksR)bn=12π∫2π0−(ν(θ)⋅uinc(R,θ))e−inθdθ,
 inRH(1)|n|(kpR)an−ksH(1)′|n|(ksR)bn=12π∫2π0−(τ(θ)⋅uinc(R,θ))e−inθdθ.

Solve and to obtain

 an = ksH(1)′|n|(ksR)2πC∫2π0(ν(θ)⋅uinc(R,θ))e−inθdθ (4.3) +inH(1)|n|(ksR)2πRC∫2π0(τ(θ)⋅uinc(R,θ))e−inθdθ,
 bn = inH(1)|n|(kpR)2πRC∫2π0(ν(θ)⋅uinc(R,θ))e−inθdθ (4.4) −kpH(1)′|n|(kpR)2πC∫2π0(τ(θ)⋅uinc(R,θ))e−inθdθ,

where

 C:=n2R2H(1)|n|(kpR)H(1)|n|(ksR)−kpksH(1)′|n|(kpR)H(1)′|n|(ksR).

From Theorem 3.1, we know that and . Based on the asymptotic behavior of the Hankel functions, (4.1) and (4.2), the compressional part and shear part of the far field pattern are

 u∞,0p(θ)=ikpϕ∞=ikpe−iπ4√2πkp∞∑n=−∞ani−|n|einθ, (4.5)
 u∞,0s(θ)=iksψ∞=ikpe−iπ4√2πks∞∑n=−∞bni−|n|einθ, (4.6)

where and are given by (4.3) and (4.4), respectively.

The plane incident wave can be written in the form

 uinc(x;d,ap,as)=apdeikpx⋅d+asd⊥eiksx⋅d,   ap,as∈C. (4.7)

Let be a rigid disk centered at with radius . Then is a translation of with respect to . Denote the corresponding far field patterns for and by

 uB∞(^x;d,ap,as):=(u∞,0p(^x);u∞,0s(^x))  and  uBz∞(^x;d,ap,as):=(u∞,zp(^x);u∞,zs(^x)),

respectively. Then the far field pattern for can be computed using the following translation relations (see (2.13)-(2.16) in [16]):

 uBz∞(^x;d,1,0)=T1 uB∞(^x;d,1,0), (4.8)

with , and

 uBz∞(^x;d,0,1)=T2 uB∞(^x;d,0,1), (4.9)

with .

### 4.2 Extended sampling method

Let . For any , we define an inner product for the Hilbert space (see [3]):

 ⟨g,h⟩:=ωkp∫Sgp(^x)¯¯¯¯¯¯¯¯¯¯¯¯¯hp(^x)ds(^x)+ωks∫Sgs(^x)¯¯¯¯¯¯¯¯¯¯¯¯¯hs(^x)ds(^x),   g,h∈L2. (4.10)

Given , the elastic Herglotz wave function with density is defined as

 vg(x):=∫S{√kpwdeikpx⋅dg