Abstract
We proposed numerical methods for solving the direct and inverse scattering problems for domains with multiple corners. Both the near field and far field cases are considered. For the forward problem, the challenges of logarithmic singularity from Green’s functions and corner singularity are both taken care of. For the inverse problem, an efficient and robust direct imaging method is proposed. Multiple frequency data are combined to capture details while not losing robustness.
1. Introduction
Direct and inverse scattering problems have important applications in radar, sonar, and geophysical exploration, in medical imaging, and in nondestructive testing [1]. However, many theoretical and numerical challenges are associated with these problems, especially when the boundary has multiple corners. For direct scattering problem, the singularity of the Green’s function for the Helmholtz equation and the corner singularity of the boundary both require special treatment. For inverse scattering problem, the nonlinearity and ill-posedness challenge the design of an accurate, efficient, and robust numerical algorithm. There are two types of methods for solving the problem. The direct methods [2–6] are efficient but less accurate; the iterative methods [7–15] are accurate but more expensive. Typically the forward and adjoint problems have to be solved in each iteration. For the direct and inverse scattering problems, the near field and far field cases should be treated differently, although the main procedure remains to be the same.
For the direct scattering problem, we extended the work in [1] from the case with one corner to dealing with domain with multiple corners. We managed to make the system well conditioned. For the inverse scattering problem, we used a method similar to the MUSIC algorithm in [3, 4]. However, unlike the MUSIC algorithm which is a projection, we keep the phase information so that multiple frequency data can be combined to obtain both accurate and robust result. Efficiency is also obtained since it is a direct imaging algorithm with no iteration needed. We study both the near field and far field cases.
The organization of the paper is as follows. In Section 2, we will discuss the forward scattering problem for domains with multiple corners. In Section 3, a direct imaging algorithm is presented to solve the inverse problem using response matrix data. In Section 4, numerical examples are presented to demonstrate the efficiency, robustness, and accuracy of our methods. We will conclude in Section 5.
2. The Forward Problem
Consider a time-harmonic plane wave, , incident on a scatterer with multiple corners, where is the wavenumber and is the incident direction. Let be the boundary of the scatterer. set . We consider the obstacle scattering problem. The total field satisfies the Helmholtz equation: The total field consists of the incident field and the scattered field : The incident field satisfies the homogeneous equation: It follows from (1), (2), and (3) that the scattered field satisfies In addition, the scattered field is required to satisfy the following Sommerfeld radiation condition: uniformly in . The uniqueness of the solution to the obstacle scattering problem is discussed in [1].
We use the combined single and double layer potential approach [16] so that the integral equation is uniquely solvable. For simplicity of notation, we first assume there is one corner at ; then By using the jump relation [1], we have the integral equation
We use change of variable and trapezoidal rule as follows [1]. Roughly speaking, about half of the points are equally distributed while the other half are accumulated near the corner. Consider
These lead to a linear system: where , , , , are given in [1].
In [1], the algorithm using graded mesh for solving direct scattering problem for sound-soft obstacle with one corner is presented. However, if we follow the same steps to treat the problem with multiple corners, large condition number for the linear system is observed; see Section 4. We propose a method to reduce the condition number. In [1], a notationally advantageous modification is made to the integral equation by inserting terms involving the fundamental solution to the Laplace equation. We observed that if this particular modification is not made, then the condition number is reduced significantly. In other words, although the notationally advantageous modification in [1] is helpful in theoretical discussion, in practice for domain with one corner it works well; however, for domain with multiple corners, it would lead to a linear system with large condition number. In this case we get around by removing this modification. See Section 4 for numerical comparison.
By solving the integral equation, we could compute both the near field and far field data. There are four cases:(1)plane incident wave, far field data,(2)plane incident wave, near field data,(3)point source, far field data,(4)point source, near field data.
The formula for far field data is
The formula for near field data is
The following formulas are implementation details for near field: Another issue is that for point source incident instead of plan wave we replace with where We have the following response matrix relations: For . For .
The corresponding singular values have the relations
We will verify these relations numerically in Section 4.
3. The Inverse Problem
Shape reconstruction has important applications in radar, sonar, and geophysical exploration, in medical imaging, and in nondestructive testing [1]. The nonlinearity and ill-posedness make it a challenging problem. There are two types of methods for solving the problem. The direct methods [2–6] are efficient but less accurate; the iterative methods [7–15] are accurate but more expensive. Typically the forward and adjoint problems have to be solved in each iteration.
Figure 1 shows a typical configuration for such a problem. The background medium is assumed to be homogeneous.
The response matrix is a collection of the scattered field data received at the th transducer which originated at the th transducer. There are two ways to obtain data for the response matrix. One is to do physical measurements. The other is to solve the forward problem given the target shapes. One way to solve for the scattered field is to truncate the unbounded domain to a bounded domain using the perfectly matched layer (PML) technique [17, 18]. This layer is shown in Figure 1. In this paper, we use our forward solver in Section 2 using a boundary integral formulation. As we discussed earlier, we expect our forward solver to have slower than exponential convergence, but accurate enough to use as input data for inverse problems.
We next review some properties of the singular value decomposition of the response matrix [3, 4].
Depending on the target size compared with the array resolution, the singular value decomposition of the response matrix can have the following three patterns.
For point targets with sizes much smaller than the array resolution, the number of significant singular values equals the number of targets. In this case, the response matrix only contains location information. It is unrealistic to expect to recover shape information.
For small targets whose sizes are smaller than, but comparable to, the array resolution, the pattern of singular values becomes more complicated [19]. The response matrix contains location and some size information.
For extended targets whose sizes are larger than the array resolution, the response matrix contains both location and geometry information of the target. It is no longer clear in the singular value plot how many singular values correspond to one target. In [3] a direct imaging algorithm is developed for extended target. The key idea in the imaging algorithm is to determine the illumination vector based on a physical factorization of the scattered field and the signal space as well as its dimension using resolution analysis.
In this paper, we consider both near field and far field data. Therefore the illumination vector should take the form in [3, 4]. Here we outline the procedure for far field data. For near field it is similar; we consider sound-soft targets. For simplicity we assume here that the outgoing directions we measure are the same as the incoming directions, . The scattered far field is then [1] where is the boundary of the targets, is a unit vector that defines the far field direction, is the total field, and is the outer normal direction on the boundary of the targets. In our setup, the element of the response matrix corresponds to the far field pattern of the scattered field in the th direction due an incident wave coming from the th direction: where the total field is due to incident plane wave coming from the direction . In matrix form where and is the vector of total fields corresponding to the incident plane waves from . Equation (18) gives a physical factorization of the scattered field into known and unknown parts. The far field pattern is a superposition of the far field patterns of point sources located on the boundary of the target; however, we do not know the weight function which depends on the total field. In other words, the scattering at the target boundary acts as sources for the scattered field. In this far field setup, it is natural to use as the illumination vector. The signal space of the response matrix should be well approximated by the span of the illumination vectors with on the well-illuminated part of the boundary of the targets.
The next step is to determine the signal space, which is spanned by appropriate singular vectors of the response matrix. It has been shown in [3] (for near field data) and [4] (for far field data) that, by using a resolution analysis based thresholding, we could determine a threshold and use the first singular vectors to image the shape of the targets.
In this paper, we propose to use the singular values as the natural weight to get around the above thresholding procedure.
One major drawback of the MUSIC-type algorithms [3, 4] is that such algorithms are projection algorithms that remove the phase information. It is not meaningful to combine multiple frequency projection results. For full aperture data, the MUSIC-type algorithms work so well that the drawback is disguised. However, for synthetic aperture data, which is more realistic in some applications, the results from the MUSIC algorithms degenerate. It is crucial to take advantage of phase coherence to overcome the challenge of lack of data [5]. We keep the phase information and use singular values as natural weight. Thus multiple frequency data can be combined. We will show in Section 4 that numerical results using low frequency data are more robust but less accurate; numerical results using high frequency data are more accurate but less robust. By combining multiple frequency results while keeping phase information, we could generate accurate and robust results. Further, efficiency is guaranteed since there is no need for iteration. The evaluation at different grid points is also independent, making it easy to be parallelized.
4. Numerical Experiments
In this section, we present some numerical examples to demonstrate the effectiveness of our method.
We first present results for a smooth target, reproducing the work from [1]. Figure 2 shows the geometry of a kite with grid points on its boundary. Table 1 shows the error of where . Exponential convergence is exhibited since the number of correct digits doubles when doubles. Note that for the fine grid compared with , the error reaches machine precision.
Next, we present an example with three corners. In [1], an example with one corner was presented. However, if we follow the procedure to treat three corners, large condition number is observed. By using the method in Section 2, we reduced the condition number. Figure 3 shows the geometry of three corners with grid points on its boundary. Note that we used graded mesh. About half of the grid points are near the corners. Table 2 shows the error of for . Due to the corner singularities, we no longer have exponential convergence. Still, high order convergence is exhibited.
Next, we presented a more complicated example with the geometry of a butterfly with many corners. The equations for each arc of the butterfly are We first denote to be the number of points on the shortest arc. By using the ratio between arc lengths, we could assign the number of points on each arc. Then we have an equal partition for a parameter that ranges between . Finally, we map the points to the graded mesh in Section 2. Figure 4 shows the butterfly geometry with graded mesh. Table 3 shows the error of for , where is the number of points on the shortest arc. Again, high order convergence is exhibited.
Now we demonstrate the relation among the response matrices for near field and far field data. We compare the following four cases:(1)plane incident wave, far field data,(2)plane incident wave, near field data with ,(3)point source with , far field data,(4)point source with , near field data with .
In Section 2, we derived the relation among the above matrices. Now we use a numerical example to verify the results. We choose the butterfly shape as the obstacle. Figure 5 is the plots for the singular values of the four response matrices. They clearly share the same pattern. Next we justify the relations (14)–(15) numerically. We use the butterfly shape. We find , , , , and . Note that and . The reason is that, with a perturbation of matrix elements without bias, the singular values are more likely to get larger.
Now we present numerical results to justify our claim in Section 2 for condition number. For the shape with three corners with , the condition numbers for the two methods are 14 and 47135. For the butterfly shape, with , they are 90 and 6353. Clearly, by modifying the method in [1], the condition number is significantly reduced. This justifies our claim that the method in [1] for one corner needs to be slightly modified for the case with multiple corners.
We now use some numerical examples to demonstrate the accuracy efficiency and robustness of our method in Section 3 for solving the inverse obstacle problem.
Our first example is the butterfly shape. The response matrix is generated by the forward solver in Section 2.
Figure 6(a) shows the result of our imaging function for low frequency data. The result is robust but not accurate. Figure 6(b) shows the result of our imaging function for high frequency data. It is accurate but not robust. By combining multiple frequency data, Figure 6(c) shows the improved result that is both accurate and robust. Efficiency is also guaranteed since no forward solver is repeatedly needed for iteration.
(a)
(b)
(c)
Next we demonstrate robustness of our method with respect to noisy data. We perturb the response matrix as follows: both the real and imaginary parts of the response matrix are perturbed by 100% with a uniform distribution. Figure 7 shows that our shape reconstruction for the butterfly is still good. The reason behind this is the robustness of the first few singular values with respect to large perturbation.
(a)
(b)
(c)
We next present an example with multiple targets. Figure 8(a) shows the result using low frequency data. It allows us to localize the 5 shapes, but with no details. Figure 8(b) shows the high frequency result with details, but less robustness. With the help from low frequency result, we would zoom the 5 regions one by one. For example, Figure 8(c) shows the zoomed result for the butterfly shape in the middle using high frequency data.
(a)
(b)
(c)
5. Conclusion
We proposed numerical methods for solving the direct and inverse scattering problems for domains with multiple corners. For the forward problem, we extended the method in [1]. The resulting method is well conditioned. For the inverse problem, we proposed a method similar to the MUSIC algorithm in [3, 4] but we keep the phase information so that multiple frequency data can be combined. Numerical results showed that our method is efficient, accurate, and robust.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
Songming Hou’s research is supported by NSF Grant DMS-1317994.