Kirchhoff Migration for Identifying Unknown Targets Surrounded by Random Scatterers

: In this paper, we take into account a two-dimensional inverse scattering problem for localizing small electromagnetic anomalies when they are surrounded by small, randomly distributed electromagnetic scatterers. Generally, subspace migration is considered to be an improved version of Kirchhoff migration; however, for the problem considered here, simulation results have conﬁrmed that Kirchhoff migration is better than subspace migration, though the reasons for this have not been investigated theoretically. In order to explain theoretical reason, we explored that the imaging function of Kirchhoff migration can be expressed by the size, permittivity, permeability of anomalies and random scatterers, and the Bessel function of the ﬁrst kind of order zero and one. Considered approach is based on the fact that the far-ﬁeld pattern can be represented using an asymptotic expansion formula in the presence of such anomalies and random scatterers. We also present results of numerical simulations to validate the discovered imaging function structures.

Kirchhoff and subspace migration are also well-known non-iterative techniques for finding location/shape of unknown inhomogeneities, and they have been applied to a variety of problems (see [38][39][40][41][42][43], for instance). Several studies have confirmed that they are fast, stable, and effective methods for finding various kinds of defects without a priori information of unknown targets. However, most of these researches was performed for homogeneous background media, and further research is still needed on the imaging performance of Kirchhoff and subspace migrations when the unknown targets are surrounded by random scatterers or embedded in an inhomogeneous medium. Notice that for subspace migration, some studies have considered inverse scattering problems in random media [44][45][46][47], random scatterers [48,49], and mathematical theory for detecting point-like scatterers with random noise [50]. However, Kirchhoff migration-based imaging in random media and related mathematical theory has not been considered yet.
In this paper, we consider Kirchhoff migration for localizing small electromagnetic anomalies surrounded by small random scatterers. We carefully analyze its imaging functions by discovering a relationship with the Bessel function of the first kind of order zero and one. This enables us to discuss its various properties and compare the imaging performance with subspace migration. The analysis is based on the structure of the singular vectors associated with non-zero singular values of so-called multi-static response (MSR) matrix and asymptotic expansion formula in the presence of small electromagnetic defects [51].
This paper is organized as follows: In Section 2, we introduce the two-dimensional direct scattering problem and asymptotic expansion formula that holds in the presence of small inhomogeneities surrounded by random scatterers. In Section 3, we introduce the imaging function of Kirchhoff migration and its mathematical structure by establishing a relationship with Bessel function of order zero and one. On the basis of analyzed structure, we discuss its intrinsic properties for several cases of anomalies and random scatterers, and compare the imaging performance with subspace migration. In Section 4, we present some numerical simulation results to support the analytical discussion. We apply various non-iterative techniques for imaging of several extended targets completely hidden in an inhomogeneous medium and discuss their imaging performances in Section 5. Finally, we present our conclusion in Section 6.

Two-Dimensional Direct Scattering Problems
In this section, we briefly survey two-dimensional direct scattering problems and introduce an asymptotic expansion formula that holds true in the presence of small anomalies. For a more detailed discussion of this topic, we recommend [51]. Let A m be anomalies with small diameters α m , m = 1, 2, · · · , M, expressed as ANO denotes the location of A m and the B m are simply connected smooth domains containing the origin. For simplicity, we let Λ be the set of all A m . Analogously, let R s be random scatterers with small radii σ s , s = 1, 2, · · · , S, which is given by RND denotes the location of R s and Υ is the set of all R s . In this paper, we assume that the A m and R s are characterized by their dielectric permittivities and magnetic permeabilities at a given positive angular frequency ω = 2π f . For simplicity, we define the piecewise-constant electric permittivity ε(x) and magnetic permeability µ(x) such that respectively. With this, we set the wavenumber k = ω √ ε 0 µ 0 = 2π/λ, where λ denotes the wavelength that satisfies λ α m and λ σ s for all m and s. Throughout this paper, we consider plane-wave illumination. For a given fixed wavenumber k, u inc (x, θ) = e ikθ·x denotes a plane-wave incident field with propagation direction θ ∈ S 1 , where S 1 denotes the two-dimensional unit circle. Let u(x, θ) be the time-harmonic total field that satisfies the Helmholtz equation with transmission conditions on the boundaries of A m and R s . With this, we let u scat (x, θ) = u(x, θ) − u inc (x, θ) be the scattered field that satisfies the Sommerfeld radiation condition uniformly in all directions ϑ = x/|x| ∈ S 1 . The far-field pattern u ∞ (ϑ, θ) of the scattered field u scat (x, θ) is defined on S 1 : Based on [51], an asymptotic expansion formula of far-field pattern u ∞ (ϑ, θ) can be written as follows. This formula plays a key role of the analysis of mathematical structure in the next section.

Lemma 1.
For sufficiently large k, the far-field pattern u ∞ (ϑ, θ) can be represented as follows.
where |B| denotes the area of B.

Introduction to Kirchhoff Migration and Its Mathematical Structure
The main purpose of this problem is to identify unknown locations x (m) ANO from measured far-field pattern data without any a priori information of targets, e.g., permittivity, permeability, size, shape, etc. For this, we present the Kirchhoff migration technique for a real-time identification of the anomalies A m . For simplicity, we will ignore the term k 2 (1 + i)/4 √ kπ in (1). To introduce this topic, we first consider the MSR matrix K: where {θ n : n = 1, 2, · · · , N} and {ϑ n : n = 1, 2, · · · , N} denote the set of incident and observation directions, respectively. For the sake, we set ϑ n = −θ n for all n = 1, 2, · · · , N, and Now, let us define a test vector and corresponding unit vector where c n ∈ R 1×3 \ {0}, n = 1, 2, · · · , N. With this, we can introduce the following imaging function adopted by the Kirchhoff migration: Then, the map of F KIR (x) will contain peak of large magnitude at the location x ∈ Λ ∪ Υ, so that it will be possible to identify unknown anomalies. A detailed description of Kirchhoff migration is discussed in Appendix A.

Remark 1 (Selection of test vector).
It is worth mentioning that the selection of c n of (2) is highly rely on the shape of A m and R s . Unfortunately, we have no a priori information on targets A m and R s so that it is impossible to select optimal vectors c n . Thus, with motivation from several researches [14,15,41,42], we adopt following unit vector W(x) instead of (3) such that The feasibility of the Kirchhoff migration can be explained based on the discussion in Appendix A. However, sometimes it is impossible to obtain good results via the Kirchhoff migration. Furthermore, appearance of various unexpected results via Kirchhoff migration cannot be explained. Hence, a careful investigation of the mathematical structure of the imaging function F KIR (x) must be considered. For this purpose, we establish the mathematical structure of imaging function. The result is following. Theorem 1. For sufficiently large N and k, F KIR (x) can then be represented as follows: where J n denotes the Bessel function of order n of the first kind.
Proof. See Appendix B.

Various Properties of Kirchhoff Migration
Identified structure (6) allows us to examine how the locations of the anomalies A m are identified using Kirchhoff migration. It is worth mentioning that unlike to the several researches, the applied wave number k, i.e., angular frequency ω must be sufficiently large: for small values of k, the results in Theorem 1 no longer hold and cannot be used to identify A m , we refer to (Figure 7, [49]) for an example of low-frequency imaging. In the following discussion, we consider following three possible cases.
R for all m and s. Then, the terms associated with x (s) RND are dominated by the ones associated with x (m) On the basis of [49], the first M singular vectors, U 1 , U 2 , · · · , U M , will be associated with the vectors W(x (m) ANO ), m = 1, 2, · · · , M, and remaining singular vectors will be associated with the vectors W(x (s) RND ), s = 1, 2, · · · , S. Thus, the singular vectors will satisfy so that it will be possible to discriminate the first M singular values. In this case, the imaging function of subspace migration F SUB (x) becomes (see [52] for instance) Hence, subspace migration will provide better results to the Kirchhoff migration. However, if the wrong singular values are selected, the subspace migration results will actually be worse than those of the Kirchhoff migration.
R for all m and s. Then, even though the first M singular vectors, U 1 , U 2 , · · · , U M will be associated with the vectors W(x (m) ANO ), m = 1, 2, · · · , M, it will still be very hard to discriminate the first M singular values. In this case, if we select the first N 0 singular vectors, This means that a map of F SUB (x) will identify not just the anomalies A m but also some random scatterers R s . Note that a map of F KIR (x) will identify almost all the anomalies A m and random scatterers R s , but the anomaly locations will show peaks of magnitude τ m , m = 1, 2, · · · , M, while the scatterer locations will only show relatively small peaks τ s , s = 1, 2, · · · , N 0 − M. Otherwise, if N 0 < M, location of some anomalies A m , m = N 0 + 1, N 0 + 2, · · · , M, cannot be identified via subspace migration. Therefore, in this case, Kirchhoff migration will provide better results than subspace migration. R for any m and s, both Kirchhoff and subspace migrations will identify a mix of anomalies and random scatterers. Hence, the detection performance of both subspace and Kirchhoff migration will be poor, even if we apply a sufficiently large k.
Based on the discussion, we can conclude that the imaging performance of the subspace migration is better than the one of the Kirchhoff migration if the background is homogeneous and discrimination of nonzero singular values of MSR matrix is clear, i.e., for the imaging of small inhomogeneities (see Figure 2, [14]) for the distribution of singular values corresponding to the target shape). In contrast, the imaging performance of the Kirchhoff migration will be better than the one of the subspace migration if the background is inhomogeneous and discrimination of nonzero singular values of MSR matrix is vague.

Simulation Results
In order to validate the results derived from Theorem 1, we now present a set of simulation results. In this section, we only consider the dielectric permittivity contrast case. For the simulation, and radius α m ≡ 0.1 were placed at the following locations: x In addition, S = 100 small scatterers were randomly distributed over where the η p (a, b), p = 1, 2, are arbitrary real values between a and b. Figure 1 shows the distribution of the three anomalies A m and random scatterers R s . The data u ∞ (ϑ j , θ l ) for the MSR matrix K was generated by solving the Foldy-Lax formulation to avoid an inverse crime (see [53] for instance). We used a total of N = 64 incident and observation directions and wavelengths of λ = 0.7 and λ = 0.2 as low and high frequencies, respectively. After obtaining the far-field data, 20 dB Gaussian random noise is added to the unperturbed data through the MATLAB command awgn included in the signal processing package.  Example 1 (Case 1). In this Example, we consider the imaging results for examining Case 1. For this, we set the permittivities and radii of the random scatterers as ε (s) R = η 3 (2, 3) and σ s = η 5 (0.03, 0.06). Figure 2 shows the distribution of K's normalized singular values and maps of F SUB (x) and F KIR (x) for both wavelengths. Note that although the 3 singular values were successfully discriminated, it is very hard to identify the anomaly locations from the F SUB (x) map for the low applied frequency (λ = 0.7) due to the presence of huge numbers of artifacts. Fortunately, the three anomalies could be extracted satisfactorily from the F KIR (x) map at this frequency. At the higher frequency (λ = 0.2), selecting the first 3 singular values allows the anomaly locations to be identified more clearly from the F SUB (x) map, supporting the discussion of Case 1. That said, the anomaly locations can also be extracted via Kirchhoff migration at this frequency. Example 2 (Case 2). In this Example, we consider the imaging results for Case 2. For this, we set the permittivities and radii of the random scatterers as ε (s) R = η 6 (3.5, 4.5) and σ s = η 8 (0.05, 0.1). Based on the Figure 3, it is far from clear that the first 3 singular values can be identified as the non-zero values of K so the exact threshold of the non-zero singular values cannot be determined. Hence, the F SUB (x) imaging results are poor at the lower frequency. Although there are some artifacts in the F KIR (x) map, all the anomaly locations can be accurately identified. In contrast, at the higher frequency, the F SUB (x) map allows all the anomaly locations to be retrieved very accurately although one random scatterer, located at [−0.58, −0.42] T , is also identified. The permittivity and radius of this random scatterer were 2.5596 and 0.0998, respectively. This scatterer's location was also identified by the F KIR (x) map, but its magnitude was small enough that the anomaly locations could still be accurately identified via Kirchhoff migration.

Further Result: Imaging of Extended Dielectric Targets in an Inhomogeneous Medium
To examine the effectiveness of Kirchhoff migration, additional numerical simulation about the imaging of extended targets completely embedded in an inhomogeneous medium were performed. Following the test configuration in [54], we consider only the permittivity contrast case and set the inhomogeneous domain Ω to be a unit circle with a permittivity of ε 0 = η(0.5, 1.5). For the sake, we set the value of permeability µ(x) = µ 0 ≡ 1 for x ∈ Ω. We assumed the existence of four extended inhomogeneities A m ⊂ Ω with smooth boundaries ∂A m satisfying ∂A m ∩ ∂Ω = ∅ and permittivities ε m , m = 1, 2, 3, 4. The shapes of A m are shown in the Figure 6, and their boundaries are expressed as : 76.5625(s + 0.6) 2 + 81.1898(s + 0.6)t + 29.6875t 2 = 1 , With this configuration, we denote the time-harmonic total field as u (n) (x), which satisfies the following boundary value problem and with transmission condition on the boundary ∂A m , m = 1, 2, 3, 4. Here, ν(x) is the unit normal to ∂Ω at x, and θ n , n = 1, 2, · · · , N(= 128) denotes a two-dimensional vector on the ∂Ω, such that Similarly, u (n) Based on [12,51], u (n) (x) − u (n) 0 (x) can be represented as where φ (n) (y) is an unknown density function, E (ε 0 ) denotes an error term that is highly depending on the value of ε 0 , and N (x, y) satisfies Here, δ denotes the dirac delta function. On the basis of (9), we cannot use the boundary 0 (x) directly to design the Kirchhoff and subspace migrations. So, motivated by [12], we consider the following normalized boundary measurement data B meas (n, n ) (10) where v (n ) (x) = e −ikθ n ·x . Now, let us consider the following MSR matrix Then, K can be decomposed as (see [14] for instance) where E(y) and D(y) denote the illumination and resulting density vectors E(y) = e −ikθ 1 ·y , e −ikθ 2 ·y , · · · , e −ikθ N ·y T and D(y) = φ (1) (y), φ (2) (y), · · · , φ (N) (y) , respectively. Based on [14], it should be noted that the range of K is determined by the span of the E(y) corresponding to A m . This means that the signal subspace can be determined by selecting the singular vectors associated with the nonzero singular values of K. Thus, by taking the test vector W(x) of (5) as imaging functions F KIR (x) and F SUB (x) can be defined similarly with (4) and (7), respectively.
Example 5 (Imaging performances of Kirchhoff and subspace migrations). Figure 6 shows the imaging results via F KIR (x) and F SUB (x) when λ = 0.3. We also mention that 20dB Gaussian random noise is added to the unperturbed data B meas (n, n ), n, n = 1, 2, · · · , N. For performing subspace migration, first 15 singular values are selected. In contrast to the traditional results, it is very hard to identify the outline shapes of every inhomogeneity via F SUB (x) while it is possible to identify the outline shape of every inhomogeneity. Although the size and permittivity of A 4 were both very small, the existence of A 4 was successfully recognized. On the basis of simulation results, we conclude that the result via Kirchhoff migration is better than the one via subspace migration for addressing problems of this type.
Example 6 (Influence of total number of incident fields). Based on the condition in Theorem 1, let us examine the effect of total number N of incident fields. Figure 7 exhibits the maps of F KIR (x) and F SUB (x) for N = 32 and 64. Similar to the traditional researches, it appears that small value of N might be a reason of poor result while large value of N will guarantee good imaging performance of Kirchhoff migration. Unfortunately, it is impossible to identify outline shape of all A m through the subspace migration with small and large N.
Example 7 (Comparison of imaging performances). For the final example, let us apply various non-iterative techniques such as MUSIC, direct sampling method, and factorization method for imaging extended targets and compare the imaging performances. In Figure 8, imaging results via subspace migration (map of F SUB (x)), MUSIC (map of F MUSIC (x)), direct sampling method (map of F DSM (x)), and factorization method (map of F FM (x)) are exhibited with the same simulation configuration of Figure 6. First, let us compare the imaging results via subspace migration and MUSIC algorithm with different selection of nonzero singular values. For this, different from the Figure 6e, N 0 = 10 and N 0 = 37 nonzero singular values are selected to define F SUB (x). However, for any selection of N 0 , an outline shape of inhomogeneities A m cannot be retrieved. Thus, we can examine that the result via Kirchhoff migration is better than the one via subspace migration and MUSIC for imaging of arbitrary shaped inhomogeneities embedded in a random medium.
Next, let us consider the imaging result via direct sampling method with a fixed incident direction θ 1 . Based on the Figure 8c, we can easily examine that the result via direct sampling method is poorer than the one via Kirchhoff migration. Notice that Kirchhoff migration and direct sampling method are identical when the number of incident fields becomes sufficiently large (see Section 5, [24]) so that total number of incident fields must be large enough to obtain a good result via direct sampling method. Figure 8c shows the imaging result via factorization method. It is interesting to observe that similar to the imaging result via Kirchhoff migration, the outline shapes of inhomogeneities are also successfully retrieved. Hence, factorization method can be regarded as an appropriate imaging technique in random medium.

Conclusions
In this paper, we have considered Kirchhoff migration techniques for locating small anomalies when they are surrounded by small random scatterers. In order to investigate the mathematical structure of imaging function and examine the imaging performance, we carefully established a relationship with Bessel functions of the first kind of order zero and one based on the asymptotic expansion formula in the presence of anomalies and random scatterers. Based on the established relationship, we showed that Kirchhoff migration can be regarded as an effective non-iterative technique for identifying the locations of small anomalies when the sizes and permittivities of the random scatterers are smaller than those of the anomalies and the applied frequency is high enough. However, when one of these conditions is not satisfied, the imaging results are somehow poor but Kirchhoff migration yields better results than various non-iterative techniques such as subspace migration, MUSIC, direct sampling method, etc. We have also presented simulation results for several different cases to validate these relationships and observed properties. Although we have discovered various properties of Kirchhoff migration, some phenomenon cannot yet be explained. Further development of mathematical theory for explaining such phenomenon will be a remarkable research topic. Here, we considered the two-dimensional problem, the analysis could be extended to the three-dimensional problem.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Basic Idea of Kirchhoff Migration
Here, a detailed description of Kirchhoff migration is discussed. For this, let us recall the MSR matrix K: Then, since each element u ∞ (ϑ j , θ l ) of K is K can be represented as e ikθ 1 ·x (θ 1 · e 1 )e ikθ 1 ·x (θ 1 · e 2 )e ikθ 1 ·x e ikθ 2 ·x (θ 2 · e 1 )e ikθ 2 ·x (θ 2 · e 2 )e ikθ 2 ·x . . . . . . . . . e ikθ N ·x (θ N · e 1 )e ikθ N ·x (θ N · e 2 )e ikθ N ·x with e 1 = [1, 0] T and e 2 = [0, 1] T . Based on this representation, we can introduce the basic concept of Kirchhoff migration as follows. The singular value decomposition (SVD) of K can be represented as where the superscript * represents the Hermitian operator, U n and V n are the left-and right-singular vectors of K, respectively, and τ n denotes singular value of K such that τ 1 ≥ τ 2 ≥ · · · ≥ τ N 0 > 0 and τ n ≈ ρ ≈ 0 for n > N 0 .
Note that on the basis of (A1), the value of τ n is significantly depending on the size, shape, permittivity, and permeability of A m and R s . Then, by comparing (A1)-(A3), orthonormal property of the singular vectors, and relationship [38], we can observe that for some y n ∈ Λ ∪ Υ, U n ≈ e iγ (1) n W(y n ), V n ≈ e −iγ (2) n W(y n ), and γ (1) n + γ (2) n = arg(τ n ), W(x), U n ≈ 1 and W(x), V n ≈ 1 if x ∈ Λ ∪ Υ, W(x), U n ≈ 0 and W(x), V n ≈ 0 if x ∈ R 2 \(Λ ∪ Υ), where U, V = U · V. With this, we can introduce the following imaging function adopted by the Kirchhoff migration: Following (A4) we can easily observe that the map of F KIR (x) will contain peak of magnitude τ n at the location x ∈ Λ ∪ Υ.

Appendix B. Derivation of Theorem 1
In order to explore the mathematical structure of of F KIR (x), we recall a useful relation derived in [41].
Since N is sufficiently large, applying Lemma A1 yields Similarly, we can get