We consider an inhomogeneous medium that fills up a bounded Lipschitz domain D ⊂ R n (n = 2 or 3). Suppose that the medium is characterized by bounded function η(y) satisfying η = 0 in R n \ D. Consider the incident plane wave
where k > 0 is the wave number and d is the direction vector of propagation. The scattering of u in by the inhomogeneous medium is described by the following model problem ∆u + k 2 (1 + η(x))u = 0, x ∈ R n , (1.1)
∂u sc ∂r -iku sc = 0, r = |x|, (1.3) where u is the total field, u sc is the scattered field, and the Sommerfeld radiation condition (1.3) holds uniformly for all directions x/|x| ∈ S n-1 . If R n \D is connected and Re (η) ≥ 0, this scattering problem is known to have a unique weak solution u sc ∈ H 1 loc (R n ), see [11]. Inverse problem. Let Ω be a bounded Lipschitz domain such that D ⊂ Ω and denote by ν(x) the outward normal unit vector to ∂Ω at x. Given u sc and ∂u sc /∂ν on ∂Ω determine D.
In this paper we develop a sampling type method and combine this method with deep learning to solve the inverse problem above. There has been a large body of literature on both theoretical and numerical studies on this inverse problem and its variants, see [10,11] and references therein. To numerically solve the inverse problem sampling methods are known to be fast, non-iterative and do not require advanced a priori information about the unknown target. More details about sampling methods can be found in [28,23,6] and references therein.
The sampling method studied in this paper is inspired by the orthogonality sampling method (OSM) by Potthast [29]. The OSM is also studied in other works under the name of direct sampling method. The OSM is attractive and promising thanks to its computational efficiency. For instance, the OSM is a regularization free method that is very robust against noise in the data, and its implementation only involves an evaluation of vector products. However, the theoretical analysis of the OSM is far less developed compared with that of the linear sampling method and the factorization method. Also most of the published results deal with the case of far-field measurements, see e.g., [29,13,17,1,18,27,15,24]. There have been only a few results on the OSM concerning the case of near-field measurements. The near-field OSM studied in [3] is applicable to the 2D case with circular measurement boundaries. The 3D case was studied in [20] under the small volume hypothesis of well-separated inhomogeneities. The sampling method proposed in this paper can be applied to near-field data (or far-field data) and it is not limited to those conditions. However, it requires Cauchy data instead of only scattered field data. Using the Lippmann-Schwinger equation and the Helmholtz integral representation we show that the imaging functional is equal to a functional that involves Bessel functions J 0 (k|y -z|) or j 0 (k|y -z|) where z is a sampling point and y is inside the unknown scatterer. These Bessel functions are crucial for justifying the behavior of the imaging functional.
It is known that in the case of scattering data generated by one incident wave the OSM can only provide reasonable reconstructions for small scattering objects [29]. It is also the case of other reconstruction methods. This motivates us to combine our orthogonality sampling type method with deep learning to obtain better reconstruction results. The orthogonality sampling type method is combined with deep learning in the following way. This combined method can be understood as a network that uses the image computed by the sampling method for the first layer, followed by the U-net structure for the remaining layers. The combination of deep learning to physics-based reconstruction methods to solve inverse problems has been studied, see, e.g., [19,14,26,4,21,32,7]. In particular, we refer to the review paper [7] and references therein for recent advances in using deep learning for solving inverse scattering problems. To our knowledge (see also [7]) most of the studies on deep learning for inverse scattering consider data generated by multiple incident fields while the the data for inverse scattering problem considered in this paper is associated with only one incident field at a fixed wave number. Our inverse problem has a minimal amount of data. In addition, the imaging functional for the proposed sampling method is a new functional and it has an advantage over the imaging functional of the OSM that it can be applied to near field data or far field data. Furthermore, to our knowledge the combination of this sampling method (or orthogonality sampling methods) and deep learning for solving inverse scattering problems has not been studied before.
In this paper we exploit a deep neural network of U-Net-Xception style [8]. We train this network using simulated data sets of scattering objects with elliptical geometry. The fast computation and the information from the results of the sampling method help speed up the training of the network. Our numerical study shows that the network is able to predict some scattering objects with geometries different from those of the training and testing data sets. The combination leads to a significant improvement in the reconstruction results initially obtained by the sampling method. We also demonstrate that the combined method is also able to invert some experimental data from the Fresnel Institute ( [5]) without any additional transfer learning.
The paper is organized as follows. The orthogonality sampling type method and its analysis are presented in Section 2. Section 3 is dedicated to a description of the deep neural network used to combine with the sampling method. The numerical study for simulated data and real data is presented in Sections 4 and 5, respectively.
In this section we develop the an orthogonality sampling type method for the inverse scattering problem. For x, y ∈ R 3 and x = y, we denote by Φ(x, y) the free-space Green's function of the scattering problem (1.1)-(1.3) which is given by
It is well known that problem (1.1)-(1.3) is equivalent to the Lippmann-Schwinger equation (see, e.g. [11])
The following Helmholtz integral representation is important to the analysis of the sampling method.
(2.2)
Proof. We refer to [11,Chapter 2] for a proof of the lemma.
We also need to the Funk-Hecke formula for the analysis in this section.
Lemma 2.2. Let Y m be the spherical harmonics of order m and j m be the spherical Bessel functions of first kind and order m. We have
Proof. A proof of this lemma can be found in [11,Chapter 2]. Now we define the imaging functional as
This functional I(z) aims to determine D in Ω and in the numerical implementation we evaluate I(z) for a finite set of sampling points z. We expect that I(z) is relatively large as z ∈ D and that it is small as z is outside D. The behavior of I(z) is analyzed in the following theorem.
Theorem 2.3. The imaging functional I satisfies
Proof. From the Lippmann-Schwinger equation we obtain that
Therefore, a substitution leads to
Since ∆Im Φ(z, y) + k 2 Im Φ(z, y) = 0 for all z, y ∈ R n and Im Φ(z, y) is regular function, the Helmholtz integral representation (2.2) implies that
Therefore, substituting this identity in (2.5) implies
Now the proof can be completed by substituting this equation in (2.4).
We know that
k 4π j 0 (k|z -y|), n = 3. Since J 0 (k|y-z|) and j 0 (k|y-z|) peak when k|y-z| = 0, we expect from Theorem 1 that I(z) takes larger values for z ∈ D and much smaller values outside D. Further, from the asymptotic behavior of J 0 (k|y -z|) and j 0 (k|y -z|) as |y -z| → ∞ we can easily estimate, for z / ∈ D, that
This incomplete analysis is common for imaging functionals of orthogonality sampling methods (see, e.g., [29,17,18,3,27,1]). To our knowledge a complete analysis for orthogonality sampling methods is still an open problem.
We now show a relation between I(z) and the imaging functional suggested by Potthast in [29]. To this end, we recall that the scattered field u sc has the asymptotic behavior
for all x ∈ S n-1 . The function u ∞ (x) is called the scattering amplitude or the far-field pattern of u sc (x). The original imaging functional of the orthogonality sampling method suggested in [29] is given by
Theorem 2.4. The two imaging functionals are related by
Proof. From the Lippmann-Schwinger equation for u sc it is known that its scattering amplitude is given by
e ikx•(z-y) ds(x)η(y)u(y)dy . Now using (2.6) and the Funk-Hecke formula (2.3) we obtain that (2.8)
k Im Φ(z, y), n = 3. Substituting (2.8) in (2.7) and comparing with Theorem 1 we complete the proof.
Theorem 2 also implies that I(z) is stable against noise in the data following the stability of I OSM (z) [29].
In this section, we use the sampling method along with a deep neural network of U-Net-Xception style to form a combined method to solve the inverse problem. Numerical simulation was done for both simulated and experimental scattering data.
3.1. U-Net Xception style. Xception is an efficient architecture which relies on two main points: Depthwise Separable Convolution and Shortcuts between Convolution blocks. Xception architecture is an alternative to classical convolution layers, which performs more efficiently in terms of computation time and accuracy. See details of Xception architecture in [8].
Firstly, Depthwise Separable Convolution includes two subnetworks: Depthwise Convolution and Pointwise Convolution. Compared to conventional convolution layers, Depthwise Separable Convolution does not need to compute convolution operations across all channels, which will eliminate the number of parameters in the model and makes it simpler. The overall architecture of Depthwise Separable Convolution is described in Figure 1.
Secondly, there are residual (shortcut/skip) connections in Xception networks. Many state-of-the-art Deep neural networks (Resnet, YOLOv3, MobileNetV2, Cascading Residual Network) have shown that accuracy is much higher when residual connections are used, see [16,30,31,2]. The pattern of residual connection in U-Net Xception architecture is shown in Figure 2.
Training the network. Given some scattering object and an incident wave we use the spectral method studied in [25] to solve the direct problem and generate the associated scattering data. The image of the scattering object is called a true image. Applying the sampling method to the scattering data we obtain a computed image of the scattering object. We call this image a preliminary image.
The training data set includes 60, 000 pairs of images, each pair consists of a true image and a preliminary image. The true images are those of either one or two ellipses whose center and radii are randomly generated within the square domain (-2, 2) 2 . More precisely, if the object consists of only one ellipse, the center is randomly generated within the square [-0.8, 0.8] 2 and the radii are randomly generated within the interval [0.1, 1]; if the object consists of two ellipses, the first one is generated in the same manner as before, while the second one has its center randomly generated within the square [-1, 1] 2 and radii randomly generated within the interval [0.1, 0.5]. Note that the two ellipses are allowed to overlap each other. The wave number associated with the data set is k = 6, which means the wavelength is about 1. As for the incident wave, we use u in (x) = e ik(x1 cos θ+x2 sin θ) where θ = 90 o for 30,000 images and θ = 45 o for the other 30,000 images. The Cauchy scattering data are computed on the circle of radius 100. We did not add any artificial noise to the scattering data for the training and testing of the network. See Figure 3 for some examples of training data pairs. All images are uniformly partitioned into 160×160 pixels. Let I = {1, 2, . . . , 160 2 }. The neural network is initialized with random parameters. Using a preliminary image as input, the network produces the corresponding output ŷ = {ŷ i } i∈I that will then be used along with the true image y = {y i } i∈I to calculate the loss function
After that, the parameters within the convolution layers are tuned using Adam method (see [22]) in order to minimize the loss function. This minimization process is done for one batch of size 32 at a time, with the loss function for each batch being the sum of those of individual training pairs in the batch. The initial learning rate is 10 -2 , and will be divided by 10 if the validation loss does not decrease or validation accuracy does not increase after 5 epochs until it reaches 10 -6 . The neural network's architecture can be briefly described by the diagram in Figure 4.
Out of 60,000 pairs of images, 50,000 were used to train the network, 5,000 were used for validation and 5,000 are used for testing. In Figure 5, we show the result of training loss versus validation loss and training accuracy versus validation accuracy. The loss value at each epoch is the average value of L over all paired images, whereas the accuracy value is the ratio of the number of pixels in the predicted image that match their counterparts in the true image to the total number of pixels. Those values for validation were computed over 5, 000 images in the validation set, and those values for training were computed over 50, 000 images in the training set.
In order to enrich the training data set, we used data augmentation. More specifically, we produced five copies of each training pair by applying five types of transformations: horizontal flip, vertical flip, 90 o rotation, -90 o rotation and zoom. Zoom transformations were done by first upscaling the image by 5%, then using a random frame of the same size as the original image to crop it. See Figure 6 for an example of such augmentations.
Our choice for some hyperparameters was based on direct comparison to other choices in terms of validation loss and accuracy. For example, we trained the model with 10, 000 pairs of images through 5 epochs with different input image sizes, then we chose 160 × 160 × 3 because it gave the best values for loss and accuracy. We did a similar experiment for two different optimization methods, which are Adam method and stochastic gradient descent with momentum (SGD), and Adam method performed better for this problem. Detailed results for the comparisons are given in Table 1 and2. As for the other hyperparameters such as the number of hidden layers and units, we followed the same choice as in [9]. With this choice, we were able to get very good loss and accuracy values, see Figure 5. Therefore, no adjustment was needed.
We also compared the performance of the U-Net Xception style with a modified U-Net in [12]. The result in Table 3 shows that U-Net Xception style works better for our problem. 4.1. Implementation of the sampling method. In our numerical study the sampling domain is chosen as [-2, 2] 2 . The sampling domain contains the unknown scatterer D and is uniformly discretized by 64 sampling points in each direction. The imaging functional I(z) in (2.4) is evaluated at each sampling point z of the sampling domain. As mentioned in the section of training the neural network the data (u sc , ∂u sc /∂ν) are given at 32 uniformly distributed points on the circle of radius 100. With the data given, the evaluation of I(z) is simple since after a change of variables using polar coordinates we numerically evaluate a single integral using the rectangular rule. In the numerical simulation the imaging functional is normalized by dividing by its (numerical) maximal value. Note that we make sure that the sampling point z and the data point x are away from each other so that the integrands are all regular functions. Recall that the scattering data are measured on ∂Ω meaning that we should look for the unknown scatterer D in a sampling domain that is strictly contained in Ω.
In this part, we present some reconstruction results of the combined method for objects in the testing data set. These objects were generated in a way that is similar to that of the training data sets, meaning the objects consist of one or two ellipses. The wave number k and incident direction θ were those the neural network was trained with. Moreover, the scattering data generated by solving the direct problem was not perturbed by noise.
In Figures 78, four test objects are presented: one large disk, two ellipses of similar sizes, two ellipses of different sizes and two overlapping ellipses. We can see that the method was able to reconstruct very well the geometry of these test objects. As it was done for the validation and training data sets in Figure 5, the accuracy for 5000 pairs of images in the testing data set is also computed as 99.4%. This once again confirms that the training process was effective as indicated by the training and validation accuracy in Figure 5.
Reconstruction of objects with no elliptical shapes. In this part, we test the method for other types of objects with noisy scattering data. We added 5%, 7%, 10% and 15% of artificial noise to the scattering data obtained from solving the direct problem with wave number k = 6. Recall that the inverse problem with one incident wave is extremely ill-posed, and there were no artificial noise added to the scattering data in the training process. The tests with noisy data for nonelliptical objects are to show that the network can generalize its training to more see Figure 9. The combined method was able to provide reasonable reconstructions for the L-shaped object with different levels of noise.
Next, in Figures 101112, we test the method with different non-elliptical shapes and noisy scattering data associated with an incident direction different from which was used in the training process. This is to further assert its flexibility. Recall that we trained the neural network with θ = 90 o and θ = 45 o . We reconstructed a peanut-shaped object, a T-shaped object and an object consisting of a disk and a rectangle using incident direction θ = 220 o , see Figure 10-12. Again, we were able to get satisfactory results regardless of the circumstances.
For scattering data from a single incident direction and a fixed wave number, it is almost impossible not only for sampling methods but also for many other inversion methods to reconstruct the shape and location of extended objects. Therefore, these results show that the deep neural network has significantly improved the reconstruction capability of the imaging functional under lack of data. Moreover, they also show that the combined method is flexible, meaning it can work well in some situations that are unseen from the training process.
The results from simulated scattering data show that the proposed deep learning approach for OSM is capable of reconstructing different types of geometries in different setups. Therefore in this section, we demonstrate the effectiveness of our method when it comes to experimental data. We used the data sets for two-dimensional homogeneous objects provided by Fresnel Institute [5]. Unlike the simulated data that are computed on an entire .exp deals with a de-centered circular cross section of radius 15 mm, and the second one is rectTM cent.exp concerning a centered rectangular cross section of dimensions 25.4 × 12.7 mm 2 . We refer to [5] for detailed descriptions of the experimental setup. We rescale 40 millimeters to be 1 unit of length in our MATLAB simulations. No other adaptation to the experimental system was necessary. The two scatterers were both tested at wave frequency 8 GHz (wave number k is about 6.7) and with incident direction θ = 90 o . See reconstructions using our combined method in Figure 13 and Figure 14.
We see that the proposed method can handle limited-aperture experimental data very well without any additional transfer learning. For the first object, we are able to achieve a reasonable reconstruction of both shape and location of the target. The second object is outside the training process, but the inverted result still can capture highly accurate information for its location and provide a pretty good estimate for its shape.