A Direct Sampling Method and Its Integration with Deep Learning for Inverse Scattering Problems with Phaseless Data

Jianfeng Ning School of Mathematics and Statistics, Wuhan University, Wuhan, China. ([email protected]). Fuqun Han Department of Mathematics, University of Califonia, Los Angeles, USA. ([email protected]). Jun Zou Department of Mathematics, The Chinese University of Hong Kong, Shatin, N.T., Hong Kong. The work of this author was substantially supported by the Hong Kong RGC General Research Fund (projects 14306623, 14308322 and 14306921). ([email protected]).

Abstract

We consider in this work an inverse acoustic scattering problem when only phaseless data is available. The inverse problem is highly nonlinear and ill-posed due to the lack of the phase information. Solving inverse scattering problems with phaseless data is important in applications as the collection of physically acceptable phased data is usually difficult and expensive. A novel direct sampling method (DSM) will be developed to effectively estimate the locations and geometric shapes of the unknown scatterers from phaseless data generated by a very limited number of incident waves. With a careful theoretical analysis of the behavior of the index function and some representative numerical examples, the new DSM is shown to be computationally efficient, easy to implement, robust to large noise, and does not require any prior knowledge of the unknown scatterers. Furthermore, to fully exploit the index functions obtained from the DSM, we also propose to integrate the DSM with a deep learning technique (DSM-DL) to achieve high-quality reconstructions. Several challenging and representative numerical experiments are carried out to demonstrate the accuracy and robustness of reconstructions by DSM-DL. The DSM-DL networks trained by phased data are further theoretically and numerically shown to be able to solve problems with phaseless data. Additionally, our numerical experiments also show the DSM-DL can solve inverse scattering problems with mixed types of scatterers, which renders its applications in many important practical scenarios.

1 Introduction

Inverse scattering problems involve the recovery of the locations, shapes, and physical properties of unknown scatterers with some measurement data. These problems have wide applications in various fields such as identification of oil reservoirs, medical imaging, geophysical prospection, and microwave remote sensing [36, 6]. Numerous numerical methods have been developed in the literature over the past few decades to address inverse scattering problems with phased data. These methods include contrast source inversion method[41], recursive linearization methods [4], reverse time migration methods [8, 9], subspace optimization methods [10], linear sampling methods [7], direct sampling methods [24, 25], etc.

However, phased data is usually difficult and expensive to obtain in real applications [20, 26], especially at frequencies beyond tens of gigahertz. In these scenarios, since collecting high-accuracy phaseless data is much easier and more cost-effective, phaseless reconstruction is usually preferred in many important practical applications even the phaseless reconstruction is more ill-posed and non-linear compared to phased reconstruction. To tackle the inverse scattering problems with phaseless data, several optimization and iterative methods have been developed in the literature (see, e.g.,[1, 5, 35, 47]). While optimization and iterative methods can generally provide accurate reconstructions of unknown scatterers, they are mostly very expensive, require reasonably good initial guesses as well as some prior knowledge of the scatterers, which may not be easy to acquire in many scenarios. To achieve fast reconstruction, several non-iterative methods have gained popularity recently in inverse scattering problems with phaseless data. In [12], the reverse time migration for phased reconstruction was extended to reconstruct unknown scatterers from phaseless total field data. An approximate factorization method was proposed in [49] to recover the shapes and locations of scatterers from the phaseless total field data induced by incident plane waves. In [48], by superpositions of plane waves at a fixed frequency, a direct imaging approach was proposed to recover scattering obstacles from phaseless far-field data without suffering the translation-invariance obstruction. Moreover, introducing a reference ball into the scattering system also offers an alternative method to break the translation-invariance for phaseless far-field data [19]. By introducing the concept of scattering coefficients, the algorithms for phased and phaseless reconstructions in the linearized case were proposed and analyzed in [2]. For some theoretical results including the uniqueness and stability concerning inverse scattering problems with phaseless data, we refer to [28, 34, 29, 45].

In addition to the severely non-linear and ill-posed challenges of phaseless reconstruction, another challenge in some applications is that measurement data is available only from one or a few incident waves. The direct sampling methods (DSMs) introduced in [24, 25, 30] have been developed for phased reconstruction and can provide reasonable estimations of the locations and shapes of unknown scatterers using only limited incident waves, without requiring any prior knowledge of the scatterers. In addition to their data efficiency, the DSMs are also computationally efficient, involving only scalar products, completely parallel, and highly robust to large noise. Therefore, it is desirable to extend the DSMs to phaseless reconstruction. We also refer to some developments of the DSMs for other inverse problems in mathematical imaging [16, 13, 15, 14].

The main aims of this work are fourfold. The first one is to develop a novel DSM for inverse scattering problems with only one set of phaseless total field data. The implementation of the new DSM is as simple as that of the phased DSMs while it is still capable of providing a robust estimation of shapes and locations of scatterers. On the other hand, it is known that the DSMs have a limitation that they can not generate very accurate reconstructions in some important scenarios, e.g., for scatterers with complex geometrical shapes. So the second aim of this work is to alleviate this main limitation, by combining the DSMs with some special machine learning techniques. Furthermore, we will show by theoretical analysis and numerical experiments that the DSM-DL network trained by phased data can also be applied to solve problems with only phaseless data available. Finally, we further exploit the advantages of our novel algorithm to demonstrate that without any prior knowledge of the scatterers and with only one trained neural network, the DSM-DL with phaseless data can also solve inverse problems with inhomogeneous media or obstacles in the computational domain, as well as their mixed problems.

Machine learning techniques, especially deep learning ones, have shown promising potential recently for solving various PDEs [38, 31, 23] and inverse problems of PDEs [3, 40]. In particular, there are several advantages of solving inverse problems by employing neural networks to learn the inverse maps. Firstly, once the network is trained, the forward pass of the network is very fast and a real-time reconstruction is usually achievable. Secondly, in classical iterative methods, a prior knowledge of the unknown targets is incorporated into a regularization term to make the inversion more stable, whereas it is often challenging to choose a proper regularization and its balancing parameter. On the other hand, neural networks can automatically learn the distribution of unknown targets from training data, which helps improve the reconstruction quality.

Over the last few years, we have witnessed many successes of deep learning methods in solving inverse scattering problems [33, 43, 27, 21], and we refer to the review papers [11, 22] for more comprehensive studies. Some deep-learning approaches have also been proposed for inverse scattering problems with phaseless data. For instance, a rough image is initially obtained from the phaseless total field using the contrast source inversion scheme with a few iterations [44], followed by refinement using neural networks to produce a high-quality image. Another two-stage method presented in [32] involves a cascaded complex U-Net compromising a phase retrieval net and an image reconstruction net. In [46], a two-layer sequence-to-sequence neural network is proposed to recover obstacles with limited phaseless far-field data, where parameters representing the boundary curve of the impenetrable obstacles are utilized. Given the success of deep learning techniques in inverse scattering problems, we will first derive a novel DSM for phaseless reconstruction in this work, and then combine the new DSM with deep learning to produce more desirable reconstructions. This combination is a natural extension of our previous work DSM-DL [33], where the inverse medium scattering problem with phased data was considered. The DSM-DL in [33] was shown to be very robust to noise, easy to implement, and computationally efficient.

The remaining part of this paper is organized as follows. In Section 2, we introduce the inverse scattering problems addressed in this work. Section 3 is dedicated to a review of the DSMs and the development of a novel DSM for phaseless reconstruction. The combination of the DSM with deep learning is discussed in Section 4. In Section 5, we present numerical experiments for phaseless reconstruction using the DSM approach with a limited number of incident fields, as well as the imaging results achieved through DSM-DL. Finally, Section 6 concludes the paper with some closing remarks.

2 Problem formulations

We shall consider the inverse scattering problems to recover the unknown obstacles or inhomogeneous media in the domain of interest with phaseless total fields. Let $u^{i}(x,d)=e^{\mathbf{i}kx\cdot d}$ be the incident plane wave, with wavenumber $k$ and incident direction $d\in\mathbb{S}^{N-1}$ , or $u^{i}(x,x_{s})=G(x,x_{s}),x_{s}\in\mathbb{R}^{N}$ be the incident source wave where $G(x,y)$ is the free-space Green’s function associated with the Helmholtz equation. Assuming that a bounded domain $D\subset\mathbb{R}^{N}$ is the support of obstacles or inhomogeneous media, then the total field $u=u^{i}+u^{s}$ induced by impenetrable obstacles satisfies the following Helmholtz equation

\Delta u+k^{2}u=0\quad\text{in}\quad\mathbb{R}^{N}\backslash D,

(2.1)

\lim_{r\rightarrow\infty}r^{(N-1)/2}(\frac{\partial u^{s}}{\partial r}-\mathbf{i}ku^{s})=0,\quad r=|x|,

(2.2)

with a boundary condition

\left\{\begin{aligned} &u=0\quad\qquad\qquad\text{on}\quad\partial D\quad\text{for sound-soft obstacles},\\ &\frac{\partial u}{\partial\nu}=0\qquad\qquad\text{on}\quad\partial D\quad\text{for sound-hard obstacles },\\ &\frac{\partial u}{\partial\nu}+\mathbf{i}k\lambda u=0\quad\text{on}\quad\partial D\quad\text{with}\quad\lambda\in C(\partial D)\quad\text{for impedance obstacles}.\end{aligned}\right.

(2.3)

For the scattering by inhomogeneous media, the total field satisfies

\Delta u+k^{2}n(x)u=0\quad\text{in}\quad\mathbb{R}^{N}\,,

(2.4)

with the same radiation (2.2), where $n(x)$ is the refractive index. For the scattered field $u^{s}$ , it holds that [18]

u^{s}(x)=\frac{\exp(\mathbf{i}k|x|)}{|x|^{(N-1)/2}}(u^{\infty}(\hat{x})+\mathcal{O}(1/|x|)),\quad|x|\rightarrow\infty\,,

(2.5)

which holds uniformly for all $\hat{x}=x/|x|\in\mathbb{S}^{N-1}$ , where $u^{\infty}$ is called the far-field pattern of $u^{s}$ .

The inverse scattering problem of our interest is to recover the boundary of impenetrable obstacles or the refractive index $n(x)$ of the inhomogeneous media from some noisy measurements of the phaseless total field $|u|=|u^{s}+u^{i}|$ , corresponding to the incident field $u^{i}$ .

3 A direct sampling method for phaseless imaging

3.1 Direct sampling methods for phased data

Consider a sampling domain $\Omega$ in $\mathbb{R}^{N}$ , and a measurement surface $\Gamma_{r}$ , which is assumed to be the circle or ball of radius $R_{r}$ in the sequel. We now first recall the direct sampling method (DSM) for the reconstruction of inhomogeneous medium inclusions and then extend it for the reconstruction of impenetrable scatterers.

Reconstruction of inhomogeneous medium inclusions. The DSM was proposed in [24] for the inverse medium scattering problem with phased data from one incident wave:

\begin{split}I_{\text{DSM}}(z):=|\langle G(z,\cdot),u^{s}\rangle_{L^{2}(\Gamma_{r})}|=\bigg{|}\int_{\Gamma_{r}}G(z,x_{r})\overline{u^{s}(x_{r})}ds(x_{r})\bigg{|},\quad z\in\Omega.\end{split}

(3.1)

The inhomogeneous medium is recovered based on the values of the index function: if it attains an extreme value at a point $z$ , the point lies most likely within the support of the inhomogeneous inclusion, whereas if the index function takes small values, the point $z$ likely lies outside the support. The DSM is shown to be computationally efficient, easy to implement, and does not require any a priori knowledge of the unknown inclusions. For the reconstruction using the far-field pattern, the index function in DSM is defined as [30]

I_{DSM}^{\infty}(z):=|\langle G^{\infty}(z,\cdot),u^{\infty}\rangle_{L^{2}(\mathbb{S}^{N-1})}|,\quad z\in\Omega\,,

(3.2)

where $G^{\infty}(z,\cdot)$ is the far-field pattern associated with the fundamental solution $G(z,x_{r})$ , given by

G^{\infty}(z,\hat{x}_{r})=\left\{\begin{aligned} &\frac{\exp(\mathbf{i}\pi/4)}{\sqrt{8k\pi}}\exp(-\mathbf{i}k\hat{x}_{r}\cdot z),\quad&N=2,\\ &\frac{1}{4\pi}\exp(-\mathbf{i}k\hat{x}_{r}\cdot z),\quad&N=3.\end{aligned}\right.

(3.3)

And we have

G(z,x_{r})=\dfrac{e^{\mathbf{i}kR_{r}}}{R_{r}^{(N-1)/2}}\bigg{\{}G^{\infty}(z,\hat{x}_{r})+\mathcal{O}(R_{r}^{-1})\bigg{\}}.

(3.4)

By combining (2.5) (3.4) and (3.1), we see the asymptotic behavior of the index function $I_{DSM}(z)$ :

I_{DSM}(z)=C_{N}\bigg{|}\int_{\mathbb{S}^{N-1}}G^{\infty}(z,\hat{x})\overline{u^{\infty}(\hat{x},d)}ds(\hat{x})\bigg{|}+\mathcal{O}(R_{r}^{-1})=C_{N}I_{DSM}^{\infty}(z)+\mathcal{O}(R_{r}^{-1})\,,

(3.5)

where $I_{DSM}^{\infty}(z)$ is independent of $R_{r}$ and $C_{N}$ depends only on the dimension $N$ .

Reconstruction of impenetrable scatterers. We now extend the above index function to help reconstruct impenetrable scatterers. This may be derived by using the fact that an impenetrable scatterer can be considered as the limit of the medium scatterer with vanishing or singular material properties [30]. Alternatively, we provide a simple and new derivation below by using the representation formula of the scattered wave.

For the sound-soft obstacle $D$ , the scattered field satisfies the Huygens’ principle [18]

u^{s}(x)=-\mathcal{S}_{D}\bigg{(}\frac{\partial u}{\partial\nu}\bigg{)}(x):=-\int_{\partial D}\frac{\partial u}{\partial\nu}(y)G(x,y)dy\approx\sum_{j}\omega_{j}G(x,y_{j}),\quad x\in\mathbb{R}^{N}\backslash D\,,

(3.6)

where $\mathcal{S}_{D}$ refers to the single-layer potential, $\{y_{j}\}$ is a general set of discrete integration points on $\partial D$ . When the radius of the measurement surface $\Gamma_{r}$ is large, we have the following approximation [24]

\int_{\Gamma_{r}}G(x_{r},z)\overline{G(x_{r},y_{j})}ds(x_{r})\approx k^{-1}\Im(G(z,y_{j}))=k^{-1}\begin{cases}\frac{1}{4}J_{0}(k|z-y_{j}|)\,,\quad N=2\,;\\ \frac{1}{4\pi}\frac{\sin(k|z-y_{j}|)}{|z-y_{j}|}\,,\quad N=3\,\end{cases}

(3.7)

for $z\in\Omega$ , $y_{j}\in\partial D.$ Then by combining the above equation with (3.1) and (3.6), we derive an approximate relation for the sound-soft obstacles:

\int_{\Gamma_{r}}\overline{G(z,x_{r})}u^{s}(x_{r})ds(x_{r})\approx k^{-1}\sum_{j}\omega_{j}\Im(G(y_{j},z))\,.

(3.8)

We notice that $\Im(G(z,y_{j}))$ ) takes a relatively large value if a point $z$ is close to some boundary point $y_{j}$ and decays quickly if $z$ moves away from the scatterers. This observation provides a heuristic justification of the DSM for sound-soft obstacles.

For sound-hard and impedance obstacles, similar arguments to the sound-soft obstacles can be made by writing the scattered field as a single layer potential [17]:

u^{s}(x)=\int_{\partial D}\phi(y)G(x,y)dy,\quad x\in\mathbb{R}^{N}\backslash D,

(3.9)

where the density $\phi\in C(\partial D)$ is a solution to the integral equation

\phi-\mathcal{K}^{T}_{D}\phi-\mathbf{i}k\lambda\mathcal{S}_{D}\phi=2\frac{\partial u^{i}}{\partial\nu}+2\mathbf{i}k\lambda u^{i},

(3.10)

with $\lambda=0$ for sound-hard obstacles, where $\mathcal{K}_{D}$ is the double-layer potential and $\mathcal{K}^{T}_{D}$ is the the dual of $\mathcal{K}_{D}$ with respect to the dual form $\langle\cdot,\cdot\rangle_{\partial D}$ . By conducting the same approximation as the case for sound-soft obstacles, we can verify the DSM for sound-hard and impedance obstacles.

3.2 Formulation and verification of DSM for phaseless imaging

We are now ready to develop a new DSM for reconstruction with phaseless data. Motivated by the corrected data adopted in [12], we propose the following novel index function to reconstruct obstacles or inhomogeneous media with phaseless total field data

I_{\text{DSM}}^{\text{phaseless}}(z):=\bigg{|}\int_{\Gamma_{r}}G(z,x_{r})\Delta(x_{r},d)ds(x_{r})\bigg{|},\quad\forall\,z\in\Omega,

(3.11)

where $\Delta(x_{r},d)$ is given by

\Delta(x_{r},d)=\frac{|u(x_{r},d)|^{2}-|u^{i}(x_{r},d)|^{2}}{u^{i}(x_{r},d)}.

(3.12)

The field $\Delta(x_{r},d)$ above is called the corrected phaseless data. We remark that the work in [12] concerns the phaseless imaging by making use of the reverse time migration with a large number of point source incident fields in $\mathbb{R}^{2}$ . Very different from [12], we focus in this work on some special important physical scenarios, namely when highly limited incident fields are available, e.g., either a few incident source waves or a few incident plane waves for $\mathbb{R}^{N},N=2,3$ . More importantly, we shall combine the proposed DSM and the deep learning strategy to achieve high-quality imaging reconstructions, even for scatterers with complicated geometric shapes and mixed types (i.e., penetrable and impenetrable inhomogeneous inclusions coexist).

The motivation for using the corrected phaseless data is that the leading order term in $I_{\text{DSM}}^{\text{phaseless}}(z)$ will be the same as $I_{\text{DSM}}(z)$ when the measurement surface is far away from $D$ , as shown below.

We now make a resolution analysis on the index function $I_{\text{DSM}}^{\text{phaseless}}$ defined in (3.11) and will show that

I_{\text{DSM}}^{\text{phaseless}}(z)=I_{\text{DSM}}(z)+\mathcal{O}(R_{r}^{(1-N)/2})\quad\text{in}\quad\mathbb{R}^{N}\,.

(3.13)

To do so, we first rewrite $\Delta(x_{r},d)$ as

\Delta(x_{r},d)=\overline{u^{s}(x_{r},d)}+\dfrac{|u^{s}(x_{r},d)|^{2}}{u^{i}(x_{r},d)}+\dfrac{u^{s}(x_{r},d)\overline{u^{i}(x_{r},d)}}{u^{i}(x_{r},d)}\,,

(3.14)

using this relation, we can rewrite $I_{\text{DSM}}^{\text{phaseless}}(z)$ for any $z\in\Omega$ as

\begin{split}I_{\text{DSM}}^{\text{phaseless}}(z)=&\bigg{|}\int_{\Gamma_{r}}G(z,x_{r})\overline{u^{s}(x_{r},d)}ds(x_{r})\\ &+\int_{\Gamma_{r}}G(z,x_{r})\dfrac{|u^{s}(x_{r},d)|^{2}}{u^{i}(x_{r},d)}ds(x_{r})+\int_{\Gamma_{r}}G(z,x_{r})\dfrac{u^{s}(x_{r},d)\overline{u^{i}(x_{r},d)}}{u^{i}(x_{r},d)}ds(x_{r})\bigg{|}\,.\end{split}

(3.15)

Planar incident waves. We first consider the case of the incident plane wave $u^{i}(x_{r},d)=e^{\mathbf{i}kx\cdot d}$ . In the following Lemmas 3.1 to 3.3, we show that the last two terms at the right-hand side of the equation (3.15) are small when $R_{r}\gg 1$ . For simplification, we use $C$ to denote a generic constant independent of $R_{r}$ while its exact values may vary at different occasions. Then we first have an estimate of the second term at the right-hand side of (3.15), which can be derived directly by using the asymptotic formulas (2.5), (3.4) and the fact that $|u^{i}(x_{r},d)|=1$ .

Lemma 3.1.

For any $z\in\Omega$ , we have

\bigg{|}\int_{\Gamma_{r}}G(z,x_{r})\dfrac{|u^{s}(x_{r},d)|^{2}}{u^{i}(x_{r},d)}ds(x_{r})\bigg{|}\leq CR_{r}^{(1-N)/2}\,.

(3.16)

Now we turn to estimating the third term at the right-hand side of (3.15). For this purpose, we need the following estimates of some oscillatory integrals from [12, Lemma 3.9].

Lemma 3.2.

For any $-\infty<a<b<\infty$ and real-valued function $u\in C^{2}[a,b]$ thats satisfies $|u^{\prime}(t)|\geq 1$ for $t\in(a,b)$ . Assume that $a=x_{0}<x_{1}<\cdots<x_{M}=b$ is a division of $(a,b)$ such that $u^{\prime}$ is monotone in each interval $(x_{i-1},x_{i}),i=1,...,M$ . Then it holds for any function $\phi$ defined on $(a,b)$ with integrable derivative and for any $\lambda>0$ that

\bigg{|}\int_{a}^{b}e^{\mathbf{i}\lambda u(t)}\phi(t)dt\bigg{|}\leq(2M+2)\lambda^{-1}\bigg{[}|\phi(b)|+\int_{a}^{b}|\phi^{\prime}(t)|dt\bigg{]}.

(3.17)

With the above Lemma, we can have the following estimate for the last term in (3.15).

Lemma 3.3.

For any $z\in\Omega$ , we have

\bigg{|}\int_{\Gamma_{r}}G(z,x_{r})\dfrac{u^{s}(x_{r},d)\overline{u^{i}(x_{r},d)}}{u^{i}(x_{r},d)}ds(x_{r})\bigg{|}\leq CR_{r}^{(1-N)/2}.

(3.18)

Proof.

By the asymptotic relations in (2.5), (3.4) and change of variables, we obtain

\bigg{|}\int_{\Gamma_{r}}G(z,x_{r})\dfrac{u^{s}(x_{r},d)\overline{u^{i}(x_{r},d)}}{u^{i}(x_{r},d)}ds(x_{r})\bigg{|}=C\bigg{|}\int_{\mathbb{S}^{N-1}}G^{\infty}(z,\hat{x})u^{\infty}(\hat{x},d)e^{-2\mathbf{i}kR_{r}\hat{x}\cdot d}ds(\hat{x})\bigg{|}+\mathcal{O}(R_{r}^{-1}).

(3.19)

We now estimate the first term in (3.19) for the cases of $N=2$ and $N=3$ separately. For simplicity, we write $\phi(\hat{x})=G^{\infty}(z,\hat{x})u^{\infty}(\hat{x},d)$ . First for $N=2$ , without loss of generality we take $d=(1,0)$ and let $\hat{x}=(\cos\theta,\sin\theta)$ , then we can write

\int_{\mathbb{S}^{1}}\phi(\hat{x})e^{-2\mathbf{i}kR_{r}\hat{x}\cdot d}ds(\hat{x})=\int_{0}^{2\pi}\phi(\hat{x})e^{-2\mathbf{i}kR_{r}\cos\theta}d\theta\,.

(3.20)

Denoting $\delta=(kR_{r})^{-1/2}$ and $Q_{\delta}=\{\theta\in(0,2\pi):|\theta-m\pi|<\delta,m=0,1,2\}$ , it is straightforward to derive

\bigg{|}\int_{Q_{\delta}}\phi(\theta)e^{-2\mathbf{i}kR_{r}\cos\theta}d\theta\bigg{|}\leq CR_{r}^{-1/2}\,.

(3.21)

For $\theta\in(0,2\pi)\backslash Q_{\delta}$ , we can easily see

2(kR_{r})^{1/2}|(\cos\theta)^{\prime}|=2(kR_{r})^{1/2}|\sin\theta|\geq 2(kR_{r})^{1/2}|\sin\delta|\geq 1\,,

(3.22)

and obviously $(\cos\theta)^{\prime}$ is piecewise monotone in $(0,2\pi)$ . Then it follows from Lemma 3.17 that

\bigg{|}\int_{(0,2\pi)\backslash Q_{\delta}}\phi(\theta)e^{-2\mathbf{i}kR_{r}\cos\theta}d\theta\bigg{|}\leq CR^{-1/2}\,.

(3.23)

Next, for $N=3$ , without loss of generality, we assume $d=(0,0,1)$ and $\hat{x}=(\sin\theta\cos\varphi,\sin\theta\sin\varphi,\cos\theta)$ . Then we can write

\int_{\mathbb{S}^{2}}\phi(\hat{x})e^{-2\mathbf{i}kR_{r}\hat{x}\cdot d}ds(\hat{x})=\int_{0}^{2\pi}\int_{0}^{\pi}\phi(\theta,\varphi)e^{-2\mathbf{i}kR_{r}\cos\theta}\sin\theta d\theta d\varphi.

(3.24)

Since $(\frac{1}{2\mathbf{i}kR_{r}}e^{-2\mathbf{i}kR_{r}\cos\theta})^{\prime}=-e^{-2\mathbf{i}kR_{r}\cos\theta}\sin\theta$ , by integration by parts, we obtain that

\int_{0}^{\pi}\phi(\theta,\varphi)e^{-2\mathbf{i}kR_{r}\cos\theta}\sin\theta d\theta=-\dfrac{1}{2\mathbf{i}kR_{r}}\bigg{[}e^{-2\mathbf{i}kR_{r}\cos\theta}\phi(\theta,\varphi)\bigg{]}_{0}^{\pi}+\dfrac{1}{2\mathbf{i}kR_{r}}\int_{0}^{\pi}\frac{\partial\phi(\theta,\varphi)}{\partial\theta}e^{-2\mathbf{i}kR_{r}\cos\theta}d\theta.

(3.25)

The proof is then completed by observing that $\phi$ and $\frac{\partial\phi}{\partial\theta}$ are bounded and independent of $R_{r}$ . ∎

Finally, we summarize the main results in this section by combining the above lemmas and (3.15). These results indicate that when $R_{r}$ is large enough, the index function with phaseless data can reconstruct the scatterers accurately and stably as $I_{DSM}$ for the inverse scattering problems with phased data.

Theorem 1.

With the incident plane wave $u^{i}(x,d)=e^{\mathbf{i}kx\cdot d}$ and $z\in\Omega$ , it holds that

I_{\text{DSM}}(z)=C_{N}I_{DSM}^{\infty}(z)+\mathcal{O}(R_{r}^{-1})\,,

(3.26)

I_{\text{DSM}}^{\text{phaseless}}(z)=I_{\text{DSM}}(z)+\mathcal{O}(R_{r}^{(1-N)/2})=C_{N}I_{\text{DSM}}^{\infty}(z)+\mathcal{O}(R_{r}^{(1-N)/2})\,.

(3.27)

Incident point source waves. We now consider the case of incident point source waves $u^{i}(x,x_{s})=G(x,x_{s})$ for $x_{s}\in\mathbb{R}^{N}$ with $R_{s}=|x_{s}|$ as the radius of the surface where point sources are located. Without loss of generality, we assume $R_{r}=\tau R_{s}$ for $\tau$ $\geq$ 1 where $R_{r}$ is the radius of the measurement surface. For the forward scattering problems, it is known that the solution to the scattering problem and its derivatives depend continuously on the incident field [18, Theorems 8.7, 3.11, 3.12, 3.16]:

$\displaystyle\\|u^{\infty}(\hat{x},x_{s})\\|_{\infty,\mathbb{S}^{N-1}}$	$\displaystyle\leq C\\|u^{i}\\|_{\infty,\Omega}\leq CR_{s}^{(1-N)/2}\,,$	(3.28)
$\displaystyle u^{s}(x_{r},x_{s})$	$\displaystyle=\dfrac{e^{\mathbf{i}kR_{r}}}{R_{r}^{(N-1)/2}}\bigg{\{}u^{\infty}(\hat{x}_{r},x_{s})+\mathcal{O}(R_{s}^{(1-N)/2}R_{r}^{-1})\bigg{\}}\,,$	(3.29)
$\displaystyle\|u^{s}(x_{r},x_{s})\|$	$\displaystyle\leq C(R_{s}R_{r})^{(1-N)/2}\,.$	(3.30)

The above estimates will be used repeatedly in our subsequent analysis. To further the resolution analysis for the phaseless index function with incident source waves, we first show a helpful estimate.

Lemma 3.4.

The following estimate holds for the source incident wave $u^{i}(x,x_{s})=G(x,x_{s})$ :

\bigg{|}\int_{\Gamma_{r}}G(z,x_{r})\dfrac{|u^{s}(x_{r},x_{s})|^{2}}{u^{i}(x_{r},x_{s})}ds(x_{r})\bigg{|}\leq CR_{s}^{1-N},\quad z\in\Omega\,.

(3.31)

Proof.

By (3.4) and (3.30), we have

\bigg{|}\int_{\Gamma_{r}}G(z,x_{r})\dfrac{|u^{s}(x_{r},x_{s})|^{2}}{u^{i}(x_{r},x_{s})}ds(x_{r})\bigg{|}\leq CR_{s}^{1-N}R_{r}^{3(1-N)/2}\int_{\Gamma_{r}}\dfrac{1}{|u^{i}(x_{r},x_{s})|}ds(x_{r})\,.

(3.32)

Next, we consider the case $N=2$ and $N=3$ separately. Firstly, when $N=2$ , since $|H_{0}^{1}(t)|^{-1}$ is bound in $(0,1/2)$ and $t^{-\frac{1}{2}}|H_{0}^{1}(t)|^{-1}$ is decreasing for $t>0$ , there exist constants $C_{1}$ and $C_{2}$ such that

|H_{0}^{(1)}(t)|^{-1}\leq C_{1}t^{1/2}+C_{2}\,.

(3.33)

Hence, for sufficiently large $R_{r}$ , we have

\int_{\Gamma_{r}}\dfrac{1}{|u^{i}(x_{r},x_{s})|}ds(x_{r})\leq CR_{r}^{3/2}\,

(3.34)

and

\bigg{|}\int_{\Gamma_{r}}G(z,x_{r})\dfrac{|u^{s}(x_{r},x_{s})|^{2}}{u^{i}(x_{r},x_{s})}ds(x_{r})\bigg{|}\leq CR_{s}^{-1}\,.

(3.35)

Finally, for $N=3$ , we can see that $|u^{i}(x_{r},x_{s})|^{-1}=4\pi|x_{r}-x_{s}|\leq CR_{r}$ . Then we can similarly obtain

\bigg{|}\int_{\Gamma_{r}}G(z,x_{r})\dfrac{|u^{s}(x_{r},x_{s})|^{2}}{u^{i}(x_{r},x_{s})}ds(x_{r})\bigg{|}\leq CR_{s}^{-2}.

(3.36)

∎

Next, we derive some results that are similar to those in Lemma 3.3, but for the incident source waves. For this purpose, we first recall the following important mixed reciprocity relation[18, 37]

Lemma 3.5.

Denote $u^{\infty}(d,x_{s})$ as the far field pattern of the scattering problem with incident wave $u^{i}(x)=G(x,x_{s})$ , and $w^{s}(x,-d)$ as the scattering solution with incident wave $w^{i}(x)=e^{-\mathbf{i}kx\cdot d}$ , then we have

u^{\infty}(d,x_{s})=\gamma_{N}w^{s}(x_{s},-d),\quad x_{s}\in\mathbb{R}^{N}\backslash D,d\in\mathbb{S}^{N-1},

(3.37)

where $\gamma_{2}=\frac{1}{4\pi}$ and $\gamma_{3}=\frac{e^{\mathbf{i}\pi/4}}{\sqrt{8k\pi}}$ .

With Lemma 3.5, equations (2.5) and (3.29), we then have the following estimate for $u^{s}(x_{r},x_{s})$ :

u^{s}(x_{r},x_{s})=\frac{\gamma_{N}e^{\mathbf{i}k(R_{r}+R_{s})}}{R_{r}^{(N-1)/2}R_{s}^{(N-1)/2}}w^{\infty}(\hat{x}_{s},-\hat{x}_{r})+\mathcal{O}(R_{s}^{(1-N)/2}R_{r}^{(N+1)/2}+R_{s}^{(1+N)/2}R_{r}^{(N-1)/2}),

(3.38)

where $w^{\infty}(\hat{x}_{s},-\hat{x}_{r})$ is the far field pattern corresponding to the incident wave $w^{i}(x)=e^{-\mathbf{i}kx\cdot\hat{x}_{r}}$ .

Lemma 3.6.

With source incident wave $u^{i}(x,x_{s})=G(x,x_{s})$ , we have

\bigg{|}\int_{\Gamma_{r}}G(z,x_{r})\dfrac{u^{s}(x_{r},x_{s})\overline{u^{i}(x_{r},x_{s})}}{u^{i}(x_{r},x_{s})}ds(x_{r})\bigg{|}\leq CR_{s}^{(1-N)}\,.

(3.39)

Proof.

We carry out the proof in two separate cases, $N=2$ and $N=3$ .

First for $N=2$ , without loss of generality, we assume $x_{s}=(R_{s},0)$ . Denoting $\delta=(kR_{s})^{-1/2}$ , $\Theta_{\delta}=\{\theta\in(0,2\pi):|\theta\pm m\pi|\leq\delta,m=0,1,2\}$ and $Q_{\delta}=\{x_{r}\in\Gamma_{r}:|\theta_{r}\in\Theta_{\delta}\}$ , then it follows easily from (3.4), (3.30) and $|Q_{\delta}|\leq CR_{r}R_{s}^{-1/2}$ that

\bigg{|}\int_{Q_{\delta}}G(z,x_{r})\dfrac{u^{s}(x_{r},x_{s})\overline{u^{i}(x_{r},x_{s})}}{u^{i}(x_{r},x_{s})}ds(x_{r})\bigg{|}\leq CR_{s}^{-1}\,.

(3.40)

For $x_{r}\in\Gamma_{r}\backslash Q_{\delta}$ , since $\sin t\geq t/2$ for $t\in(0,\pi/2)$ , we have

k|x_{r}-x_{s}|\geq 2kR_{s}\bigg{|}\sin\frac{\theta_{r}-\theta_{s}}{2}\bigg{|}\geq 2kR_{s}\sin\frac{\delta}{2}\geq\frac{1}{2}(kR_{s})^{1/2}\,.

(3.41)

Thus, by the asymptotic behavior of the Hankel functions [18]

H_{n}^{(1)}(t)=\bigg{(}\frac{2}{\pi t}\bigg{)}^{1/2}e^{\mathbf{i}(t-\frac{n\pi}{2}-\frac{\pi}{4})}+R_{n}(t),\quad|R_{n}(t)|\leq Ct^{-3/2},\quad t>0\,,

(3.42)

we can write

\frac{\overline{u^{i}(x_{r},x_{s})}}{u^{i}(x_{r},x_{s})}=e^{-2\mathbf{i}k|x_{r}-x_{s}|+\mathbf{i}\frac{\pi}{2}}+\rho_{1}(x_{r},x_{s})\,,

(3.43)

where $|\rho_{1}(x_{r},x_{s})|\leq CR_{s}^{-1/2}$ . Then, by combining (3.4), (3.43), and (3.38), we obtain

\begin{split}&\bigg{|}\int_{\Gamma_{r}\backslash Q_{\delta}}G(z,x_{r})\dfrac{u^{s}(x_{r},x_{s})\overline{u^{i}(x_{r},x_{s})}}{u^{i}(x_{r},x_{s})}ds(x_{r})\bigg{|}\\ \leq&CR_{s}^{-1/2}\bigg{|}\int_{(0,2\pi)\backslash\Theta_{\delta}}G^{\infty}(\hat{x}_{r},z)w^{\infty}(\hat{x}_{s},-\hat{x}_{r})e^{-2\mathbf{i}k|x_{r}-x_{s}|}d\theta_{r}\bigg{|}+\mathcal{O}(R_{s}^{-1})\\ =&CR_{s}^{-1/2}\bigg{|}\int_{(0,2\pi)\backslash\Theta_{\delta}}G^{\infty}(\hat{x}_{r},z)w^{\infty}(\hat{x}_{s},-\hat{x}_{r})e^{2\mathbf{i}kR_{s}v(\theta_{r})}d\theta_{r}\bigg{|}+\mathcal{O}(R_{s}^{-1})\,,\end{split}

(3.44)

where $v(\theta_{r})=-\sqrt{1+\tau^{2}-2\tau\cos\theta_{r}}$ . For small $\delta$ , The derivative of $v$ can be computed explicitly as

|v^{\prime}(\theta_{r})|\geq\tau|\sin\delta|/|v(\theta_{r})|\geq\frac{\tau}{1+\tau}\frac{\delta}{2}\geq\frac{1}{4}(kR_{s})^{-1/2}\,.

(3.45)

Employing the above two estimates and utilizing Lemma (3.17), we have

\bigg{|}\int_{\Gamma_{r}\backslash Q_{\delta}}G(z,x_{r})\dfrac{u^{s}(x_{r},x_{s})\overline{u^{i}(x_{r},x_{s})}}{u^{i}(x_{r},x_{s})}ds(x_{r})\bigg{|}\leq CR_{s}^{-1}\,.

(3.46)

Now we analyze for the case $N=3$ . We first see from (3.4) and (3.38) that

\begin{split}&\bigg{|}\int_{\Gamma_{r}}G(z,x_{r})\dfrac{u^{s}(x_{r},x_{s})\overline{u^{i}(x_{r},x_{s})}}{u^{i}(x_{r},x_{s})}ds(x_{r})\bigg{|}\\ =&CR_{s}^{-1}\bigg{|}\int_{\Gamma_{r}}G^{\infty}(\hat{x}_{r},z)w^{\infty}(\hat{x}_{s},-\hat{x}_{r})e^{-2\mathbf{i}k|x_{r}-x_{s}|}d(x_{r})\bigg{|}+\mathcal{O}(R_{s}^{-2}).\end{split}

(3.47)

Without loss of generality, we assume $x_{s}=(0,0,R_{s})$ and let $x_{r}=R_{r}(\sin\theta\cos\varphi,\sin\theta\sin\varphi,\cos\theta)$ , $\phi(\hat{x}_{r})=G^{\infty}(\hat{x}_{r},z)w^{\infty}(\hat{x}_{s},-\hat{x}_{r})$ , we can then rewrite the leading order term in (3.47) as

\int_{\Gamma_{r}}\phi(\hat{x}_{r})e^{-2\mathbf{i}k|x_{r}-x_{s}|}d(x_{r})=\int_{0}^{2\pi}\int_{0}^{\pi}\phi(\hat{x}_{r})e^{-2\mathbf{i}kR_{s}\sqrt{\tau^{2}+1-2\tau\cos\theta}}\sin\theta d\theta d\varphi\,.

(3.48)

Denoting $v(\theta)=\sqrt{\tau^{2}+1-2\tau\cos\theta}$ and knowing $v^{\prime}(\theta)=\tau\sin\theta/v(\theta)$ , by integration by parts, we obtain that

\int_{0}^{\pi}\phi(\hat{x}_{r})e^{-2\mathbf{i}kR_{s}v(\theta)}\sin\theta d\theta=-\dfrac{1}{2\mathbf{i}k\tau R_{s}}\bigg{(}\bigg{[}\phi(\hat{x}_{r})v(\theta)e^{-2\mathbf{i}kR_{s}v(\theta)}\bigg{]}_{\theta=0}^{\theta=\pi}-\int_{0}^{\pi}e^{-2\mathbf{i}kR_{s}v(\theta)}\frac{\partial(\phi(\hat{x}_{r})v(\theta))}{\partial\theta}d\theta\bigg{)}\,.

We then have

\bigg{|}\int_{\Gamma_{r}}G^{\infty}(\hat{x}_{r},z)w^{\infty}(\hat{x}_{s},-\hat{x}_{r})e^{-2\mathbf{i}k|x_{r}-x_{s}|}d(x_{r})\bigg{|}\leq CR_{s}^{-1}\,.

(3.49)

Then the desired estimate follows from this and (3.47). ∎

By substituting the estimates in Lemmas 3.4 and 3.6 into the decomposition of the index function in (3.15), we come to the following conclusions, where the second result comes from replacing $u^{i}(x,d)$ with $u^{i}(x,x_{s})$ in (3.11) and (3.12).

Theorem 2.

With the incident source wave $u^{i}(x,x_{s})=G(x,x_{s})$ , and assuming $R_{r}=\tau R_{s}$ for $\tau$ $\geq$ 1, we then have for $z\in\Omega$ :

I_{\text{DSM}}(z)=C_{N}I_{DSM}^{\infty}(z)+\mathcal{O}(R_{s}^{(1-N)/2}R_{r}^{-1})=\mathcal{O}(R_{s}^{(1-N)/2})\,.

(3.50)

and

I_{\text{DSM}}^{\text{phaseless}}(z)=I_{\text{DSM}}(z)+\mathcal{O}(R_{s}^{1-N})=C_{N}I_{DSM}^{\infty}(z)+\mathcal{O}(R_{s}^{1-N})\,.

(3.51)

4 The direct sampling-based deep learning approach

In this section, we present a direct sampling-based deep learning approach (DSM-DL) that combines DSM for phaseless data, as described in section 3, with deep learning techniques to perform phaseless reconstruction for both penetrable and impenetrable scatterers, as well as their mixed problems. We apply deep learning techniques to address a common trade-off in classical methods for solving inverse scattering problems, which often involves a compromise between computation time and accuracy. Deep learning has the potential to overcome this limitation by leveraging large datasets and optimization algorithms during the training process. A well-trained neural network can offer rapid, stable, and highly accurate image reconstructions. Additionally, a key advantage of the deep learning approach is its ability to alleviate the need for prior knowledge about unknown scatterers, as it can learn the distribution information from training data.

The DSM-DL consists of two steps: we first train a neural network, denoted as $\mathcal{G}_{\Theta}$ , to learn the mapping between the index functions $\{I_{\text{DSM},i}^{\text{phaseless}}\}_{i=1}^{N_{i}}$ and the medium profile $n(x)$ , where $N_{i}$ is the number of incidences. This mapping is relatively easier to learn compared to the highly nonlinear inverse mapping from phaseless scattered data to $n(x)$ because $\{I_{\text{DSM},i}^{\text{phaseless}}\}_{i=1}^{N_{i}}$ provides a stable and rough estimate for $n(x)$ . Furthermore, both the index function and $n(x)$ may be regarded to be in the same function space, e.g., $L^{2}(\mathbb{R}^{n})$ .

Subsequently, to solve inverse problems, we will first compute $\{I_{\text{DSM},i}^{\text{phaseless}}\}_{i=1}^{N_{i}}$ by (3.11) and feed them into the neural network $\mathcal{G}_{\Theta}$ that we obtained above to have an image $\tilde{n}(x)=\mathcal{G}_{\Theta}(\{I_{\text{DSM},i}^{\text{phaseless}}\}_{i=1}^{N_{i}})$ . The detailed algorithm for training is summarized in Algorithms 1 .

Input:

\bullet

Given fixed

N_{i}

incident fields

\{u_{p}^{i}\}_{p=1}^{N_{i}}

\bullet

N

true refractive functions

\{n_{j}(x)\}_{j=1}^{N}

(for impenetrable scatterer

D

, we set

n(x)=0

for point

x

inside the scatterer

D

) defined in a selected sampling domain

\Omega

and their corresponding phaseless total-field data

\{|u|_{p,j},p=1,\cdots,N_{i}\}_{j=1}^{N}

measured on a surface

\Gamma_{r}

\bullet

The free-space Green’s function

\{G(x,y),x\in\Omega,y\in\Gamma_{r}\}

\bullet

Number of epochs:

L

; Batch number:

Q

; Learning rates:

\{\tau_{q}\}_{q=1}^{L}

1 Compute the index functions

I_{p,j}^{Phaseless}(z)=\bigg{|}\int_{\Gamma_{r}}G(z,x_{r})\Delta(x_{r},d)ds(x_{r})\bigg{|},\quad\forall\,z\in\Omega;j=1,\cdots,N;p=1,\cdots,N_{i}

. where

\Delta(x_{r},d)

is given by

\Delta(x_{r},d)=\frac{|u(x_{r},d)|^{2}_{p,j}-|u^{i}(x_{r},d)|^{2}}{u^{i}(x_{r},d)}.

2Introduce a loss function

Loss

such as (4.1).

3 Construct and initialize a network

\mathcal{G}_{\Theta}

, where

\Theta

denotes the parameters in the neural network.

4 for $q=1,2,\cdots,L$ do

5 for $l=1,2,\cdots,Q$ do

6 Update

\Theta\leftarrow\Theta-\frac{1}{|\mathcal{I}_{l}|}\sum_{j\in\mathcal{I}_{l}}\tau_{q}\nabla_{\Theta}Loss\big{(}\mathcal{G}_{\Theta}(\{I_{p,j}^{Phaseless}\}_{p=1}^{N_{i}}),n_{r,j}\big{)}

\mathcal{I}_{l}

refers to the index set for the

l_{th}

batch,

\tau_{q}

refers to learning rate for the

q_{th}

epoch. */

8 end for

10 end for

Algorithm 1 DSM-DL for Phaseless Data (Training)

Refer to caption — Figure 1: The architecture of the U-Net used in the numerical experiments.

In our DSM-DL approach, we employ the U-Net architecture [39] as the neural network, as shown in Fig. 1. The U-Net architecture is suitable for our task because it can handle multi-channel input, which in our case corresponds to the index functions. The U-Net follows a symmetric U-shaped structure. The left-hand side, known as the contracting path, involves repeated convolutions, batch normalization, rectified linear unit (ReLU), and max pooling operations. And the right-hand side of the U-Net forms the expansive path, which is similar to the contracting path but employs $3\times 3$ up-convolutions instead of $2\times 2$ max pooling. Skip connections are used to recover lost spatial information during down-sampling by transferring features from the contracting path to the expansive path. Since the index functions have larger values near the scatterers and relatively smaller values away from the scatterers, the relationship between the index functions and the true contrast images is likely to be mainly local which is rather consistent with the local nature of CNNs.

While it is possible to employ extended versions of the classical U-Net to potentially enhance the results, it is important to note this is an endless pursuit of improvement and is not the primary focus of this work. Therefore, we chose the standard U-Net architecture for our neural networks, which proved to be able to provide quite satisfactory reconstructions, as shown by our numerical experiments in the sequel.

We propose the following loss function to generate an effective network for numerical reconstructions:

Loss=\frac{1}{T_{n}}\sum_{i=1}^{T_{n}}\bigg{(}\|X_{i}-Y_{i}\|^{2}_{2}+\alpha_{1}\text{TV}(X_{i})+\alpha_{2}(1-\text{SSIM}(X_{i},Y_{i})\bigg{)}\,,

(4.1)

where $T_{n}$ is the total number of the training data, $X_{i}$ and $Y_{i}$ are the output of the neural network and the true image, respectively, while $\alpha_{1}$ and $\alpha_{2}$ are the regularization parameter and weight coefficient. The first term in (4.1) is a mean square error to ensure the accuracy of reconstruction, the second term is a total variation regularization term, with $\text{TV}(u):=\int_{\Omega}|\nabla u|dx$ , that is added to detect edges of images more efficiently which is critical for many practical situations as the boundary information of scatters is usually most demanding. The third term is a structural similarity function (SSIM) [42] to help enhance the accuracy of the reconstruction regarding its luminance and contrast of true and reconstructed images.

We now list several attractive advantages of the DSM-DL approach for solving inverse scattering problems.

•

Firstly, once the neural network is trained, the DSM-DL method requires only inexpensive computations of the index functions and a forward pass through the network, enabling real-time estimation of the unknown scatterers.
•

The DSM-DL method is very robust to large noise, as the noise is generally significantly smoothed out during the DSM stage.
•

Moreover, the DSM-DL method does not require a prior knowledge of the physical properties of the scatterers. This feature enables the applications of the DSM-DL to the mixed problem, where both medium scatterers and impenetrable scatterers exist within the domain of interest.
•

Importantly, Theorems 3.27 and 2 demonstrate that with sufficient large $R_{r}$ a DSM-DL network trained using phased data can be also applied to problems with only phaseless data available which will be illustrated in numerical experiments in the next section.
•

Additionally, numerical demonstrations confirm that the DSM-DL method possesses a high potential for generalization, allowing it to solve problems that differ significantly from the training data, which is crucial in the applications of deep learning to general inverse scattering problems.
•

The DSM-DL can generally do a much better reconstruction than the classical DSM , namely the latter can only recover the locations and shapes of the scatterers, while the former is also able to recover the refractive index values of the scatterers.

5 Numerical Experiments

We present several representative numerical examples to illustrate the performance of the DSM and DSM-DL for reconstructing the unknown scatterers from noisy phaseless data. The wavelength is chosen as $\lambda=0.75$ and the wavenumber $k$ will be $2\pi/\lambda\approx 8.37$ . The noisy data $|u_{\delta}|$ are generated point-wisely by the formula

|u_{\delta}(x)|=|u(x)|+\delta\zeta(x)\|u\|_{2},\quad x\in\Gamma_{r},

(5.1)

where $\delta$ refers to the relative noise level, and $\zeta(x)$ follows the standard normal distribution.

We take a sampling domain to be $\Omega=[-1,1]^{2}$ , which may contain both inhomogeneous media and sound-soft scatterers. Moreover, to reconstruct impenetrable scatterers with DSM-DL, a fundamental question is a proper representation of scatterers using pixel values. Since the scattering property of impenetrable scatterers depends solely on their boundaries, we can assign a pixel value of $0$ to points inside the scatterers and $1$ to points outside the scatterers to represent them. Remarkably, our numerical experiments demonstrate that a single trained DSM-DL network can simultaneously solve both inverse medium and obstacle scattering problems, as well as their mixed problems (see section 5.2.3). We first present the numerical results for the new DSM with phaseless data in Section 5.1, then the numerical results for DSM-DL in Section 5.2.

5.1 Numerical results for phaseless reconstruction by DSM

For phaseless reconstruction by DSM, we are interested in the important case when very limited data is available, e.g., only one or a few incident fields. Unless otherwise specified, one incident plane wave with incident direction $d=(\cos(\frac{\pi}{4}),\sin(\frac{\pi}{4}))$ is employed, and the data $|u|$ is measured at 100 points uniformly distributed on a circle of radius 4. We normalize the index function $I(z)$ by $\hat{I}(z)=I(z)/\max_{z}I(z)$ so that its maximum value is 1.

Example 1. In this example, we consider one square scatterer of width 0.15 located at the origin. The coefficient $n(x)$ of the scatterer is 3.

The numerical reconstructions are presented in Fig. 2. We observe that for sampling points near the scatterer, the index function is indeed relatively larger than points outside the scatterer. From this example, we can see the index function can provide a reliable and visible indicator for locating the scatterer, even when the noise level is $\delta=5\%$ or $\delta=10\%$ in the phaseless data.

Example 2. In this example, we consider two small square sound-soft scatterers of width 0.15. The two scatterers are located at $(-0.5,0.5)$ and $(0.5,-0.5)$ which are well separated, respectively; see Fig. 3.

From the plots in Fig. 3, we can see that the index function can provide satisfactory reconstructions, with only one incident field, and is quite robust under the noise level $\delta=5\%$ or and $10\%$ .

Example 3. This example considers two scatterers with width $0.15$ located at $(0.1,0.1)$ and $(-0.1,-0.1)$ , respectively. The distance between them is less than one-half of the wavelength ( $\lambda/2$ ), which is known to be very challenging for imaging reconstruction. Nonetheless, the index function can still allow us to distinguish two close scatterers (cf. Fig. 4).

Example 4. In this example, the scatterer is a ring located at the origin, with the outer and the inner boundaries of the ring very close, with the radius being 0.3 and 0.25, respectively.

The ring-shaped scatterer is known to be one of the most challenging objects to recover in inverse scattering problems, especially in the case when the thickness of the ring is very small. From Fig. 5, we notice that employing one single incident wave is only able to recover some parts of the ring structure which is reasonable as the lack of measurement data. We then try to use more incident waves $\{e^{ikx\cdot d_{i}}\}_{i=1}^{N_{i}}$ , where $d_{i}=(\cos(\theta_{i}),\sin(\theta_{i})),\theta_{i}=\frac{2\pi(i-1)}{N_{i}}+\frac{\pi}{4}$ , with the following average index function

I(z)=\dfrac{1}{N_{i}}\sum_{i=1}^{N_{i}}I_{i}(z),

(5.2)

where $I_{i}(z)$ denotes the index function corresponding to the $i$ th incident wave. We investigate $N_{i}=1,3$ and $5$ and observe that with more incident waves the method can produce more accurate and stable reconstructions for this challenging case, as shown in Fig. 5.

5.2 Numerical results of phaseless reconstruction by DSM-DL

In this subsection, we evaluate the performance of the DSM-DL scheme by training the neural networks with three very different datasets: a polygon dataset, a MNIST dataset, and a mixed circle dataset. We discretize the index functions and the true images by $64\times 64$ pixels. We scale all the inputs (index functions) of the neural networks by multiplying a constant $2/W$ in the numerical experiments, where $W$ is the maximum value of the index functions among all training data. The Adam optimizer is used to minimize the loss function (4.1) and $\alpha_{1}=\alpha_{2}=0.5$ which is determined by grid-search as described in section 4. We add $1\%$ Gaussian noise to the phaseless data in the training stage, while in the testing stage, much higher noise levels are considered. To evaluate the quality of the recovered images, we employ the relative L2 error:

R_{e}(X,Y)=\frac{\|X-Y\|_{2}}{\|Y\|_{2}},

(5.3)

where $X$ is the output of the network and $Y$ is the true image. We also employ the following accuracy metric:

Acc(X,Y):=\frac{1}{N_{d_{x}}N_{d_{y}}}\sum_{i=1}^{N_{d_{x}}}\sum_{j=1}^{N_{d_{y}}}a_{i,j},\quad\text{where}\quad a_{i,j}=1\quad\text{if}\quad X_{i,j}=Y_{i,j},\quad\text{else}\quad a_{i,j}=0.

(5.4)

The metric $Acc(X,Y)$ is suitable for impenetrable scatterers as it is independent of the pixel values we choose, and is effective when we are concerned with only the locations of the unknown scatterers.

5.2.1 Training with polygon dataset

In the first example, we use a regular polygon with a random location to simulate either a medium scatterer or a sound-soft obstacle in each image. The number of the sides of each polygon is randomly set between 3 and 6, and the circumradius of each polygon is taken from the uniform distribution $U(0.3,0.5)$ . For medium scatterers, the coefficient $n(x)$ inside the scatterer is set to 3. For the sound-soft obstacles, the pixel values of the points inside the obstacles are set to 0. In addition, we randomly rotate the polygon to enhance the diversity of the data. In this example, We employ 7000 images as the training data and 200 images as the testing data. The batch size is set to 10 and a total of 30 epochs is used in the training. For each batch of training data, there are 5 examples concerning medium scatterers and 5 examples concerning sound-soft obstacles. The learning rate starts at 0.001 and decreases by a factor of 0.5 every 3 epochs. We will employ $N_{i}$ incident plane waves with $N_{i}=1$ or $4$ , and employ $100$ receivers that are equally distributed on a circle of radius 4 centered at the origin.

5.2.1.1 Tests with Polygon testing data

After the training, we use the networks to test the examples from the testing data of the Polygon dataset. We further apply cutoff values 0.5 and 2.0 to the output of neural network so that the pixel values of the images only take from $\{0,1,3\}$ . In Fig. 6, we present reconstructed images with different noise levels and number of incidences for 5 examples, where 3 images are for medium scatterers and 2 images are for sound-soft obstacles. We observe that for $\delta=2\%$ , even with only one single incident wave, the DSM-DL can still exactly identify the type of the scatterer and provide a very accurate estimation for the shapes and locations of the scatterers. As we increase the noise level to $10\%$ , the method with one incidence can still identify the type and location of the scatterers, while the shapes of reconstructed scatterers are distorted, especially for medium scatterers. But when we use multiple incidents, even still a small number, e.g., 4 incident fields, higher-quality reconstructions can be obtained. The average accuracy metric $Acc$ over the testing data is listed in Table. 1. We notice that for $N_{i}=1$ , the accuracy is decreased by $1.77\%$ when the noise level is increased from $2\%$ to $10\%$ , while for $N_{i}=4$ the accuracy is only decreased by $0.61\%$ . This indicates that employing more incidences can make the reconstruction more robust against noise.

5.2.2 Training with MNIST dataset

In this subsection, we employ a modification of the well-known MNIST dataset to model the inhomogeneous medium and train the neural networks. The resolution of each image in the MNIST dataset is $28\times 28$ and the pixel values range from 0 and 1. The resolution is then rescaled to $64\times 64$ in our numerical experiments and a threshold is applied with value 0.3 so that the pixel values only take 0 and 1. To increase the diversity of the data, the digits are randomly rotated, and a circle with a radius from $U(0.2,0.4)$ is added to each image. The coefficient $n(x)$ of the digits and the circles are randomly taken from $U(1.2,1.7)$ . In this example, We employ 10000 images as the training data and 200 images as the testing data. The batch size is set to 10 and a total of 30 epochs is used in the training. The learning rate starts at 0.001 and decreases by a factor of 0.5 every 3 epochs. We will employ $N_{i}$ incident plane waves with $N_{i}=4$ or $16$ , and employ $100$ receivers that are equally distributed on a circle of radius 4 centered at the origin.

Example	$\mathbf{N_{i}=1,\delta=2\%}$	$\mathbf{N_{i}=1,\delta=10\%}$	$\mathbf{N_{i}=4,\delta=2\%}$	$\mathbf{N_{i}=4,\delta=10\%}$
Polygon	$99.49\%$	$97.72\%$	$99.77\%$	$99.16\%$

Table 1: Examples 5.2.1.1: The accuracy

Acc

for the polygon example.

5.2.2.1 Tests with MNIST testing data.

To evaluate the performance of the trained networks, we first test the examples from MNIST testing data. In Fig. 7, reconstructions of five examples from the MNIST testing data with different noise levels and number of incidences are presented. With $N_{i}=4$ , the method can produce some reasonable reconstructions and can distinguish the hand-written digit and the circle while the accuracy is not very high. However, if 16 incidences are employed, high-quality and stable reconstructions can be achieved. This shows that multiple data are important for recovering complicated scatterers and the DSM-DL can fully make use of multiple data. The average relative error for the testing data is presented in Table. 2.

5.2.2.2 Tests with Chinese characters.

To further test the generalization ability of the trained models, we consider the reconstruction of five Chinese characters as shown in Fig. 8, whose shapes are very different from those of the MNIST dataset. The coefficient value $n(x)$ of the scatterers is set to 1.5. The recovered images presented in Fig. 8 also show the ability of the DSM-DL to fully extract the hidden information in the measurement data and thus provide more accurate and stable reconstructions with more incidences used. The average relative error for the five examples is presented in Table. 2

5.2.2.3 Tests with “Austria Rings”.

Finally, we apply the trained neural networks to test two different “Austria rings” as shown in Fig. 9. The “Austria ring” is known to be a challenging profile in inverse scattering problems. For the first “Austria” example, the coefficient $n(x)$ of the scatterers is set to 1.5. While for the second “Austria” example, the coefficient $n(x)$ of the two circles is changed to 2.0, which is out of the range in the training data. When $N_{i}=4$ , the method can only provide very rough reconstructions, and the recovered images are heavily affected by high-level noises. By using $N_{i}=16$ incidences, the method can produce more accurate reconstructions and can distinguish the strong scatterers and the weak scatterer in the second example.

Example	$\mathbf{N_{i}=4,\delta=5\%}$	$\mathbf{N_{i}=4,\delta=10\%}$	$\mathbf{N_{i}=16,\delta=5\%}$	$\mathbf{N_{i}=16,\delta=10\%}$
MNIST	0.0827	$0.1043$	$0.0617$	$0.0755$
Chinese characters	$0.1096$	$0.1252$	$0.0721$	$0.0854$
Austria Ring 1	$0.1163$	$0.1258$	$0.0851$	$0.0922$
Austria Ring 2	$0.1897$	$0.1810$	$0.1260$	$0.1367$

Table 2: Examples 5.2.2.1, 5.2.2.2 and 5.2.2.3: The relative L2 errors for different examples with the networks trained by the MNIST dataset.

5.2.3 Training with a mixed circle dataset by phased data

We now consider a mixed problem in which the domain of interest may consist of both medium scatterers and sound-soft obstacles. To further show the convergence property of the index functions, the network is trained by phased data, i.e, with the index functions(3.1) computed form phased data as the input of network, and then applied to solve problems with phaseless data. This is a very important advantage of the DSM-DL as it is usually difficult to recover phased information from phaseless data. In the training dataset, $1\sim 3$ circles with varying radius drawn from the uniform distribution $U(0.2,0.3)$ are randomly added to the domain to simulate scatterers. Scatterers are not allowed to overlap and each scatterer has an equal probability of being a medium or sound-soft obstacle. The coefficient value of the medium scatterer is randomly set from $U(1.5,3.0)$ , and we still set the pixel value for points inside the sound-soft obstacles to be 0. We will employ 10 incident plane waves and employ $180$ receivers that are equally distributed on a circle of radius 8 centered at the origin to collect both phased data and phaseless data. We use 20000 datasets to train the network and 200 datasets for testing. The batch size is taken as 20 and a total of 30 epochs in the training stage. The learning rate starts at 0.001 and decreases by a factor of 0.5 every 3 epochs.

5.2.3.1 Tests with mixed circle testing data by phaseless data.

In this example, we employ the model trained by phased data to solve 200 mixed circle examples with phaseless data, and compare the results by adding different noises to the measurement data. The reconstructed images of five typical cases are presented in Fig. 10. The shown results not only prove the capability of the DSM-DL for dealing with inverse scattering problems with mixed scatterers in a unified framework, but also show that the DSM-DL network trained by phased data can also be applied to phaseless data. We observe that the method can accurately identify the boundaries, locations, sizes, and coefficient values of the scatterers. In particular, the results are still quite satisfactory when high-level noises are added to the data.

6 Conclusions

We have studied the inverse scattering problems under highly challenging conditions, where the measurement data are phaseless total fields and available only from one or a few incident waves. To address this demanding problem, we introduce a novel technique known as the Direct Sampling Method (DSM) for phaseless data. The DSM offers a robust and reasonable means of estimating the shapes and locations of unknown scatterers. Remarkably, our proposed DSM is applicable to both inverse medium and obstacle problems and does not necessitate prior knowledge of the physical properties of the scatterers.

To enhance the precision of our reconstructions, we present the DSM-DL, a fusion of the DSM with deep learning methodologies, offering high-quality reconstructions. The DSM-DL not only recovers the geometrical attributes and positions of the scatterers but also identifies their physical properties and the coefficient values of the medium scatterers. Moreover, it is capable of addressing mixed problems involving both medium and obstacle scatterers. It is worth noting that, supported by both theoretical analysis and numerical results, the DSM-DL network trained with phased data can also effectively tackle problems with phaseless data which implies a broad application of the algorithm.

In conjunction with our present work, there are several promising directions for future exploration. Firstly, a more rigorous theoretical treatment for DSM on phaseless near-field data would be an interesting direction, where we have observed robust results. Secondly, in the context of DSM, it would be noteworthy to investigate the feasibility of employing a neural network to compute optimal probing functions under specific noise levels. Lastly, our proposed DSM-DL has the potential to address a wide range of other challenging and practically critical inverse problems, including phaseless reconstruction without information about the incident field, limited aperture scattering problems, and electromagnetic inverse problems.

References

[1] A. D. Agaltsov, T. Hohage, and R. G. Novikov. An iterative approach to monochromatic phaseless inverse scattering. Inverse Problems, 35(2):024001, 2018.
[2] H. Ammari, Y. T. Chow, and J. Zou. Phased and phaseless domain reconstructions in the inverse scattering problem via scattering coefficients. SIAM J. Appl. Math., 76(3):1000–1030, 2016.
[3] S. Arridge, P. Maass, O. Öktem, and C.-B. Schönlieb. Solving inverse problems using data-driven models. Acta Numerica, 28:1–174, 2019.
[4] G. Bao, P. Li, J. Lin, and F. Triki. Inverse scattering problems with multi-frequencies. Inverse Problems, 31(9):093001, 2015.
[5] G. Bao, P. Li, and J. Lv. Numerical solution of an inverse diffraction grating problem from phaseless data. J. Opt. Soc. Am. A, 30(3):293–299, 2013.
[6] A. Bulyshev, A. Souvorov, S. Y. Semenov, V. Posukh, and Y. Sizov. Three-dimensional vector microwave tomography: theory and computational experiments. Inverse problems, 20(4):1239, 2004.
[7] F. Cakoni, D. Colton, and P. Monk. The linear sampling method in inverse electromagnetic scattering. SIAM, 2011.
[8] J. Chen, Z. Chen, and G. Huang. Reverse time migration for extended obstacles: acoustic waves. Inverse Problems, 29(8):085005, 2013.
[9] J. Chen, Z. Chen, and G. Huang. Reverse time migration for extended obstacles: electromagnetic waves. Inverse Problems, 29(8):085006, 2013.
[10] X. Chen. Subspace-based optimization method for solving inverse-scattering problems. IEEE Trans. Geosci. Remote Sens., 48(1):42–49, 2009.
[11] X. Chen, Z. Wei, M. Li, and P. Rocca. A review of deep learning approaches for inverse scattering problems. Prog. Electromagn. Res., 167:67–81, 2020.
[12] Z. Chen and G. Huang. Phaseless imaging by reverse time migration: acoustic waves. Numer. Math. Theory Methods Appl., 10(1):1–21, 2017.
[13] Y. T. Chow, F. Han, and J. Zou. A direct sampling method for simultaneously recovering inhomogeneous inclusions of different nature. SIAM J. Sci. Comput., 43(3):A2161–A2189, 2021.
[14] Y. T. Chow, F. Han, and J. Zou. A direct sampling method for the inversion of the radon transform. SIAM J. Imaging Sci., 14(3):1004–1038, 2021.
[15] Y. T. Chow, F. Han, and J. Zou. A direct sampling method for simultaneously recovering electromagnetic inhomogeneous inclusions of different nature. J. Comput. Phys., 470:111584, 2022.
[16] Y. T. Chow, K. Ito, and J. Zou. A direct sampling method for electrical impedance tomography. Inverse Problems, 30(9):095003, 2014.
[17] D. Colton and R. Kress. Integral equation methods in scattering theory. SIAM, 2013.
[18] D. Colton and R. Kress. Inverse Acoustic and Electromagnetic Scattering Theory, volume 93. Springer, 4th edition, 2019.
[19] H. Dong, D. Zhang, and Y. Guo. A reference ball based iterative algorithm for imaging acoustic obstacle from phaseless far-field data. Inverse Probl. Imaging, 13(1):177–195, 2018.
[20] M. d’Urso, K. Belkebir, L. Crocco, T. Isernia, and A. Litman. Phaseless imaging with experimental data: facts and challenges. J. Opt. Soc. Am. A, 25(1):271–281, 2008.
[21] Y. Gao, H. Liu, X. Wang, and K. Zhang. On an artificial neural network for inverse scattering problems. J. Comput. Phys., 448:110771, 2022.
[22] R. Guo, T. Huang, M. Li, H. Zhang, and Y. C. Eldar. Physics-embedded machine learning for electromagnetic data imaging: Examining three types of data-driven imaging methods. IEEE Signal Process. Mag., 40(2):18–31, 2023.
[23] J. Han, A. Jentzen, and W. E. Solving high-dimensional partial differential equations using deep learning. Proc. Natl Acad. Sci., 115(34):8505–8510, 2018.
[24] K. Ito, B. Jin, and J. Zou. A direct sampling method to an inverse medium scattering problem. Inverse Problems, 28(2):025003, 2012.
[25] K. Ito, B. Jin, and J. Zou. A direct sampling method for inverse electromagnetic medium scattering. Inverse Problems, 29(9):095018, 2013.
[26] O. Ivanyshyn and R. Kress. Inverse scattering for surface impedance from phase-less far field data. J. Comput. Phys., 230(9):3443–3452, 2011.
[27] Y. Khoo and L. Ying. SwitchNet: a neural network model for forward and inverse scattering problems. SIAM J. Sci. Comput., 41(5):A3182–A3201, 2019.
[28] M. V. Klibanov. Phaseless inverse scattering problems in three dimensions. SIAM J. Appl. Math., 74(2):392–410, 2014.
[29] M. V. Klibanov and V. G. Romanov. Reconstruction procedures for two inverse scattering problems without the phase information. SIAM J. Appl. Math., 76(1):178–196, 2016.
[30] J. Li and J. Zou. A direct sampling method for inverse scattering using far-field data. Inverse Probl. Imaging, 7(3):757–775, 2013.
[31] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar. Fourier neural operator for parametric partial differential equations. arXiv preprint arXiv:2010.08895, 2020.
[32] F. Luo, J. Wang, J. Zeng, L. Zhang, B. Zhang, K. Xu, and X. Luo. Cascaded complex U-Net model to solve inverse scattering problems with phaseless-data in the complex domain. IEEE Trans. Antennas Propag., 70(8):6160–6170, 2021.
[33] J. Ning, F. Han, and J. Zou. A direct sampling-based deep learning approach for inverse medium scattering problems. Inverse Problems, 40(1):015005, 2023.
[34] R. Novikov. Formulas for phase recovering from phaseless scattering data at fixed frequency. Bull. Sci. Math., 139(8):923–936, 2015.
[35] L. Pan, Y. Zhong, X. Chen, and S. P. Yeo. Subspace-based optimization method for inverse scattering problems utilizing phaseless data. IEEE Trans. Geosci. Remote Sens., 49(3):981–987, 2010.
[36] R. Persico. Introduction to ground penetrating radar: inverse scattering and data processing. John Wiley & Sons, 2014.
[37] R. Potthast. Point sources and multipoles in inverse scattering theory. CRC Press, 2001.
[38] M. Raissi, P. Perdikaris, and G. E. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J. Comput. Phys., 378:686–707, 2019.
[39] O. Ronneberger, P. Fischer, and T. Brox. U-Net: Convolutional networks for biomedical image segmentation. In Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18, pages 234–241. Springer, 2015.
[40] D. N. Tanyu, J. Ning, T. Freudenberg, N. Heilenkötter, A. Rademacher, U. Iben, and P. Maass. Deep learning methods for partial differential equations and related parameter identification problems. Inverse Problems, 39(10):103001, 2023.
[41] P. M. van den Berg. Forward and inverse scattering algorithms based on contrast source integral equations. John Wiley & Sons, 2021.
[42] Z. Wang, E. P. Simoncelli, and A. C. Bovik. Multiscale structural similarity for image quality assessment. In The Thrity-Seventh Asilomar Conference on Signals, Systems & Computers, 2003, volume 2, pages 1398–1402. IEEE, 2003.
[43] Z. Wei and X. Chen. Deep-learning schemes for full-wave nonlinear inverse scattering problems. IEEE Trans. Geosci. Remote Sens., 57(4):1849–1860, 2018.
[44] K. Xu, L. Wu, X. Ye, and X. Chen. Deep learning-based inversion methods for solving inverse scattering problems with phaseless data. IEEE Trans. Antennas Propag., 68(11):7457–7470, 2020.
[45] X. Xu, B. Zhang, and H. Zhang. Uniqueness in inverse scattering problems with phaseless far-field data at a fixed frequency. SIAM J. Appl. Math., 78(3):1737–1753, 2018.
[46] W. Yin, W. Yang, and H. Liu. A neural network scheme for recovering scattering obstacles with limited phaseless far-field data. J. Comput. Phys., 417:109594, 2020.
[47] B. Zhang and H. Zhang. Recovering scattering obstacles by multi-frequency phaseless far-field data. J. Comput. Phys., 345:58–73, 2017.
[48] B. Zhang and H. Zhang. Fast imaging of scattering obstacles from phaseless far-field measurements at a fixed frequency. Inverse Problems, 34(10):104005, 2018.
[49] B. Zhang and H. Zhang. An approximate factorization method for inverse acoustic scattering with phaseless total-field data. SIAM J. Appl. Math., 80(5):2271–2298, 2020.

Ground Truth	$N_{i}=4,\delta=5\%$	$N_{i}=4,\delta=10\%$	$N_{i}=16,,\delta=5\%$	$N_{i}=16,,\delta=10\%$

Ground
Truth
$N_{i}=10$
$\delta=5\%$
$N_{i}=10$
$\delta=10\%$

Ground
Truth
$N_{i}=1$
$\delta=2\%$
$N_{i}=1$
$\delta=10\%$
$N_{i}=4$
$\delta=2\%$
$N_{i}=4$
$\delta=10\%$