Application of neural networks to inverse elastic scattering problems with near-field measurements

1.   Introduction

Direct and inverse time-harmonic elastic wave scattering by a bounded body ^[1] has a wide range of applications, such as radar/sonar, geophysical exploration, medical imaging, seismology and materials sciences ^[2,3,4,5]. The forward problem is to determine and approximate the diffracted field given an incident wave and the elastic body. We refer the reader to, e.g., ^[6,7] for the existence and uniqueness of forward solutions via the variational approach. The inverse scattering problem (see e.g., ^[8,9,10]) is to determine the physical or geometric information of the scatterer based on knowledge of the scattered field. This paper is concerned with both direct and inverse problems.

In recent years, numerical approaches for position and shape reconstruction in many elastic scattering problems have drawn the attention of a large number of scholars. In ^[11], Arens proved that in two dimensions the shape of an elastically rigid body can be determined from the information of compressional and shear far-field patterns for all incident elastic plane waves with distinct directions at a fixed frequency. In ^[7], Bao et al. investigated the direct and inverse problems for a penetrable anisotropic elastic body embedded in a homogeneous and isotropic background medium. They show the uniqueness and existence of weak solutions by using the Fredholm alternative together with properties of the Dirichlet-to-Neumann map in both two and three dimensions. Moreover, an iterative method was designed to recover the interface from the multi-frequency near-field data. Charalambopoulos et al. ^[12] evaluated the inverse elastic scattering problem by using penetrable isotropic bodies embedded in an isotropic host environment via the factorization method. Hu et al. ^[13] investigated some uniqueness results for a ball or a convex Lipschitz polyhedron and adapted the factorization method to recover the shape of a rigid body from the scattered S-waves or P-waves corresponding to all incident plane shear or pressure waves. Recently, a reverse time migration method was proposed by Chen and Huang ^[14]. In ^[15], Li et al. presented an iterative method to solve the inverse elastic scattering problem. The domain derivative expresses the rate of change with respect to the surface displacement. A continuation method using multiple frequency data was established for the inverse problem.

The method of fundamental solutions (MFS) is an effective method for solving both direct and inverse problems governed by partial differential equations (PDEs). Since it is very simple and easy to program, it has been attracting an increasing amount of attention in recent years. Various engineering problems, including the acoustics scattering problems ^[16,17] and the acoustics-elastic scattering problem ^[18], can be solved by the MFS. The MFS can also be used to solve the inverse acoustic scattering problems (see e.g., ^[19,20,21]). Alves et al. ^[22] used the MFS to solve non-homogeneous elastic wave equations.

In recent years, neural networks, one of the most popular schemes in machine learning, have been applied in many research studies, including inverse scattering problems ^[23,24]. The seq2seq neural network has been applied for natural language processing and machine translation in the literature ^[25,26]. In ^[27], a new SwitchNet with a sparse neural network structure was proposed to solve the wave equation-based backscattering problems. Neural networks are especially suitable for dealing with inaccurate and vague information processing problems involving many parameters and constraints, and they can fully approximate complex nonlinear maps. Yin et al. ^[23] proposed a two-layer neural network method to recover the location and shape of an obstacle by using phaseless far-field data.

In this paper, we consider inverse elastic scattering problems related to recovering the shape of a rigid body and a cavity by applying a two-layer gated controlled recurrent unit (GRU) neural network. The convergence of the presented method will be verified. The contributions of this paper are summarized as follows.

(1) We propose a fast numerical scheme for solving the direct problem;

(2) We prove the convergence of a two-layer GRU neural network;

(3) We test the accuracy of the reconstruction method by numerical examples.

The rest of the paper is arranged as follows. In Section 2, we formulate the elastic scattering problem and describe the numerical scheme for forward solutions. In Section 3 we introduce the architecture of the neural network. In Section 4, we perform theoretical analysis for the convergence of the network. In Section 5, numerical experiments are presented to verify the efficiency and effectiveness of the neural network.

2.   Formulation and solution method

Let $D \subset \Omega$ be a two-dimensional elastic body with the smooth boundary $\partial D$, where $\Omega\subset \mathbb{R}^2$ is a subdomain. The exterior domain ${D^c:} = \mathbb{R}^{2}\setminus{\overline{D}}$ is supposed to be connected and filled with a homogeneous and isotropic elastic medium. The displacement ${{\mathit{\boldsymbol{u}}}}(\textbf{x}, t)$ satisfies

In general, the displacement ${{\mathit{\boldsymbol{u}}}}(\textbf{x}, t)$ is time harmonic and given by ${{\mathit{\boldsymbol{u}}}}(\textbf{x}, t) = e^{{\rm i}\omega t}\mathit{\boldsymbol{v}}(\textbf{x})$ (where ${\rm i}$ is the imaginary unit); thus, the displacement $\mathit{\boldsymbol{v}}(\textbf{x})$ satisfies the Navier equation

where $\omega > 0$ is the angular frequency, $\lambda, \mu$ are the Lamé constants satisfying that $\mu > 0$ and $\lambda +\mu > 0$, and $\rho$ is the density function, which is assumed to be constant.

The forward problem considered in this paper is to give an incident elastic plane wave ${{\mathit{\boldsymbol{u}}}^{incident}} = {\mathit{\boldsymbol{d}}}e^{{{\rm i}}k_{p}{\mathit{\boldsymbol{d}}}\cdot \textbf{x}}+ {\mathit{\boldsymbol{d}}}^{\bot}e^{{{\rm i}}\kappa_{s}{\mathit{\boldsymbol{d}}}\cdot \textbf{x}}$ with ${\mathit{\boldsymbol{d}}} = (\cos\vartheta, \sin\vartheta)$ and ${\mathit{\boldsymbol{d}}}^{\bot} = {(-\sin\vartheta, \cos\vartheta)}, \vartheta\in [0, 2\pi)$, as well as to find the scattered field ${\mathit{\boldsymbol{v}}} = (v_1, v_2)^\top$ that satisfies the Navier equation.

By the well-known Hodge decomposition, any solution $\mathit{\boldsymbol{v}}$ to (2.1) can be decomposed into

where $\mathit{\boldsymbol{v}}_p: = (-1/{k_p^2}){\rm grad}\; {\rm div} \mathit{\boldsymbol{v}}$ is known as the compressional part of $\mathit{\boldsymbol{v}}$ and $\mathit{\boldsymbol{v}}_s: = (-1/{k_s^2}){\rm grad}^\bot\; {\rm div}^\bot \mathit{\boldsymbol{v}}$, with ${\rm grad}^\bot = \left(-{\partial}_2, {\partial}_1\right)^\bot$ and ${\rm div}^\bot{\mathit{\boldsymbol{v}}} = -{\partial}_2v_1+{\partial}_1v_2$ as the shear part. Here, $k_s = \omega\sqrt{\frac{\rho}{\mu}}$ is the shear wave number and $k_p = \omega\sqrt{\frac{\rho}{\lambda+2\mu}}$ is the compressional wave number.

For a rigid body $D$, the scattered field $\mathit{\boldsymbol{v}}$ satisfies the first (Dirichlet) kind of boundary condition, i.e., it holds that

so that the total displacement is zero.

If $D$ is a cavity, the following second (Neumann) kind of boundary condition holds:

with ${\mathit{\boldsymbol{t}}}: = {\mathit{\boldsymbol{t}}}_n{{\mathit{\boldsymbol{u}}}^{incident}}$. Here, the traction operator ${\mathit{\boldsymbol{t}}}_n$ is defined by

where $\hat{\mathit{\boldsymbol{e}}}_i$, $i = 1, 2,$ designates the unit vector along the $x_i$ axis and ${\mathit{\boldsymbol{n}}} = (n_1, n_2)$ denotes the unit normal direction point into $D^c$.

To ensure uniqueness, we require the scattered field $\mathit{\boldsymbol{v}}$ to satisfy the two-dimensional Kupradze-Sommerfeld radiation condition at infinity, that is,

In this paper, we will employ the MFS (see ^[28] for the acoustic scattering problem) to calculate the forward solution $\mathit{\boldsymbol{v}}$. The main idea of the MFS is to approximate the scattered field in the solution domain $D^c$ by using a linear combination of fundamental solutions ^[18,29] to the Navier equation (2.1) with $N$ source points $\textbf{y}_j\in D$, i.e.,

Here, $a_j, b_j$ are the coefficients to be determined by the boundary condition on $\partial D$ and $\mathit{\boldsymbol{U}}(\textbf{x}, \textbf{y})$ is Green's tensor of the two-dimensional Lam$\acute{e}$ system of the form

where ${{{\mathit{\boldsymbol{I}}}}}$ is the 2-by-2 identity matrix and

is the fundamental solution of the Helmholtz equation associated with the wavenumber $k$. Here, $H_0^{(1)}$ is the Hankel function of the first kind of order zero.

Remark 2.1. Following the idea in ^[29], the solution ${\mathit{\boldsymbol{v}}}(\textbf x)$ of this problem can be approximated by the single-layer potential function, with $\boldsymbol{\psi}\in (L^2(\partial B))^2$, as follows:

where $\overline{B}\subset D$. If we discretize the integral operator $S$, we can get a linear combination of the fundamental solutions with respect to some source points located on $\partial B$ with the unknown density function ${\boldsymbol{\psi}}(\textbf y_j)$. It can be seen that (2.5) is a certain discrete form of (2.6) with some weight coefficients.

We respectively define the compressional and shear parts of $\mathit{\boldsymbol{U}}(\textbf{x}, \textbf{y})$ as follows:

and

Then $\mathit{\boldsymbol{U}}(\textbf{x}, \textbf{y}) = \mathit{\boldsymbol{U}}_p(\textbf{x}, \textbf{y})+\mathit{\boldsymbol{U}}_s(\textbf{x}, \textbf{y})$. Thus we have

and

Recalling the definitions of the traction vector (2.4) and the fundamental solution to the Navier equation, one can approximate the traction vector on the boundary $\partial D$ as follows:

Here, the matrix $\textbf{T}(\textbf{x}, \textbf{y}_j)$ is defined by

with the entries

To determine the coefficients $a_j, b_j, j = 1, 2, ...N,$ we employ the collocation method by choosing $M$ collocation points $\left\{\textbf x_j\right\}^M_{j = 1}\subset \partial D$ and impose the boundary condition on $\partial D$. Set ${\textbf{d}'}: = (a_1, \cdot\cdot\cdot, a_N, b_1, \cdot\cdot\cdot, b_N)^\top$ to be the vector of $2N$ unknown parameters. It can be calculated by solving a linear system with $2M$ equations of the form

where the matrix $\mathcal {A}\in \mathbb{C}^{2M \times 2N}$ on the left-hand side and the vector $\mathit{\boldsymbol{h'}}\in \mathbb{C}^{2M\times 1}$ on the right-hand side are defined as follows depending on the boundary condition on $\partial D$.

● Dirichlet boundary condition (rigid body)

● Neumann boundary condition (cavity)

In our numerics, we choose the number $M = N$ in all examples. The algebraic linear system (2.10) is solved directly by using Matlab with the command ${\textbf{d}'} = \mathcal {A} \backslash\mathit{\boldsymbol{h'}}$. Then we can get the scattered near-field data $\mathit{\boldsymbol{v}}$ on $\partial \Omega$ from (2.5). One can also get the shear part $\mathit{\boldsymbol{v}}_{s}|_{\partial \Omega}$ and compressional part $\mathit{\boldsymbol{v}}_{p}|_{\partial \Omega}$ of $\mathit{\boldsymbol{v}}$ by replacing Green's tensor ${\mathit{\boldsymbol{U}}}$ by its shear and compressional parts in (2.7) and (2.8), respectively. These synthetic data will be used as the measurement data in our inverse problem. In our examples below, we suppose, for simplicity that the boundary $\partial \Omega$ can be parameterized as follows:

This means that the measurement data are taken on the circle centered at $(a, b)\in \mathbb{R}^2$ with the radius $R > 0$.

Definition 2.1. Let $\hat{\mathit{\boldsymbol{d}}}_{m}: = R(\cos\theta_m, \sin\theta_m)+(a, b), m = 1, 2, \cdots, \mathcal{N}_{1};$ indicate the measurement positions where $\mathcal{N}_{1}\in \mathbb{N}$ and ${\theta_m\in(0, 2\pi]}$. Given the obstacle $D$ and the compressional and shear wave-numbers $k_{p}, k_{s}$, respectively, we define the measurement set $\mathit{\boldsymbol{I}} = (I(1), I(2), \cdots, I(T))$ with $T = 2\mathcal{N}_{1}$ as follows:

The inverse problem in this paper is to recover the shape of the obstacle or the cavity, namely $\partial D$, by applying knowledge of the near-field data on $\partial\Omega$ associated with a fixed incident field:

Assumption 2.1. The boundary of the obstacle $D$ is assumed to take the form

where $p(\theta)$ permits the following (truncated) Fourier representation:

where $Q\in \mathbb{N}$. Below we define the set of the Fourier coefficients $a_0, a_{m}$ and $b_{m}{m = 1}, 2, \cdots, Q,$ by $\mathit{\boldsymbol{y}} = \{y(1), y(2), \cdots, y(N)\}$ with $N = 2Q+1.$ In the remaining part of this paper, we shall consider the inverse problem of recovering $\mathit{\boldsymbol{y}}$ from our measurement data set.

3.   Architecture of neural networks

In this section, we introduce the main components of the neural network first. Then how to construct the neural network to deal with the inverse elastic scattering problem is given.

Figure 1 gives a simple map of the exhibition of its working status. We introduce the GRU cells (see e.g., ^[23,30]), i.e., the principal element of the update network. This network consists of many GRU neurons.

$(1)$ The cell receives a stimulus, i.e., an input $I(t)(t = 1, 2, \cdots, T)$. Then the cell will respond to the stimulus. We describe this response as an activation function acting on the previous output $h(t-1)$ and the present stimulus $I(t)$ to produce a reset gate ${r(t)}$, which is to control candidate status $\tilde{h}(t)$ based on whether it depends on the prior status of learning.

$(2)$ It preserves historical output $h(t-1)$ and presents input $I(t)$ to obtain update gate $z(t)$, which can be used to choose and forget memory. In other words, it can forget some information passed down in $h(t-1)$ and select and memorize the information in $\tilde{h}(t)$.

In order to see the notations clearly, we list the notations used in the neural network in Table 1.

3.1.   Construction of the network

The whole model consists of three parts, namely, the encoding end, the middle status and the decoding end. The role of the encoding end is to extract data features, and the intermediate state serves to connect the encoding and decoding ends. The function of the decoder is to decode the information obtained by the encoder to obtain the output sequence. Next, the relevant details of the network structure are illustrated as below.

Step 1: The scattered data features are extracted by using GRU neurons on the encoding end, i.e.,

where $\zeta_{t}$ signifies the t$^{th}$ GRU cells in the course of the work. That is, it is equivalent to extracting the input data features by using (3.1)–(3.5). Then the first layer extracts the scattered data as $\mathit{\boldsymbol{h}}^{1} = (h^{1}(1), h^{1}(2), \cdots, h^{1}(T))$.

Step 2: Extracting features again.

Take the features extracted for the first time, i.e., $\mathit{\boldsymbol{h}}^{1} = (h^{1}(1), h^{1}(2), \cdots, h^{1}(T))$ as input to extract the features again by means of (3.1)–(3.5):

Step 3: First decoding process on the decoding terminal.

The decoder side will use $N$ GRU neurons to decode the second layer output information $h^2(T)$ and the obstacle boundary Fourier expansion coefficients ${\mathit{\boldsymbol{y}}} = (y(1), y(2), \cdots, y(N))$:

where $\zeta_{n}$ signifies operation of the $n^{th}$ GRU. It plays the same role as $\zeta_{t}$ in Step $1$.

Step 4: Second decoding process at the decoder end.

Get the output $\mathit{\boldsymbol{h}}^{4}$ for the fourth layer by decoding the output $\mathit{\boldsymbol{h}}^{3} = (h^{3}(1), h^{3}(2), \cdots, h^{3}(N))$ of the third layer. That is, this layer will preserve profitable information and discard unprofitable information

Finally, the latest output from the decoding side will be converted into a predictive output by the fully connected layer. That is to say, we will calculate $\hat{\mathit{\boldsymbol{y}}} = (\hat{y}(1), \hat{y}(2), \cdots, \hat{y}(N))$, i.e., the predicted output of the network by using the following formula

4.   Theoretical analysis of convergence

This section is devoted to the convergence analysis of the neural network.

Definition 4.1. Let $\{y(n)\}^{N}_{n = 1}$ be the initial values of the boundary curve and $\{\hat{y}(n)\}^{N}_{n = 1}$ be the network prediction output. Define the common mean absolute value loss (error) function as

When $y(n)-\hat{y}(n)\neq0$, the derivative of the loss function (4.1) with respect to the weight $\mathit{\boldsymbol{w}}$ in (3.1)–(3.5) is given by

For the sake of convenience, we introduce some symbols. The input weights and the internal weights of every layer are gathered in the weight matrices $\mathit{\boldsymbol{W}}_{I}$ and $\mathit{\boldsymbol{W}}_{L}$. We denote $\mathit{\boldsymbol{W}} = (\mathit{\boldsymbol{W}}_{I}, \mathit{\boldsymbol{W}}_{L})$. The specific forms are as follows

Here, $\mathit{\boldsymbol{W}}_{I}$ has the same form as $\mathit{\boldsymbol{W}}_{L}$.

Next, we write the input of every layer as

The input vectors for each layer of the encoding and decoding ends are given as follows:

Definition 4.2. Given a vector $\mathit{\boldsymbol{\alpha}} = (\alpha_{1}, \alpha_{2}, \cdots, \alpha_{n})$, the activation function acts on a vector as follows:

The derivatives of the activation function acts similarly

In other words, the activation function acts on a vector componentwise.

The weight gradient of the fourth-layer neural network update gate can be calculated directly:

From the initial condition constant $\hat{y}(0)\equiv 0$, we have

For any given initial weight $\mathit{\boldsymbol{w}}^{0}$, the weight sequence $\{ \mathit{\boldsymbol{w}}^{k} \}$ can be generated from an iterative formula:

where

and $\gamma$ is the learning rate.

Assumption 4.1. (^[31,32]) $|\sigma(a)|, |\sigma'(a)|, |\sigma''(a)|$ and $|\tanh(a)|, |\tanh'(a)|, |\tanh''(a)|$ are uniformly bounded for $a\in \mathbb{R}$.

Assumption 4.2. (^[31,32]) $\|\mathit{\boldsymbol{W}}^{k}_{L}\|_{(k = 1, 2, \cdots)}$ is bounded in the learning course described by (4.5).

In order to check the convergence of the neural network, we claim the following theorem.

Theorem 4.1. The loss function is defined as the formula (4.1), and the weight sequence $\{\mathit{\boldsymbol{w}}^{k} \}$ can be generated from the iterative formula (4.5). If Assumptions 4.1 and 4.2 are valid, then we have the following:

(1) $L(\mathit{\boldsymbol{w}}^{k+1}) \leq L(\mathit{\boldsymbol{w}}^{k})$;

(2) $\exists \ L^{*} \geq 0$ s.t. $\lim_{k\rightarrow \infty} L(\mathit{\boldsymbol{w}}^{k}) = L^{*}$.

Proof. First, we can also use the method employed in ^[32], which adopts the induction argument to calculate $\Delta \mathit{\boldsymbol{v}}^{4}(\mathit{\boldsymbol{w}}^{k}, n)$ as follows:

where $\delta (\mathit{\boldsymbol{w}}^{k}, n)$ will be given in the following remark. Then, apply the Taylor expansion to $\mathit{\boldsymbol{v}}^{4}(\mathit{\boldsymbol{w}}^{k}, n)$; we have

Next, substituting (4.7) into the first order term of $\Delta \mathit{\boldsymbol{v}}^{4}(\mathit{\boldsymbol{w}}^{k}, n)$ of the above equation, we have

Then, combining (4.2), (4.3) and (4.6), we will get

By Assumptions 4.1 and 4.2, together with the Cauchy-Schwarz inequality, we obtain that $|\varphi (\mathit{\boldsymbol{w}}^{k})| \leq C\| \Delta \mathit{\boldsymbol{w}}^{k}\|^{2}.$ Hence $L(\mathit{\boldsymbol{w}}^{k+1})-L(\mathit{\boldsymbol{w}}^{k}) \leq (- \frac{1}{\gamma}+C)\| \Delta \mathit{\boldsymbol{w}}^{k}\|^{2}.$ If the learning rate is sufficiently small and satisfies that $0 < \gamma < 1/C$ ($C$ is a constant), then the first conclusion of Theorem 4.1 is established.

Since $L(\mathit{\boldsymbol{w}}^{k}) \geq 0$ is monotonically bounded, there exists an $L^{*}$ such that $\lim_{k\rightarrow \infty}L(\mathit{\boldsymbol{w}}^{k}) = L^{*}.$ The proof is completed. □

Remark 4.1. In (4.7), we have

with

when $d = 1$. Here, $\tau (n-d)$ is the intermediate quantity for $\mathit{\boldsymbol{v}}(\mathit{\boldsymbol{w}}^{k+1}, n)$ and $\mathit{\boldsymbol{v}}(\mathit{\boldsymbol{w}}^{k}, n).$ From (4.3) and (4.4), we have

When $d = 0$, $\prod\limits_{l = n-d+1}^{n}\mathit{\boldsymbol{W}}^{4}_{L}{((1-z_{4}'(l))} \odot \ h^{4}(l-1) \oplus \ z_{4}'(l) \odot \tilde {h}(l))\sigma'(v^{4}(l)) = 1.$

5.   Numerical examples

In this section, we will give some numerical examples to show the effectiveness of the presented method.

The scattered data described by ${\mathit{\boldsymbol{I}}}_{D} = (I(1), I(2), \cdots, I(T))$ corresponding to the obstacle $D$ are given by Definition 2.1. The shape parameters ${\mathit{\boldsymbol{y}}}_{D} = (y(1), y(2), \cdots, y(N))$ are defined as in Assumption 2.1. We denote the training data set by

where $\mathcal{D}$ denotes the set of shape parameters. For the elements in $T'$, ${\mathit{\boldsymbol{I}}}_{D}$ is also updated as the shape $D$ changes, and we randomly selected $2\%$ of the scattered data as the testing set to get the final prediction results.

The inverse problem can be solved in two steps. First, for a given shape $D$, the elastic wave equation (2.1) is solved by using the MFS described in Section 2. Second, the network was trained to identify the obstacle from the scattered data, i.e., the network was trained to approximate the Fourier coefficients of the obstacle.

To test the sensitivity to the noise, we suppose that noisy data is generated by adding Gaussian random noise pointwise, i.e.,

where $\delta\geq 0$ denotes the noise level and the random variable $\xi(i)$ is the standard Gaussian distribution.

In all numerical examples, the parameters are given as follows.

$(1)$ We shall use the Lamé constants $\lambda = 2$ and $\mu = 1$ and the frequency $\omega = 3$, $\rho = 1.$

$(2)$ The symbol "-" signifies the exact boundary of an obstacle, and "--" signifies the forecast boundary of an obstacle.

$(3)$ The number of neurons of the GRU neural network is generally set as $2^{p}(p \in \mathbb{R}^{+}),$ and the network parameters were determined through a large number of numerical experiments. The number of GRU neurons was set as 128. Dropout was set as 0.5. Batch processing was set as 500. And the number of iterations is 40.

Remark 5.1. The scattering data are complex numbers and the current neural network method can only deal with real numbers; thus, we took out the real and imaginary parts of the scattering data as the scattering data information in the training data set. In addition, we need to use a random function to generate random numbers as the shape parameters in the training data set. This is the preparation of the training set. After the training set is prepared, the network model is trained by using the data in the training set through the two-layer network which has been built. Finally the output of the network, i.e., the predicted values of the shape parameters, is obtained by using the trained model to use the test data set.

Example 5.1 (Dirichlet boundary condition). First, the scattered data used in the training network are noise free, and the scattered data are collected on a circle with the radius $R = 3$ centered at the origin. We chose to investigate the numerical recovery results by changing the number of measurements $T$. The ${\mathit{\boldsymbol{d}}} = (1, 0)$ is the incident direction. The boundaries of the scatterers used in this example are parameterized as follows.

The first case is a heart-shaped domain parameterized by

The second case is a round-square domain parameterized by

The polar diameter has the following (truncated) Fourier representation: $\frac{a_{0}}{2}+\sum_{m = 1}^{Q}a_{m}\cos(m\theta).$ We take $Q = 6$, and the idea is to reconstruct the Fourier coefficients $a_m$. The coefficients used in the training process, i.e., $a_0\in[0.5, 4]$, $a_1, a_2\in[-0.5, 0.5]$ and $a_m\in[-0.05, 0.05], m = 3, 4, 5, 6,$ are uniformly distributed random numbers.

Figure 2 presents the recovery results for a fixed number $T\in\{4000$, 20,000, 40,000$\}$. From the numerical recovery results in Figure 2, we can see that the accuracy in the numerical results will be improved by increasing the number $T$.

Second, in addition to the usual rounding errors, we also contaminated the forward scattered data by adding random noise. We chose to investigate the stability of the inversion schemes by changing the noise level $\delta$. Figure 3 presents the recovery results for the noise level $\delta\in\{5\%, 10\%, 15\%\}$. From the fourth column in Figure 3, we can see that the errors are becoming smaller as the noise level decreases. We can also see that the numerical reconstruction results are poorest where the body has high curvature. This meets our expectations, since more expansion coefficients are required at the corner.

Example 5.2 (Neumann boundary condition). For comparison, in this example we consider the scattering problem by cavities, i.e., the Neumann boundary condition is imposed. In this example, we will consider a heart-shaped domain parameterized by

The training sets are the same as in Example 5.1.

We chose the direction ${\mathit{\boldsymbol{d}}} = (1, 0)$. The scattered data were collected on a circle with the radius $R = 3$ centered at the origin. The numerical recovery results for the noise level $\delta\in\{0\%, 5\%, 10\%, 15\%\}$ with a fixed number $T = 40,000$ are shown in Figure 4. From this, we observe that the locations and shapes are still well captured. From the fourth row in Figure 4, we can see that the errors are becoming smaller as the noise level decreases. We can also see that the numerical reconstruction results are poorest where the body has high curvature. This further shows that our methods are independent of the physical properties of the underlying scatterers.

Example 5.3 (Neumann boundary condition). In this example, we will consider the case of an apple-shaped domain parameterized by

The polar diameter has the following (truncated) Fourier representation $\frac{a_{0}}{2}+\sum_{m = 1}^{Q}b_{m}\sin(m\theta).$ We take $Q = 6$, and the idea is to reconstruct the Fourier coefficients $a_0$ and $b_m$. The coefficients used in the training process, i.e., $a_0\in[0.5, 3]$, $b_1, b_2, b_3, \in[-0.5, 0.5]$ and $b_4, b_5, b_6, \in[-0.1, 0.1]$, are uniformly distributed random numbers.

We chose the direction ${\mathit{\boldsymbol{d}}} = (1, 0)$. The scattered data were collected on a circle with the radius $R = 3$ centered at the origin. The numerical recovery results for the noise level $\delta\in\{5\%, 10\%, 15\%\}$ with a fixed number $T = 40,000$ are shown in Figure 5. The error for all testing sets obtained with different noise levels $\delta\in\{5\%, 10\%, 15\%\}$ and $T = 40,000$ is demonstrated by the fourth column in Figure 5. From this, we observe that the locations and shapes are still well captured.

The above examples give the effectiveness of the presented method for the case of a single obstacle or cavity. The following example will give the case of two disjoint components. Meanwhile, we should know that $D$ is the case of two disjoint components as the a priori information. The Fourier coefficients were tracked as two independent groups.

Example 5.4 (Mixed boundary conditions). As shown in Figure 6, we consider a scatterer with two disjoint components. We set the scatterer to be a combination of a triangle-shaped domain and a peanut-shaped domain. Suppose that $D_1$ is a big triangle-shaped domain of the Dirichlet kind, parameterized by

$D_2$ is a mini peanut-shaped cavity. The boundary $\partial D_2$ is parameterized by

We chose the direction ${\mathit{\boldsymbol{d}}} = (\frac{\sqrt{3}}{2}, \frac{1}{2})$ in this example. The scattered data were collected on a circle with the radius $R = 3$ centered at the origin. Figure 6 shows the numerical results obtained from the knowledge of the compressional near-field measurements, the shear near-field measurements and the full near-field measurements with $\delta = 10\%$. From these subfigures, we observe that the locations and shapes can be captured. This further shows that the presented method is independent of the physical properties of the underlying scatterers, even for two disconnected obstacles with different boundary conditions. We also note that the numerical error is worse than that in the above examples.

Example 5.5 (Two obstacles). As shown in Figures 7 and 8, the underling scatterer $D$ is given as the union of two disjoint obstacles $D = D_1\cup D_2$. We consider two cases in this example.

● Dirichlet boundary condition: We chose the direction ${\mathit{\boldsymbol{d}}} = (1, 0)$. $D_1$ is a peanut-shaped domain, and the boundary $\partial D_1$ is parameterized by

$D_2$ is a heart-shaped domain, and the boundary $\partial D_2$ is parameterized by

● Neumann boundary condition: We chose the direction ${\mathit{\boldsymbol{d}}} = (\frac{\sqrt{2}}{2}, \frac{\sqrt{2}}{2})$. $D_1$ is a round-square domain, and the boundary $\partial D_1$ is parameterized by

$D_2$ is a kite-shaped domain, and the boundary $\partial D_2$ is parameterized by

The numerical recovery results for the noise level $\delta\in\{0\%, 10\%\}$ with a fixed number $T = 40,000$ are shown in Figures 7 and 8. From this, we observe that the locations and shapes are still well captured.

Example 5.6 (Dirichlet boundary condition). The last example considered here is a different option of the physical parameters. The Lamé constants are given as $\lambda = 3.88, \mu = 2.56$ and the frequency $\omega = \pi$. The scattered data were collected on a circle with the radius $R = 3$ centered at the origin. The incident direction is ${\mathit{\boldsymbol{d}}} = (-\frac{1}{2}, \frac{3}{2})$. The boundary of the scatterer is a heart-shaped domain parameterized by

Figure 9 presents the recovery results for a fixed number $T = 40,000$ with the noise level $\delta\in\{5\%, 10\%, 15\%\}$. It can be seen that the recovery results are independent of the physical parameters.

6.   Conclusions

In this paper, we applied a neural network trained on solutions from the MFS to solve inverse elastic wave scattering problems. Both convex and non-convex geometries with smooth boundaries have been investigated. The numerical results confirm the effectiveness of the presented method. From the presented examples, we can see that the presented method is independent of the physical properties of the underlying scatterers. When the boundary consists of two parts, this method also gives a good approximation of the exact boundary curves.

Use of AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgements

We would like to thank the reviewers for their valuable comments that helped to improve the quality of the submitted manuscript. We also thank professor G. Hu for the valuable comments on this paper. The research was supported by the Natural Science Foundation of China (Nos. 12101603, 11501566), the Fundamental Research Funds for the Central Universities (No. 3122019159) and Tianjin Education Commission Research Project (No. 2022KJ072).

Conflict of interest

The authors declare that there are no conflicts of interest.