An inverse acoustic-elastic interaction problem with phased or phaseless far-field data

TL;DR

This paper develops a mathematical and numerical framework for solving inverse acoustic-elastic scattering problems using phased or phaseless far-field data, including techniques to address translation invariance in phaseless cases.

Contribution

It introduces a coupled boundary value problem approach, boundary integral equations, and a reference ball technique for the first time to solve inverse elastic obstacle problems with phaseless data.

Findings

The method achieves high accuracy in locating and shaping elastic obstacles.

The reference ball technique effectively breaks translation invariance in phaseless data.

Numerical experiments confirm the robustness and effectiveness of the proposed algorithms.

Abstract

Consider the scattering of a time-harmonic acoustic plane wave by a bounded elastic obstacle which is immersed in a homogeneous acoustic medium. This paper concerns an inverse acoustic-elastic interaction problem, which is to determine the location and shape of the elastic obstacle by using either the phased or phaseless far-field data. By introducing the Helmholtz decomposition, the model problem is reduced to a coupled boundary value problem of the Helmholtz equations. The jump relations are studied for the second derivatives of the single-layer potential in order to establish the corresponding boundary integral equations. The well-posedness is discussed for the solution of the coupled boundary integral equations. An efficient and high order Nystr\"{o}m-type discretization method is proposed for the integral system. A numerical method of nonlinear integral equations is developed for the inverse problem. For the case of phaseless data, we show that the modulus of the far-field pattern is invariant under a translation of the obstacle. To break the translation invariance, an elastic reference ball technique is introduced. We prove that the inverse problem with phaseless far-field pattern has a unique solution under certain conditions. In addition, a numerical method of the reference ball technique based nonlinear integral equations is also proposed for the phaseless inverse problem. Numerical experiments are provided to demonstrate the effectiveness and robustness of the proposed methods.

Tables

Algorithm I: Iterative algorithm for the phased IAEIP
Step 1	Send an incident plane wave $u^{\mathrm{inc}}$ with a fixed wave number ${\kappa }_a$ and a fixed incident direction $d\in \Omega$, and then collect the corresponding far-field data $u_{\infty }$ for the scatterer $D$;
Step 2	Select an initial star-like curve ${\Gamma }^{\left(0\right)}$ for the boundary ${\Gamma }_D$ and the error tolerance $ϵ$. Set $k=0$;
Step 3	For the curve ${\Gamma }^{\left(k\right)}$, compute the densities ${\phi }_l$, $l=1,2,3$ from (5.1);
Step 4	Solve (6.3) to obtain the updated approximation ${\Gamma }^{\left(k+1\right)}:={\Gamma }^{\left(k\right)}+q$ and evaluate the error $E_{k+1}$ defined in (6.4);
Step 5	If $E_{k+1}\ge ϵ$, then set $k=k+1$ and go to Step 3. Otherwise, the current approximation ${\Gamma }^{\left(k+1\right)}$ is taken to be the final reconstruction of ${\Gamma }_D$.

Table 2. Table 1. Parametrization of the exact boundary curves.

Type	Parametrization
Apple-shaped	$p_D\left(t\right)=\frac{0.55\left(1+0.9\cos t+0.1\sin 2t\right)}{1+0.75\cos t}(\cos t,\sin t),t\in [0,2\pi ]$
Peanut-shaped	$p_D\left(t\right)=0.65\sqrt{0.25{\cos }^{2}t+{\sin }^{2}t}(\cos t,\sin t),t\in [0,2\pi ]$

Algorithm II: Iterative algorithm for the phaseless IAEIP
Step 1	Send an incident plane wave $u^{\mathrm{inc}}$ with a fixed wave number ${\kappa }_a$ and a fixed incident direction $d\in \Omega$, and then collect the corresponding far-field data $u_{\infty }$ for the scatterer $D\cup B$;
Step 2	Select an initial star-like curve ${\Gamma }^{\left(0\right)}$ for the boundary ${\Gamma }_D$ and the error tolerance $ϵ$. Set $k=0$;
Step 3	For the curve ${\Gamma }^{\left(k\right)}$, compute the densities ${\phi }_{1,\sigma }$, ${\phi }_{2,\sigma }$ and ${\phi }_{3,\sigma }$ from field equations;
Step 4	Solve (6.16) to obtain the updated approximation ${\Gamma }^{\left(k+1\right)}:={\Gamma }^{\left(k\right)}+q$ and evaluate the error $E_{k+1}$ defined in (6.17);
Step 5	If $E_{k+1}\ge ϵ$, then set $k=k+1$ and go to Step 3. Otherwise, the current approximation ${\Gamma }^{\left(k+1\right)}$ is served as the final reconstruction of ${\Gamma }_D$.

Equations

$$\mathrm{curl}\mathbf{U}={\partial }_{x_1}{U}_2-{\partial }_{x_2}{U}_1,\mathbf{curl}u=({\partial }_{x_2}u,-{\partial }_{x_1}u{)}^{\top }.$$

$$\mu \Delta \mathbf{U}+(\lambda +\mu )\nabla \nabla \cdot \mathbf{U}+{\omega }^2{\rho }_{\mathrm{e}}\mathbf{U}$$

$$\Delta u+{\kappa }_{\mathrm{a}}^{2}u$$

$$T(\mathbf{U})=-u\nu ,\mathbf{U}\cdot \nu =\frac{1}{{\omega }^2{\rho }_{\mathrm{a}}}{\partial }_{\nu }u\mathrm{on}{\Gamma }_D.$$

$$\underset{r\to \infty }{\lim }{r}^{\frac{1}{2}}({\partial }_r{u}^{\mathrm{s}}-\mathrm{i}{\kappa }_{\mathrm{a}}{u}^{\mathrm{s}})=0,r=|x|.$$

$$T(\mathbf{U}):=\mu {\partial }_{\nu }\mathbf{U}+(\lambda +\mu )(\nabla \cdot \mathbf{U})\nu .$$

$$\mathbf{U}=\nabla \varphi +\mathrm{curl}\psi ,$$

$$\nabla [(\lambda +2\mu )\Delta \varphi +{\omega }^2{\rho }_{\mathrm{e}}\varphi ]+\mathrm{curl}(\mu \Delta \psi +{\omega }^2{\rho }_{\mathrm{e}}\psi )=0,$$

$$\Delta \varphi +{\kappa }_{\mathrm{p}}^2\varphi =0,\Delta \psi +{\kappa }_{\mathrm{s}}^2\psi =0.$$

$${\kappa }_{\mathrm{p}}=\omega {\left(\frac{{\rho }_{\mathrm{e}}}{\lambda +2\mu }\right)}^{1/2},{\kappa }_{\mathrm{s}}=\omega {\left(\frac{{\rho }_{\mathrm{e}}}{\mu }\right)}^{1/2},$$

$$\mu \nu \cdot {\partial }_{\nu }\nabla \varphi +\mu \nu \cdot {\partial }_{\nu }\mathrm{curl}\psi -(\lambda +\mu ){\kappa }_{\mathrm{p}}^2\varphi +u^{\mathrm{s}}$$

$$\mu \tau \cdot {\partial }_{\nu }\nabla \varphi +\mu \tau \cdot {\partial }_{\nu }\mathrm{curl}\psi$$

$${\partial }_{\nu }\varphi +{\partial }_{\tau }\psi -{\partial }_{\nu }{u}^{\mathrm{s}}/({\omega }^2{\rho }_{\mathrm{a}})$$

$$f_1=-u^{\mathrm{inc}},f_2=0,f_3={\partial }_{\nu }{u}^{\mathrm{inc}}/({\omega }^2{\rho }_{\mathrm{a}}).$$

$$\{\begin{array}{ll}\Delta \varphi +{\kappa }_{\mathrm{p}}^2\varphi =0,\Delta \psi +{\kappa }_{\mathrm{s}}^2\psi =0,& \mathrm{in}D,\\ \Delta u^{\mathrm{s}}+{\kappa }_{\mathrm{a}}^2{u}^{\mathrm{s}}=0,& \mathrm{in}{\mathbb{R}}^2\setminus \overline{D},\\ \mu \nu \cdot {\partial }_{\nu }\nabla \varphi +\mu \nu \cdot {\partial }_{\nu }\mathrm{curl}\psi -(\lambda +\mu ){\kappa }_{\mathrm{p}}^2\varphi +u^{\mathrm{s}}=f_1,& \mathrm{on}{\Gamma }_D,\\ \tau \cdot {\partial }_{\nu }\nabla \varphi +\tau \cdot {\partial }_{\nu }\mathrm{curl}\psi =f_2,& \mathrm{on}{\Gamma }_D,\\ {\partial }_{\nu }\varphi +{\partial }_{\tau }\psi -{\partial }_{\nu }{u}^{\mathrm{s}}/({\omega }^2{\rho }_{\mathrm{a}})=f_3,& \mathrm{on}{\Gamma }_D,\\ \underset{r\to \infty }{\lim }{r}^{\frac{1}{2}}({\partial }_r{u}^{\mathrm{s}}-\mathrm{i}{\kappa }_{\mathrm{a}}{u}^{\mathrm{s}})=0,& r=|x|.\end{array}$$

$$-{\int }_{{\Gamma }_D}{u}^{\mathrm{s}}\overline{{\partial }_{\nu }{u}^{\mathrm{s}}}\mathrm{d}s=$$

$$=$$

$$+{\omega }^2{\rho }_{\mathrm{a}}{\int }_{{\Gamma }_D}\left(\mu {\partial }_{\nu }(\nabla \varphi +\mathrm{curl}\psi )\cdot \tau \right)\left((\nabla \overline{\varphi }+\mathrm{curl}\overline{\psi })\cdot \tau \right)\mathrm{d}s$$

$$=$$

$$\cdot \left((\nabla \overline{\varphi }+\mathrm{curl}\overline{\psi })\cdot \nu \nu +(\nabla \overline{\varphi }+\mathrm{curl}\overline{\psi })\cdot \tau \tau \right)\mathrm{d}s$$

$$=$$

$$\nabla \cdot (\nabla \overline{\varphi }+\mathrm{curl}\overline{\psi })-{\omega }^2{\rho }_{\mathrm{e}}(\nabla \varphi +\mathrm{curl}\psi )\cdot (\nabla \overline{\varphi }+\mathrm{curl}\overline{\psi }))\mathrm{d}x=0,$$

$$\Im {\int }_{{\Gamma }_D}{u}^{\mathrm{s}}\overline{{\partial }_{\nu }{u}^{\mathrm{s}}}\mathrm{d}s=0.$$

$$\nabla \varphi =-\mathbf{curl}\psi \text{mbox}inD,$$

$$u^{\rm s}(x)=\frac{\mathrm{e}^{\mathrm{i}\kappa_{\rm a}|x|}}{\sqrt{|x|}}\Big{\{}u_{\infty}(\hat{x})+\mathcal{O}\left(\frac{1}{|x|}\right)\Big{\}}\quad\text{as}~{}|x|\to\infty,$$

$$⟨a,b⟩=a{b}^{\top },\nabla a=(\nabla a_1,\nabla a_2{)}^{\top },\nabla a^{\top }=(\nabla a_1,\nabla a_2)=(\nabla a{)}^{\top }.$$

$$\nabla_{y}(\nabla_{x}f)=\langle\nabla_{x},\nabla_{y}\rangle f=\nabla_{x}\nabla_{y}^{\top}f=\left[\begin{array}[]{cc}\partial^{2}_{x_{1}y_{1}}f&\partial^{2}_{x_{1}y_{2}}f\\[2.0pt] \partial^{2}_{x_{2}y_{1}}f&\partial^{2}_{x_{2}y_{2}}f\\ \end{array}\right]$$

$$\displaystyle(\nabla_{x}\nabla_{y}^{\top}f,\nu)=\left[\begin{array}[]{cc}\partial^{2}_{x_{1}y_{1}}f&\partial^{2}_{x_{1}y_{2}}f\\[2.0pt] \partial^{2}_{x_{2}y_{1}}f&\partial^{2}_{x_{2}y_{2}}f\\ \end{array}\right]\left[\begin{array}[]{c}\nu_{1}\\ \nu_{2}\\ \end{array}\right],\quad(\boldsymbol{\rm curl}_{x}\nabla_{y}^{\top}f,\nu)=\left[\begin{array}[]{cc}\partial^{2}_{x_{2}y_{1}}f&\partial^{2}_{x_{2}y_{2}}f\\[2.0pt] -\partial^{2}_{x_{1}y_{1}}f&-\partial^{2}_{x_{1}y_{2}}f\\ \end{array}\right]\left[\begin{array}[]{c}\nu_{1}\\ \nu_{2}\\ \end{array}\right].$$

$$\Phi (x,y;\kappa )=\frac{\mathrm{i}}{4}{H}_0^{(1)}(\kappa |x-y|),x\ne y,$$

$$\omega (x)={\int }_{{\Gamma }_D}\Phi (x,y;\kappa )g(y)\mathrm{d}s(y),\chi (x)={\int }_{{\Gamma }_D}\frac{\partial \Phi (x,y;\kappa )}{\partial \nu (y)}g(y)\mathrm{d}s(y),x\in {\mathbb{R}}^2\setminus {\Gamma }_D.$$

$$\zeta (x)={\int }_{{\Gamma }_D}\frac{\partial \Phi (x,y;\kappa )}{\partial \tau (y)}g(y)\mathrm{d}s(y),x\in {\mathbb{R}}^2\setminus {\Gamma }_D.$$

$$(\nabla \omega {)}_{\pm }(x)={\int }_{{\Gamma }_D}{\nabla }_x\Phi (x,y;\kappa )g(y)\mathrm{d}s(y)\mp \frac{1}{2}\nu (x)g(x),x\in {\Gamma }_D,$$

$$(\nabla \omega {)}_{\pm }(x):=\underset{h\to 0^{+}}{\lim }(\nabla \omega )(x\pm h\nu (x)).$$

$$(\mathbf{curl}\omega {)}_{\pm }(x)={\int }_{{\Gamma }_D}{\mathbf{curl}}_x\Phi (x,y;\kappa )g(y)\mathrm{d}s(y)\pm \frac{1}{2}\tau (x)g(x),x\in {\Gamma }_D.$$

Full text

1. Introduction

Consider the scattering of a time-harmonic acoustic plane wave by a bounded penetrable obstacle, which is immersed in an open space occupied by a homogeneous acoustic medium such as some compressible inviscid air or fluid. The obstacle is assumed to be a homogeneous and isotropic elastic medium. When the incident wave impinges the obstacle, a scattered acoustic wave will be generated in the open space and an elastic wave is induced simultaneously inside the obstacle. This scattering phenomenon leads to an acoustic-elastic interaction problem (AEIP). Given the incident wave and the obstacle, the direct acoustic-elastic interaction problem (DAEIP) is to determine the pressure of the acoustic wave field and the displacement of the elastic wave field in the open space and in the obstacle, respectively; the inverse acoustic-elastic interaction problem (IAEIP) is to determine the elastic obstacle from the far-field pattern of the acoustic wave field. The AEIPs have received ever-increasing attention due to their significant applications in seismology and geophysics [29]. Despite many work done so far for both of the DAEIP and IAEIP, they still present many challenging mathematical and computational problems due to the complex of the model equations and the associated Green tensor, as well as the nonlinearity and ill-posedness.

The phased IAEIP referres to the IAEIP that determines the location and shape of the elastic obstacle from the phased far-field data, which contains both the phase and amplitude information. It has been extensively studied in the recent decades. In [10, 11], an optimization based variational method and a decomposition method were proposed to the IAEIP. The direct imaging methods, such as the linear sampling method [34, 35] and the factorization method [23, 40, 13], were also developed to the corresponding inverse problems with far-field and near-field data. For the theoretical analysis, the uniqueness results may be found in [34, 36] for the phased IAEIP.

The phaseless IAEIP is to determine the location and shape of the elastic obstacle from the modulus of the far-field acoustic scattering data, which contains only the amplitude information. Due to the translation invariance property of the phaseless far-field field, it is impossible to uniquely determine the location of the unknown object by plane incident wave, which makes the phaseless inverse problem much more challenging than the phased counterpart. Various numerical methods have been proposed to solve the phaseless inverse obstacle scattering problems, especially for the acoustic waves which are governed by the scalar Helmholtz equation. For the shape reconstruction with one incident plane wave, we refer to the Newton iterative method [27], the nonlinear integral equation method [14, 15], the fundamental solution method [22], and the hybrid method [30]. In particular, the nonlinear integral equation method, which was proposed by Johansson and Sleeman [20], was extended to reconstruct the shape of a sound-soft crack by using phaseless far-field data from a single incident plane wave [12]. To reconstruct the location and shape simultaneously, Zhang et al. [45, 46] proposed an iterative method by using the superposition of two plane waves with different incident directions to reconstruct the unknown object. In [17], a phase retrieval technique combined with the direct sampling method was proposed to reconstruct the location and shape of an obstacle from phaseless far-field data. The method was extended to the phaseless inverse elastic scattering problem and phaseless IAEIP [16]. We refer to [38, 41, 37, 44, 43] for the uniqueness results on the inverse scattering problems by using phaseless data. Related phaseless inverse scattering problems as well as numerical methods can be found in [1, 31, 18, 5, 3, 4, 42, 24]. Recently, a reference ball technique based nonlinear integral equations method was proposed in [9] to break the translation invariance from phaseless far-field data by one incident plane wave. In our recent work [8], we extended this method to the inverse elastic scattering problem with phaseless far-field data by using a single incident plane wave to recover both the location and shape of a rigid elastic obstacle.

In this paper, we consider both the DAEIP and IAEIP. In particular, we study the IAEIP of determining the location and shape of an elastic obstacle from the phased or phaseless far-field data with a single incident plane wave. The goal of this work is fivefold:

1. (1)

deduce the jump relations for the second derivatives of the single-layer potential and the coupled system of boundary integral equations; 2. (2)

prove the well-posedness of the solution for the coupled system and develop a Nyström-type discretization for the boundary integral equations; 3. (3)

show the translation invariance of the phaseless far-field pattern and present a uniqueness result for the phaseless IAEIP; 4. (4)

propose a numerical method of nonlinear integral equations to reconstruct the obstacle’s location and shape by using the phased far-field data from a single plane incident wave; 5. (5)

develop a reference ball based method to reconstruct both the obstacle’s location and shape by using phaseless far-field data from a single plane incident wave.

For the direct problem, instead of considering directly the coupled acoustic and elastic wave equations, we make use of the Helmholtz decomposition and reduce the model problem into a coupled boundary value problem of the Helmholtz equations. The method of boundary integral equations is adopted to solve the coupled Helmholtz system. However, the boundary conditions are more complicated, since the second derivatives of surface potentials are involved due to the traction operator. Therefore, we investigate carefully the jump relations for the second derivatives of the single-layer potential and establish coupled boundary integral equations. Moreover, we prove the existence and uniqueness for the solution of the coupled boundary integral equations, and develop a Nyström-type discretization to efficiently and accurately solve the direct acoustic-elastic interaction problem. The proposed method is extremely efficient for the direct scattering problem since we only need to solve the scalar Helmholtz equations instead of solving the vector Navier equations. Related work on the direct acoustic-elastic interaction problems and time-domain acoustic-elastic interaction problem can be found in [2, 19, 33, 39].

For the inverse problem, motivated by the reference ball technique [32, 41] and the recent work [8, 9], we give a uniqueness result for the phaseless IAEIP by introducing an elastic reference ball, and also propose a nonlinear integral equations based iterative numerical scheme to solve the phased and phaseless IAEIP. Since the location of reference ball is known, the method breaks the translation invariance and is able to recover the location information of the obstacle with negligible additional computational costs. Numerical results show that the method is effective and robust to reconstruct the obstacle with either the phased or phaseless far-field data.

The paper is organized as follows. In Section 2, we introduce the coupled acoustic-elastic interaction problem and show the uniqueness for the coupled boundary value problem by using the Helmholtz decomposition. In Section 3, we study the jump properties for the second derivatives of the single-layer potential and establish the coupled boundary integral equations. The existence and uniqueness of the solution for the coupled boundary integral equations are given. Section 4 is devoted to the translation invariance and the uniqueness for the phaseless IAEIP. Section 5 presents a high order Nyström-type discretization to solve the coupled boundary value problem. In Section 6, a method of nonlinear integral equations and a reference ball based method are developed to solve the phased and phaseless inverse problems, respectively. Numerical experiments are provided to demonstrate the effectiveness of the proposed methods in Section 7. The paper is concluded with some general remarks and directions for future work in Section 8.

2. Problem formulation

Consider the scattering problem of a time-harmonic acoustic plane wave by a two-dimensional elastic obstacle $D$ with $\mathcal{C}^{2}$ boundary $\Gamma_{D}$. The elastic obstacle $D$ is assumed to be homogeneous and isotropic with a mass density $\rho_{\rm e}$; the exterior domain $\mathbb{R}^{2}\setminus\overline{D}$ is assumed to be filled with a homogeneous and compressible inviscid air or fluid with a mass density $\rho_{\rm a}>0$. Denote by $\nu=(\nu_{1},\nu_{2})^{\top}$ and $\tau=(-\nu_{2},\nu_{1})^{\top}$ the unit normal vector and the tangential vector on $\Gamma_{D}$, respectively. Let $\nu_{\perp}=(\nu_{2},-\nu_{1})^{\top}=-\tau$. Given a vector function $\boldsymbol{U}=(U_{1},U_{2})^{\top}$ and a scalar function $u$, we introduce the scalar and vector curl operators

[TABLE]

Specifically, the time-harmonic acoustic plane wave is given by $u^{\rm inc}(x)=e^{{\rm i}\kappa_{\rm a}x\cdot d}$, where $d=(\cos\theta,\sin\theta)^{\top}$ is the propagation direction vector, $\theta\in[0,2\pi)$ is the incident angle. Given the incident field $u^{\rm inc}$, the direct problem is to find the elastic wave displacement $\boldsymbol{U}\in(C^{2}(D)\cap C(\overline{D}))^{2}$ and the acoustic wave pressure $u\in C^{2}(\mathbb{R}^{2}\setminus\overline{D})\cap C(\mathbb{R}^{2}\setminus D)$, which satisfy the Navier equation and the Helmholtz equation, respectively:

[TABLE]

Moreover, $\boldsymbol{U}$ and $u$ are required to satisfy the transmission conditions

[TABLE]

The scattered acoustic wave pressure $u^{\rm s}:=u-u^{\rm inc}$ is required to satisfy the Sommerfeld radiation condition

[TABLE]

Here $\omega>0$ is the angular frequency, $\kappa_{\rm a}=\omega/c$ is the wavenumber in the air/fluid with the sound speed $c$, and $\lambda,\mu$ are the Lamé parameters satisfying $\mu>0,\lambda+\mu>0$. The traction operator $T$ is defined by

[TABLE]

It can be shown (cf. [33, 39]) that the scattering problem (2.1)–(2.4) admits a unique solution $(\boldsymbol{U},u)$ for all but some particular frequencies $\omega$, which are called the Jones frequencies [21]. At the Jones frequency, the acoustic wave field $u$ is unique, but the elastic field $\boldsymbol{U}$ is not unique. Since the Jones frequency happens only for some special geometries [21], for simplicity, we assume that $D$ does not admit any Jones mode in this work.

For any solution $\boldsymbol{U}$ of the elastic wave equation (2.1), we introduce the Helmholtz decomposition

[TABLE]

where $\phi,\psi$ are two scalar potential functions. Substituting (2.5) into (2.1) yields

[TABLE]

which is fulfilled if $\phi$ and $\psi$ satisfy the Helmholtz equation with a different wavenumber, respectively:

[TABLE]

Here

[TABLE]

are the compressional wavenumber and the shear wavenumber, respectively.

Substituting the Helmholtz decomposition into (2.3) and taking the dot product with $\nu$ and $\tau$, respectively, we obtain

[TABLE]

where

[TABLE]

In summary, the scalar potential functions $\phi,\psi$ and the scattered acoustic wave $u^{\rm s}$ satisfy the following coupled boundary value problem

[TABLE]

The following result concerns the uniqueness of the boundary value problem (2.6).

Theorem 2.1**.**

The coupled boundary value problem (2.6) has at most one solution for $\kappa_{\rm p}>0,\kappa_{\rm s}>0,\kappa_{\rm a}>0$.

Proof.

It suffices to show that $\phi=\psi=u^{\rm s}=0$ when $f_{1}=f_{2}=f_{3}=0$. It follows from straightforward calculations that

[TABLE]

where $A:B={\rm tr}(AB^{\top})$ is the Frobenius inner product of square matrices $A$ and $B$. The last two identities follow from Green’s formula and the Navier equation (2.1). Taking the imaginary part of the above equation yields

[TABLE]

which gives that $u^{\rm s}=0$ in $\mathbb{R}^{2}\backslash\overline{D}$ by Rellich’s lemma. Using the continuity conditions (2.3), we conclude that $\boldsymbol{U}$ is identically zero in $D$ provided that there is no Jones mode in $D$. Hence,

[TABLE]

which implies $\Delta\phi=0$ and $\Delta\psi=0$. The proof is completed by noting that $\Delta\phi=-\kappa^{2}_{\rm p}\phi=0$ and $\Delta\psi=-\kappa^{2}_{\rm s}\psi=0$ in $D$. ∎

It is known that a radiating solution of the Helmholtz equation (2.2) has the asymptotic behaviour of the form

[TABLE]

uniformly in all directions $\hat{x}:=x/|x|$. The function $u_{\infty}$, defined on the unit circle $\Omega=\{x\in\mathbb{R}^{2}:|x|=1\}$, is known as the far-field pattern of $u^{\rm s}$. Let $B=\left\{x\in\mathbb{R}^{2}:|x-x_{0}|\leq R\right\}\subset\mathbb{R}^{2}$ be an artificially added elastic ball centered at $x_{0}$ such that $D\cap B=\emptyset$. The problem geometry is shown in Figure 1. For brevity, we denote the boundary of $D$ and $B$ by $\Gamma_{D}$ and $\Gamma_{B}$, respectively. The phased and phaseless IAEIP can be stated as follows:

Problem 1** (Phased IAEIP).**

Given an incident plane wave $u^{\rm inc}$ with a single incident direction $d$ and the corresponding far-field pattern $u_{\infty}(\hat{x}),~{}\forall\hat{x}\in\Omega$ due to the unknown obstacle $D$, the inverse problem is to determine the location and shape of the boundary $\Gamma_{D}$.

Problem 2** (Phaseless IAEIP).**

Given an incident plane wave ${u}^{\rm inc}$ with a single incident direction $d$ and the corresponding phaseless far-field pattern $|u_{\infty}(\hat{x})|,~{}\forall\hat{x}\in\Omega$ due to the scatterer $D\cup B$, the inverse problem is to determine the location and shape of the boundary $\Gamma_{D}$.

3. Boundary integral equations

In this section, we derive the boundary integral equations for the coupled boundary value problem (2.6) and discuss their well-posedness.

3.1. Jump relations

We begin with investigating the jump relations for the surface potentials at the boundary $\Gamma_{D}$.

For given vectors $a=(a_{1},a_{2})^{\top}\in\mathbb{R}^{2}$ and $b=(b_{1},b_{2})^{\top}\in\mathbb{R}^{2}$, denote

[TABLE]

For a given scalar function $f(x,y)$, define

[TABLE]

and

[TABLE]

Denote the fundamental solution of the two-dimensional Helmholtz equation by

[TABLE]

where $H_{0}^{(1)}$ is the Hankel function of the first kind of order zero. The single- and double-layer potentials with density $g$ are defined by

[TABLE]

In addition, we define the tangential-layer potential by

[TABLE]

The jump relations can be found in [6] for the single- and double-layer potentials as $x\to\Gamma_{D}$. It is necessary to study the jump properties for the derivatives of those layer potentials in order to derive the boundary integral equations for the coupled boundary value problem (2.6).

Lemma 3.1**.**

The first derivatives of the single-layer potential $\omega$ with density $g\in C^{0,\alpha}(\Gamma_{D})$, $0<\alpha<1$, can be uniformly extended in a Hölder continuous fashion from $\mathbb{R}^{2}\setminus\overline{D}$ into $\mathbb{R}^{2}\setminus{D}$ and from $D$ into $\overline{D}$ with the limiting values

[TABLE]

where

[TABLE]

Moreover, for the single-layer potential $\omega$ with density $g\in C^{0,\alpha}(\Gamma_{D})$, $0<\alpha<1$, we have

[TABLE]

Proof.

Noting

[TABLE]

and

[TABLE]

we may similarly show (3.1) by following the proof of Theorem 2.17 in [6]. It is clear to note (3.2) by combining the fact that ${\bf curl}~{}\omega=(\nabla\omega)_{\perp}$, $\nu_{\perp}=-\tau$ and the jump relation (3.1). ∎

Lemma 3.2**.**

The first derivatives of the double-layer potential $\chi$ with density $g\in C^{1,\alpha}(\Gamma_{D})$, $0<\alpha<1$, can be uniformly extended in a Hölder continuous fashion from $\mathbb{R}^{2}\setminus\overline{D}$ into $\mathbb{R}^{2}\setminus{D}$ and from $D$ into $\overline{D}$ with the limiting values

[TABLE]

Proof.

Using the jump relation (3.2) and the identities

[TABLE]

and

[TABLE]

we may easily show (3.2) by following the proof of Theorem 7.32 in [26] and Theorem 2.23 in [6]. ∎

Theorem 3.3**.**

For the tangential-layer potential $\zeta$ with density $g\in C^{1,\alpha}(\Gamma_{D})$, $0<\alpha<1$, we have

[TABLE]

Proof.

Using the integration by parts, we have

[TABLE]

which implies (3.6) by noting the jump relation (3.1). ∎

Theorem 3.4**.**

The second derivatives of the single-layer potential $\omega$ with density $g\in C^{1,\alpha}(\Gamma_{D})$, $0<\alpha<1$, can be uniformly extended in a Hölder continuous fashion from $\mathbb{R}^{2}\setminus\overline{D}$ into $\mathbb{R}^{2}\setminus{D}$ and from $D$ into $\overline{D}$ with the limiting values

[TABLE]

and

[TABLE]

Proof.

Using (3.3)–(3.4), we have from taking the second derivatives of the single-layer potential that

[TABLE]

Combining the above equation and the jump relations (3.2)–(3.6) gives (3.4).

Analogously, noting ${\bf curl}~{}\omega=(\nabla\omega)_{\perp}$, $\nu_{\perp}=-\tau$ and $\tau_{\perp}=\nu$, we have

[TABLE]

Applying the jump relation (3.2)–(3.6) again yields (3.4). ∎

In view of $\partial_{\nu}(\nabla\omega)=(\nabla\nabla^{\top}\omega,\nu)$, $\partial_{\nu}({\bf curl}~{}\omega)=({\bf curl}\nabla^{\top}\omega,\nu)$ and Theorem 3.4, we have the following result.

Corollary 3.5**.**

For the single-layer potential $\omega$ with density $g\in C^{1,\alpha}(\Gamma_{D})$, $0<\alpha<1$, we have on $\Gamma_{D}$ that

[TABLE]

and

[TABLE]

where

[TABLE]

3.2. Boundary integral equations

We introduce the single-layer integral operator and the corresponding far-field integral operator

[TABLE]

the normal derivative integral operator

[TABLE]

and the tangential derivative integral operator

[TABLE]

Let the solution of (2.6) be given in the form of single-layer potentials, i.e.,

[TABLE]

where the densities $g_{1}\in C^{1,\alpha}(\Gamma_{D})$, $g_{2}\in C^{1,\alpha}(\Gamma_{D})$, and $g_{3}\in C^{1,\alpha}(\Gamma_{D})$.

Letting $x\in{D}$ tend to boundary $\Gamma_{D}$ in (3.9)–(3.10) and $x\in\mathbb{R}^{2}\setminus\overline{D}$ tend to boundary $\Gamma_{D}$ in (3.11), using the jump relations of the single-layer potentials, Lemmas 3.1–3.2, Corollary 3.5, and the boundary conditions of (2.6), we obtain on $\Gamma_{D}$ that

[TABLE]

We point out that $\nu$ and $\tau$ inside of $[\cdot]$ or $\langle\cdot\rangle$ are with respect to the variable $y$; otherwise $\nu$ and $\tau$ are taken with respect to the variable $x$. For brevity, we shall adopt the same notations in the rest of the paper but they should be clear from the context. The far-field pattern is

[TABLE]

where $\gamma_{\rm a}=e^{\mathrm{i}\pi/4}/{\sqrt{8\kappa_{\rm a}\pi}}$.

Now we discuss the uniqueness and existence of the solution for the integral equations (3.12).

Theorem 3.6**.**

There exists at most one solution to the boundary integral equations (3.12) if $\kappa_{\rm a}$ is not the eigenvalue of the interior Dirichlet problem of the Helmholtz equation in $D$.

Proof.

It suffices to show that $g_{1}=g_{2}=g_{3}=0$ if $f_{1}=f_{2}=f_{3}=0$ for equations in (3.12). For $x\in\mathbb{R}^{2}\setminus\Gamma_{D}$, we define single-layer potentials

[TABLE]

Let

[TABLE]

and

[TABLE]

Since $\phi_{i}$, $\psi_{i}$ and $u^{\rm s}_{e}$ satisfy the boundary value problem (2.6), they are identically zero by Theorem 2.1. Using the jump condition of single layer potentials, we have on $\partial D$ that

[TABLE]

Combining (3.14) and the fact $\phi_{i}=\psi_{i}=0$ in $D$, we derive that $\phi_{e}$ and $\psi_{e}$ satisfy the zero boundary condition on $\partial D$. By the uniqueness of the exterior problem for the Helmholtz equation, it holds that $\phi_{e}=\psi_{e}=0$ in $\mathbb{R}^{2}\backslash\overline{D}$. We conclude that $g_{1}=g_{2}=0$ by (3.15). Similarly, we have on $\partial D$ that

[TABLE]

By (3.16), we see that $u^{\rm s}_{i}$ satisfies the zero Dirichlet boundary condition. Since $\kappa_{\rm a}$ is not the eigenvalue of the interior Dirichlet problem, we conclude that $u^{\rm s}_{i}$ is identically zero in $D$, which implies $g_{3}=0$ by (3.17). ∎

Theorem 3.7**.**

There exists a unique solution to the boundary integral equations (3.12) if none of $\kappa_{\rm p}$, $\kappa_{\rm s}$ and $\kappa_{\rm a}$ is the eigenvalue of the interior Dirichlet problem of the Helmholtz equation in $D$.

Proof.

Since the original coupled equations (2.1)–(2.4) admit a unique solution $(u,\boldsymbol{U})$, by the Helmholtz decomposition $\boldsymbol{U}=\nabla\phi+\mathbf{curl}~{}\psi$, we have on $\partial D$ that

[TABLE]

Plugging the single layer representations (3.9)–(3.10) for $\phi$ and $\psi$ into the above equations and using the jump property of boundary integral operators, we have

[TABLE]

where $I$ is the identity operator. Following the idea in [28], we can show that the boundary integral equation (3.18) admits a solution $(g_{1},g_{2})$ when neither of $\kappa_{\rm p}$ or $\kappa_{\rm s}$ is the eigenvalue of the interior Dirichlet problem of the Helmholtz equation in $D$. For $g_{3}$, since $\kappa_{\rm a}$ is not an interior Dirichlet eigenvalue either, the single layer operator $S_{\kappa_{\rm a}}$ is invertible. Therefore

[TABLE]

Based on the construction, one can easily see that $g_{1}$, $g_{2}$ and $g_{3}$ satisfy the boundary integral equations (3.12). ∎

4. translation invariance and a uniqueness result

In this section, we prove the translation invariance of the phaseless far-field pattern and present a uniqueness result for the phaseless IAEIP.

Theorem 4.1**.**

Under the assumption of Theorem 3.7, let $u_{\infty}$ be the far field pattern of the scattered waves $u^{\rm s}$ with the incident plane wave ${u}^{\rm inc}(x)=e^{{\rm i}\kappa_{\rm a}d\cdot x}$. For the shifted domain $D_{h}:=\{x+h:x\in D\}$ with a fixed vector $h\in\mathbb{R}^{2}$, the far-field pattern $u_{\infty}^{h}$ satisfies

[TABLE]

Proof.

We assume that the densities $g_{1}^{h},g_{2}^{h}$ and $g_{3}^{h}$ solve the boundary integral equations (3.12) with $\Gamma_{D}$ replaced by $\Gamma_{D_{h}}$. We claim that if $g_{1},g_{2}$ and $g_{3}$ solve the equations (3.12), then

[TABLE]

In fact, by substituting the above equations into the right hand side of (3.12) with $\Gamma_{D}$ replaced by $\Gamma_{D_{h}}$ and noting that

[TABLE]

where $\tilde{x}=x-h$, $\tilde{y}=y-h$, the first equation in (3.12) becomes

[TABLE]

where the boundary integral operators defined on $\Gamma_{D_{h}}$ are denoted by $S^{h},K^{h}$ and $H^{h}$. Similarly, the second and third equations of (3.12) can be verified in the same way. Thus, (4.1) follows from the fact that the system of boundary integral equations (3.12) for $D_{h}$ has a unique solution (cf. Theorems 3.6 and 3.7).

Combining (3.13) and (4.1), we obtain

[TABLE]

which completes the proof. ∎

Remark 4.2**.**

The assumption on the wavenumber $\kappa_{\rm p}$, $\kappa_{\rm s}$ and $\kappa_{\rm a}$ in Theorem 4.1 can be removed by using combined boundary layer potentials to represent the solution. We refer to [28] for a related discussion.

Theorem 4.1 implies that the location of the obstacle can not be uniquely recovered by the modules of far-field pattern when the plane wave is used as an incident field. To overcome this difficulty, motivated by [41], we may introduce an elastic reference ball $B=B(x_{0},R)=\{x\in\mathbb{R}^{2}:|x-x_{0}|<R\}$ to the scattering system in order to break the translation invariance.

Assume that $P$ is a disk (with positive radius) such that $P\subset\mathbb{R}^{2}\setminus(\overline{D}\cup\overline{B})$ and $\kappa_{\rm a}^{2}$ is not a Dirichlet eigenvalue of $-\Delta$ in $P$. Denote the boundary of $P$ by $\partial P$. Consider that the incident wave is given by a plane wave $u^{\rm inc}(x,d)$ and a point source $v^{\rm inc}(x,z)$, i.e., $u^{\rm inc}(x,d)=e^{{\rm i}\kappa_{\rm a}x\cdot d}$ and $v^{\rm inc}(x,z)=\Phi(x,z;\kappa_{\rm a})$, where $z\in\partial P$ is the source location. Assume further that $\{u_{D\cup B}^{\rm s}(x)(x,d),u_{D\cup B}^{\infty}(\hat{x},z)\}$ and $\{v_{D\cup B}^{\rm s}(x,z),v_{D\cup B}^{\infty}(\hat{x},z)\}$ are the scattered field and the far-field pattern generated by $D\cup B$ corresponding to the incident field $u^{\rm inc}(x,d)$ and $v^{\rm inc}(x,z)$, respectively.

Now we present a uniqueness result for the phaseless inverse scattering problem. A similar uniqueness result may be found in [41, Theorem 4.1] for the phaseless inverse medium scattering problem.

Theorem 4.3**.**

Let $D_{1}$ and $D_{2}$ be two elastic obstacles with $\mathcal{C}^{2}$ boundaries, and $\omega$ is not a Jones frequency either for $D_{1}$ or $D_{2}$. Suppose that the far-field patterns satisfy the following conditions:

[TABLE]

for a fixed $d_{0}\in\Omega$, then $D_{1}=D_{2}$.

Proof.

By (4.2)–(4.4), we have

[TABLE]

In view of (4.2) and (4.3), we assume that

[TABLE]

where $r(\hat{x},d_{0})=|u_{D_{j}\cup B}^{\infty}(\hat{x},d_{0})|$, $s(\hat{x},z)=|v_{D_{j}\cup B}^{\infty}(\hat{x},z)|$, $\alpha_{j}(\hat{x},d_{0})$ and $\beta_{j}(\hat{x},z)$ are real-valued functions, $j=1,2$. Following the proof of [7, Theorem 3.16], we can show the mixed reciprocity relation

[TABLE]

Using (4.5) and similar arguments in [41, Theorem 3.1], we obtain that

[TABLE]

where $\gamma(\hat{x}):=\alpha_{1}(\hat{x},d_{0})-\alpha_{2}(\hat{x},d_{0})-2m\pi$, $\hat{x}\in S$, $m\in\mathbb{Z}$, and $S\subset\Omega$ is an open arc. Furthermore, for $x\in\Gamma_{B}$, $-d\in S$, we get

[TABLE]

Noting that the total fields $(u_{D_{1}\cup B},\boldsymbol{U}_{D_{1}\cup B})$ and $(u_{D_{2}\cup B},\boldsymbol{U}_{D_{2}\cup B})$ are the solutions of (2.1)–(2.4) corresponding to the scatterers $D_{1}\cup B$ and $D_{2}\cup B$, respectively, we find that

[TABLE]

satisfy the Navier equation and the Helmholtz equation

[TABLE]

and the transmission conditions on $\Gamma_{B}$

[TABLE]

Suppose that $\left(w(x,d_{0}),\boldsymbol{W}(x,d_{0})\right)$ is the solution of (2.1)–(2.4) corresponding to the single reference ball $B$ with incident plane wave $u^{\rm inc}(x,d_{0})$, then the far-field $w^{\infty}(\hat{x},d_{0})\not\equiv 0$, $\hat{x}\in\Omega$. Using the Betti formula and the transmission condition (4.7), and noting the identity

[TABLE]

we have

[TABLE]

We claim $|w^{\infty}(-d,d_{0})|\not\equiv 0$, $\forall-d\in S$. Otherwise, we obtain by using the analytic continuation that $w^{\infty}(\hat{x},d_{0})=0$, $\forall\hat{x}\in\Omega$. This is a contradiction. By continuity, there exists an open curve $\tilde{S}\subset S$, such that $|w^{\infty}(-d,d_{0})|\not=0$, $\forall-d\in\tilde{S}$, which implies that $e^{{\rm i}\gamma(-d)}=1$ for $-d\in\tilde{S}$. From (4.6), we have

[TABLE]

Again, using the reciprocity relation and the analyticity of $u_{D_{j}\cup B}^{\infty}(\hat{x},d)$ for $j=1,2$, we obtain that the far-field patterns $u_{D_{1}\cup B}^{\infty}$ and $u_{D_{2}\cup B}^{\infty}$ coincide for all observation and incident directions $\hat{x},d\in\Omega$. We conclude from [34, Theorem 4.1] that $D_{1}=D_{2}$. ∎

Remark 4.4**.**

In view of the proof of Theorem 4.3, we can also assume that the scatterers $D_{1},D_{2}$ and the reference ball $B$ possess different mass densities and Lamé parameters.

Remark 4.5**.**

Theorem 4.3 only gives a sufficient condition to uniquely reconstruct the the unknown obstacle $D$ with phaseless data. We expect the uniqueness result also holds with much less data. In particular, for our numerical practice, we do not use the phaseless far-field data generated by the point source $v^{\rm inc}(x,z)$ with $z\in\partial P$, although in this case, the uniqueness result is still under investigation.

5. Nyström-type discretization for boundary integral equations

In this section, we introduce a Nyström-type discretization for the boundary integral equations and present some effective numerical quadrature to handle the singular integrals.

5.1. Parametrization

For simplicity, the boundary $\Gamma_{D}$ is assumed to be a starlike curve with the parametrized form

[TABLE]

where $\Omega=\{\hat{x}(t)=(\cos t,\sin t)^{\top};~{}0\leq t<2\pi\}$. We introduce the parametrized integral operators which are still represented by $S_{\kappa}$, $S_{\kappa}^{\infty}$, $K_{\kappa}$ and $H_{\kappa}$ for convenience, i.e.,

[TABLE]

where $\vartheta(\varsigma)=G(\varsigma)g(p(\varsigma))$, $G(\varsigma):=|p^{\prime}(\varsigma)|=\sqrt{(r^{\prime}(\varsigma))^{2}+r^{2}(\varsigma)}$ is the Jacobian of the transformation,

[TABLE]

and

[TABLE]

Thus, (3.12) can be reformulated as the parametrized integral equations

[TABLE]

where $w_{j}=2(f_{j}\circ p)$, $\varphi_{j}=(g_{j}\circ p)$, ${\varphi}^{\prime}_{j}=(g_{j}\circ p)^{\prime}$, $j=1,2,3$, and $\tilde{\tau}^{\prime}:=(\tilde{\tau}^{\prime}_{1},\tilde{\tau}^{\prime}_{2})^{\top}$, $\tilde{\nu}^{\prime}:=(\tilde{\nu}^{\prime}_{1},\tilde{\nu}^{\prime}_{2})^{\top}$.

To avoid calculating the derivative of the Jacobi $G$ in numerical discretization, we transform the parametrized integral equations (5.1) to

[TABLE]

where $\tilde{\varphi}_{l}=\varphi_{l}/G$, $l=1,2,3$, $\mathsf{n}^{\prime}=(p_{2}^{\prime\prime},-p_{1}^{\prime\prime})^{\top}$, and ${\mathsf{n}^{\perp}}^{\prime}=(p_{1}^{\prime\prime},p_{2}^{\prime\prime})^{\top}$.

5.2. Discretization

The kernel $\widetilde{M}$ and $\widetilde{K}$ of the parametrized single-layer and normal derivative integral operators can be written in form of

[TABLE]

where

[TABLE]

and the diagonal terms are given as

[TABLE]

with the Euler constant $E_{\rm c}=0.57721\cdots$.

For the kernel $\widetilde{H}$ of parametrized tangential derivative integral operator, analogously to [8], we split the kernel in the form

[TABLE]

where

[TABLE]

turn out to be analytic with the diagonal terms

[TABLE]

Let $\varsigma_{j}^{(n)}:=\pi j/n$, $j=0,\cdots,2n-1$ be an equidistant set of quadrature nodes. For the integral of weakly singular part, by making use of quadrature rule in our previous work [8, eqn. (4.6)], we employ the following quadrature rules

[TABLE]

where the function $Q$ is required to be continuous, and the quadrature weights are given by

[TABLE]

We also refer to [26] for details of (5.3). On the other hand, with the help of trapezoidal rule

[TABLE]

and Lagrange bases

[TABLE]

for the trigonometric interpolation, we derive the following quadrature rules for the integration with derivative involved

[TABLE]

where we have set $d_{m-j}^{(n)}=\mathcal{L}^{\prime}_{j}(\varsigma_{m}^{(n)})$, and the quadrature weights can be given by

[TABLE]

Now, in view of the quadrature rules (5.3)–(5.8), we employ following quadrature operators

[TABLE]

[TABLE]

as the approximation of integral operators $S_{\kappa}[Q\vartheta;p]$, $K_{\kappa}[Q\vartheta;p]$, $K_{\kappa}[Q\vartheta^{\prime};p]$, $H_{\kappa}[Q\vartheta;p]$ and $H_{\kappa}[Q\vartheta^{\prime};p]$.

To obtain a Nyström-type of discretization, we express the following combination

[TABLE]

with unknowns $\varUpsilon^{(l)}_{j}:=\tilde{\varphi}_{l}(\varsigma_{j}^{(n)})$ as finite dimensional approximation of the densities $\tilde{\varphi}_{l}$, $l=1,2,3$. Then, the derivative $\tilde{\varphi}^{\prime}_{l}$ can be approximate by

[TABLE]

Hence, the full discretization of (5.2) can be deduced as the form

[TABLE]

where

[TABLE]

[TABLE]

with

[TABLE]

and $w_{l,i}^{(n)}=w_{l}(\varsigma_{i}^{(n)})$ for $i=0,\cdots,2n-1$, $l=1,2,3$.

6. Reconstruction methods

In this section, we introduce the iterative methods and the algorithms for the phased and phaseless IAEIP.

6.1. Iterative method for the phased IAEIP

We assume that the field equations are (3.12), and the data equation is given by

[TABLE]

Thus, the field equations and data equation can be reformulated as the parametrized integral equations (5.1) and

[TABLE]

where $w_{\infty}=u_{\infty}\circ p$.

In the reconstruction process, when an approximation of the boundary $\Gamma_{D}$ is available, the field equations (5.1) are solved for the densities $\varphi_{l}$, $l=1,2,3$. Once the approximated densities $\varphi_{l}$ are computed, the update of the boundary $\Gamma_{D}$ can be obtained by solving the linearized data equation (6.1) with respect to $\Gamma_{D}$.

6.1.1. Iterative scheme

The linearization of (6.1) with respect to a given $p$ requires the Fréchet derivative of the parameterized integral operator $S_{\kappa}^{\infty}$, which can be easily computed and is given by

[TABLE]

where $q(\varsigma)=(\Delta c_{1},\Delta c_{2})+\Delta r(\varsigma)(\cos\varsigma,\sin\varsigma)$ denoted as the update of the boundary $\Gamma_{D}$. Then, the linearization of (6.1) leads to

[TABLE]

where

[TABLE]

As usual for iterative algorithms, the stopping criteria is necessary to justify the convergence numerically. With regard to our iterative procedure, the relative error estimator is chosen as follows

[TABLE]

for some sufficiently small parameter $\epsilon>0$ depending on the noise level, where $p^{(k)}$ is the $k$th approximation of the boundary $\Gamma_{D}$.

We are now in a position to present the iterative algorithm for the inverse obstacle scattering problem with phased far-field data as Algorithm I.

6.1.2. Discretization

We use the Nyström-type method which is described in Section 5 for the full discretizations of (5.1). Now we discuss the discretization of the linearized equation (6.3) and obtain the update by using the least squares with Tikhonov regularization [25]. As for a finite dimensional space to approximate the radial function $r$ and its update $\Delta r$, we choose the space of trigonometric polynomials of the form

[TABLE]

where the integer $M>1$ is the truncation number. For simplicity, we reformulate the equation (6.3) by introducing the following definitions

[TABLE]

Then, by combining (6.2) and (6.3) together and using trapezoidal rule (5.5), we get the discretized linear system

[TABLE]

to determine the real coefficients $\Delta c_{1}$, $\Delta c_{2}$, $\alpha_{m}$ and $\beta_{m}$, where $\varsigma_{i}^{(\tilde{n})}:=\pi i/\tilde{n},~{}i=0,1,\cdots,2\tilde{n}-1$ are the far-field observation points in $[0,2\pi]$,

[TABLE]

for $l=1,2$, and

[TABLE]

In general, $2M+1\ll 2\tilde{n}$, and due to the ill-posedness, the overdetermined system (6.5) is solved via the Tikhonov regularization. Hence the linear system (6.5) is reformulated into minimizing the following function

[TABLE]

with $H^{2}$ penalty term, where $\lambda>0$ is a regularization parameter. It is easy to show that the minimizer of (6.6) is the solution of the system

[TABLE]

where

[TABLE]

Thus, we obtain the new approximation

[TABLE]

6.2. Iterative method for the phaseless IAEIP

To incorporate the reference ball, we find the solution of (2.6) with $D$ replaced by $D\cup B$ in the form of single-layer potentials with densities $g_{1,\sigma}$, $g_{2,\sigma}$ and $g_{3,\sigma}$:

[TABLE]

for $x\in\mathbb{R}^{2}\setminus\Gamma_{D\cup B}$, where $\sigma=D,B$.

Furthermore, we introduce following integral operators, i.e.,

[TABLE]

where $\varrho=D,B$. Then, letting $x\in\mathbb{R}^{2}\setminus\overline{D\cup B}$ tend to boundaries $\Gamma_{D}$ and $\Gamma_{B}$ respectively in (6.8)–(6.10), and making use of the jump relation of the single-layer potentials and the boundary condition of (2.6) for $\Gamma_{D\cup B}$, we can readily deduce the following field equations in the operator form on $\Gamma_{D}$:

[TABLE]

and on $\Gamma_{B}$ the field equations are the same as above with superscript/subscript $D$ and $B$ interchanged. The phaseless data equation is given by

[TABLE]

In the reconstruction process, the field equations are solved for $g_{1,\sigma},g_{2,\sigma}$ and $g_{3,\sigma}$ with an approximation of the boundary $\Gamma_{D}$. Then, by keeping $g_{1,\sigma},g_{2,\sigma}$ and $g_{3,\sigma}$ fixed, the update of the boundary $\Gamma_{D}$ can be obtained by linearizing (6.14) with respect to $\Gamma_{D}$.

6.2.1. Parametrization and iterative scheme

For simplicity, the boundary $\Gamma_{D}$ and $\Gamma_{B}$ are assumed to be starlike curves with the parametrized form

[TABLE]

where $\Omega=\{\hat{x}(t)=(\cos t,\sin t)^{\top};~{}0\leq t<2\pi\}$. We assume $G_{D}(\varsigma):=|p^{\prime}(\varsigma)|=\sqrt{(r^{\prime}(\varsigma))^{2}+r^{2}(\varsigma)}$ and $G_{B}=R$ denoted as the Jacobian of the transformation.

Now, we reformulate the phaseless data equation (6.14) as the parametrized integral equations

[TABLE]

where $\varphi_{3,\sigma}=g_{3}\circ p_{\sigma}$, $\sigma=D,B$. By recalling the Fréchet derivative operator ${S_{\kappa}^{\infty}}^{\prime}[p;\varphi]q$ in (6.2), the linearization of (6.15) leads to

[TABLE]

where

[TABLE]

Again, with regard to our iterative procedure, the relative error estimator is chosen as following

[TABLE]

for some sufficiently small parameter $\epsilon>0$ depending on the noise level, where $p_{D}^{(k)}$ is the $k$th approximation of the boundary $\Gamma_{D}$.

The iterative algorithm for the phaseless IAEIP is given by Algorithm II.

6.2.2. Discretization

Noting that the kernels of $S^{\sigma,\varrho}_{\kappa}$, $K^{\sigma,\varrho}_{\kappa}$ and $H^{\sigma,\varrho}_{\kappa}$ are weakly singular when $\sigma=\varrho$. With the help of quadrature rules (5.3)–(5.8), the full discretization of (6.11)–(6.13) can be handled the same as those described in Section 5.

In addition, we introduce the following definition

[TABLE]

Then, we get the discretized linear system

[TABLE]

to determine the real coefficients $\Delta c_{1}$, $\Delta c_{2}$, $\alpha_{m}$ and $\beta_{m}$, where

[TABLE]

for $l=1,2$, and

[TABLE]

Similarly, the overdetermined system (6.18) is also solved via the Tikhonov regularization with $H^{2}$ penalty term which is introduced in Section 6.1.2.

7. Numerical experiments

In this section, we present some numerical examples to illustrate the feasibility of the iterative reconstruction methods. We use a single plane wave to illuminate the obstacle. The synthetic far-field data and phaseless far-field data are numerically generated by the Nyström-type method described in Section 5. In order to avoid the inverse crime, the number of quadrature nodes used in the inverse solver ($n=64$) is chosen to be different from that of the forward solver ($n=100$). The noisy data $u_{\infty,\delta}$ and $|u_{\infty,\delta}|^{2}$ are generated in the following way,

[TABLE]

where $\breve{\eta}=\breve{\eta}_{1}+\mathrm{i}\breve{\eta}_{2}$, $\breve{\eta}_{1}$, $\breve{\eta}_{2}$ and $\eta$ are normally distributed random numbers ranging in $[-1,1]$, $\delta>0$ is the relative noise level. In addition, we denote the $L^{2}$ relative error between the reconstructed and exact boundaries by

[TABLE]

In the iteration process, we obtain the update $\xi$ from a scaled Newton step with the Tikhonov regularization and $H^{2}$ penalty term, i.e.,

[TABLE]

where the scaling factor $\rho\geq 0$ is fixed throughout the iterations. Analogously to [9], the regularization parameter $\lambda$ in (6.7) is chosen as

[TABLE]

In all of the following figures, the exact boundary curves are displayed in solid lines, the reconstructed boundary curves are shown in dashed lines $--$, and all the initial guesses are taken to be a circle which is indicated in the dash-dotted lines $\cdot-$. The incident directions are denoted by directed line segments with arrows. Throughout all the numerical examples, we take $\lambda=3.88,\mu=2.56$, the angular frequency $\omega=0.7\pi$, the scaling factor $\rho=0.9$, and the truncation $M=6$. We present the results for two commonly used examples: an apple-shaped obstacle and a peanut-shaped obstacle. The parametrization of the exact boundary curves for these two obstacles are given in Table 1.

Example 1: The IAEIP with far-field data. We consider the inverse problem of reconstructing an elastic obstacle from far-field data by using Algorithm I. The synthetic far-field data is numerically generated at 128 points, i.e. $\tilde{n}=64$. In Fig. 2 and Fig. 3, the reconstructions of an apple-shaped and a peanut-shaped obstacles with $1\%$ and $5\%$ noise are shown, respectively. Moreover, the relative $L^{2}$ error $Err_{k}$ between the reconstructed and exact boundaries and the error $E_{k}$ defined in (6.4) are also presented with respect to the number of iterations. As we can see from the figures, the trend of two error curves is basically the same for larger number of iteration. Therefore, the choice of the stopping criteria is reasonable. The reconstructions with different initial guesses for the two curves are given in Fig. 4 and Fig. 5, and the reconstructions with different directions of incident waves are presented in Fig. 6 and Fig. 7. As shown in these results, the location and shape of the obstacle could be simultaneously and satisfactorily reconstructed for a single incident plane wave.

Example 2: The IAEIP with phaseless far-field data and a reference ball. By adding a reference ball to the inverse scattering system, we consider the inverse problem of reconstructing an elastic obstacle from phaseless far-field data based on the algorithm of Table 2. The synthetic phaseless far-field data is numerically generated at 64 points, i.e. $\bar{n}=32$. The reconstructions with $1\%$ noise and $5\%$ noise are shown in Fig. 8 and Fig. 9, respectively. Again, the relative $L^{2}$ error $Err_{k}$ and the error $E_{k}$ are presented in the figures. The reconstructions with different initial guesses for the two curves are given in Fig. 10 and Fig. 11. The reconstructions with different reference balls are shown in Fig. 12 and Fig. 13. From this example, we found that the translation invariance property of the phaseless far-field pattern can be broken down by introducing a reference ball. Based on this algorithm, both the location and shape of the obstacle can be satisfactorily reconstructed from the phaseless far-field data for a single incident plane wave.

8. Conclusions

In this paper, we have studied the two-dimensional inverse acoustic scattering problem by an elastic obstacle with the phased and phaseless far-field data for a single incident plane wave. Based on the Helmholtz decomposition, the coupled acoustic-elastic wave equation is reformulated into a coupled boundary value problem of the Hemholtz equations, and the uniqueness of the solution for this boundary problem is proved. We investigate the jump relations for the second derivatives of single-layer potential and establish coupled boundary integral equations. We prove the well-posedness of the solution for the coupled boundary integral equations, and develop an efficient and accurate Nyström-type discretization to solve the coupled system. The method of nonlinear integral equations is developed for the inverse problem. In addition, we show that the phaseless far-field pattern is invariant under translation of the obstacle. To locate the obstacle, an elastic reference ball is introduced to the scattering system in order to break the translation invariance. We establish the uniqueness for the IAEIP with phaseless far-field pattern. A reference ball technique based nonlinear integral equations method is proposed for the inverse problem. Numerical results show that the location and shape of the obstacle can be satisfactorily reconstructed. Future work includes the uniqueness for the phaseless inverse scattering with one incident plane wave and the extension of the method to the three-dimensional inverse scattering problem.