Reservoir prescriptive management combining electric resistivity tomography and machine learning

1.   Introduction

Inverse problems for wave equations arise widely in non-destructive testing ^[1,2], underwater acoustics ^[3,4], antenna synthesis ^[5,6] and other application areas. Being two prototypical models in the field of inverse problems, the time-harmonic inverse obstacle scattering and inverse source problem have received significant investigations in the literature. For the nonlinear inverse obstacle scattering problem, the emanating source (incident wave) is given and the inhomogeneous scatterer serves as the target. Reciprocally, in the linear inverse source problem, the source is unknown while the passive scatterering inhomogeneity is vacant. Classical algorithms for inverse obstacle scattering problems involve the optimization methods ^[7,8,9], the singular source method ^[10] and the recursive linearization algorithm ^[11]. In addition, there are some sampling-type methods with qualitative features, for instance, the factorization method ^[12], the linear and direct sampling methods ^[13,14], the single-shot schemes ^[15,16] and the reverse time migration approach ^[17,18]. Recently, a variety of reconstruction methods have also been developed to deal with the inverse source problems. Bao et al. ^[19] applied a recursive algorithm to recover unknown sources of acoustic field with multi-frequency measurement data. Zhang et al. ^[20] used the direct sampling method to locate multipolar acoustic sources for the Helmholtz equation from near-field Cauchy data. Bousba et al. ^[21] utilized the direct sampling method to identify point sources for the Helmholtz equation from far-field data. The algebraic method has also been applied to the Helmholtz equation for determining monopolar sources ^[22].

In this paper, we are concerned with an inverse problem to simultaneously reconstruct the acoustic obstacle and its excitation source points from the near-field measurements. This co-inversion problem can be mathematically formulated as an intrinsically coupled problem comprising the aforementioned inverse obstacle scattering and inverse source problem. In comparison with the vast studies on the conventional inverse obstacle/source problems, focuses on identification of both unknown source and scatterer are relatively rare. Over recent years, some interesting attempts have been made in dealing with such co-inversion problems. In ^[23], Liu et al. established sufficient conditions to recover both the sound speed of the medium and the source in thermo- and photo-acoustic tomography. The simultaneous recovery of an embedded obstacle and its surrounding medium by formally-determined scattering data was studied in ^[24]. In ^[25], the authors proved the uniqueness in simultaneously recovering the source functions and the obstacle for the wave equation. An adjoint state method for recovering both the source and attenuation in the attenuated X-ray transform was proposed in ^[26]. In addition, we also refer to ^[27,28,29] for some recent studies on the determination of both potential and source in the Schrödinger equations.

In contrast to the standard problem of reconstructing either the source or obstacle solely, it is substantially more challenging to identify these unknown targets concurrently. It is worthwhile noting that only the total field is available and the contributions from inaccessible source and scatterer are superimposed in the total field measurements in the scenario of co-reconstruction. Hence, decoupling these incident and scattered components would inevitably suffer from inherent difficulties. As a result, the usual nonlinearity, ill-posedness as well as the additional lack of information constitute the mathematical and numerical challenges of our current study.

In this work, we propose an optimization method to recover an obstacle and its illuminating point sources simultaneously from the measured total field. The salient features of our method can be summarized as follows. First, the proposed method is capable of simultaneously and quantitatively determining both the obstacle and source points since the reconstruction does not rely on any alternating iteration between the source and scatterer. Second, by the idea of the decomposition method and reformulating the inverse problem as a nonlinear optimization problem, the solution process for the direct problem is unnecessary. Third, by designing proper indicator functions, we provide an easy-to-implement and effective strategy for choosing the initial guesses for the optimization method, which avoids the a prior information of the unknown targets. Finally, the proposed method is insensitive to the noise-contaminated data and works well even with limited-aperture measurement data despite the severe ill-posedness.

We now introduce the precise formulation of the co-inversion problem under consideration. Let $D\subset \mathbb{R}^2$ be an open and simply connected domain with $C^2$ boundary $\partial D.$ Given the source point located at $z\in\mathbb{R}^2\backslash\overline{D},$ the incident field $u^i$ due to $z$ is given by

where $k > 0$ is the wave number, $H_0^{(1)}$ is the Hankel function of the first kind of order zero. Then the forward scattering problem for the sound-soft obstacle $D$ is to find the scattered field $u^s\in H_{\rm loc}^1(\mathbb{R}^2\backslash\overline{D})$ such that

where $u = u^i+u^s$ is the total field. The Sommerfeld radiation condition (1.4) holds uniformly with respect to all directions $\hat{x} = x/|x|\in\mathbb{S}: = \{x\in\mathbb{R}^2:|x| = 1\}$, and (1.1) is also known as the fundamental solution to the Helmholtz equation (1.2). It is well known (see, e.g., ^[30]) that the direct problem (1.2)–(1.4) admits a unique solution $u^s\in H^1_{\rm loc}\left(\mathbb{R}^2\backslash\overline{D}\right)$.

In this paper, we assume that $D\subset B_R: = \{x\in\mathbb{R}^2:|x| < R\}$ and $B_R\backslash\overline{D}$ is connected. Let $S: = \cup_{j = 1}^N\{z_j\}\subset B_R\backslash\overline{D}$ be a set of distinct source points with $N\in\mathbb{N}_+$ the number of source points. Denote by $u^s(x; z_j)$ the scattered field corresponding to the incident field $u^i(x; z_j)$ of the form (1.1), and we collect the total field $u(x; z_j) = u^i(x; z_j)+u^s(x; z_j)$ on the measurement curve $\Gamma_R: = \partial B_R = \{x\in\mathbb{R}^2: |x| = R\}$. Define the measured total field data by

then the inverse problem under consideration is to determine the obstacle-source pair $(\partial D, S)$ from the measurements $\mathbb{U}$, namely,

We refer to Figure 1 for an illustration of the geometry setting of the inverse obstacle and source problem (1.5). In addition, it deserves noting that, both the single source case ($N = 1$) and the multiple source case ($N\ge 2$) are admissible scenarios. The proposed algorithms in the subsequent sections can be applicable to both scenarios. We also want to point out that, the discrete total field measurement $\bf{U}$ is a complex multistatic response matrix whose $(i, j)$-th entry is the total field data due to the $j$-th source point and the $i$-th receiver. In other words, the total fields are collected separately for each source point one by one in our model and thus not directly related to the incident wave produced by the superposition of source points. Under the above configuration, the number $N$ of the point sources is known as a priori, which is exactly the number of columns of the measurement matrix $\bf{U}$.

The rest of this paper is organized as follows: In Section 2, an optimization method is proposed to reconstruct the source-obstacle pair and the corresponding convergence theory is established. In Section 3, an effective strategy based on the sampling-type method is provided to select the initial guess for the optimization method. Numerical experiments are given in Section 4 to verify the performance of the method. Finally, some conclusions are given in Section 5.

2.   Optimization methods

Based on the decomposition method for inverse obstacle scattering problem ^[9], we propose a novel optimization method in this section by incorporating the source locations as new variables in the cost functional in order to reconstruct the obstacle-source pair $(\partial D, S)$. The basic idea of the optimization method is to seek for the positions of source points and the boundary of obstacle as the locations where the Dirichlet boundary condition holds, namely, the total field $u^i+u^s$ vanishes. We would like to emphasize that, since the admissible source locations and obstacles are radically distinct, the new optimization problem will thereby encounter nontrivial difficulties in the analysis and implementation.

2.1.   Optimization method for the co-inversion problem

In this subsection, we present the optimization method for the inverse obstacle-source problem. According to the initial guess by the sampling results (see the next section for more details) about the obstacle-source pair $(\partial D, S),$ we can place a closed $C^2$ curve $\Lambda$ inside $D$ such that $k^2$ is not a Dirichlet eigenvalue for the negative Laplacian in the interior of $\Lambda$. Choosing some proper admissible set $U\times V$ of the obstacle-source pair $(\Gamma, z')$ such that $\Gamma\in U,$ $z' = (z'_1, z'_2, \cdots, z'_N)\in V,$ the scattered field $u^s(x; z'_j)$ corresponding to $z'_j$ can be represented by a single-layer potential

with the unknown density $\varphi_j\in L^2(\Lambda).$ Given the measured total field $u(x; z_j),$ the density $\varphi_j$ are supposed to satisfy the integral equation of the first kind

where $u^i(x; z'_j)$ is the incident wave due to source point $z'_j$ and the compact integral operator $\mathcal{S}:L^2(\Lambda)\to L^2(\Gamma_R)$ is defined by

The fact that the kernel of $\mathcal{S}$ is analytic leads to the severe ill-posedness of (2.2), which motivates us to seek for a stable numerical solution by adopting the Tikhonov regularization strategy to replace (2.2) by

where $\alpha > 0$ is a regularization parameter and the adjoint operator $\mathcal{S}^*:L^2(\Gamma_R)\to L^2(\Lambda)$ of $\mathcal{S}$ is given by

where the overline signifies the complex conjugate.

After obtaining the regularized density $\varphi_j^\alpha$ by solving (2.3), we can calculate the approximation $u^s_\alpha(x; z'_j)$ for the scattered field $u^s(x; z'_j)$ by inserting the solution $\varphi_j^\alpha$ of (2.3) into the single-layer potential representation (2.1). Then, the obstacle-source pair $(\partial D, z)$ can be sought for by minimizing the defect

over some suitable class $U\times V$ of admissible obstacle-source pairs $(\Gamma, z').$ In (2.4), $u^i(x; z'_j)$ is the incident field due to the source point $z'_j.$

Remark 2.1. We would like to point out that, as $\alpha$ tends to zero, the solution $\varphi_j^\alpha$ of (2.3) converges only if (2.2) is solvable. Following [7,Section 5.5] for the classic decomposition method, the solvability of (2.2) is related to regularity properties of the scattered field which, in general, cannot be known in advance for an unknown obstacle $D$. In other words, $u^s_\alpha(x; z'_j)$ may not converge to the $u(x; z_j)-u^i(x; z'_j)$ since (2.2) may not be solvable.

For the minimization problem, the admissible class $U$ of $\Gamma$ is chosen to be a compact set (with respect to the $C^{1, \beta}$ norm, $0 < \beta < 1$) of all star-like closed $C^2$ curves, described by

where $r\in C^2(\mathbb{S})$ is the radial function satisfying the prior information

with given constants $a$ and $b.$ The admissible class $V$ of $z' = (z'_1, z'_2, \cdots, z'_N)$ is chosen to be a compact set (with respect to the Euclidean norm in $\mathbb{R}^2$) of all $N$ ordered points such that $z'_j\in B_R\backslash \overline{D_\Gamma}, j = 1, \cdots, N$, where $D_\Gamma$ is the domain enclosed by $\Gamma$ such that $\Gamma = \partial D_\Gamma$.

Remark 2.2. It deserves noting that, the purpose for choosing the admissible curve to be star-like is to provide convenience for the subsequent analysis and calculation, while the exact boundary curve $\partial D$ does not necessarily require to be star-like. In fact, even if $\partial D$ is not star-like, the obstacle-source pair may also be well reconstructed, which will be illustrated by the numerical experiments in Section 4.

2.2.   Convergence analysis

In this subsection, we will investigate the convergence properties of the optimization method with an analogous analysis to that of the decomposition methods ^[7,8].

For later use, we define the acoustic single-layer operator $\mathcal{S}_1:L^2(\Lambda)\to L^2(\Gamma)$ by

where $\Gamma$ is a closed $C^2$ curve containing $\Lambda$ in its interior. To start with, we state some properties of the operators $\mathcal{S}$ and $\mathcal{S}_1$ (cf.^[7]).

Lemma 2.1. The single-layer operators $\mathcal{S}$ and $\mathcal{S}_1$, defined by (2.1) and (2.6), respectively, are injective and have dense range provided that $k^2$ is not a Dirichlet eigenvalue for the negative Laplacian in the interior of $\Lambda.$

Combining the minimization of the Tikhonov functional for (2.4) and the defect minimization (2.2), we define the cost functional $\mu:(L^2(\Lambda))^N\times U\times V\to \mathbb{R}$ by

where $\varphi = (\varphi_1, \cdots, \varphi_N)\in(L^2(\Lambda))^N,$ $z' = (z'_1, z'_2, \cdots, z'_N)\in V,$ $\alpha > 0$ is the regularization parameter in (2.3).

Now, we turn to analyze the convergence properties of the optimization method. To this end, we first introduce the following definition of optimal obstacle-source pair.

Definition 2.1. Given the measured total field $u(x; z_j), j = 1, \cdots, N$ and a regularization parameter $\alpha > 0,$ an obstacle-source pair $(\Gamma^*, z^*)\in U\times V$ is called optimal if there exists $\varphi^*\in (L^2(\Lambda))^N$ such that $\varphi^*$ and $(\Gamma^*, z^*)$ minimize the cost functional (2.7) simultaneously over all $\varphi\in(L^2(\Lambda))^N$ and $(\Gamma, z')\in U\times V,$ i.e.,

where

In terms of Definition 2.1, the following theorem establishes the existence of the optimal solution.

Theorem 2.1. For each $\alpha > 0,$ there exists an optimal obstacle-source pair $(\Gamma^*, z^*)\in U\times V.$

Proof. Let $(\varphi^n, \Gamma^n, z^n)$ be a minimizing sequence in $\Big(L^2(\Lambda)\Big)^N\times U\times V$, i.e.,

Since $U, V$ are compact, we can assume that $\Gamma^n\to\Gamma^*,$ $z^n\to z^*,$ $n\to\infty.$ In this paper, the convergence $\Gamma^n\to\Gamma^*$ is understood in the sense that $\|r^n-r^*\|_{C^{1, \beta}(\mathbb{S})}\to0,$ $n\to\infty,$ with the functions $r^n$ and $r^*$ representing $\Gamma^n$ and $\Gamma^*$ via (2.5), respectively. From

with $\alpha > 0,$ we conclude that the sequence $\{\varphi^n\}$ is bounded, i.e., there exists some constant $C$ such that $\|\varphi^n\|_{(L^2(\Lambda))^N}\le C, \forall n.$ Further, we can assume that $\varphi^n\rightharpoonup \varphi^*\in(L^2(\Lambda))^N, n\to\infty.$ Since $\mathcal{S}$ is compact, it follows that $\mathcal{S}\varphi_j^n\to\mathcal{S}\varphi_j^*, n\to\infty, j = 1, \cdots, N.$ The fact that the Hankel function $H_0^{(1)}$ is continuous leads to the following convergence

Hence,

By Taylor's formula, one can estimate that for all $\hat{x}\in\mathbb{S},$

with

From the Schwarz inequality, it holds that

for all $\hat{x}\in\mathbb{S}.$ Denote the dependence of $\mathcal{S}_1:L^2(\Lambda)\to L^2(\Gamma^n)$ on $n$ by writing $\mathcal{S}_1^n,$ and the uniform convergence $\mathcal{S}_1^n\varphi_j^n\to\mathcal{S}_1\varphi^*_j(n\to\infty$) holds. Therefore, we can deduce that

Further,

Since $\varphi_j^n\rightharpoonup \varphi_j^*, n\to\infty, j = 1, \cdots, N,$ it follows that

Thus, we have the norm convergence $\left\|\varphi_j^n-\varphi_j^*\right\|\to0, n\to\infty$ for $j = 1, \cdots, N.$ Finally, it holds that

and the proof is complete.

Before heading for the convergence analysis of the optimization method, we shall state the convergence property of the cost functional first.

Theorem 2.2. Let $z_j,$ $j = 1, \cdots, N$ be the exact source points, $u(x; z_j)$ be the exact total fields corresponding to the sound-soft obstacle $D$ such that $\partial D$ belongs to $U.$ Then the cost functional $m(\alpha)$ is convergent in the sense that

Proof. By Lemma 2.1, given $\varepsilon_j > 0,$ there exists $\varphi_j\in L^2(\Lambda)$ such that

Since the scattered field of a radiating solution to the Helmholtz equation depends continuously on the boundary (cf. ^[7]), there exists some constant $C_1$ such that

i.e.,

where we have used the fact that $u(x; z_j) = 0, x\in\partial D$. The continuity of the Hankel function $H_0^{(1)}$ implies that there exists $\delta_j > 0$ such that for $z'_j$ satisfying $|z_j-z'_j|\le\delta_j$, it holds that

Based on the above arguments, it holds that

Further, we deduce that

The theorem is accomplished due to the arbitrariness of each $\varepsilon_j, j = 1, 2, \cdots, N.$

We are now in a position to state the main convergence result.

Theorem 2.3. Let $\{\alpha_n\}$ be a null sequence, $(\Gamma^n, z^n)$ be the corresponding sequence of optimal obstacle-source pairs for the regularization parameter $\alpha_n.$ Then there exists a convergent subsequence of $\left\{(\Gamma^n, z^n)\right\}.$ Assume that $u$ is the exact total field of a domain $D$ such that $\partial D$ is contained in $U.$ Then every limit point $(\Gamma^*, z^*)$ of $(\Gamma^n, z^n)$ represents an optimal obstacle-source pair. In other words, if the incident wave due to the source point $z^*$ is scattered by the sound-soft obstacle $D^*$ bounded by $\Gamma^*$, then the total field vanishes on $\Gamma^*$.

Proof. The existence of convergent subsequences of $(\Gamma^n, z^n)$ follows from the compactness of $U\times V$. Let $(\Gamma^*, z^*)$ be the limit point. Without loss of generality, we assume that $(\Gamma^n, z^n)\to(\Gamma^*, z^*), n\to\infty,$ i.e., $\Gamma^n\to\Gamma^*,$ $z^n\to z^*,$ $n\to\infty.$ Let $u_j^*$ be the solution to the direct scattering problem with the incident wave emitted by the source point $z_j^*$ for the obstacle with boundary $\Gamma^*,$ i.e., the boundary condition reads

Since $(\Gamma^n, z^n)$ is optimal for the parameter $\alpha_n$, there exists $\varphi^n\in (L^2(\Lambda))^N$ such that

for $n = 1, 2, \cdots$. From Theorem 2.2, we conclude that the boundary data satisfies

From the condition that $\Gamma^n\to\Gamma^*, z_j^n\to z_j^*, n\to\infty,$ we obtain that $u_j^s = u_j^*$ with $u_j^s = \lim\limits_{n\to\infty}\mathcal{S}_1^n\varphi_j^n,$ and the total field $u^s_j(\cdot)+u^i(\cdot; z^*_j)$ must vanish on $\Gamma^*$ as a result of (2.10).

Remark 2.3. As far as we know, even for the standard decomposition method for the inverse obstacle problem, a uniqueness result for the optimal solution is not available. It is commented in ^[7] that

"Since we do not have uniqueness either for the inverse scattering problem or for the optimization problem, in the above convergence analysis we cannot expect more than convergent subsequences."

For the more complicated co-inversion problem considered here, the uniqueness for the optimal solution in Theorem 2.1 is still open. Therefore, the above comments from ^[7] are also applicable to our convergence analysis.

Based on the previous discussions, the rudiments of our optimization procedure is summarized in Algorithm 1.

Remark 2.4. For simplicity of implementation, the regularization parameter $\alpha$ is chosen by the trial-and-error strategy and the error tolerance is independent of $\alpha$ in this paper. To improve the performance of the algorithm, $\alpha$ could be chosen by, for instance, the Morozov's discrepancy principle according to the noise level. Meanwhile, the error tolerance should dependent on the value of $\alpha$, and may not be reached if $\alpha$ is too large.

At the end of this section, we would like to remark that, the nonlinear optimization problems generally suffer from the notorious difficulty of local minimums. Specifically, performance of any iterative solver for the optimization problem (2.8) would be unavoidably affected by the initial guess $(\Gamma^0, z^0)$. Therefore, choosing an appropriate initial value is highly crucial for the success of our optimization method. Nevertheless, the a priori information of the rough $(\Gamma^0, z^0)$ is not straightforwardly available. And it is thus vital to develop a fast and cheap strategy for resolving this crux, which will be the main concern of the next section.

3.   Direct sampling for choosing the initial guess

In the previous section, we develop an optimization method to reconstruct the obstacle-source pair $(\partial D, S)$ from the total field $\mathbb{U}$ and it has been pointed out that the initial guess $(\Gamma^0, z^0)$ should play an important role in the quality of reconstruction. Noticing that the direct sampling methods (DSM) are independent of any a priori information of the unknown objects and easy to implement, we wish to utilize the sampling-type methods to seek the initial guess for the source $S$ and the obstacle $D$, respectively.

3.1.   Determine the initial guess of source points

Motivated by the direct sampling schemes for inverse obstacle/source scattering problems, we first propose the total field version of indicator functions to reconstruct the source points by knowledge of $\mathbb{U}$. Let $\Omega_1\subset B_R$ be a bounded sampling domain such that $D\cup S\subset\Omega_1.$ For any sampling point $y\in\Omega_1$, we define $N$ indicator functions as follows:

Note that the indicators rely on the inner product of the fundamental solution with respect to the total field data rather than the usual incident or scattered data. We expect to locate the peak of $I_j(y)$ for each $j$ and take the corresponding maximum point $\tilde{z}_j$ as an approximation of the exact source $z_j$.

Before analyzing the characteristics of $I_j(y),$ we first claim two lemmas which will play a key role in the subsequent analysis.

Lemma 3.1. [17,Lemma 3.1] Let $g \in H^{1 / 2}\left(\partial{D}\right).$ Then the scattering problem

admits a unique solution $w \in H_{\mathrm{loc}}^{1}\left(\mathbb{R}^2 \backslash \overline{D}\right)$. Moreover, there exists a constant $C > 0$ such that

with $\nu$ the unit outer normal to the boundary $\partial D$.

For any $g \in H^{1 / 2}\left(\partial{D}\right)$, the Dirichlet-to-Neumann mapping $\mathbb{T}: H^{1 / 2}\left(\partial{D}\right) \rightarrow H^{-1 / 2}\left(\partial{D}\right)$ for the scattering problems (3.2)–(3.4) is defined as

By Lemma 3.1, $\mathbb{T}$ is a bounded linear operator and we will denote $\|\mathbb{T}\|$ its operator norm in this paper.

Analogous to [17,Lemma 3.2], one can immediately obtain the following estimate concerning the scattered field.

Lemma 3.2. Let $L: = \min_{1\le j\le N}{\rm dist}\left(z_j, \overline{D}\right)$. Then for each $j = 1, \cdots, N$, the following estimate holds

where $x\in\Gamma_R$ is the receiver and the constant $C_{js}$ is dependent of $k, R$ and $L$.

The following theorem sheds light on the indicating properties of $I_j(y), j = 1, \cdots, N$.

Theorem 3.1. Let $I_j(y), j = 1, \cdots, N,$ be defined by (3.1) and assume that

Then there exists a constant $C_j > 0$ such that for each $y\in \Omega_1, \ \ y\ne z_j,$

Proof. We first recall the following estimate for the Hankel function [17,(3.8)]:

Then from the representation of incident field $u^i(x; z_j),$ it can be easily derived that

Further, by Lemma 3.2, we deduce that there is a positive constant $C_{js}$ such that

Combining (3.8) and (3.9), together with the aid of $|x-y|\ge\left||x-z_j|-|y-z_j|\right|$, we have

where

Assumption (3.5) implies that there exists a positive constant $C_{j2}$ such that $L_j\ge (1+C_{j2})\sup_{y\in\Omega_1}|y-z_j|$. Then

where

and this completes the proof.

In virtue of the above theorem, each function $I_j(y)$ should decay as the sampling point $y$ recedes from the corresponding source point $z_j$. Hence, the initial detection can be conveniently chosen by locating the point $\tilde{z}_j$ where $I_j(y)$ reaches its apex. Then in the optimization process, the exact source location $z_j$ could be effortlessly and accurately found in the vicinity of $\tilde{z}_j$. This indicating behavior will be numerically verified by the experiments in the next section.

3.2.   Determine the initial guess of the obstacle

The indicator function proposed in the previous subsection has the capability of locating the approximate source points, but it usually fails to image the shape of obstacle. Thus a further sampling step for identifying the rough obstacle should be taken into account. By integrating the approximate source points into the reverse time migration (RTM) approach ^[18], we propose an approximate reverse time migration method to reconstruct the obstacle in this subsection.

To reconstruct the obstacle, information about the scattered field is indispensable in the classical RTM. However, the scattered field can not be measured directly due to the existence of unknown source points. So we are going to approximate the scattered field by subtracting the incident field $u^i(x; \tilde{z}_j)$ due to the numerical source point $\tilde{z}_j$ from the measured total field $u(x; z_j),$ i.e.,

and develop an imaging functional $I_D(y)$ using the approximate scattered field $u_j^s(x)$. To this end, confining the imaging domain to $\Omega_2\subset\Omega_1$ such that $\overline{D}\subset \Omega_2,$ $S\subset\Omega_1\backslash\overline{\Omega}_2,$ and we define the imaging function for $y\in\Omega_2$ by

Note that, in (3.11), instead of the exact scattered field $u^s(x; z_j)$ corresponding to the $j$-th source point $z_j,$ we use an approximate scattered field $u^s_j(x)$ for $z_j$. That is why this method is called the approximate reverse time migration method. Once the source points are retrieved, the indicating behavior of the approximate RTM immediately follows the original RTM ^[17,18]. The imaging function $I_D(y)$ provides a qualitative information of the obstacle $D$, which will be used as a useful clue to discover the obstacle in the optimization method. We emphasize that, if this sampling procedure for $D$ is not implemented, then it would be difficult to find a good initial guess for $D$ or a suitable internal curve $\Lambda$. In this case, such an apriori information about the obstacle is usually necessary to facilitate the realization of Algorithm 1.

We end this section with a brief description of Algorithm 2, which can be used to determine the initial guess for the optimization method proposed in Section 2.

4.   Numerical experiments

In this section, we will present several numerical experiments to demonstrate the performance of the proposed method. In our experiments, synthetic total field data are generated by reformulating the direct problems (1.2)–(1.4) as a boundary integral equation and the integral equation is solved by the Nyström method ^[7]. Here we use a numerical quadrature rule with $64$ equidistant grid points on $[0, 2\pi]$. The receivers are set to be $x_r = 4(\cos\theta_r, \sin\theta_r),$ $\theta_r = {r\theta}/{N_R},$ $r = 1, \cdots, N_R$ with $\theta\in (0, 2\pi]$ the observation aperture. For the full-aperture case, $\theta = 2\pi,$ $N_R = 120$. For the limited-aperture case, we consider two cases with $\theta = 3\pi/2,$ $N_R = 90$ and $\theta = \pi,$ $N_R = 60$, respectively. The forward solver provides us with the total field data

To test the stability of the numerical method, random noise of different levels is added to the total field data $u(x; z_j)$ to produce the noisy total field data

where $r_1,$ $r_2$ are two uniformly distributed random numbers ranging from $-1$ to $1$ and $\varepsilon > 0$ is the noise level. We use $\varepsilon = 10\%$ in Examples 1–3 and $\varepsilon = 5\%$ in Example 4.

In all the subsequent examples, the sampling domains $\Omega_1$ are chosen as squares centered at the origin with $\mathcal{T}_1$ being a $200\times200$ uniform grid. The imaging domains $\Omega_2$ are individualized case by case, according to the sampling results for the point sources. The imaging points $\mathcal{T}_2$ are selected as $100\times100$ uniform grids distributed over $\Omega_2$. For the optimization method, we choose the auxiliary curve $\Lambda$ to be a circle centered at the origin with radius $r_\Lambda = 0.6.$ The integral over the unit circle $\mathbb{S}$ is numerically approximated by the trapezoidal rule with $64$ grid points. To solve the nonlinear least-squares problem, the Levenberg-Marquardt algorithm is adopted with both functional value stopping criterion and successive iterate stopping criterion chosen to be $10^{-6}$ for the full-aperture case and $10^{-4}$ for the limited-aperture problems. The starting curve $\Gamma_0$ for the Levenberg-Marquardt algorithm is a circle centered at the origin with radius $r_{\Gamma_0} = \left|\arg\max_{y\in\Omega_2}I_D(y)\right|$. For $j = 1, \cdots, N$, the starting guesses for the source points $\{z_j^0\}_{j = 1}^N$ are chosen to be $\{\tilde{z}_j\}_{j = 1}^N$. The regularization parameter $\alpha$ is selected by trial and error. The ansatz function $r(t)$ is chosen in the $(2M+1)$-dimensional subspace $U_M\subset U$ and can be expanded by trigonometric polynomials of degree less than or equal to $M,$ i.e.,

with some $M\in\mathbb{N}_+$. Unless otherwise specified, $M = 8$ is used in what follows.

To evaluate the reconstruction quantitatively, we introduce an equidistant set of knots on $[0, 2\pi]$ by $t_i: = 2\pi i/N_1, i = 0, 1, \cdots, N_1-1,$ with $N_1$ the number of knots chosen to be 256 in our numerical experiments. The discrete relative $L^2$ error of the boundary curve is defined by

with $r^*$ standing for the exact boundary curve and $\tau_i = \tau(t_i)$ being $t_i$ or some transformation of $t_i$ depending on whether the exact boundary is star-like or not, respectively.

In the following figures regarding geometrical settings of the problem, the red and blue points denote respectively the exact source locations and receivers, the green solid line denotes the boundary of the obstacle, the black and purple dashed lines denote the boundary of the sampling domains $\Omega_1$ and $\Omega_2,$ respectively. In the figures about the results of the numerical experiments, the small black '+' markers are the exact source points, the small red circles stand for the reconstructed source points, the green solid line denotes the exact boundary of the obstacle, the blue dashed line denotes the reconstructed boundary. In addition, for ease of comparison, the indicator functions $I_j(y),$ $j = 1, \cdots, N,$ and the imaging function $I_D(y)$ are normalized so that the maximum value of each function is $1.$

Example 1. In the first example, we consider the reconstruction of a starfish-shaped obstacle and the influence of source quantities. The exact boundary of the starfish-shaped obstacle is parameterized by

We first consider the case where the five source points are respectively located at

The sampling domains are chosen to be $\Omega_1 = [-2.5, 2.5]\times[-2.5, 2.5]$, $\Omega_2 = [-1.5, 1.5]\times[-1.5, 1.5]$. The regularization parameter and the wave number are set to be $\alpha = 10^{-8}$ and $k = 5$, respectively.

Figure 2(a) exhibits the geometry setting. Figure 2(b) shows the reconstructed locations of $S_1$ by the DSM and we can find that the reconstructed locations of the source points are very close to the exact locations. Figure 2(c) shows the image of the imaging function $I_D(y)$ over the sampling domain $\Omega_2.$ We find that in Figure 2(c), the shape of the obstacle can be roughly captured.In Figure 2(d), reconstruction of the optimization method is depicted, where the initial guess for the source points are chosen to be the numerical source points obtained by the DSM and the initial curve is represented by the black dashed circle in Figure 2(c). Further, we compute the error of boundary curve by (4.2) with $\tau_i = t_i, i = 1, \cdots, N$ and in this case $E_D = 5.73\%$. In Table 1, we list the exact and reconstructed locations of source points $S_1$ by the DSM and the optimization method, respectively. According to Table 1 and Figure 2, we want to point out that both the DSM and optimization method could produce satisfactory reconstructions of the source points. A promising advantage of the optimization method is that the obstacle can also be reconstructed quantitatively, i.e., a significantly more sharper profile can be identified by the optimization method.

Next we consider the reconstruction from a different set of source points, which is given by $S_2 = \{(1.3, 0), (1.35, 0.1), (0, 2), (0, -2), (-2, 0)\}$. In this case, the sampling domains are adapted to be $\Omega_1 = [-2.2, 2.2]\times[-2.2, 2.2],$ $\Omega_2 = [-1.2, 1.2]\times[-1.2, 1.2]$. The other parameters remain unchanged. We refer to Figure 3(a) for the geometry setting.

Figure 3(b) shows the indicator $I_1(y)$ over the sampling domain $\Omega_1$. One can observe that the maximizer of $I_1(y)$ coincides well with the exact source location. Figure 3(c) shows that all the reconstructed source locations match well with the exact positions. In Figure 3(d), the imaging function $I_D(y)$ is plotted. From Figure 3(d), we find that the maximizer of $I_D(y)$ incurs on the boundary of the obstacle.

Compared with Figure 2(c), it seems that the current result deteriorates slightly. A reason accounting for this fact is the specific distribution of the source points where two points are almost adjacent to each other.In Figure 3(e), we show the reconstruction of the starfish by the optimization method. Figure 3(e) demonstrates that the quality of reconstruction for the left side of the obstacle is better than that of the right side. Physically speaking, this is due to the fact that the right side of the obstacle is affected by the nearby source points 1 and 2. Moreover, the adjacent source points 1 and 2 could be individually identified as long as the sampling grid is sufficiently fine. To see the local reconstruction more clearly, we intercept a part $[1, 1.4]\times[-0.1, 0.2]$ of Figure 3(e) and it is zoomed in as Figure 3(f).

In Table 2, we list the exact and reconstructed source locations. To evaluate the reconstruction of the obstacle quantitatively, we compute the errors between the exact boundary curve and the reconstructed boundary curve as (4.2) with $\tau_i = t_i$. In this case, the error of the boundary curve is $E_D = 7.90\%$.

Comparing the above two cases, it can be found that our method could simultaneously produce a rather accurate reconstruction of the obstacle and sources, no matter whether the source points are relatively sparse or dense, or even quite close to the obstacle.

In Figure 4, we consider the reconstruction with $N( = 2, 4, 6, 8)$ source points. The wave number and the regularization parameter is chosen to be $k = 8$ and $\alpha = 10^{-6}$, respectively.The sampling domains are chosen to be $\Omega_1 = [-2.5, 2.5]\times[-2.5, 2.5],$ $\Omega_2 = [-1.5, 1.5]\times[-1.5, 1.5].$ Further, we take $k = 5, 8,$ respectively, and list the relative $L^2$ errors for the reconstruction of the starfish-shaped obstacle in Table 3. Table 3 shows that the error of reconstruction decreases as the number of source points increases.

Example 2. In this example, we aim to investigate the accuracy of reconstruction by considering a circular obstacle as well as different number of source points. The parametric boundary of the obstacle is given by $x_C(t) = (\cos t, \sin t), \ 0\leq t < 2\pi,$ and $\alpha = 10^{-8},$ $\Omega_1 = [-2.5, 2.5]\times[-2.5, 2.5],$ $\Omega_2 = [-1.2, 1.2]\times[-1.2, 1.2].$

To analyze the reconstruction quantitatively, the relative errors of the boundary curves are computed as (4.2) with $\tau_i = t_i$. Table 4 shows the reconstructions in different cases by varying the wave numbers and the number of source points. In particular, a suitable choice of the parameters (for instance, $k = 8$ and $N = 6$) would produce a stunningly accurate reconstruction. Figure 5 exhibits the reconstruction of the circle and $N( = 1, 2, 4, 6, 8)$ source points with wave number $k = 5$. In addition, it deserves noting that, even if we take $N = 1,$ both the source and the obstacle can be reconstructed, though the error is slightly large (see Figure 5(a) and the second column in Table 4).

Further, we test the effectiveness of the DSM proposed in Section 3 for choosing the initial guess. The parameters are chosen to be $\alpha = 10^{-9}, k = 8, N = 8$. The exact source points are given by

For comparison purpose, we consider the following two sets of artificially chosen initial guess

and

respectively. Taking $S_4$ and $S_5$ as the initial guesses for the source points respectively, we then compare the numerical results with the case where the initial guess is determined by the DSM. The relative $L^2$ errors for $k = 5$ are listed in Table 5. Furthermore, we present the reconstructions of the circle as well as the source points in Figure 6. From Table 5 and Figure 6, one can observe that the DSM plays a significant role in the reconstruction by providing a good initial guess for the optimization process, which is a key step for achieving a satisfactory reconstruction. If the initial source points somewhat deviate from the true ones, then the reconstruction would drastically deteriorate, see Figure 6(b). Therefore, results of the DSM indeed exert an immense impact on the effectiveness of the optimization. And a good initial knowledge generated by a high-resolution sampling grid $\mathcal{T}_1$ is the essential prerequisite for achieving the high-quality reconstruction.

Example 3. This example is designed to check the validity of our method when the boundary of obstacle is not star-like. Now we consider the reconstruction of two kite-shaped obstacles with boundaries $\partial D$ described by the parametric representations

respectively. In this example, the sampling domains are taken to be $\Omega_1 = [-2.5, 2.5]\times[-2.5, 2.5],$ $\Omega_2 = [-1.5, 1.5]\times[-1.5, 1.5]$. To compute the accuracy of reconstruction, the relative errors of the boundary curves are computed by (4.2) with

In Figure 7, we take the ansatz function space of different dimensions and present the reconstruction of the kite-I with $k = 5$ and $\alpha = 10^{-6}.$ By taking different $M( = 2, 4, 6, 8, 12, 20)$ in (4.1), we compute the reconstruction errors defined by (4.2). From the errors in Table 6, we observe that the accuracy of inversion is insensitive to moderately large $M$. Except for the cases $M = 2$ and $M = 4$, the reconstruction is almost negligibly influenced by the dimension of $U_M$.

Next, we fix $M = 8$ and reconstruct the two kites. In Table 7, we provide the errors for reconstructing the kites with different number of source points. These results illustrate that the reconstructions of kite-I are superior to those of kite-II, which is probably due to the fact that the boundary fluctuation of kite-I is relatively mild in comparison to that of kite-II.

In Figure 8, we plot the reconstruction of kite-II with $N(N = 2, 6, 8, 10)$ source points. It can be found that both the source points and the obstacle are reconstructed well. So the proposed method also works for the recovery of non-starlike scatterers.

Example 4. In the last example, we plan to reconstruct the obstacle-source pair with limited-aperture measurements. Motivated by realistic applications where only partial data is available, we test the proposed method when the total field is only accessed on a portion of the measurement curve. In this example, we collect the total field recorded by the receivers located on a part of $\Gamma_R$ with the aperture chosen to be $\theta = 3\pi/2$ and $\pi,$ respectively. The regularization parameter is chosen to be $\alpha = 10^{-6}$. The sampling domains are $\Omega_1 = [-2.5, 2.5]\times[-2.5, 2.5]$ and $\Omega_2 = [-1.5, 1.5]\times[-1.5, 1.5]$.

In Table 8, we list the relative $L^2$ errors of the boundary curves with respect to different parameters $k$, $\theta$ and $N$. Figure 9 shows the reconstruction for the starfish together with $N( = 4, 8)$ source points in the limited-aperture case with $k = 8$. Figure 9 illustrates that the part of boundary curve encompassed by the receivers can always be well reconstructed. Typically, due to the lack of information, recovery of the opposite side oriented towards the open aperture is less accurate. Anyhow, all the numerical results show that, even with limited-aperture data, the reconstruction method still reasonably works.

5.   Conclusions

In this work, we develop a novel optimization method to simultaneously reconstruct a sound-soft obstacle and multiple source points from near field scattering data. To obtain a good initial guess for the optimization method, a two-step sampling method is also developed to respectively image the rough source points and obstacle. Theoretically, we analyze the convergence of the optimization method and the indicating behaviors of the sampling scheme. Numerical results are provided to demonstrate the feasibility and promising features of the proposed method. Our future work consists of the theoretical justification of the uniqueness issue, which may be established by exploiting the nonlinearity of the DtN map and the Complex Geometrical Optics (CGO) solutions as in ^[29]. We will also consider the extension to three-dimensional model and the other boundary conditions, as well as the applicability in the scenarios of electromagnetic or elastic waves.

Acknowledgments

This work was supported by the NSFC grant under NO. 11971133 and the Fundamental Research Funds for the Central Universities. The authors would like to thank Professor Xianchao Wang for valuable discussions. We would also like to thank the anonymous referees for their insightful and constructive suggestions on improving the quality of this paper.

Conflict of interest

The authors declare there is no conflicts of interest.