Direct sampling methods for inverse elastic scattering problems

# Direct sampling methods for inverse elastic scattering problems

Xia Ji LSEC, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China (jixia@lsec.cc.ac.cn).    Xiaodong Liu Institute of Applied Mathematics, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, 100190 Beijing, China(xdliu@amt.ac.cn).    Yingxia Xi LSEC, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China (yxiaxi@lsec.cc.ac.cn).
###### Abstract

We consider the inverse elastic scattering of incident plane compressional and shear waves from the knowledge of the far field patterns. Specifically, three direct sampling methods for location and shape reconstruction are proposed using the different component of the far field patterns. Only inner products are involved in the computation, thus the novel sampling methods are very simple and fast to be implemented. With the help of the factorization of the far field operator, we give a lower bound of the proposed indicator functionals for sampling points inside the scatterers. While for the sampling points outside the scatterers, we show that the indicator functionals decay like the Bessel functions as the sampling point goes away from the boundary of the scatterers. We also show that the proposed indicator functionals continuously dependent on the far field patterns, which further implies that the novel sampling methods are extremely stable with respect to data error. For the case when the observation directions are restricted into the limited aperture, we firstly introduce some data retrieval techniques to obtain those data that can not be measured directly and then use the proposed direct sampling methods for location and shape reconstructions. Finally, some numerical simulations in two dimensions are conducted with noisy data, and the results further verify the effectiveness and robustness of the proposed sampling methods, even for multiple multiscale cases and limited-aperture problems.

Key words. Elastic scattering; sampling method; far field pattern; limited-aperture problem

AMS subject classifications. 35P25, 45Q05, 78A46, 74B05

## 1 Introduction

Scattering of elastic waves plays an important role in such different areas as seismic imaging, nondestructive testing, geophysical exploration, and medical diagnosis. In the last twenty years, non-iterative sampling methods for location and shape reconstruction in inverse elastic scattering problems have attracted a lot of interest. Typical examples include the Linear Sampling Method by Arens [2], the Factorization Method by Alves and Kress [1], Arens [2], Charalambopoulos & Kirsch et al. [6], and recently by Hu, Kirsch & Sini [13], and the method based on topological derivative by Guizina & Chikivhev [11]. The basic idea is to design an indicator which is big inside the underlying scatterer and relatively small outside. Recently, in [7], Chen & Huang introduced the Reverse Time Migration where the indicator functional is defined as the imaginary part of the cross-correlation of the weighted elastic Green function and the weighted back-propagated elastic wave field. Different to the previous sampling methods, the indicator functional decays as the sampling points go away from the boundary. We also refer to Li, Wang et al. [18] and Bao & Hu et al. [3] for iterative methods by using multiple frequency data. For the readers interested in a more comprehensive treatment of the inverse elastic scattering problems, we suggest consulting Bonnet and Constantinescu [5] on this subject.

The purpose of this paper is to introduce some novel indicator functionals for inverse elastic scattering problems with the far field patterns. These indicator functionals have a lower bound for sampling points inside the scatterers and decay like the Bessel functions as the sampling point goes away from the boundary of the scatterers. Besides, some data retrieval techniques will be introduced for limited aperture problems. We will focus on the case of two-dimensional scattering and present both the theoretical and numerical results for this case, illustrating the fast, effective and robust reconstructions, even for the limited aperture problems.

We begin with the notations used throughout this paper. All vectors will be denoted in bold script. For a vector , we introduce the two unit vectors and obtained by rotating anticlockwise by . For simplicity, we write for the usual partial derivative . Then, in addition to the usual differential operators and , we define two auxiliary differential operators and by and , respectively. It is easy to deduce the differential indentities .

Let be an open and bounded Lipschitz domain such that the exterior of is connected. Here and throughout the paper we denote the closure of the set . A confusion with the complex conjugate of is not expected. The propagation of time-harmonic elastic wave equation in an isotropic homogeneous media outside is governed by the reduced Navier equation

where denotes the total displacement field and is the circular frequency, and are Lam constants satisfying . It is well known that any solution of the Navier equation (LABEL:elastic) has a decomposition in the form

 u=up+us,

where

are known as the compressional (longitudinal) and shear (transversal) parts of , respectively. Here, and denote the compressional wave number and the shear wave number, respectively. It is clear that

 Δup+k2pup=0,div⊥up=0,Δus+k2sus=0,divus=0inR2∖¯¯¯¯Ω.

For a rigid body , the total displacement field satisfies the first (Dirichlet) boundary condition

 u=0on∂Ω. \hb@xt@.01(1.2)

For a cavity, we impose the second (Neumann) boundary condition

 Tνu=0on∂Ω, \hb@xt@.01(1.3)

where denotes the surface traction operator defined by

in terms of the exterior unit normal vector on .

Let denote the unit circle in . In elastic scattering, an important case is the scattering of a plane wave with incident direction which takes the form

 uin=apuinp+asuins,ap,as∈C,

where is a plane compressional wave and is a plane shear wave, respectively. The scatterer gives rise to a scattered field . The scattered field is a solution to (LABEL:elastic) and has a decomposition with and . The scattered field satisfies the Kupradze’s radiation conditions

 ∂uscp∂r−ikpuscp=o(r−12),∂uscs∂r−iksuscs=o(r−12), \hb@xt@.01(1.4)

uniformly in all directions as . For the unique solvability of the scattering problems (LABEL:elastic)-(LABEL:KupradzeRC) in the space we refer to Kupradze [17] and Li, Wang et al. [18].

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

 usc(x) = kpωeiπ/4√8πωeikp|x|√|x|u∞p(^x)^x \hb@xt@.01(1.6) +ksωeiπ/4√8πωeiks|x|√|x|u∞s(^x)x⊥+O(|x|−3/2),x→∞

uniformly in all direction . The functions and are known as the compressional and shear far field pattern of , and are analytic functions on . We want to remark here that, to simplify the subsequent representations, the coefficients in (LABEL:usasymptotic) are slightly different to those given in [2, 21]. We also have the asymptotic behaviour [2]

 T^xusc(x) = iωeiπ/4√8πωeikp|x|√|x|u∞p(^x)^x \hb@xt@.01(1.8) +iωeiπ/4√8πωeiks|x|√|x|u∞s(^x)x⊥+O(|x|−1),x→∞,

where denotes the surface traction operator on a circle centered on the origin of radius .

Throughout this paper, we will denote the pair of far field patterns

 [u∞p(^x,d,ap,as);u∞s(^x,d,ap,as)]

of the corresponding scattered field by , indicating the dependence on the observation direction , the incident direction and the pair of coefficients . We want to remark that the dependence on the coefficients is linear. For convenience, we use the following notations

 u∞pp(^x,d):=u∞p(^x,d,1,0),u∞sp(^x,d):=u∞p(^x,d,0,1), u∞ps(^x,d):=u∞s(^x,d,1,0),u∞ss(^x,d):=u∞s(^x,d,0,1).

Let be a subset of with nonempty interior. Then, we are interested in the following three inverse elastic scattering problems.

• IP-FF: Determine the location and shape of the scatterer from the knowledge of the full far field pattern for all observation directions , all incident directions and for and .

• IP-PP: Determine the location and shape of the scatterer from the knowledge of the compressional far field pattern for all observation directions , and all incident directions .

• IP-SS: Determine the location and shape of the scatterer from the knowledge of the shear far field pattern for all observation directions , and all incident directions .

Clearly, the measurement data in IP-PP and IP-SS is much less than those in IP-FF. For uniqueness of IP-FF, we refer to Hähner and Hsiao [12], while the corresponding result of IP-PP and IP-SS can be found in the recent works by Gintides & Sini [10] and Hu, Kirsch & Sini [13]. Due to analyticity, for uniqueness it suffices to know the far field pattern on the subset . However, as one would expect, the quality of the reconstructions decreases drastically for this so called limited-aperture problem. In particular, the traditional sampling type methods studied in Alves and Kress [1], Arens [2], Charalambopoulos & Kirsch et al. [6], Hu, Kirsch & Sini [13], and Sevroglou [21] fail to work. In this paper, we will introduce some data retrieval techniques to compute those data that can not be measured directly. Combining this and using the proposed sampling methods, the limited-aperture problems are desired to be solved partially.

This paper is organized as follows. In the next section, the theoretical analysis of the proposed reconstruction scheme will be established. A lower bound of the indicators for the sampling points inside the scatterers is obtained with the help of an inf-criterion characterization. The decay behavior of the indicators will then be studied for sampling points away from the boundary of the scatterers. A stability statement will also be established to reflect the important feature of the reconstruction scheme. For the limited-aperture problems, a data retrieval scheme will be introduced to numerically obtain the far field patterns that can not be measured directly. One then may combine the data retrieval technique and the previous sampling methods for inverse elastic scattering problems. Some numerical simulations in two dimensions will be presented in the last section to indicate the efficiency and robustness of the proposed methods.

## 2 Novel direct sampling methods and their mathematical basis

For any and any polarization define the two functions and by

 ϕpz(\bbtheta):=e−ikpz⋅\bbtheta(q⋅\bbtheta)andϕsz(\bbtheta):=e−iksz⋅\bbtheta(q⋅⊥),\bbtheta∈S, \hb@xt@.01(2.1)

respectively. Then for IP-FF, IP-PP and IP-SS we introduce the indicator functionals

 \hb@xt@.01(2.2) \hb@xt@.01(2.3) \hb@xt@.01(2.4)

respectively.

### 2.1 Lower bound estimate of ImissingFF(z), ImissingPP(z) and ImissingSS(z) for z∈Ω

We will use the short-hand . For any , we adopt the notation with

 gp(d):=g(d)⋅dandgs(d):=g(d)⋅d⊥,d∈S.

The Hilbert space will be equipped with the inner product

 (g,h)S:=∫Sgp(d)¯¯¯¯¯¯¯¯¯¯¯¯¯hp(d)ds(d)+∫Sgs(d)¯¯¯¯¯¯¯¯¯¯¯¯¯hs(d)ds(d),g,h∈L2. \hb@xt@.01(2.5)

It is clear that there holds the decomposition

 g(d)=gp(d)+gs(d),gp(d):=gp(d)d,gs(d):=gs(d)d⊥,

where and . For later use, we introduce the orthogonal projection operators and , i.e., for all , and . Their adjoint operators and are just the inclusions from and , respectively, to .

Consider the elastic Herglotz wave function

 vg(x):=∫S{eikpx⋅dgp(d)d+eiksx⋅dgs(d)d⊥}ds(d),x∈R2, \hb@xt@.01(2.6)

with a vector Herglotz kernel . We now introduce the far field operator defined by

 (Fg)(^x) := ∫Su∞(^x,d,gp(d),gs(d))ds(d) \hb@xt@.01(2.7) = ∫S(u∞pp(^x,d)u∞sp(^x,d)u∞ps(^x,d)u∞ss(^x,d))(gp(d)gs(d))ds(d), \hb@xt@.01(2.8)

i.e., is the far field pattern for the scattering of elastic Herglotz wave function with kernel . The far field operator plays an essential role in the investigations of the sampling type methods for inverse problems, we refer to [2] for a survey on the state of the art of its properties and applications.

For any , recall the two test functions and , we define by

 z(\bbtheta):=(ϕpz;ϕsz)=e−ikpz⋅\bbtheta(q⋅\bbtheta)\bbtheta+e−iksz⋅\bbtheta(q⋅⊥)⊥,\bbtheta∈S. \hb@xt@.01(2.9)

By interchanging orders of integration, we may rewrite our indicator given by (LABEL:Indicatorff) in a very simple form

 ImissingFF(z)=|(Fz,z)S|,z∈R2, \hb@xt@.01(2.10)

where is the inner product of the space given in (LABEL:innerproduct). Similarly, by defining and , we found that the indicators and given by (LABEL:Indicatorpp) and (LABEL:Indicatorss) can be written as a very simple form,

 ImissingPP(z):=∣∣(Fppz,pz)S∣∣andImissingSS(z):=∣∣(Fssz,sz)S∣∣,z∈R2, \hb@xt@.01(2.11)

with and , respectively.

Recall the Green’s tensor

of the Navier equation in in terms of the identity matrix and the Hankel function of the first kind of order zero. For any and any polarization , an elastic point source with source point and polarization is given by . From the asymptotics for the Hankel function it follows that the far field pattern of the point source is given by

 Φ∞p(^x,y;q) = e−ikp^x⋅y(q⋅^x)=ϕpz(^x),^x∈S, \hb@xt@.01(2.12) Φ∞s(^x,y;q) = e−iks^x⋅y(q⋅^x⊥)=ϕsz(^x),^x∈S. \hb@xt@.01(2.13)

Now we introduce the elastic single-layer operator , given by

 (S\bbphi)(x):=∫∂ΩΦ(x,y)\bbphi(y)ds(y),x∈∂Ω. \hb@xt@.01(2.14)

Then the following property of the single-layer operator is important for our subsequent analysis.

###### Lemma 2.1

Assume that is not a Dirichlet eigenvalue of in . Then there exist with

 |⟨\bbphi,S\bbphi⟩|≥c∥ϕ∥2[H−1/2(∂Ω)]2,\bbphi∈[H−1/2(∂Ω)]2,

where denotes the duality pairing in .

Proof. This property follows immediately from Lemma 1.17 of [15] in combination with Lemma 4.2 of [2].

Define the data-to-pattern operator by

 Gh:=u∞, \hb@xt@.01(2.15)

where is the far field pattern of the solution to the Dirichlet boundary value problem with boundary data . By a standard argument, we have that the data-to-pattern operator is compact, injective with dense range in . To characterize the scatterer by the corresponding date-to-pattern operator, we have the following lemma.

###### Lemma 2.2

For any and , define the function by (LABEL:phiz). Then the following holds.

• if and only if  .

• if and only if  .

• if and only if  .

Proof. The first result has been proved by Arens in Theorem 4.7 of [2]. The other two can be done analogously in combination with the definitions of the projection operators.

We introduce the Hergoltz wave operator by setting

 (Hg)(x):=∫S{eikpx⋅dgp(d)d+eiksx⋅dgs(d)d⊥}ds(d),x∈∂Ω.

is the trace on of the elastic Herglotz wave function (LABEL:eHerglotz) with vector Herglotz kernel . Its adjoint is given by

 (H∗\bbphi)(d):=(∫∂Ωe−ikpx⋅dd⋅\bbphi(x)ds(x);∫∂Ωe−iksx⋅dd⊥⋅\bbphi(x)ds(x)),d∈S.

We note that is exactly the far field pattern of the elastic single-layer potential , hence

 H∗=GSand thereforeH=S∗G∗,

where and denote the adjoints of and , respectively. On the other hand, we observe that is the far field pattern of the solution of the exterior Dirichlet problem with boundary data , which implies . Combining the previous operator equality, we deduce the factorization of the far field operator

 F=−GS∗G∗, \hb@xt@.01(2.16)

Consequently, we have the factorizations of the operators and ,

 Fp=−(PpG)S∗(PpG)∗,Fs=−(PsG)S∗(PsG)∗. \hb@xt@.01(2.17)

For all and polarization , define a subspace by

 Az:={\bbpsi∈L2:(\bbpsi,z)S=1},

where is the test function given in (LABEL:phiz) and is the inner product of defined by (LABEL:innerproduct). Now we are in a position to state the main result of this subsection.

###### Theorem 2.3

Consider the inverse elastic scattering by a rigid body . Assume that is not a Dirichlet eigenvalue of in . Let be the polarization. Then , if and only if,

 inf{|(F\bbpsi,\bbpsi)S|:\bbpsi∈Az}>0. \hb@xt@.01(2.18)

Furthermore, we have the lower bound estimate

 ImissingFF(z)≥c∥Φ(⋅,z)q∥2[H1/2(∂Ω)]2,z∈Ω, \hb@xt@.01(2.19)

for some constant which is independent of . Similar, we have

 ImissingPP(z)≥c∥Φ(⋅,z)q∥2[H1/2(∂Ω)]2andImissingSS(z)≥c∥Φ(⋅,z)q∥2[H1/2(∂Ω)]2,z∈Ω, \hb@xt@.01(2.20)

for some constant which is independent of .

Proof. First, (LABEL:infcreterior) follows directly by applying Theorem 1.16 in [15] to the factorization (LABEL:ffactorization), Lemma LABEL:Gproperty and noting the fact that the operator is coercive by Lemma LABEL:Scoercive. Furthermore, note that , using Theorem 1.16 of [15] again we deduce that

 inf{|(F\bbpsi,\bbpsi)S|:\bbpsi∈Az}≥c∥Φ(⋅,z)q∥2[H1/2(∂Ω)]2

for some constant which is independent of . Second, we observe that

 (z,z)S=2π, \hb@xt@.01(2.21)

which implies . From this we derive that

 ImissingFF(z) = |(Fz,z)S| = 4π2|(F0,0)S| ≥ 4π2inf{|(F\bbpsi,\bbpsi)S|:\bbpsi∈Az} ≥ c∥Φ(⋅,z)q∥2[H1/2(∂Ω)]2

for some constant which is independent of . Thus, the lower bound estimate (LABEL:estimateff) for the indicator in has been proved. The other two lower bound estimates in (LABEL:estimateppss) can be shown analogously using the factorizations (LABEL:ppssfactorization).

We conclude this subsection by a remark that the analogous result of Theorem LABEL:theorem1 holds for the Neumann boundary condition. Our analysis rely on the factorization of the far field operator. The data-to-pattern operator is now defined to map into the far field pattern of the exterior Neumann boundary value problem with boundary data . We introduce the hypersingular integral operator by

 (N\bbphi)(x):=Tν(x)∫∂Ω(Tν(y)Φ(x,y))T\bbphi(y)ds(y),x∈∂Ω,

for . Then one can derive the factorization of the far field operator for Neumann problem in the form

 F=−GN∗G∗.

Based on this, the analogous results of Lemmas LABEL:Scoercive, LABEL:Gproperty and Theorem LABEL:theorem1 can be derived.

Finally, we want to remark that the assumption on interior eigenvalues in Lemma LABEL:Scoercive and Theorem LABEL:theorem1 has only its theory interest. It is well-known that the classical sampling type methods for solving inverse scattering problems fail if the wave number is an eigenvalue of a corresponding interior eigenvalue problem. We refer to a modification proposed by Kirsch & Liu[16] to avoid the eigenvalues. However, our sampling methods are independent of the interior eigenvalues from the numerical point of view.

### 2.2 Indicator behavior for the sampling points away from the boundary ∂Ω

We have known from the previous subsection that the values of the indicator functionals can not be small arbitrarily for sampling points inside. In this subsection, we study the behavior of the indicator functionals for sampling points away from the boundary.

To simplify the subsequent representations, we introduce

 J(w):=wwT|w|2,w∈R2.

Straightforward calculations show that . For the scattering problem (LABEL:elastic)-(LABEL:KupradzeRC), the far-field pattern has the following representation [2]

 u∞p(^x,d)^x = \hb@xt@.01(2.23) −J(^x)e−ikp^x⋅yTν(y)usc(y,d)}ds(y), u∞s(^x,d)^x⊥ = \hb@xt@.01(2.25) −J(^x⊥)e−iks^x⋅yTν(y)usc(y,d)}ds(y).

By straightforward calculations, it can be seen that, for all ,

 [Tν(y)(J(^x)e−ikp^x⋅y)]^x=Tν(y)(^xe−ikp^x⋅y)

and

 [Tν(y)(J(^x⊥)e−ikp^x⋅y)]^x⊥=Tν(y)(^x⊥e−ikp^x⋅y).

Inserting this into (LABEL:upinfty)-(LABEL:usinfty), we find that

 u∞p(^x,d) = ∫∂Ω{uinp(y,−^x)⋅Tν(y)usc(y,d) \hb@xt@.01(2.27) −[Tν(y)uinp(y,−^x)]⋅usc(y,d)}ds(y), u∞s(^x,d) = ∫∂Ω{uins(y,−^x)⋅Tν(y)usc(y,d) \hb@xt@.01(2.29) −[Tν(y)uins(y,−^x)]⋅usc(y,d)}ds(y).

Now we introduce the following auxiliary functions

 Gpp(z,d):=∫Su∞pp(^x,d)¯¯¯¯¯¯¯¯¯¯¯¯¯ϕpz(^x)ds(^x),Gsp(z,d):=∫Su∞sp(^x,d)¯¯¯¯¯¯¯¯¯¯¯¯¯ϕpz(^x)ds(^x), Gps(z,d):=∫Su∞ps(^x,d)¯¯¯¯¯¯¯¯¯¯¯¯¯ϕsz(^x)ds(^x),Gss(z,d):=∫Su∞ss(^x,d)¯¯¯¯¯¯¯¯¯¯¯¯¯ϕsz(^x)ds(^x).

Then the indicators , and involve the terms

 ∫SGpp(z,d)ϕpz(d)ds(d),∫SGps(z,d)ϕpz(d)ds(d), \hb@xt@.01(2.30) ∫SGsp(z,d)ϕsz(d)ds(d),∫SGss(z,d)ϕsz(d)ds(d). \hb@xt@.01(2.31)

Noting that and , by the well known Riemann-Lebesgue Lemma, we obtain that all the expressions in (LABEL:G1)-(LABEL:G2) go to as , so we have

 ImissingFF(z)→0,ImissingPP(z)→0,andImissingSS(z)