A novel quantitative inverse scattering scheme using interior resonant modes

In this article, we are concerned with the inverse acoustic scattering problem of reconstructing an impenetrable obstacle by the associated far-field measurement. To begin with, we present the mathematical setup of the inverse scattering problem for our study.

Let $k\in\mathbb{R}_{+}$ be the wavenumber of a time harmonic wave and $D\subset\mathbb{R}^{m},m=2,3$ be a bounded domain with a Lipschitz-boundary $\partial D$ and a connected complement $\mathbb{R}^{m}\setminus\overline{D}$, which signifies the target obstacle in our study. We take the incident field $u^{i}$ to be a time harmonic plane wave of the form

where $\mathrm{i}:=\sqrt{-1}$ is the imaginary unit, $d\in\mathbb{S}^{m-1}$ is the direction of propagation and $\mathbb{S}^{m-1}:=\{x\in\mathbb{R}^{m}:|x|=1\}$ is the unit sphere in $\mathbb{R}^{m}$. Let $u^{s}$ and $u:=u^{i}+u^{s}$ signify the scattered and total wave fields, respectively. The forward scattering problem is described by the following Helmholtz system:

where the last limit is known as the Sommerfeld radiation condition which holds uniformly in $\hat{x}:=x/|x|\in\mathbb{S}^{m-1}$ and characterizes the outgoing nature of the scattered field. In (1.1), $\mathcal{B}(u)=u$ or $\mathcal{B}(u)=\partial u/\partial\nu$, respectively, signify the physical scenarios that the obstacle $D$ is sound-soft or sound-hard. Here and also in what follows, $\nu\in\mathbb{S}^{m-1}$ denotes the exterior unit normal vector to $\partial D$. The well-posedness of the scattering system (1.1) can be conveniently found in [25]. There exists a unique solution $u\in H_{loc}^{1}(\mathbb{R}^{m}\backslash\overline{D})$, and it admits the following asymptotic expansion [11]:

which holds uniformly for all directions $\hat{x}:=x/|x|\in\mathbb{S}^{m-1}$. The complex-valued function $u^{\infty}$ defined on the unit sphere $\mathbb{S}^{m-1}$ is known as the far-field pattern of the scatter field $u^{s}$, which encodes the information of the scattering obstacle $D$. The inverse scattering problem of our concern is to recover $D$ by knowledge of $u^{\infty}$; that is,

where $V$ is an open interval in $\mathbb{R}_{+}$, and $\mathcal{F}$ is an abstract operator defined by (1.1). Though the forward scattering problem (1.1) is linear, it can be straightforwardly verified that the inverse problem (1.2) is nonlinear.

The inverse scattering problem (1.2) is prototypical model of fundamental importance for many scientific and industrial applications including radar/sonar, medical imaging, geophysical exploration and nondestructive testing. There are rich results in the literature for (1.2), both theoretically and numerically. It is known that one has the unique identifiability for (1.2), namely the correspondence between $u^{\infty}$ and $D$ is one-to-one. Many numerical reconstruction algorithms have been developed in solving (1.2) as well. In fact, it is impossible for us to list and discuss all of the existing numerical studies for (1.2) in the literature. In what follows, we only discuss a few selective ones to motivate our current study. Roughly speaking, the existing numerical approaches can be classified into two categories: quantitative ones and qualitative ones. Quantitative approaches employ Newton’s linearization and/or optimization strategies to tackle (1.2) directly. We refer the reader to the Newton’s iterative method [11, 15, 36], the decomposition method [11, 14, 16] and the recursive linearization method [2, 3] for the relevant results. The other class of approaches resorts to establishing a criterion to distinguish the interior and exterior of the obstacle $D$, and hence can qualitatively reconstruct the boundary of $D$. These include the linear sampling method [5, 7, 10, 22], the factorization method [17], the direct sampling method [13, 21, 28] and the probe method [31, 32].

In this article, we develop a novel quantitative reconstruction scheme for (1.2). A key ingredient of our method is to connect the exterior scattering problem (1.1) to its interior counterpart:

where $v\in H^{1}(D)$. (1.3) is the classical Dirichlet/Neumann Laplacian eigenvalue problem and $v$ represents an interior resonant mode. It has been well understood, say e.g. there exist infinitely many discrete eigenvalues accumulating only at $\infty$, and the eigenspace associated with each eigenvalue is finite-dimensional. It turns out that the the exterior problem (1.1) and the interior problem (1.3) are dual to each other, in particular in the sense that one can read off the spectral system of (1.3) from the far-field data of (1.1). This is done by employing a certain indicating behaviour of the linear sampling method and by solving a constrained optimization problem, both only involving the far-field data in a simple manner, which form the first phase of the proposed method. In the second phase, we derive an iteration formula of Newton’s type which can be used to quantitatively reconstruct the boundary of the obstacle. By utilizing the interior eigenfunctions determined in the first phase, the shape derivative involved in the Newton’s iteration can be calculated in a very efficient way without the need to solve any forward scattering problem. To highlight the novel contributions of the current study, we present three remarks as follows.

First, determining the interior eigenvalues from the associated far-field data via the linear sampling method has been studied in [6, 19, 23, 29]. In [23], the determination of the interior Dirichlet Laplacian eigenfunctions has been further investigated. In fact, it is proposed in [23] that the Dirichlet eigenfunctions can be used to recover $\partial D$ (in a qualitative manner). This is natural since the Dirichlet eigenfunctions vanish on $\partial D$. However, such an idea cannot be extended to the sound-hard case since identifying the place where the Neumann data vanish requires knowing the normal vector of $\partial D$ which is equivalent to knowing $\partial D$. It is our aim of developing an imaging scheme that can recover sound-hard obstacles by making use of the interior resonant modes. It turns out that the scheme we develop can not only recover sound-hard obstacles but can also yield highly-accurate quantitative reconstructions. The key idea is to establish an iteration formula of Newton’s type by using the homogeneous Neumann condition on $\partial D$, which can then be used for identifying the unknown $\partial D$. A salient feature of this Newton’s iteration approach is that it is essentially derivative-free in the sense that the involved shape derivatives can be explicitly calculated by using the Neumann eigenfunctions that have been already determined, and there is no need to calculate any direct scattering problem. Moreover, the idea can be directly extended to the sound-soft case to produce quantitative reconstructions that outperform the qualitative reconstructions in [23]. However, it can be seen even at this point that the sound-hard case is technically more challenging than the sound-soft case. Hence, we shall mainly stick our study to the sound-hard case and only remark the sound-soft case throughout the rest of the paper.

Second, the far-field data used in (1.2) is significantly overdetermined. In fact, most of the existing methods in the literature, in particular those qualitative ones, make use of far-field data corresponding to a fixed $k$. Moreover, it is widely conjectured that one can determine an obstacle by at most a few or even a single far-field measurements, namely a few $k$ and $d$ or even fixed $k$ and $d$ in (1.2); see e.g. [8, 9, 24, 27, 33] for related theoretical uniqueness and stability results and [21] for related numerical reconstructions if a-priori geometrical knowledge is available on $D$. The overdetermination in (1.2), especially on $k$, is mainly needed to compute the interior eigenvalues located inside $V$, which is the prerequisite to the eigenfunction determination. As mentioned in the first remark, the eigenfunctions are critical to avoid computing shape derivatives for the Newton’s iteration in our method. Hence, the overdetermined data are an unobjectionable cost for achieving both high accuracy and high efficiency. It is interesting to note that all of the aforementioned qualitative methods inevitably involve shape derivatives and hence the calculation of a large amount of forward scattering problems, and moreover suffer from the local minima issue. In fact, to our best knowledge, no 3D quantitative reconstruction of a sound-hard obstacle was ever conducted in the literature due to the highly complicated computational nature. In Section 4, we present both 2D and 3D reconstructions and it can be seen that our method is computationally straightforward. Moreover, it can be seen that our method is robust and insensitive to the initial guess, and this is physically justifiable as the interior resonant modes carry the geometrical information of $D$ in a sensible way (though implicitly in the sound-hard case). It is also interesting to note the similarity shared by our method and the machine learning approaches for inverse obstacle problems. In [12, 35], machine learning approaches were developed that can yield a highly-accurate reconstruction of a target obstacle by using only a few far-field measurements. However, a large amount of data as well as computations are needed to train the neural networks therein, and hence are computationally more costly.

Third, we would like to mention two practical scenarios where our method might find applications to corroborate our viewpoint in the above remark. In many practical applications, say e.g. photo-acoustic tomography, time-dependent data are collected which can yield multiple frequency data as needed in (1.2) via temporal Fourier transform [26]. The other application is the bionic approach of generating human body shape [20], where overdetermined data are not a practical drawback, but highly accurate reconstructions are needed. We shall consider the application of our method in those practical setups in our future work.

The rest of the paper is organized as follows. Section 2 introduces the linear sampling method to determine the interior eigenvalues from the multi-frequency far-field data. Then by the Herglotz wave approximation, we present an efficient optimization algorithm to reconstruct the interior eigenfunctions. In Section 3, a novel Newton iterative formula is proposed to reconstruct the sound-hard obstacle via the use of the interior eigenfunctions. Numerical experiments are conducted in 4 to verify the promising features of our method.

In this section, we aim to determine the interior eigenvalues and eigenfunctions to (1.3) from the far-field data in (1.2) corresponding to the unknown $D$. As discussed in the previous section, we shall mainly focus on the sound-hard case, namely $\mathcal{B}(u)=\partial u/\partial\nu$ and only remark the extension to the sound-soft case.

To begin with, we introduce the linear sampling method for reconstructing the Neumann eigenvalues. The linear sampling method is a widely used method to recover the shape of the scatterer without a priori information of the boundary condition of the scatterer. The basic idea of the linear sampling method is choosing an approximate indicator function and distinguishing whether the sampling point is inside or outside the scatterer. To that end, we present the indicator function of the linear sampling method. Define the test function by

where $z\in\mathbb{R}^{m}$ denotes the sampling point. The key ingredient of the linear sampling method is to find the kernel $g\in L^{2}(\mathbb{S}^{m-1})$ as a solution to the following integral equation

where $F_{k}$ is the far field operator from $L^{2}(\mathbb{S}^{m-1})$ to $L^{2}(\mathbb{S}^{m-1})$ and it is defined by

We usually adopt Tikhonov regularization to solve the Fredholm integral equation (2.1) since this equation is ill posed. Due to the existence of noisy data in practice, we suppose that $F_{k}^{\delta}$ is the corresponding operator to the noisy measurements $u^{\infty,\delta}$ and then instead seek the unique minimizer $g_{z,k}^{\delta}\in L^{2}(\mathbb{S}^{m-1})$ of the following functional

where $\epsilon$ is a regularization parameter.

Next, we discuss how to recover the Neumann eigenvalues by using the linear sampling method. According to the rigorous justification of the rationale behind the linear sampling method (see [1, 29]), one can use the following lemma to distinguish whether $k$ is a Neumann eigenvalue or not.

Next, we proceed to recover the corresponding Neumann eigenfunctions. Recall that the Herglotz wave function is defined by

where $g\in L^{2}(\mathbb{S}^{m-1})$ is called the Herglotz kernel of $v_{g,k}$. We first show that the Herglotz wave function can be used to approximate the solution of the Helmholtz equation by the following lemma.

The following theorem plays an important role to determine the Neumann eigenfunctions.

By Theorem 2.1 and normalization if necessary, we can say that the following optimization problem:

exists at least one solution $g\in L^{2}(\mathbb{S}^{m-1})$ when $k$ is a Neumann eigenvalue in $D$. Since $D$ is unknown, it is unpractical to solve the optimization problem (2.12) with the constraint term $\|v_{g,k}\|_{L^{2}(D)}=\mathcal{O}(1)$. Thus, we consider the following optimization problem instead:

where $B$ is bounded domain such that $D\subset B$. In practice, one can simply choose $B$ to be a large ball containing $D$. In summarizing our discussion above, we can formulate the following scheme (Scheme I) to determine the Neumann eigenvalues and the corresponding eigenfunctions. It is emphasized that all the results in this section hold equally for the sound-soft case.

In this section, we develop a novel Newton-type method to reconstruct the shape of a sound-hard obstacle based on the interior resonant modes determined in the previous section.

Let $\gamma$ be a closed curve (in $\mathbb{R}^{2}$) or a closed surface (in $\mathbb{R}^{3}$), which is parameterized by

Here, $\gamma$ signifies the boundary of a bounded domain which shall be involved in the Newton’s iteration in what follows. For a fixed total field $u$, we define the operator $G$ that maps the boundary contour $\gamma$ onto the trace of $\partial u/\partial\nu$ on $\gamma$, that is,

In terms of the above operator $G$, we seek the parameterization $z$ of the boundary contour $\gamma$ such that the total field $u$ satisfies the Neumann boundary condition, i.e.,

Next, we introduce the Newton-type method to recover the boundary of the obstacle. Let $\gamma_{n}$ be the approximation to the boundary $\partial D$ at the $n$-iteration, $n=0,1,2,\cdots$. Following the idea of Newton method, we replace the previous nonlinear equation (3.2) by the linearized equation

where $h$ denotes the shift. Actually, the key point for solving the last linearized equation is to determine the Fréchet derivative of the operator $G$. Inspired by the work [18], we seek to improve and present a simpler Fréchet derivative in two and three dimensions, respectively.

For the two-dimensional case, we assume that $z=(z^{(1)},z^{(2)})^{\top}$. Then the exterior normal vector can be represented by

where $z_{\phi}^{\bot}=\left({\partial z^{(2)}}/{\partial\phi},\,-{\partial z^{(1)}}/{\partial\phi}\right)^{\top}$. Now we characterize the Fréchet derivative in the following theorem.

For the three-dimensional case, we define two orthogonal unit tangential vector fields $\tau_{1}$ and $\tau_{2}$ on the surface $\gamma$, that is,

where $z_{\theta}=\partial z/\partial\theta$ and $z_{\phi}=\partial z/\partial\phi$. Therefore, the normal vector can be defined by

In the following theorem, we characterize the Fréchet derivative of $G$ in three dimension.

To avoid calculating the direct scattering at each iteration, we use the eigenfunction $v_{g,k}$ to approximate the wave field $u$, namely,

Then the gradient and the Hessian Matrix of $u$ defined in Theorem 3.1 and 3.2 can be represented by

where the Herglotz kernel $g$ is determined by solving the optimization problem (2.13).

We would like to emphasize that it is non-trivial to show that $G^{\prime}$ is injective, namely, the operator $G^{\prime}$ may be not invertible. Therefore, we employ the standard Tikhonov regularization scheme with a regularization parameter $\alpha$ for solving linearized equation (3.3) at each iteration. As discussed above, we propose the following Newton-type scheme.

Noting that the Newton iteration method usually produces local minima for solving the inverse obstacle problems. To overcome local minimum and extend the convergence range, we propose a multifrequency Newton-type iterative algorithm.

Finally, we would like to remark the extension to the sound-soft case. Indeed, by modifying the definition of the operator $G$ in (3.1) via replacing $\partial u/\partial\nu$ by $u$, all the results in this section can be readily extended to the sound-soft case.

In this section, several numerical experiments are presented to verify the effectiveness and efficiency of the proposed methods. All the numerical results are conducted for recovering sound-hard obstacles, which present more challenges than the sound-soft case. To avoid inverse crime, the artificial far-field data are calculated by the finite element method, which is written as

Here $\hat{x}_{i}$ denotes the observation direction, $d_{j}$ denotes the incident direction and $k_{s}$ denotes the wavenumber. The observation and incident directions are chosen as the equidistantly distributed points on the unit circle (2D) or unit sphere (3D). Moreover, the wavenumbers are also equidistantly distributed in the open interval $V\subset\mathbb{R}_{+}$. To test the stability of the method, for any fixed wavenumber, we add some random noise to the measurement matrix:

where $\delta>0$ is the noise level, and $R_{1}$ and $R_{2}$ are two uniform random matrices that ranging from $-1$ to $1$.