Multi-frequency iterative methods for the inverse medium scattering problems in elasticity

TL;DR

This paper develops multi-frequency iterative algorithms to accurately reconstruct elastic parameters of an inhomogeneous medium in elasticity using scattering data, with numerical validation showing effectiveness.

Contribution

It introduces two Landweber iterative algorithms leveraging multi-frequency data and phaseless measurements for elastic parameter reconstruction.

Findings

Algorithms successfully reconstruct elastic parameters with high accuracy.

Plane pressure waves improve reconstruction quality.

Numerical examples validate the proposed methods.

Abstract

This paper concerns the reconstruction of multiple elastic parameters (Lam\'e parameters and density) of an inhomogeneous medium embedded in an infinite homogeneous isotropic background in $\mathbb{R}^2$. The direct scattering problem is reduced to an equivalent system on a bounded domain by introducing an exact transparent boundary condition and the wellposedness of the corresponding variational problem is established. The Fr\'{e}chet differentiability of the near-field scattering map is studied with respect to the elastic parameters. Based on the multi-frequency measurement data and its phaseless term, two Landweber iterative algorithms are developed for the reconstruction of the multiple elastic parameters. Numerical examples, indicating that plane pressure incident wave is a better choice, are presented to show the validity and accuracy of our methods.

Figures

Tables

Table 1. Table 1: Final reconstruction errors of q λ subscript 𝑞 𝜆 q_{\lambda} , q μ subscript 𝑞 𝜇 q_{\mu} and q ρ subscript 𝑞 𝜌 q_{\rho} for Examples 1-4.

Figure	Error of $q_{\lambda }$	Error of $q_{\mu }$	Error of $q_{\rho }$
2	0.74	0.36	0.39
3	0.40	0.44	0.39
4(a,b,c)	0.43	0.47	0.42
4(d,e,f)	0.44	0.47	0.42
5	0.47	0.51	0.48
6	0.43	0.46	0.37
7	0.66	0.73	0.61

Table 2. Table 2: Final reconstruction errors of q ρ subscript 𝑞 𝜌 q_{\rho} for Examples 5-6.

Figure	9	10	11(a)	11(b)	12(a)	12(b)	12(c)	13
Error of $q_{\rho }$	0.035	0.24	0.80	0.78	0.14	0.65	0.59	0.46

Equations

$$\text{mbox}supp\{q\}\subset B_R,q\in \mathcal{K}:=\{q\in L^{\infty }(B_R{)}^3:-1<C_1\le q_{\lambda },q_{\mu },q_{\rho }\le C_2<\infty \}.$$

$$\nabla \cdot {\sigma }_{0,0}(u^{in})+{\rho }_0{\omega }^2{u}^{in}=0\text{mbox}in{\mathbb{R}}^2,$$

$${\sigma }_{q_{\lambda },q_{\mu }}(u):={\lambda }_0(1+q_{\lambda })(\text{mbox}divu)\mathbf{I}+2{\mu }_0(1+q_{\mu })\mathcal{E}(u),\mathcal{E}(u):=\frac{1}{2}\left(\nabla u+\nabla u^{\top }\right),$$

$${\Delta }_{{\lambda }_0,{\mu }_0}^{\ast }{u}^{in}+{\rho }_0{\omega }^2{u}^{in}=0\text{mbox}in{\mathbb{R}}^2,$$

$${\Delta }_{\lambda ,\mu }^{\ast }=\mu \text{mbox}div\text{mbox}grad+(\lambda +\mu )\text{mbox}grad\text{mbox}div.$$

$$\displaystyle u^{in}=u_{s}^{in}:=d^{\perp}e^{ik_{s}x\cdot d},\quad d=(\cos\theta^{in},\sin\theta^{in})^{\top}\in{\color[rgb]{0.000,0.000,0.000}\Gamma_{1}},{\color[rgb]{0.000,0.000,0.000}d^{\perp}=(-\sin\theta^{in},\cos\theta^{in})^{\top},}$$

$$\displaystyle u^{in}=u_{p}^{in}:=de^{ik_{p}x\cdot d},\quad d\in{\color[rgb]{0.000,0.000,0.000}\Gamma_{1}},$$

$$k_s=\omega \sqrt{\frac{{\rho }_0}{{\mu }_0}},k_p=\omega \sqrt{\frac{{\rho }_0}{{\lambda }_0+2{\mu }_0}},$$

$$\nabla \cdot {\sigma }_{q_{\lambda },q_{\mu }}(u)+{\omega }^2\rho u=0\text{mbox}in{\mathbb{R}}^2.$$

$$u^{sc}=u_p^{sc}+u_s^{sc},u_p^{sc}=-\frac{1}{k_p^2}\text{mbox}grad\text{mbox}div{u}^{sc},u_s^{sc}=\frac{1}{k_s^2}\overrightarrow{\text{mbox}curl}\text{mbox}curl{u}^{sc},$$

$$\text{mbox}curlv={\partial }_1{v}_2-{\partial }_2{v}_1,v=(v_1,v_2{)}^{\top },\overrightarrow{\text{mbox}curl}f:=({\partial }_{2}f,-{\partial }_{1}f{)}^{\top }.$$

$$\displaystyle(\Delta+k_{\alpha}^{2})\,u_{\alpha}^{sc}=0,\qquad\alpha=p,s,\qquad{\rm div\,}u_{s}^{sc}=0,\quad{\color[rgb]{0.000,0.000,0.000}\mbox{curl}}\,u_{p}^{sc}=0.$$

$$\underset{r\to \infty }{\lim }{r}^{\frac{1}{2}}\left(\frac{\partial u_t^{sc}}{\partial r}-i{k}_t{u}_t^{sc}\right)=0,r=|x|,t=p,s,$$

$${\Delta }_{{\lambda }_0,{\mu }_0}^{\ast }{v}^{sc}+{\rho }_0{\omega }^2{v}^{sc}$$

$$v^{sc}$$

$$T_{{\lambda }_0,{\mu }_0}u:=\nu \cdot {\sigma }_{0,0}(u)=2{\mu }_0{\partial }_{\nu }u+{\lambda }_0\nu \mathrm{div}u-{\mu }_0{\nu }^{\perp }\mathrm{curl}\mathrm{u},$$

$${\alpha }_n(t_{\xi }):=\frac{{H_n^{(1)}}^{\prime }(t_{\xi })}{H_n^{(1)}(t_{\xi })},{\beta }_n(t_{\xi }):=\frac{{H_n^{(1)}}^{\prime \prime }(t_{\xi })}{H_n^{(1)}(t_{\xi })},\xi =p,s,$$

$$t_{\xi }=k_{\xi }R.$$

$$\mathcal{B}w=\sum _{n\in \mathbb{Z}}\frac{1}{2\pi R}{M}_{\theta }^{\top }{W}_n{\int }_0^{2\pi }{M}_{\varphi }w(R,\varphi )e^{in(\theta -\varphi )}d\varphi ,$$

$$M_{\theta }:=\left[\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right]$$

$$W_n=B_n{A}_n^{-1},$$

$$A_n:=\left[\begin{array}{cc}{t}_p{\alpha }_n(t_p)& in\\ in& -t_s{\alpha }_n(t_s)\end{array}\right],$$

$$B_n:=\left[\begin{array}{cc}2{\mu }_0{t}_p^2{\beta }_n(t_p)-{\lambda }_0{t}_p^2& 2i{\mu }_{0}n\left(t_s{\alpha }_n(t_s)-1\right)\\ 2i{\mu }_{0}n\left(t_p{\alpha }_n(t_p)-1\right)& -2{\mu }_0{t}_s^2{\beta }_n(t_s)-{\mu }_0{t}_s^2\end{array}\right].$$

$$\mathcal{B}\phi :=\sum _{n\in \mathbb{Z}}\frac{1}{2\pi R}{M}_{\theta }^{\top }[W_n{]}^{\top }{\int }_0^{2\pi }{M}_{\varphi }\phi e^{in(\varphi -\theta )}d\varphi ,$$

$${\mathcal{B}}^{\ast }\phi =\overline{\mathcal{B}\overline{\phi }}.$$

$$⟨{\mathcal{B}}^{\ast }\phi ,\psi {⟩}_{{\Gamma }_R}$$

$$\mathcal{B}w$$

$$W_{-n}=B_{-n}{A}_{-n}^{-1},$$

$$A_{-n}:=\left[\begin{array}{cc}{t}_p{\alpha }_n(t_p)& -in\\ -in& -t_s{\alpha }_n(t_s)\end{array}\right],$$

$$B_{-n}:=\left[\begin{array}{cc}2{\mu }_0{t}_p^2{\beta }_n(t_p)-{\lambda }_0{t}_p^2& -2i{\mu }_{0}n\left(t_s{\alpha }_n(t_s)-1\right)\\ -2i{\mu }_{0}n\left(t_p{\alpha }_n(t_p)-1\right)& -2{\mu }_0{t}_s^2{\beta }_n(t_s)-{\mu }_0{t}_s^2\end{array}\right].$$

$${H_n^{(1)}}^{\prime \prime }(z)=\left(\frac{n^2}{z^2}-1\right)H_n^{(1)}(z)-\frac{1}{z}{H_n^{(1)}}^{\prime }(z),$$

$${\beta }_n(t_p)=\frac{n^2}{t_p^2}-1-\frac{1}{t_p}{\alpha }_n(t_p),{\beta }_n(t_s)=\frac{n^2}{t_s^2}-1-\frac{1}{t_s}{\alpha }_n(t_s).$$

Full text

Keywords: Elastic wave, inverse medium problem, iterative method, multi-frequency

1 Introduction

Time-harmonic elastic scattering problems play important roles in many fields of applications and the linear elasticity theory provides an essential tool for analysis and design of mechanic systems and engineering structures ([3, 21]). In this paper, we consider several inhomogeneous isotropic elastic bodies embedded in an infinite homogeneous isotropic background medium in ${\mathbb{R}}^{2}$. Denote by $\lambda$, $\mu$ the Lamé parameters with $\mu>0$, $\lambda>0$, and by $\rho$ the density of the elastic medium. Suppose that $\lambda=\lambda_{0}(1+q_{\lambda})$, $\mu=\mu_{0}(1+q_{\mu})$ and $\rho=\rho_{0}(1+q_{\rho})$ where $\lambda_{0}$, $\mu_{0}$ and $\rho_{0}$ are constants representing the Lamé parameters and density of the background elastic medium. Denote $B_{R}:=\{x\in{\mathbb{R}}^{2}:|x|<R\}$ and $\Gamma_{R}:=\partial B_{R}$. Set $q:=(q_{\lambda},q_{\mu},q_{\rho})^{\top}$. Throughout, we make the following assumption:

Assumption:

there exists some $R>0$ and constants $C_{1},C_{2}>0$ such that

[TABLE]

Let $u^{in}$ be a plane incident field satisfying

[TABLE]

where $\omega>0$ is the frequency and the stress tensor $\sigma_{q_{\lambda},q_{\mu}}(u)$ is defined as

[TABLE]

and ${\bf I}$ stands for the $2\times 2$ identity matrix. The elliptic equation (1.1) can be restated as

[TABLE]

where the Lamé operator $\Delta_{\lambda,\mu}^{*}$ is defined as

[TABLE]

In this paper, the incident wave is allowed to be either a plane shear wave taking the form

[TABLE]

or a plane pressure wave taking the form

[TABLE]

where

[TABLE]

are the wave numbers of pressure wave and shear wave, respectively and $d$ and $\theta^{in}$ are referred as the direction and angle of the incidence, respectively.

The total displacement field $u=(u_{1},u_{2})^{\top}$ can be modeled by the reduced Navier equation

[TABLE]

Since the background medium is unbounded, an appropriate radiation condition at infinity must be imposed on the scattered field $u^{sc}:=u-u^{in}$ to ensure well-posedness of the scattering problem. The scattered field in ${\mathbb{R}}^{2}\backslash\overline{B_{R}}$ can be decomposed into the sum of the compressional (longitudinal) part ${u}^{sc}_{p}$ and the shear (transversal) part ${u}^{sc}_{s}$ as follows :

[TABLE]

where the two-dimensional operators curl and $\overrightarrow{\mbox{curl}}$ are defined respectively by

[TABLE]

It then follows from the decompositions in (1.3) that

[TABLE]

The scattered field is required to satisfy the Kupradze radiation condition (see e.g. [31])

[TABLE]

uniformly with respect to all \hat{x}=x/|x|\in{\color[rgb]{0.000,0.000,0.000}\Gamma_{1}}.

Given the incident field $u^{in}$, the direct problem is to determine the scattered field $u^{sc}$ for the known elastic parameters $\lambda,\mu,\rho$. In practice, the original boundary value problem for the Lamé system can be reduced to an equivalent system on a bounded domain via introducing an exact transparent boundary condition (TBC) on an artificial boundary enclosing the inhomogeneous bodies. The TBC can be formulated by the so-called Dirichlet-to-Neumann (DtN) map taking in the form of a Fourier series ([8, 23, 32, 33]). Based on the properties of the DtN map and Fredholm alternative theorem, the uniqueness and existence of weak solutions of the equivalent system can be derived(Section 2, Theorem 2.4).

The main purpose of this paper is to study numerical algorithm for the inverse medium problem in the elastic scattering, that is, to determine the unknown elastic parameters $q_{\lambda},q_{\mu},q_{\rho}$ from the measurements of near-field data $u|_{\Gamma_{R}}$, given the incident field $u^{in}$. For static case, i.e., $\omega=0$, uniqueness of the inverse medium problem has been investigated under appropriate assumptions on the elastic parameters in [2, 18, 22, 28, 35, 37, 36]. For the time-harmonic case, Beretta et al ([19]) proves uniqueness when the Lamé parameters and the density are assumed to be piecewise constant on a given domain partition. For the mathematical analysis of the stability for the inverse medium problems in elasticity, we refer to [1, 18, 27, 35, 36, 37].

Recently, it has been realized that the use of multi-frequency data is an effective approach to overcome the major difficulties associated with the inverse medium problems: the ill-posedness and the presence of many local minima. Based on the multi-frequency measurements and the Fréchet derivative of the solution operator, a stable recursive linearization method is proposed in [10] for the inverse medium problems in acoustics, see also in [12, 11, 13, 14]. The idea to use multi-frequency data has also been widely developed in proving uniqueness and increasing stability for inverse source problems in acoustics, elastodynamics and electromagnetics ([9, 5, 15, 6, 16, 20, 26, 29, 34, 33]). It still remains open whether or not the multi-frequency measurements uniquely determine the Lamé parameters and the density? This paper is designed to study the capability of iterative methods for inverse medium problems in elasticity using multi-frequency measurements. Compared with the acoustic and electromagnetic case, the elasticity problem appears to be more complicated because of the coexistence of pressure and shear waves that propagate at different speeds and the aim to reconstruct multiple parameters. In particular, if the Lamé parameters are constants, the inverse medium problem in elasticity is consistent with that in acoustics, see Section 4. Relying on the variational arguments for the direct problems, we derive the Fréchet derivative of the solution operator with respect to the elastic parameters and investigate the adjoint of the Fréchet derivative. Then we employ the Landweber iterative method based on the multi-frequency measurements to find the unknown elastic parameters. At each iteration step, the forward problem and an adjoint one need to be solved and the correctness of the parameters needs to be evaluated.

The outline of the paper is as follows. In Section 2, we derive the well-posedness of the direct scattering problem using variational approach and investigate the Fréchet differentiability of the near-field scattering map. We develop the Landweber iterative methods for solving the inverse medium problem in Section 3. In Section 4, we give a brief discussion about a special case that $q_{\lambda}=q_{\mu}=0$. Numerical examples are presented in Section 5.

2 Direct scattering problem

In this section, we discuss the well-posedness of the direct elastic scattering problem and investigate the Fréchet differentiability of the near-field scattering map.

2.1 Variational formulation

We first reduce the original scattering problem described in section 1 on a bounded domain via introducing a TBC on an artificial boundary $\Gamma_{R}$ enclosing the inhomogeneity inside. The TBC is formulated by the so-called DtN map defined as follows.

Definition 2.1**.**

For any $w\in H^{1/2}(\Gamma_{R})^{2}$, The DtN map $\mathcal{B}$ applied to $w$ is defined as $T_{\lambda_{0},\mu_{0}}v^{sc}|_{\Gamma_{R}}$, where $v^{sc}$ satisfies

[TABLE]

and the Kupradze radiation condition. Here, $T_{\lambda_{0},\mu_{0}}$ is the traction operator defined by

[TABLE]

where $\nu=(\nu_{1},\nu_{2})^{\top}$ denotes the exterior unit normal vector to $\Gamma_{R}$ and the corresponding tangential vector is given by $\nu^{\perp}:=(-\nu_{2},\nu_{1})^{\top}$.

The DtN map $\mathcal{B}$ is well-defined since the Dirichlet-kind boundary value problem (2.1)-(2.2) is uniquely solvable in $H^{1}_{loc}({\mathbb{R}}^{2}\backslash\overline{B_{R}})^{2}$, see Corollary 2.3 in [8]. For all $n\in{\mathbb{Z}}$, denote

[TABLE]

where

[TABLE]

Following the procedure described in [8], it can be derived that

[TABLE]

where the matrix $M_{\theta}$ dependent on the angle $\theta\in[0,2\pi)$ is defined as

[TABLE]

and the coefficient matrix $W_{n}$ is given by

[TABLE]

with

[TABLE]

For the invertibility of matrix $A_{n}$, we refer to Lemma 2.11 in [8]. Denote by $\Lambda_{n}$ the determinant of $A_{n}$.

Lemma 2.2**.**

For all $\varphi\in(H^{1/2}(\Gamma_{R}))^{2}$, the DtN mapping $\mathcal{B}$ can be expressed equivalently as

[TABLE]

In addition, it holds that

[TABLE]

Proof.

We first prove (2.6). We can derive that for all $\varphi,\psi\in(H^{1/2}(\Gamma_{R}))^{2}$,

[TABLE]

where $\langle\cdot,\cdot\rangle_{\Gamma_{R}}$ is the $L^{2}$ duality pairing between $(H^{-1/2}(\Gamma_{R}))^{2}$ and $(H^{1/2}(\Gamma_{R}))^{2}$. It remains to prove (2.5). It follows from (2.3) that

[TABLE]

where

[TABLE]

with

[TABLE]

It follows from the properties of Hankel function that

[TABLE]

giving rise to the identities

[TABLE]

From the expressions of $A_{n}^{-1}$ and $B_{n}$ we get the entries $W_{n}^{(i,j)}$ of $W_{n}$, given by (see also Lemma 2.13 in [8])

[TABLE]

Then we can obtain that $W_{-n}^{(1,1)}=W_{n}^{(1,1)}$, $W_{-n}^{(2,2)}=W_{n}^{(2,2)}$ and $W_{-n}^{(1,2)}=W_{n}^{(2,1)}$. These further imply that $W_{-n}=[W_{n}]^{\top}$ which completes the proof. ∎

The following property of the coefficient matrix $W_{n}$ can be obtained directly from Lemma 2.13 in [8].

Lemma 2.3**.**

The matrix $\widetilde{W}_{n}=-(W_{n}+W_{n}^{*})/2$ is positive definite for sufficiently large $|n|$.

Using the DtN map, we can impose the following TBC for the scattered field

[TABLE]

Note that

[TABLE]

Then the original scattering problem is equivalently reduced to the following nonlocal boundary value problem

[TABLE]

where $g:=T_{\lambda_{0},\mu_{0}}u^{in}-\mathcal{B}u^{in}$. Then the variational formulation of (2.7)-(2.8) reads as follows: find $u=(u_{1},u_{2})^{\top}\in(H^{1}(B_{R}))^{2}$ such that

[TABLE]

where the sesquilinear form $a_{q}(\cdot,\cdot):(H^{1}(B_{R}))^{2}\times(H^{1}(B_{R}))^{2}\rightarrow{\mathbb{C}}$ is defined by

[TABLE]

and

[TABLE]

The double dot notation appeared in $B_{q_{\mu}}$ is understood in the following way. If tensors ${\bf A}$ and ${\bf B}$ have rectangular Cartesian components $a_{ij}$ and $b_{ij}$, $i,j=1,2$, respectively, then the double contraction of ${\bf A}$ and ${\bf B}$ is

[TABLE]

The well-posedness of the variational formulation (2.9) is a consequence of the following theorem.

Theorem 2.4**.**

For any $f\in L^{2}(B_{R})^{2}$ and $g\in H^{-1/2}(\Gamma_{R})^{2}$, there exists a unique weak solution $u\in H^{1}(B_{R})^{2}$ to the boundary value problem

[TABLE]

and

[TABLE]

where $\gamma_{q}>c_{0}>0$ is a constant independent of $u$.

Proof.

We know from Lemma 2.3 that $-\mbox{Re}(W_{n})$ is positive definite for large $|n|$. The operator $-\mathcal{B}$ can be decomposed into the sum of an operator $\mathcal{B}_{1}$ whose real part is positive definite and a finite rank operator $\mathcal{B}_{2}$ from $H^{1/2}(\Gamma_{R})^{2}$ to $H^{-1/2}(\Gamma_{R})^{2}$. Define the sesquilinear form $A(\cdot,\cdot):H^{1}(B_{R})^{2}\times H^{1}(B_{R})^{2}\rightarrow{\mathbb{C}}$ as

[TABLE]

Then the variational equation corresponding to (2.11)-(2.12) reads:

[TABLE]

We split the sequilinear form $A$ into the sum $A=A_{1}+A_{2}$, where

[TABLE]

Recalling the Korn’s inequality (see e.g. [25]), we have

[TABLE]

with some constant $c_{1}>0$. Moreover, applying the Cauchy-Schwarz inequality yields

[TABLE]

for some constant $c_{2}>0$. From the compact imbedding $H^{1}(B_{R})\hookrightarrow L^{2}(B_{R})$ and the compactness of $\mathcal{B}_{2}$, we conclude that the sesquilinear form $A$ is strongly elliptic (see Definition 2.5) over $H^{1}(B_{R})^{2}\times H^{1}(B_{R})^{2}$. The sesquilinear form $A$ obviously generates a continuous linear operator $\mathcal{A}:H^{1}(B_{R})^{2}\rightarrow(H^{1}(B_{R})^{2})^{\prime}$ such that

[TABLE]

Here $(H^{1}(B_{R})^{2})^{\prime}$ denotes the dual space of $H^{1}(B_{R})^{2}$ with respect to the duality $\langle\cdot,\cdot\rangle$ extending the $L^{2}$ scalar product in $L^{2}(B_{R})^{2}$. We know from the Rellich’s lemma in elasticity (see Lemma 2.14 in [8]) that the homogeneous operator equation $\mathcal{A}u=0$ has only the trivial solution $u=0$. Then it follows from the Fredholm alternative that the variational formulation (2.9) is uniquely solvable. Finally, the inf-sup condition

[TABLE]

with some constant $\gamma_{q}>0$ generated from the general theory in Babuška and Aziz[4] implies the estimate (2.13). In fact, for $F\in(H^{1}(B_{R})^{2})^{\prime}$, considering the operator equation $\mathcal{A}u=F$ related to the variational equation (2.14), it follows from the Babuška’s theory that the operator equation is well-posed if and only if the conditions

[TABLE]

and

[TABLE]

hold and they are equivalent to the uniqueness and existence of solution of the operator equation, respectively. Then it follows from the variational equation (2.14) and trace theorem that

[TABLE]

which completes the proof. ∎

Definition 2.5**.**

A bounded sesquilinear form $a(\cdot,\cdot)$ on some Hilbert space $X$ is called strongly elliptic if there exists a compact form $q(\cdot,\cdot)$ such that

[TABLE]

2.2 Near-field scattering map

For given perturbed parameters $q\in\mathcal{K}$, we define the scattering operator $S:\mathcal{K}\rightarrow H^{1}(B_{R})^{2}$ by $S(q)=u$, where $u\in H^{1}(B_{R})^{2}$ is the unique weak solution of (2.7)-(2.8). It is easily seen that the map $S$ is nonlinear with respect to $q$. A direct application of Theorem 2.4 gives the following result.

Lemma 2.6**.**

The map $S$ is bounded with the estimate

[TABLE]

where $c>0$ is a constant.

Lemma 2.7**.**

Given the perturbed parameters $q_{l}\in\mathcal{K}$, $l=1,2$, we have the estimate

[TABLE]

where $u_{2}$ is the the unique weak solution of the problem (2.7)-(2.8) with perturbed parameters $q_{2}$ and $c>0$ is a constant.

Proof.

Let $u_{1},u_{2}\in H^{1}(B_{R})^{2}$ be the unique weak solutions of the problem (2.7)-(2.8) with perturbed parameters $q_{1}$ and $q_{2}$, respectively. Set $w=u_{2}-u_{1}$. Then we have,

[TABLE]

Then the the inf-sup condition (2.15) together with the boundedness of $a_{q_{2}-q_{1}-1}$ implies the desired estimate. ∎

For any $\delta q:=(\delta q_{\lambda},\delta q_{\mu},\delta q_{\rho})^{\top}\in\mathcal{K}$, assume that $w\in H^{1}(B_{R})^{2}$ is the unique weak solution of the following variational problem

[TABLE]

Let the map $\mathcal{T}_{q}:\mathcal{K}\rightarrow H^{1}(B_{R})^{2}$ be such that

[TABLE]

Lemma 2.8**.**

Given the perturbed parameters $q,\delta q\in\mathcal{K}$, we have the estimate

[TABLE]

where $C>0$ is a constant.

Proof.

Let $u_{1},u_{2}\in H^{1}(B_{R})^{2}$ be the unique weak solutions of the problem (2.7)-(2.8) with perturbed parameters $q$ and $q+\delta q$, respectively. Let $\widetilde{u}=u_{2}-u_{1}$, it follows that

[TABLE]

Then we have the estimate

[TABLE]

where $C>0$ is a constant defined as

[TABLE]

∎

Let $\gamma:H^{1}(B_{R})^{2}\rightarrow H^{1/2}(\Gamma_{R})^{2}$ be the trace operator to the boundary $\Gamma_{R}$ and define the near-field scattering map $N$ as $N(q)=\gamma S(q)$. By combining Lemmas 2.6-2.8, we arrive at the following theorem.

Theorem 2.9**.**

The near-field scattering map $N$ is Fréchet differentiable with respect to $q$ and its Fréchet derivative is $N^{\prime}_{q}=\gamma\mathcal{T}_{q}$.

3 Inverse medium problem

In this section, we consider the inverse medium scattering problem of reconstructing the unknown perturbed elastic parameters and develop a Landweber iterative method. Assume that the total-field data $u$ is available over a range of frequencies $\omega\in[\omega_{min},\omega_{max}]$ which can be divided into $\omega_{min}=\omega_{1}<\omega_{2}<\cdots<\omega_{N-1}<\omega_{N}=\omega_{max}$ and over a range of incident directions $\theta^{in}\in[\theta_{min},\theta_{max}]$ which can be divided into $\theta_{min}=\theta_{1}<\theta_{2}<\cdots<\theta_{M-1}<\theta_{M}=\theta_{max}$. Let $u^{i,j}:=u(x,\omega_{i},\theta_{j})|_{\Gamma_{R}}$ be the unique solution of the direct scattering problem with $\omega=\omega_{i}$ and $\theta^{in}=\theta_{j}$. Consider the following inverse medium problem:

(IP): Given the elastic parameters $\lambda_{0},\mu_{0}$ and $\rho_{0}$ of the background medium, reconstruct the perturbed Lamé parameters $q_{\lambda},q_{\mu}$ and perturbed density $q_{\rho}$ from the multi-frequency measurements $u^{i,j}$, $i=1,\cdots,N$, $j=1,\cdots,M$.

The inverse problem (IP) can be formulated as: Given $\lambda_{0},\mu_{0}$ and $\rho_{0}$, find $q_{\lambda},q_{\mu}$ and $q_{\rho}$ such that

[TABLE]

In particular, the nonlinearity and ill-posedness of the inverse problem cause mathematical challenges from both theoretical and computational points of view. The nonlinearity leads to a nonconvex optimization problem, and the ill-posedness requires certain form of regularization to get a reasonable approximation. Here, to solve the operator equation (3.1) we apply the Landweber iteration method taking the form ([24])

[TABLE]

where $\alpha$ is the step size parameter.

3.1 Multi-frequency iterative algorithm

For the adjoint of the operator $N_{q}^{\prime}$, we have the following result.

Theorem 3.1**.**

Let $u\in H^{1}(B_{R})^{2}$ be the unique weak solution of (2.7)-(2.8). Then for any $h\in H^{1/2}(\Gamma_{R})^{2}$,

[TABLE]

where $\varphi\in H^{1}(B_{R})^{2}$ is the unique weak solution of the following boundary value problem

[TABLE]

Proof.

For any $\delta q\in\mathcal{K}$, let $w\in(H^{1}(B_{R}))^{2}$ be the unique weak solution of the following variational problem

[TABLE]

Replacing $v$ by $\varphi$, we obtain

[TABLE]

in which we have used the relation (2.6) and $\langle\cdot,\cdot\rangle_{s,-s}$ denotes the $L^{2}$ duality pairing between $H^{s}(\Gamma)^{2}$ and $H^{-s}(\Gamma)^{2}$. Since it holds for any $\delta q\in\mathcal{K}$, we complete the proof. ∎

Denote

[TABLE]

where the index $i,j,l$ are related to the frequency $\omega_{i}$, the incident direction $\theta_{j}$ and the current inner Landweber iteration number. Given initial guesses $q_{1,1,0}=0$, we now describe a procedure that determines a better approximation $q_{i,j}:=q_{i,j,L}$ at the frequency $\omega=\omega_{i}$ with incident direction $\theta^{in}=\theta_{j}$ for $i=1,\cdots,N$, $j=1,\cdots,M$ in an increasing manner. For each $i,j$, we apply $L$ steps of Landweber iterations, i.e., $l=1,\cdots,L$. For fixed $i,j$, suppose now that an approximation of the scatterer $q_{i,j,l-1}$ has been recovered. For the recovered scatterer $q_{i,j,l-1}$, we solve at $\omega=\omega_{i}$ and $\theta^{in}=\theta_{j}$ the direct problem

[TABLE]

Then for any $h\in H^{1/2}(\Gamma_{R})^{2}$,

[TABLE]

where $\varphi^{i,j,l}\in H^{1}(B_{R})^{2}$ is the unique weak solution of the following boundary value problem

[TABLE]

Then the Landweber iteration leads to

[TABLE]

Note that the elastic parameters are all real values. Therefore, at each step of iterations, we apply a simple regularization as

[TABLE]

The multi-frequency iterative algorithm for solving the inverse medium scattering is summarized in Algorithm 3.1.

Algorithm 3.1 (Multi-frequency iterative algorithm).

• •

Collect the near-field data over all frequencies $\omega_{i}$, $i=1,\cdots,N$ and all incident directions $d_{j}$, $j=1,\cdots,M$.

• •

Set initial approximations $q_{1,1,0}=0$.

• •

Apply the following iteration:

DO $i=1,\cdots,N$

DO $j=1,\cdots,M$

DO $l=1,\cdots,L$

Update the elastic parameters by the formula

[TABLE]

ENDDO

Set $q_{i,j+1,0}=q_{i,j,L}$

ENDDO

Set $q_{i+1,1,0}=q_{i,M,L}$

ENDDO

3.2 Multifrequency iterative algorithm from phaseless data

We now consider the reconstruction from phaseless data. Define the phaseless near-field scattering map $F:\mathcal{K}\rightarrow H^{1/2}(\Gamma_{R})$ by

[TABLE]

Lemma 3.2**.**

The phaseless near-field scattering map $F$ is Fréchet differentiable with respect to $q$ and its Fréchet derivative is given by

[TABLE]

The adjoint of $F^{\prime}_{q}$ is given in the following Theorem.

Theorem 3.3**.**

Let $u\in H^{1}(B_{R})^{2}$ be the unique weak solution of (2.7)-(2.8). Then for any $\hbar\in H^{1/2}(\Gamma_{R})$,

[TABLE]

where $\psi\in H^{1}(B_{R})^{2}$ is the unique weak solution of the following boundary value problem

[TABLE]

Proof.

The proof is similar as Theorem 3.1 and we omit it here. ∎

We now describe the Landweber iterative algorithm based on phaseless data. Use the same notations in section 3.1. For any $\hbar\in H^{1/2}(\Gamma_{R})$, let $\psi^{i,j,l}\in H^{1}(B_{R})^{2}$ be the unique weak solution of the following boundary value problem

[TABLE]

Then the Landweber iteration leads to

[TABLE]

We conclude the algorithm for solving the inverse medium scattering from phaseless data is summarized in Algorithm 3.2.

Algorithm 3.2 (Multi-frequency iterative algorithm for phaseless data).

• •

Collect the near-field data over all frequencies $\omega_{i}$, $i=1,\cdots,N$ and all incident directions $d_{j}$, $j=1,\cdots,M$.

• •

Set initial approximations $q_{1,1,0}=0$.

• •

Apply the following iteration:

DO $i=1,\cdots,N$

DO $j=1,\cdots,M$

DO $l=1,\cdots,L$

Update the elastic parameters by the formula

[TABLE]

ENDDO

Set $q_{i,j+1,0}=q_{i,j,L}$

ENDDO

Set $q_{i+1,1,0}=q_{i,M,L}$

ENDDO

3.3 Convergence

In this section, we briefly discuss the convergence of the proposed algorithms based on the classical analysis of Landweber iteration method for nonlinear ill-posed problems. For more detailed analysis of the Landweber iteration and its modification form, we refer to [24, 30] and the references therein. It should be pointed out that although we apply the multi-frequency strategy in this paper, it is extremely difficult to investigate the dependence of the convergence of the algorithms on the number and interval of the selected frequencies. For the corresponding analysis of the inverse medium scattering problems in acoustics that take a similar form as the special case discussed in section 4, we refer to [17].

We consider the following operator equation

[TABLE]

where $\mathcal{D}(F)\subset X$ and $X,Y$ are all Hilbert spaces with inner products $(\cdot,\cdot)_{X},(\cdot,\cdot)_{Y}$ and norms $\|\cdot\|_{X},\|\cdot\|_{Y}$, respectively. In particular, $X=L^{2}(B_{R})^{3},Y=H^{1/2}(\Gamma_{R})^{2}$ in this paper. As mentioned above, the operator equation (3.6) is strictly nonlinear and ill-posed. For simplicity, denote $F(q)=F(q,\omega)$. Then the nonlinear Landweber iteration takes the form

[TABLE]

for exact data $y$ and the form

[TABLE]

for inexact data $y^{\delta}$ satisfying $\sup_{\omega\in[\omega_{min},\omega_{max}]}\{\|y^{\delta}-y\|_{Y}\}\leq\delta$. Here, assume that we take the frequency as the innermost loop.

Let $\mathcal{B}_{r}(q_{0})$ denote a closed ball of radius $r$ around $q_{0}$. Assume that $\alpha>0$ is sufficiently small such that $\alpha\|F_{q}^{\prime}\|<1$ in $\mathcal{B}_{2r}(q_{0})$. For simplicity, we consider $\alpha=1$ and $\|F_{q}^{\prime}\|<1$. Resulting from the regularity theory, we can obtain that $u\in H_{loc}^{2}({\mathbb{R}}^{2}\backslash\overline{B_{R-\epsilon}})^{2}$, for any small $\epsilon>0$ such that $\mbox{sup}\{q\}\subset B_{R-\epsilon}$, which implies that $u\in H^{3/2}(\Gamma_{R})^{2}\hookrightarrow H^{1/2}(\Gamma_{R})^{2}$. Thus, the operator $F$ here is compact. According to the analysis in Lemma 2.8 and inverse trace theorem, let $r>0$ and $q,\widetilde{q}\in\mathcal{B}_{2r}(q_{0})\subset\mathcal{D}(F)$ satisfying $\|q-\widetilde{q}\|_{L^{\infty}(B_{R})^{3}}\leq\epsilon$ with sufficiently small $\epsilon>0$, it holds that

[TABLE]

The convergence of the Landweber iteration is given in the following theorem, see Corollary 2.3 and Theorem 2.4 in [30].

Theorem 3.4**.**

Assume that $F(q)=y$ is solvable in $\mathcal{B}_{r}(q_{0})$. Then the nonlinear Landweber iteration (3.7) converges to a solution of $F(q)=y$. Furthermore, if $\mathcal{N}(F_{q^{\dagger}}^{\prime})\subset\mathcal{N}(F_{q}^{\prime})$ for all $q\in\mathcal{B}_{r}(q^{\dagger})$, then $q_{k}$ converges to $q^{\dagger}$ as $k\rightarrow+\infty$.

In case of inexact data, the iteration procedure should be combined with a stopping rule in order to act as a regularization method. For example, for each frequency, one can employ the discrepancy principle, i.e., the iteration is stopped after $k^{\prime}$ steps with

[TABLE]

where $\tau$ is an appropriately chosen positive number. Before stating the convergence of the iteration (3.8) for inexact data case, we need to discuss the selection of $\tau$. Note that the parameter $\eta$ must be dependent on $\omega$.

Proposition 3.5**.**

Assume $F(q)=y$ has a solution $q_{*}\in\mathcal{B}_{r}(q_{0})$ and $q_{k}^{\delta}\in\mathcal{B}_{r}(q_{*})$.

• (i).

A sufficient condition for $q_{k+1}^{\delta}$ generated from (3.8) to be a better approximation of $q_{*}$ than $q_{k}^{\delta}$ is that

[TABLE]

This further leads to the selection of $\tau$ as

[TABLE]

• (ii).

At each frequency, let $k^{\prime}$ be chosen according to the stopping rule, we have

[TABLE]

Proof.

The proof is similar to that of the Proposition 2.2 and Corollary 2.3 in [30]. In fact, from (2.13) in [30], it holds that

[TABLE]

which completes the proof of (i) under the sufficient condition. Furthermore, the special choice of $\tau$ implies that

[TABLE]

Adding up these inequalities for $k$ from 0 through $k^{\prime}-1$ leads to

[TABLE]

Then (ii) can be proved by combining this inequality and the stoping rule. ∎

Now we can have the following convergence result for inexact data, see Theorem 2.6 in [30].

Theorem 3.6**.**

Assume that $F(q)=y$ is solvable in $\mathcal{B}_{r}(q_{0})$ and let $k^{\prime}$ be chosen according to the stopping rule. The nonlinear Landweber iteration (3.8) converges to a solution of $F(q)=y$. If $\mathcal{N}(F_{q^{\dagger}}^{\prime})\subset\mathcal{N}(F_{q}^{\prime})$ for all $q\in\mathcal{B}_{r}(q^{\dagger})$, then $q_{k^{\prime}}^{\delta}$ converges to $q^{\dagger}$ as $\delta\rightarrow 0$.

Remark 3.7**.**

According to the convergence of Landweber iteration, the choice of the initial approximation may affect the efficiency of the algorithms. In fact, in this paper we reconstruct the relative factors of the elastic parameters of the inhomogeneous medium compared with the parameters of the background medium, not the elastic parameters of the inhomogeneous medium themselves, i.e., we reconstruct

[TABLE]

and they are of $\mathcal{O}(1)$. Thus, a simple choice of the initial approximation, that we applied in this paper, is $q_{\lambda}=q_{\mu}=q_{\rho}=0$ which means that there is no inhomogeneity and it can be seen from the numerical examples that this choice works well. As suggested in [5] for the inverse medium problems in acoustics, one can obtain an initial guess from the Born approximation and the scattering data with the higher frequency must be used in order to recover more Fourier modes of the true scatterer.

4 Discussion of a special case

In this section, we briefly discuss the special case that $q_{\lambda}=q_{\mu}=0$. In this case, the considered inverse medium problem in elasticity is consistent with that in acoustics and electromagnetics discussed in [10, 12, 11, 13, 14].

Alternatively, one can consider near-field scattering map $\widetilde{N}$ as \widetilde{N}(q_{\rho})={\color[rgb]{0.000,0.000,0.000}\widetilde{u}}|_{\Gamma_{R}}, where $u=u^{sc}+u^{in}\in H^{1}(B_{R})^{2}$ is the unique weak solution of (2.7)-(2.8). Following the steps discussed in sections 2 and 3 (see also in [10]), we can obtain the following result.

Theorem 4.1**.**

The near-field scattering maps $N,\widetilde{N}$ are Fréchet differentiable with respect to $q_{\rho}$ and their Fréchet derivatives are given by

[TABLE]

where $v\in H^{1}(B_{R})^{2}$ is the unique weak solution of the following boundary value problem

[TABLE]

Moreover, for any $h\in H^{1/2}(\Gamma_{R})^{2}$,

[TABLE]

where $\varphi\in H^{1}(B_{R})^{2}$ is the unique weak solution of the following boundary value problem

[TABLE]

Similar as Algorithm 3.1, we can summarize the multi-frequency iterative algorithm for the reconstruction of inhomogeneous density function $q_{\rho}$ in Algorithm 3.3.

Algorithm 3.3 (Multi-frequency iterative algorithm for $q_{\rho}$)

• •

Collect the near-field data over all frequencies $\omega_{i}$, $i=1,\cdots,N$ and all incident directions $d_{j}$, $j=1,\cdots,M$.

• •

Set an initial approximation $q_{\rho}^{1,1,0}=0$.

• •

Apply the following iteration:

DO $i=1,\cdots,N$

DO $j=1,\cdots,M$

DO $l=1,\cdots,L$

Update the elastic parameters by the formula

[TABLE]

ENDDO

Set $q_{\rho}^{i,j+1,0}=q_{\rho}^{i,j,L}$

ENDDO

Set $q_{\rho}^{i+1,1,0}=q_{\rho}^{i,M,L}$

ENDDO

5 Numerical examples

In this section, we present several numerical examples to verify the accuracy and effectiveness of the recursive algorithm. We always choose $\lambda_{0}=2$, $\mu_{0}=1$, $\rho_{0}=1$ and $R=1$. The near-field measurements are obtained by using finite element method solving the forward scattering problem. Define the relative error

[TABLE]

where $q$ and $\widetilde{q}$ are the true and reconstructed values of the parameter of the scatterer. The true values of the elastic parameters $q_{\lambda},q_{\mu},q_{\rho}$ are shown in Fig.1. Ten equally spaced frequencies are used in the construction, starting from the lowest frequency $\omega_{\min}=1$ and ending at the highest frequency $\omega_{\max}=10$. The number of incident directions is taken as $M=16$ and $\theta_{j}=2(j-1)\pi/M$ for $j=1,\cdots,M$. At each incident direction, 10 Landweber iteration steps are taken for one frequency. Corresponding to the stiffness tensor of the background elastic medium using Voigt notation, we chose the relaxation parameter $\alpha$ as a matrix

[TABLE]

We collect the final reconstruction error for the following examples in Table 1 and 2. In addition, the number of iterations shown in the following figures are from 0 to $N\times M$.

Example 1. We consider the reconstruction from multi-frequency measurements with multiple incident directions. The reconstructed elastic parameters are presented in Fig.2(a,b,c) and Fig.3(a,b,c). The relative errors shown in Fig.2(d,e,f) and Fig.3(d,e,f) indicate that the relative errors decrease as frequency and number of iteration increase. However, it can be seen that if we choose plane shear incident waves, the reconstruction of $q_{\lambda}$, which can further affect the reconstructions of $q_{\mu}$ and $q_{\rho}$, is worse than that using plane pressure incident waves. A possible explanation to this phenomenon is that, in comparison with the plane pressure waves, the plane shear incident waves only contain the information of $\mu_{0}$ and $\rho_{0}$. In the following, we consider the plane pressure incident waves only. To verify the stability of our method, the reconstructions from noised data with noise levels $\delta=3\%,5\%$ are presented in Fig.4.

Example 2. Note that one can chose $\alpha$ as a scale value. For simplicity, we chose $\alpha=0.01/\omega$ and the reconstruction results and the corresponding relative errors are presented in Fig.5. It can be seen that in this example, the iterations using the step size matrix (5.1) is more stable than those using scale step size. It should be pointed out that we can not prove that using the special choice of step size $\alpha$ which looks like the stiffness tensor of the background elastic medium, we can always have better reconstruction than using just a scale one. Our starting pointing is that for isotropic elastic medium, the two Lamé parameters have some connection since they are both determined by the Young’s modulus and Poisson’s ratio and have no direct relation with density. Thus, we take a similar form of stiffness tensor as the choice of step size for example which means that the modification for both $q_{\lambda}$, $q_{\mu}$ are determined by the combinations of the first two components of $(N_{q}^{\prime})*(h)$ and the modification for density $q_{\rho}$ is only determined by the third component of $(N_{q}^{\prime})*(h)$. There should be better choices of $\alpha$ than the special one we used.

Example 3. In this exampe, we consider the reconstruction from phaseless data. The numerical results are shown in Fig.6 which indicate the effectiveness of our method for the reconstruction from phaseless data.

Example 4. We use the measurements generated by the plane pressure incident wave with one fixed direction $d=(1,0)^{\top}$ (i.e. $M=1$). In this case the number of iterations at each frequency is set as $L=50$. The reconstruction results and relative errors are shown in Fig.7.

Example 5. In this example, we consider the special case discussed in section 4. The exact value of $q_{\rho}$ is given by

[TABLE]

see Fig. 8. Choose $\alpha=0.01$ and $\omega_{max}=11$. The reconstructed $q_{\rho}$ and relative errors from multi-frequency measurements with plane pressure incident wave are shown in Fig.9 and Fig.10.

Example 6. Finally, as a comparison, we consider the reconstruction of mass density only from data at a fixed frequency. For small frequency, it can be seen from Fig.11 that the reconstruction results are extremely bad no matter we have only one or multiple directions of incident wave. But for high frequency, see Fig.12(a), we still can have good reconstruction if we have multiple incident waves. Once we only have one fixed incident wave with direction $d=(0,1)^{\top}$, we can not obtain good reconstruction result by increasing $L$, see Fig.12(b,c). However, by increasing the number of frequency ($N=11$), we still can reconstruct some information of $q_{\rho}$, see Fig.13 and the results in Example 4. This further indicate the advantages to take multi-frequency data.

We conclude from the above numerical tests that satisfactory reconstructions are obtained through the proposed Landweber iterative algorithms.

Acknowledgments

The work of G. Bao is supported in part by an NSFC Innovative Group Fund (No.11621101), an Integrated Project of the Major Research Plan of NSFC (No. 91630309), and an NSFC A3 Project (No. 11421110002). The work of F. Zeng is supported by the NSFC grant (No. 11501063, No. 11771068), the Chongqing Research Program of Basic Research and Frontier Technology (No. CSTC2017JCYJAX0294) and the Fundamental Research Funds for the Central University (No. 106112016CDJXY100004). The authors also would like to thank Prof. Peijun Li for his suggestions on this work.