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

TL;DR

This paper introduces a simple, effective sampling method for inverse elastic scattering problems that uses data from a single incident wave to determine obstacle location and size, even with limited data.

Contribution

The paper presents a novel sampling approach that reconstructs obstacle features from limited far field data using linear integral equations and indicators.

Findings

Successfully determines obstacle location and size with limited data

Effective with limited aperture data and minimal prior information

Numerical results confirm accuracy and robustness

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.

Equations

$${\Delta }^{\ast }\mathbf{u}+\rho {\omega }^2\mathbf{u}=0,$$

$$\mathbf{u}={\mathbf{u}}_p+{\mathbf{u}}_s,$$

$${\mathbf{u}}_p:=-\frac{1}{k_p^2}\text{ grad div }\mathbf{u},{\mathbf{u}}_s:=-\frac{1}{k_s^2}{\text{ grad}}^{\perp }{\text{div}}^{\perp }\mathbf{u}$$

$$\left\{\begin{array}[]{lll}&\Delta^{*}{\bm{u}}+\omega^{2}{\bm{u}}={\bm{0}},\ \ \ \textrm{in}\ \mathbb{R}^{2}\setminus\overline{D},\\ &\lim\limits_{{r}\rightarrow\infty}\sqrt{r}(\partial{\bm{u}}_{p}/\partial r-{i}k_{p}{\bm{u}}_{p})={\bm{0}},\ \ \ r=|{\bm{x}}|,\\ &\lim\limits_{{r}\rightarrow\infty}\sqrt{r}(\partial{\bm{u}}_{s}/\partial r-{i}k_{s}{\bm{u}}_{s})={\bm{0}},\ \ \ r=|{\bm{x}}|.\end{array}\right.$$

$${\bm{u}}({\bm{x}})=\frac{e^{ik_{p}|{\bm{x}}|}}{\sqrt{|{\bm{x}}|}}{u}_{p}^{\infty}({\hat{\bm{x}}}){\hat{\bm{x}}}+\frac{e^{ik_{s}|{\bm{x}}|}}{\sqrt{|{\bm{x}}|}}{u}_{s}^{\infty}({\hat{\bm{x}}}){\hat{\bm{x}}}^{\perp}+O\big{(}|{\bm{x}}|^{-3/2}\big{)},\ \ \ |{\bm{x}}|\rightarrow\infty,$$

$$\mathbb{S}:=\{\hat{\mathbf{x}}|\hat{\mathbf{x}}\in {\mathbb{R}}^2,|\hat{\mathbf{x}}|=1\}.$$

$$\mathbf{u}=-{\mathbf{u}}^{\text{inc}}\text{on}\partial D,$$

$$T_{\mathbf{\nu }}\mathbf{u}=-T_{\mathbf{\nu }}{\mathbf{u}}^{\text{inc}}\text{on}\partial D,$$

$$T_{\mathbf{\nu }}\mathbf{u}+i\sigma \mathbf{u}=-T_{\mathbf{\nu }}{\mathbf{u}}^{\text{inc}}-i\sigma {\mathbf{u}}^{\text{inc}}\text{on}\partial D,$$

$$\left\{\begin{array}[]{lll}&\Delta{\phi}+k^{2}_{p}{\phi}=0,\ \ \ \textrm{in}\ \mathbb{R}^{2}\setminus\overline{D},\\ &\Delta{\psi}+k^{2}_{s}{\psi}=0,\ \ \ \textrm{in}\ \mathbb{R}^{2}\setminus\overline{D},\\ &\lim\limits_{{r}\rightarrow\infty}\sqrt{r}(\partial\phi/\partial r-{i}k_{p}\phi)=0,\ \ \ \ r=|{\bm{x}}|,\\ &\lim\limits_{{r}\rightarrow\infty}\sqrt{r}(\partial\psi/\partial r-{i}k_{s}\psi)=0,\ \ \ r=|{\bm{x}}|.\end{array}\right.$$

$${\varphi }_{\infty }(\hat{\mathbf{x}})=\frac{1}{i{k}_p}{u}_p^{\infty }(\hat{\mathbf{x}}),{\psi }_{\infty }(\hat{\mathbf{x}})=\frac{1}{i{k}_s}{u}_s^{\infty }(\hat{\mathbf{x}}).$$

$$\mathbf{u}(r,\theta )$$

$$\frac{\partial }{\partial r}{u}_1(r,\theta )=i{k}_p\frac{e^{i{k}_{p}r}}{\sqrt{r}}{u}_p^{\infty }(\theta )\cos \theta -i{k}_s\frac{e^{i{k}_{s}r}}{\sqrt{r}}{u}_s^{\infty }(\theta )\sin \theta +O(r^{-3/2}),$$

$$\frac{\partial }{\partial r}{u}_2(r,\theta )=i{k}_p\frac{e^{i{k}_{p}r}}{\sqrt{r}}{u}_p^{\infty }(\theta )\sin \theta +i{k}_s\frac{e^{i{k}_{s}r}}{\sqrt{r}}{u}_s^{\infty }(\theta )\cos \theta +O(r^{-3/2}).$$

$$\text{div}\mathbf{u}$$

$${\text{div}}^{\perp }\mathbf{u}$$

$$\varphi (\mathbf{x})=-\frac{1}{k_p^2}\text{div}\mathbf{u}=\frac{1}{i{k}_p}\frac{e^{i{k}_p|\mathbf{x}|}}{\sqrt{|\mathbf{x}|}}{u}_p^{\infty }(\hat{\mathbf{x}})+O(|\mathbf{x}{|}^{-3/2}),$$

$$\psi (\mathbf{x})=-\frac{1}{k_s^2}{\text{div}}^{\perp }\mathbf{u}=\frac{1}{i{k}_s}\frac{e^{i{k}_s|\mathbf{x}|}}{\sqrt{|x|}}{u}_s^{\infty }(\hat{\mathbf{x}})+O(|\mathbf{x}{|}^{-3/2}),$$

$$\phi({\bm{x}})=\frac{e^{ik_{p}|{\bm{x}}|}}{\sqrt{|{\bm{x}}|}}\bigg{\{}\phi_{\infty}(\hat{\bm{x}})+O\bigg{(}\frac{1}{|{\bm{x}}|}\bigg{)}\bigg{\}},\ \ \ |{\bm{x}}|\rightarrow\infty,$$

$$\psi({\bm{x}})=\frac{e^{ik_{s}|{\bm{x}}|}}{\sqrt{|{\bm{x}}|}}\bigg{\{}\psi_{\infty}(\hat{\bm{x}})+O\bigg{(}\frac{1}{|{\bm{x}}|}\bigg{)}\bigg{\}},\ \ \ |{\bm{x}}|\rightarrow\infty.$$

$$\mathbf{u}=\mathrm{grad}\varphi +{\mathrm{grad}}^{\perp }\psi .$$

$$\frac{\partial \varphi }{\partial \mathbf{\nu }}+\frac{\partial \psi }{\partial \mathbf{\tau }}=-\mathbf{\nu }\cdot {\mathbf{u}}^{\mathrm{inc}},\frac{\partial \varphi }{\partial \mathbf{\tau }}-\frac{\partial \psi }{\partial \mathbf{\nu }}=-\mathbf{\tau }\cdot {\mathbf{u}}^{\mathrm{inc}},\text{on}\partial D.$$

$$\left\{\begin{array}[]{lll}&\Delta{\phi}+k^{2}_{p}{\phi}=0,\ \ \ \textrm{in}\ \mathbb{R}^{2}\setminus\overline{D},\\ &\Delta{\psi}+k^{2}_{s}{\psi}=0,\ \ \ \textrm{in}\ \mathbb{R}^{2}\setminus\overline{D},\\ &\frac{\partial{\bm{\phi}}}{\partial{\bm{\nu}}}+\frac{\partial\psi}{\partial\tau}=-{\bm{\nu}}\cdot{\bm{u}}^{\textrm{inc}},\ \ \ \textrm{on}\ \partial D,\\ &\frac{\partial{\bm{\phi}}}{\partial{\bm{\tau}}}-\frac{\partial\psi}{\partial{\bm{\nu}}}=-{\bm{\tau}}\cdot{\bm{u}}^{\textrm{inc}},\ \ \ \textrm{on}\ \partial D,\\ &\lim\limits_{{r}\rightarrow\infty}\sqrt{r}(\partial\phi/\partial r-{i}k_{p}\phi)=0,\ \ \ \ r=|{\bm{x}}|,\\ &\lim\limits_{{r}\rightarrow\infty}\sqrt{r}(\partial\psi/\partial r-{i}k_{s}\psi)=0,\ \ \ r=|{\bm{x}}|,\end{array}\right.$$

$$\left\{\begin{array}[]{lll}&\Delta u+k_{t}^{2}u=0,\ \ \ \textrm{in}\ \mathbb{R}^{2}\setminus\overline{B_{\bm{z}}},\\ &u=-{e}^{{i}k_{t}{\bm{x}}\cdot{\bm{d}}},\ \ \ \textrm{on}\ \partial B_{\bm{z}}\\ &\lim\limits_{{r}\rightarrow\infty}\sqrt{r}(\partial u/\partial r-{i}k_{t}u)=0,\ \ \ r=|{\bm{x}}|,\end{array}\right.$$

$${\int }_{\mathbb{S}}{u}_{\infty }^{B_{\mathbf{z}}}(\hat{\mathbf{x}};k_t,\mathbf{d})g_{\mathbf{z}}(\mathbf{d})ds(\mathbf{d})=\frac{1}{i{k}_t}{U}_t^{\infty }(\hat{\mathbf{x}}),\hat{\mathbf{x}}\in \mathbb{S}.$$

$$\bigg{\|}\int_{\mathbb{S}}u^{B_{\bm{z}}}_{\infty}(\hat{\bm{x}};k_{t},{\bm{d}})g_{\bm{z}}^{\varepsilon}({\bm{d}})ds({\bm{d}})-\frac{1}{ik_{t}}U_{t}^{\infty}(\hat{\bm{x}})\bigg{\|}_{L^{2}(\mathbb{S})}<\varepsilon$$

$$v_{g_{\mathbf{z}}^{\epsilon }}(\mathbf{x}):={\int }_{\partial B_{\mathbf{z}}}{e}^{i{k}_t\mathbf{x}\cdot \mathbf{d}}{g}_{\mathbf{z}}^{\epsilon }(\mathbf{d})ds(\mathbf{d}),\mathbf{x}\in B_{\mathbf{z}}$$

$$\underset{\epsilon \to 0}{\lim }\Vert v_{g_{\mathbf{z}}^{\epsilon }}{\Vert }_{H^1(B_{\mathbf{z}})}=\infty .$$

$$\phi_{t}({\bm{x}}):=\left\{\begin{array}[]{ll}&-\frac{1}{k_{p}^{2}}{\textrm{div}}\ {\bm{U}},\ \ \ t=p,\\ &-\frac{1}{k_{s}^{2}}{\textrm{div}}^{\perp}{\bm{U}},\ \ \ t=s,\end{array}\right.$$

$$\left\{\begin{array}[]{lll}&\Delta{\phi}+k^{2}_{p}{\phi}=0,\ \ \ \textrm{in}\ \mathbb{R}^{2}\setminus\overline{B},\\ &\Delta{\psi}+k^{2}_{s}{\psi}=0,\ \ \ \textrm{in}\ \mathbb{R}^{2}\setminus\overline{B},\\ &\frac{\partial\phi}{\partial{\bm{\nu}}}+\frac{\partial\psi}{\partial{\bm{\tau}}}=-{\bm{\nu}}\cdot{\bm{u}}^{\textrm{inc}},\ \ \ \textrm{on}\ \partial B,\\ &\frac{\partial\phi}{\partial{\bm{\tau}}}-\frac{\partial\psi}{\partial{\bm{\nu}}}=-{\bm{\tau}}\cdot{\bm{u}}^{\textrm{inc}},\ \ \ \textrm{on}\ \partial{B},\\ &\lim\limits_{{r}\rightarrow\infty}\sqrt{r}(\partial\phi/\partial r-{i}k_{p}\phi)=0,\ \ \ \ r=|{\bm{x}}|,\\ &\lim\limits_{{r}\rightarrow\infty}\sqrt{r}(\partial\psi/\partial r-{i}k_{s}\psi)=0,\ \ \ r=|{\bm{x}}|.\end{array}\right.$$

Full text

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 ${\bm{x}}=(x_{1};x_{2}):=(x_{1},x_{2})^{T}\in\mathbb{R}^{2}$, let $\hat{{\bm{x}}}:={\bm{x}}/|{\bm{x}}|$ and $\hat{\bm{x}}^{\perp}$ be obtained by rotating ${\hat{\bm{x}}}$ $\pi/2$ anticlockwise. Denote $\partial/\partial x_{i},i=1,2$ by $\partial_{i}$ for simplicity. In addition to the usual differential operators $\textrm{grad}{u}:=({\partial_{1}u};{\partial_{2}u})$ and $\textrm{div}{\bm{u}}:={\partial_{1}u_{1}}+{\partial_{2}u_{2}}$, we will make use of $\textrm{grad}^{\perp}{u}:=(-{\partial_{2}u};{\partial_{1}u})$ and $\textrm{div}^{\perp}{\bm{u}}:={\partial_{1}u_{2}}-{\partial_{2}u_{1}}$.

The propagation of time-harmonic waves in an isotropic homogeneous medium with Lame constants $\lambda$, $\mu$ ($\mu>0$, $2\mu+\lambda>0$) and density $\rho$ is governed by the Navier equation

[TABLE]

where ${\bm{u}}$ denotes the displacement field and $\omega$ denotes the circular frequency. The differential operator $\Delta^{*}{\bm{u}}:=\mu\Delta{\bm{u}}+(\lambda+\mu)\textrm{ grad div}\,{\bm{u}}$. In this paper, we assume $\rho\equiv 1$ for simplicity.

The solution ${\bm{u}}$ of (2.1) can be decomposed as

[TABLE]

where

[TABLE]

are known as the compressional part of ${\bm{u}}$ associated with the wave number $k_{p}:={\omega}/{\sqrt{2\mu+\lambda}}$, and the shear part of ${\bm{u}}$ associated with the wave number $k_{s}:={\omega}/{\sqrt{\mu}}$.

The direct elastic scattering problem for an obstacle is as follows. Given a bounded domain $D$ of class $C^{2}$ and an incident field ${\bm{u}}^{\textrm{inc}}$ such that ${\bm{u}}^{\textrm{inc}}|_{\partial D}\in[C(\partial D)]^{2}$ and ${\bm{u}}^{\textrm{inc}}$ is a solution of (2.1) in a neighborhood of $\partial D$, find the scattered field ${\bm{u}}\in[C^{2}(\mathbb{R}^{2}\setminus\overline{D})\bigcap C^{1}(\mathbb{R}^{2}\setminus D)]^{2}$ such that

[TABLE]

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

[TABLE]

uniformly in all directions ${\hat{\bm{x}}}$, where ${u}_{p}^{\infty}$ and ${u}_{s}^{\infty}$ are analytic functions on the unit circle

[TABLE]

Throughout the paper, ${\bm{u}}_{\infty}(\hat{\bm{x}}):=({u}_{p}^{\infty}(\hat{\bm{x}});{u}_{s}^{\infty}(\hat{\bm{x}}))$ is defined as the far field pattern of ${\bm{u}}$, and ${u}_{p}^{\infty}$ and ${u}_{s}^{\infty}$ 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 $\partial D$, which depend on the physical properties of the scatterer. The scattered field ${\bm{u}}$ satisfies

the Dirichlet boundary condition

[TABLE]

for a rigid body;

• 2)

the Neumann boundary condition

[TABLE]

for a cavity, where ${T_{\bm{\nu}}}:=2\mu\frac{\partial}{\partial{\bm{\nu}}}+\lambda{\bm{\nu}}\textrm{ div }-\mu{\bm{\nu}}^{\perp}\textrm{div}^{\perp}$ denotes the surface traction operator and ${\bm{\nu}}$ is the unit outward normal to $\partial D$;

• 3)

the impedance boundary condition

[TABLE]

with some real-valued parameter $\sigma\geq 0$.

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:

• IP-P:

Determine the location and size of the scatterer $D$ from the knowledge of compressional part ${u}_{p}^{\infty}(\hat{\bm{x}}),\hat{\bm{x}}\in\mathbb{S}$ or shear part ${u}_{s}^{\infty}(\hat{\bm{x}}),\hat{\bm{x}}\in\mathbb{S}$ of the far field pattern due to one incident wave.

• IP-F:

Determine the location and size of the scatterer $D$ from the knowledge of the far field pattern ${\bm{u}}_{\infty}(\hat{\bm{x}})=({u}_{p}^{\infty}(\hat{\bm{x}});{u}_{s}^{\infty}(\hat{\bm{x}})),\hat{\bm{x}}\in\mathbb{S}$ 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 ${\bm{u}}$ is the scattered field of (2.3) and ${\bm{u}}_{\infty}(\hat{\bm{x}}):=({u}_{p}^{\infty}(\hat{\bm{x}});u_{s}^{\infty}(\hat{\bm{x}}))$ is the corresponding far field pattern. The following theorem shows that ${\bm{u}}$ and ${\bm{u}}_{\infty}$ are related with the radiating solutions and their far field patterns of some Helmholtz equations, respectively.

Theorem 3.1**.**

Let ${\bm{u}}$ be a solution of (2.3). Then $\phi:=-\frac{1}{k_{p}^{2}}{\rm div}\,{\bm{u}}$ and $\psi:=-\frac{1}{k_{s}^{2}}{\rm div}^{\perp}\,{\bm{u}}$ are the radiating solutions of

[TABLE]

Furthermore, for the far field patterns $\phi_{\infty}$ and $\psi_{\infty}$ of $\phi$ and $\psi$, respectively,

[TABLE]

Proof.

It can be easily verified that $\phi=-\frac{1}{k_{p}^{2}}{\textrm{div}}\ {\bm{u}}$ and $\psi=-\frac{1}{k_{s}^{2}}{\textrm{div}}^{\perp}{\bm{u}}$ 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

[TABLE]

Consequently, we have that

[TABLE]

[TABLE]

Using these two equations, we obtain

[TABLE]

and

[TABLE]

Thus

[TABLE]

[TABLE]

which imply that $\phi$ and $\psi$ satisfy the Sommerfeld radiation conditions in (3.1).

Since $\phi=-\frac{1}{k_{p}^{2}}{\textrm{div}}\ {\bm{u}}$ and $\psi=-\frac{1}{k_{s}^{2}}{\textrm{div}}^{\perp}{\bm{u}}$ are radiating solutions of the Helmholtz equations, the far field patterns $\phi_{\infty}(\hat{\bm{x}})$ and $\psi_{\infty}(\hat{\bm{x}})$ have the following asymptotic expansions (see [10])

[TABLE]

[TABLE]

Then from (3.5) and (3.6), (3.2) is proved. ∎

Remark 3.2**.**

If the obstacle $D$ is a rigid body, the solution of the elastic scattering problem (2.3) satisfies the Dirichlet boundary condition ${\bm{u}}=-{\bm{u}}^{\rm{inc}}$. From (2.2), $\phi=-\frac{1}{k_{p}^{2}}{\rm div}\ {\bm{u}}$ and $\psi=-\frac{1}{k_{s}^{2}}{\rm div}^{\perp}{\bm{u}}$ satisfy

[TABLE]

Denote by ${\bm{\tau}}=(\tau_{1};\tau_{2})$ the unit tangent vector and by ${\bm{\nu}}=(\nu_{1};\nu_{2})=(\tau_{2};-\tau_{1})$ the unit outward normal vector on $\partial D$. By straightforward calculation, $\phi$ and $\psi$ satisfy the coupled boundary conditions

[TABLE]

From Theorem 3.1, $(\phi,\psi)$ satisfies the following Helmholtz equations with coupled boundary conditions

[TABLE]

which has a unique solution (see [22]). On the other hand, from (2.2), if $(\phi,\psi)$ solves (3.8), direct computation shows that ${\bm{u}}={\rm grad}\ \phi+{\rm grad}^{\perp}\psi$ satisfies the elastic scattering problem (2.3) with Dirichlet boundary condition.

3.2 Extended sampling method

Denote by ${\bm{U}}({\bm{x}}),{\bm{x}}\in\mathbb{R}^{2}\setminus\overline{D}$ and ${\bm{U}}_{\infty}(\hat{\bm{x}}):=(U_{p}^{\infty};U_{s}^{\infty}),\hat{\bm{x}}\in\mathbb{S}$ the scattered field and far field pattern of the unknown scatterer $D$ 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 $U_{p}^{\infty}(\hat{\bm{x}})$ or the shear part $U_{s}^{\infty}(\hat{\bm{x}})$. In this subsection, for simplicity we use $t$ to represent either $p$ or $s$.

Let $B_{\bm{z}}$ be a sound soft disk centered at ${\bm{z}}$ with radius $R$ and let $u^{B_{\bm{z}}}({\bm{x}};k_{t},{\bm{d}})$, ${\bm{x}}\in\mathbb{R}^{2}\setminus\overline{B_{\bm{z}}}$ solve

[TABLE]

where ${\bm{d}}\in\mathbb{S}$ is the incident direction. Denote $u_{\infty}^{B_{\bm{z}}}(\hat{\bm{x}};k_{t},{\bm{d}})$ the far field pattern of $u^{B_{\bm{z}}}({\bm{x}};k_{t},{\bm{d}})$ (for its series expansion, see e.g., [10, 19]). Then for IP-P, we introduce the far field equation

[TABLE]

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

Theorem 3.3**.**

Let $B_{\bm{z}}$ be a sound soft disk centered at ${\bm{z}}$ with radius $R$. Let $D$ be an obstacle. Assume that $k_{t}R$ does not coincide with any zero of the Bessel functions $J_{n},n=0,1,2,\cdots$. Then the following results hold for the far field equation (3.9):

If $D\subset B_{\bm{z}}$, for a given $\varepsilon>0$, there exists a function $g_{\bm{z}}^{\varepsilon}\in L^{2}(\mathbb{S})$ such that

[TABLE]

and the Herglotz wave function

[TABLE]

converges to the solution $w\in H^{1}(B_{\bm{z}})$ of the Helmholtz equation with $w=-\phi_{t}$ on $\partial B_{\bm{z}}$ as $\varepsilon\rightarrow 0$.

• 2.

If $D\cap B_{\bm{z}}=\emptyset$, every $g_{\bm{z}}^{\varepsilon}\in L^{2}(\mathbb{S})$ that satisfies (3.10) for a given $\varepsilon>0$ is such that

[TABLE]

Proof.

Denote the right hand side of the far field equation (3.9) by $\phi_{t}^{\infty}(\hat{\bm{x}}):=\frac{1}{ik_{t}}U_{t}^{\infty}(\hat{\bm{x}})$. From Theorem 3.1, we know that $\phi_{t}^{\infty}(\hat{\bm{x}})$ is the far field pattern of

[TABLE]

which is a radiating solution of the Helmholtz equation outside $D$. 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 $\Omega$ be a domain containing $D$. For a sampling point ${\bm{z}}\in\Omega$, we consider the linear ill-posed integral equation (3.9). By Theorem 3.3, one expects that $\|g^{\alpha}_{\bm{z}}\|_{L^{2}(\mathbb{S})}$ is relatively large when $D$ is outside $B_{\bm{z}}$ and relatively small when $D$ is inside $B_{\bm{z}}$. Consequently, a reconstruction of the location and support of $D$ can be obtained by plotting $\|g^{\alpha}_{\bm{z}}\|_{L^{2}(\mathbb{S})}$ for all sampling points ${\bm{z}}\in\Omega$.

• The Extended Sampling Method for IP-P

Generate a set $T$ of sampling points for $\Omega$ which contains $D$.

• 2.

For each sampling point ${\bm{z}}\in T$,

• a.

compute $u^{B_{\bm{z}}}_{\infty}(\hat{\bm{x}};k_{t},{\bm{d}})$ and set up a discrete version of (3.9);

• b.

use the Tikhonov regularization to compute an approximate solution $g^{\alpha}_{\bm{z}}$ of (3.9).

• 3.

Find the global minimum point ${\bm{z}}^{*}\in T$ for $\|g^{\alpha}_{\bm{z}}\|_{L^{2}(\mathbb{S})}$.

• 4.

Choose $B_{{\bm{z}}^{*}}$ to be the reconstruction for $D$.

As in [19], one can use a multilevel technique to find a suitable radius $R$ of the sampling disks $B_{\bm{z}},{\bm{z}}\in T$.

• The Multilevel ESM for IP-P

Choose the sampling disks with a large radius $R_{0}$. Generate a proper set $T_{0}$ of sampling points. Using ESM, determine the global minimum point ${\bm{z}}_{0}\in T_{0}$ for $\|g^{\alpha}_{\bm{z}}\|_{L^{2}(\mathbb{S})}$ and an approximation $D_{0}$ for $D$.

• 2.

For $j=1,2,\cdots$

• –

Let $R_{j}=R_{j-1}/2$ and generate a proper set $T_{j}$ of sampling points.

• –

Find the minimum point ${\bm{z}}_{j}\in T_{j}$ for $\|g^{\alpha}_{\bm{z}}\|_{L^{2}(\mathbb{S})}$ and an approximation $D_{j}$ for $D$. If ${\bm{z}}_{j}\not\in D_{j-1}$, go to Step 3.

• 3.

Choose ${\bm{z}}_{j-1}$ and $D_{j-1}$ to be the location and approximate support of $D$, 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 $B\subset\mathbb{R}^{2}$ be a rigid disk centered at the origin with radius $R$. Let ${\bm{u}}^{B}({\bm{x}})$ and ${\bm{u}}^{B}_{\infty}=(u_{p}^{\infty,{\bm{0}}};u_{s}^{\infty,{\bm{0}}})$ be the radiating solution and far field pattern of the elastic scattering problem (2.3) for $B$ due to an incident wave ${\bm{u}}^{\textrm{inc}}({\bm{x}})$, respectively. From Remark 3.2, we know that ${\bm{u}}^{B}({\bm{x}})=\textrm{grad}\ \phi+\textrm{grad}^{\perp}\psi$, where $(\phi,\psi)$ is the unique solution of

[TABLE]

Denote by $H_{n}^{(1)}(\cdot)$ the Hankel function of the first kind of order $n$. In polar coordinates ${\bm{x}}=(r\cos\theta,r\sin\theta)$, $\phi$ and $\psi$ can be written as (see [10])

[TABLE]

[TABLE]

Since $\phi$ and $\psi$ satisfy boundary conditions $\frac{\partial{\bm{\phi}}}{\partial{\bm{\nu}}}+\frac{\partial\psi}{\partial{\bm{\tau}}}=-{\bm{\nu}}\cdot{\bm{u}}^{\textrm{inc}}$ and $\frac{\partial{\bm{\phi}}}{\partial{\bm{\tau}}}-\frac{\partial\psi}{\partial{\bm{\nu}}}=-{\bm{\tau}}\cdot{\bm{u}}^{\textrm{inc}}$ on $\partial B$, using the fact that $\frac{\partial}{\partial{\bm{\nu}}}=\frac{\partial}{\partial r}$ and $\frac{\partial}{\partial{\bm{\tau}}}=\frac{1}{R}\frac{\partial}{\partial\theta}$, we have

[TABLE]

[TABLE]

Multiplying these two equations by $e^{-in\theta},n=-\infty,\cdots,\infty$, and integrating with respect to $\theta$, we obtain a linear system for the coefficients $a_{n}$ and $b_{n}$,

[TABLE]

[TABLE]

Solve $a_{n}$ and $b_{n}$ to obtain

[TABLE]

[TABLE]

where

[TABLE]

From Theorem 3.1, we know that $u_{p}^{\infty,{\bm{0}}}(\theta)=ik_{p}\phi_{\infty}$ and $u_{s}^{\infty,{\bm{0}}}(\theta)=ik_{s}\psi_{\infty}$. 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 ${\bm{u}}^{B}_{\infty}=(u_{p}^{\infty,{\bm{0}}};u_{s}^{\infty,{\bm{0}}})$ are

[TABLE]

[TABLE]

where $a_{n}$ and $b_{n}$ are given by (4.3) and (4.4), respectively.

The plane incident wave ${\bm{u}}^{\textrm{inc}}$ can be written in the form

[TABLE]

Let $B_{\bm{z}}:=\{{{\bm{x}}+{\bm{z}}};{\bm{x}}\in B\}$ be a rigid disk centered at ${\bm{z}}\in\,\mathbb{R}^{2}$ with radius $R$. Then $B_{\bm{z}}$ is a translation of $B$ with respect to ${\bm{z}}$. Denote the corresponding far field patterns for $B$ and $B_{\bm{z}}$ by

[TABLE]

respectively. Then the far field pattern for $B_{\bm{z}}$ can be computed using the following translation relations (see (2.13)-(2.16) in [16]):

[TABLE]

with $T_{1}:={diag}(e^{i(k_{p}{\bm{d}}-k_{p}\hat{\bm{x}})\cdot{\bm{z}}},e^{i(k_{p}{\bm{d}}-k_{s}\hat{\bm{x}})\cdot{\bm{z}}})$, and

[TABLE]

with $T_{2}:={diag}(e^{i(k_{s}{\bm{d}}-k_{p}\hat{\bm{x}})\cdot{\bm{z}}},e^{i(k_{s}{\bm{d}}-k_{s}\hat{\bm{x}})\cdot{\bm{z}}})$.

4.2 Extended sampling method

Let $\mathbb{L}^{2}:=[L^{2}(\mathbb{S})]^{2}$. For any ${\bm{g}}:=(g_{p};g_{s}),{\bm{h}}:=(h_{p};h_{s})\in\mathbb{L}^{2}$, we define an inner product for the Hilbert space $\mathbb{L}^{2}$ (see [3]):

[TABLE]

Given ${\bm{g}}\in\mathbb{L}^{2}$, the elastic Herglotz wave function with density ${\bm{g}}$ is defined as

[TABLE]

Denote by ${\bm{U}}({\bm{x}}),{\bm{x}}\in\mathbb{R}^{2}\setminus\overline{D}$ and ${\bm{U}}_{\infty}(\hat{\bm{x}}):=(U_{p}^{\infty};U_{s}^{\infty}),\hat{\bm{x}}\in\mathbb{S}$ the scattered wave and far field pattern of the unknown scatterer $D$ due to an incident wave, respectively. The inverse problem IP-F is to determine the location and size of the scatterer $D$ from the far field pattern ${\bm{U}}_{\infty}(\hat{\bm{x}})$.

Consider the far field equation

[TABLE]

where ${\bm{g}}\in\mathbb{L}^{2}$. The far field operator $F_{\bm{z}}:\mathbb{L}^{2}\rightarrow\mathbb{L}^{2}$ is defined as

[TABLE]

where ${\bm{u}}^{B_{\bm{z}}}_{\infty}(\hat{\bm{x}};{\bm{d}},a_{p},a_{s}):=(u_{p}^{\infty,{\bm{z}}};u_{s}^{\infty,{\bm{z}}})$ is the far field pattern of the rigid disk $B_{\bm{z}}$ due to the plane incident wave (4.7). As introduced in subsection 4.1, ${\bm{u}}_{\infty}^{B_{\bm{z}}}(\hat{\bm{x}};{\bm{d}},1,0)$ and ${\bm{u}}_{\infty}^{B_{\bm{z}}}(\hat{\bm{x}};{\bm{d}},0,1)$ can be easily computed by using the series expansion and the translation properties (4.8)-(4.9). Define the operator $G:[H^{1/2}(\partial B_{z})]^{2}\rightarrow\mathbb{L}^{2}$ by

[TABLE]

where ${\bm{u}}_{\infty}$ is the far-field pattern of the solution to the elastic Dirichlet boundary value problem with boundary value ${\bm{h}}$. The uniqueness for the boundary value problem implies that $G$ is injective. The far field operator $F_{\bm{z}}$ has the following properties (see [3]).

Lemma 4.1**.**

Assume that $\omega^{2}$ is not a Dirichlet eigenvalue of $-\Delta^{*}$ in $B_{\bm{z}}$. The ranges of $G$ and $(F_{\bm{z}}^{*}F_{\bm{z}})^{1/4}$ coincide.

The following theorem is the main result for the far field equation (4.12).

Theorem 4.1**.**

Assume that $\omega^{2}$ is not a Dirichlet eigenvalue of $-\Delta^{*}$ in $B_{\bm{z}}$.

If $D\subset B_{\bm{z}}$, the far field equation (4.12) is solvable.

• 2.

If $D\cap B_{\bm{z}}=\emptyset$, the far field equation (4.12) has no solution.

Proof.

1. Assume that $D\subset B_{\bm{z}}$. Define ${\bm{h}}:={\bm{U}}|_{\partial B_{\bm{z}}}$. According to the definition of $G$ in (4.18), ${\bm{h}}$ is a solution of $G{\bm{h}}={\bm{U}}_{\infty}$. Hence ${\bm{U}}_{\infty}$ is in the range of $G$. From Lemma 4.1, ${\bm{U}}_{\infty}$ is also in the range of $(F_{\bm{z}}^{*}F_{\bm{z}})^{1/4}$.

2. Assume that $D\cap B_{\bm{z}}=\emptyset$ and the far field equation (4.12) has a solution. From Lemma 4.1, ${\bm{U}}_{\infty}$ is in the range of $G$. Assume that $G{\bm{h}}={\bm{U}}_{\infty}$ for some ${\bm{h}}\in[H^{1/2}(\partial B_{\bm{z}})]^{2}$. Let ${\bm{u}}$ denote the radiating solution of the exterior Dirichlet problem for $B_{\bm{z}}$ with boundary value ${\bm{h}}$. Then ${\bm{u}}_{\infty}=G{\bm{h}}={\bm{U}}_{\infty}$. From Rellich’s lemma, we can identify the radiating solution ${\bm{v}}:={\bm{u}}={\bm{U}}$ in $\mathbb{R}^{2}\setminus(\overline{B_{\bm{z}}}\cup\overline{D})$. Since ${\bm{u}}$ is defined in $\mathbb{R}^{2}\setminus\overline{B_{\bm{z}}}$, ${\bm{U}}$ is defined in $\mathbb{R}^{2}\setminus\overline{D}$ and ${\bm{v}}$ can be extended from $\mathbb{R}^{2}\setminus(\overline{B_{\bm{z}}}\cup\overline{D})$ into $\mathbb{R}^{2}$. That is, ${\bm{v}}$ is an entire solution to the Navier equation. Since ${\bm{v}}$ also satisfies the radiation condition, it must vanish identically in $\mathbb{R}^{2}$. Then the total field coincides with the incident field and this leads to a contradiction since the incident field ${\bm{u}}^{\textrm{inc}}$ cannot satisfy the boundary condition. ∎

In consistency with the linear sampling method and the use of partial far field data, instead of solving (4.12), we consider the equation

[TABLE]

where $\alpha>0$ is the regularization parameter.

Theorem 4.2**.**

Let $F_{\bm{z}}$ be the far field operator (4.13) and assume that $\omega^{2}$ is not a Dirichlet eigenvalue of $-\Delta^{*}$ in $B_{\bm{z}}$. Denote by ${\bm{g}}_{\bm{z}}$ the solution of $(F_{\bm{z}}^{*}F_{\bm{z}})^{1/4}{\bm{g}}={\bm{U}}_{\infty}$ if it exists. For $\alpha>0$, let ${\bm{g}}_{\bm{z}}^{\alpha}$ denote the solution of (4.19) and ${\bm{v}}_{{\bm{g}}_{\bm{z}}^{\alpha}}$ denote the Herglotz wave function (4.11) with kernel ${\bm{g}}_{\bm{z}}^{\alpha}$.

If $D\subset B_{\bm{z}}$, then $\lim_{\alpha\rightarrow 0}{\bm{v}}_{{\bm{g}}_{\bm{z}}^{\alpha}}({\bm{x}})$ exists and

[TABLE]

for some positive $c$ depending only on $B_{\bm{z}}$.

• 2.

If $D\cap B_{\bm{z}}=\emptyset$, then $\lim_{\alpha\rightarrow 0}{\bm{v}}_{{\bm{g}}_{\bm{z}}^{\alpha}}({\bm{x}})\rightarrow\infty$.

Proof.

With the property of $F_{\bm{z}}$ in Lemma 4.1, and also Theorem 4.1, the proof is analogous to the proof of Theorem 5.41 and Corollary 5.42 of [10]. ∎

Now we are ready to present the extended sampling method (ESM) for IP-F. Let $\Omega$ be a domain containing $D$. For a sampling point ${\bm{z}}\in\Omega$, by Theorem 4.2, one expects that the solution $\|{\bm{g}}^{\alpha}_{\bm{z}}\|_{L^{2}(\mathbb{S})}$ of (4.19) is relatively large when $D$ is outside $B_{\bm{z}}$ and relatively small when $D$ is inside $B_{\bm{z}}$. Consequently, an approximation of the location and size of the scatterer $D$ can be obtained by plotting $\|{\bm{g}}^{\alpha}_{\bm{z}}\|_{L^{2}(\mathbb{S})}$ for all sampling points ${\bm{z}}\in\Omega$.

• The Extended Sampling Method for IP-F

Generate a set $T$ of sampling points for $\Omega$ which contains $D$.

• 3.

For each sampling point ${\bm{z}}\in T$,

• a.

compute ${\bm{u}}^{B_{\bm{z}}}_{\infty}(\hat{\bm{x}};{\bm{d}},1,0)$ and ${\bm{u}}^{B_{\bm{z}}}_{\infty}(\hat{\bm{x}};{\bm{d}},0,1)$, $\hat{\bm{x}},{\bm{d}}\in\mathbb{S}$;

• b.

choose a Tikhonov regularization parameter $\alpha$ and compute an approximate solution ${\bm{g}}^{\alpha}_{\bm{z}}$ to (4.19).

• 4.

Find the global minimum point ${\bm{z}}^{*}\in T$ for $\|{\bm{g}}^{\alpha}_{\bm{z}}\|_{L^{2}(\mathbb{S})}$.

• 5.

Choose $B_{{\bm{z}}^{*}}$ to be the reconstruction for $D$.

Remark 4.3**.**

As the Multilevel ESM for IP-P, one can use a similar multilevel technique to choose a suitable radius $R$ of the sampling disks.

5 Numerical examples

We present some numerical examples to show the performance of the proposed method. The synthetic far field data is generated using the boundary integral equation method in [22].

We consider three obstacles: a rigid pear given by

[TABLE]

a cavity peanut given by

[TABLE]

and a kite with impedance boundary condition ($\sigma=2$) given by

[TABLE]

For all numerical examples, $\omega=\pi,\mu=1,\lambda=2$. The incident plane wave is ${\bm{u}}^{inc}({\bm{x}})={\bm{d}}e^{ik_{p}{\bm{x}}\cdot{\bm{d}}}+{\bm{d}}^{\perp}e^{ik_{s}{\bm{x}}\cdot{\bm{d}}}$ with ${\bm{d}}=(1/2;\sqrt{3}/2)$.

For the inverse problem IP-P, which uses only the compressional part or shear part of the far field pattern, the synthetic data is a $52\times 1$ vector ${\bm{f}}_{P}=U_{p}^{\infty}(\hat{\bm{x}}_{j},{\bm{d}})$ or ${\bm{f}}_{P}=U_{s}^{\infty}(\hat{\bm{x}}_{j},{\bm{d}})$ with $52$ observation directions $\hat{\bm{x}}_{j},j=1,2,\cdots,52$, uniformly distributed on the unit circle.

For the inverse problem IP-F, {\bm{f}}_{F}=\big{(}U_{p}^{\infty}(\hat{\bm{x}}_{j},{\bm{d}});U_{s}^{\infty}(\hat{\bm{x}}_{j},{\bm{d}})\big{)},j=1,2,\cdots,52, is a $104\times 1$ vector.

5.1 Examples for IP-P

Let $\Omega=[-5,5]\times[-5,5]$ and choose the samplings points to be

[TABLE]

For each mesh point ${\bm{z}}\in T$, we use the Tikhonov regularization with a fixed parameter $\alpha=10^{-5}$. Equation (3.9) leads to a linear system $A^{\bm{z}}\vec{\bm{g}}_{\bm{z}}={\bm{f}}_{P}$, where $A^{\bm{z}}$ is the matrix given by

[TABLE]

The regularized solution is given by

[TABLE]

where $I$ is the identity matrix. We plot the contours for the indicator function

[TABLE]

for all the sampling points ${\bm{z}}\in T$.

Figure 1 shows the results for the rigid pear. The asterisks are the minimum locations of $I_{\bm{z}}$. The solid lines are the reconstructions and the red dashed lines are the exact boundaries. Similar results for the cavity peanut and the kite with impedance boundary condition are shown in Figure 2 and Figure 3, respectively.

Since the size of the scatterer is not known in advance, one can determine the radius $R$ of the sampling disks using the multilevel ESM. We start with a large sampling disk ($R=2.4$) and decrease the radius until a suitable $R$ is found.

For the rigid pear, $R$ is found to be $0.6$ using either the compressional part or shear part of the far field pattern. For the cavity peanut and impedance kite, the radius $R$ is $0.3$ using the compressional part of the far field pattern. The radius $R$ is $1.2$ using the shear part of the far field pattern.

Figure 4 shows the reconstructions of the multilevel ESM for the pear with Dirichlet boundary condition, the peanut with Neumann boundary condition, and the kite with the impedance boundary condition ($\sigma=2$).

5.2 Examples for IP-F

We use the same $T$ defined in (5.1) and the Tikhonov regularization with $\alpha=10^{-5}$. Equation (4.19) leads to the following linear system

[TABLE]

where

[TABLE]

According to the definition of the inner product (4.10), we have

[TABLE]

with the diagonal matrix $D_{ps}$ given by

[TABLE]

Then the regularized solution is given by

[TABLE]

We plot the contours for the indicator function

[TABLE]

for all the sampling points ${\bm{z}}\in T$.

Figure 5 shows the contour plots of $I_{\bm{z}}$ for the rigid pear and the reconstruction result. Figure 6 and Figure 7 show the reconstructions for the cavity peanut and impedance kite, respectively.

Again, we use the multilevel ESM starting with a large sampling disk ($R=2.4$). For the rigid pear, the cavity peanut and impedance kite ($\sigma=2$), the radius of the sampling disks are all found to be $R=0.6$. Figure 8 shows the reconstructions of the multilevel ESM for the pear, the peanut, and the kite.