# Interpretation of MUSIC for Location Detecting of Small Inhomogeneities Surrounded by Random Scatterers

Authors: Won-Kwang Park
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

We consider the MUltiple SIgnal Classification (MUSIC) algorithm for identifying the locations of small electromagnetic inhomogeneities surrounded by random scatterers. For this purpose, we rigorously analyze the structure of MUSIC-type imaging function by establishing a relationship with zero-order Bessel function of the first kind. This relationship shows certain properties of the MUSIC algorithm, explains some unexplained phenomena, and provides a method for improvements.


---

## Body

## 1. Introduction

One of the purposes of the inverse scattering problem is to identify the characteristics (location, shape, material properties, etc.) of small inhomogeneities from the scattered field or far-field pattern. This problem, which arises in fields such as physics, engineering, and biomedical science, is highly relevant to human life; thus, it remains an important research area. Related works can be found in [1–5] and references therein.Attempts to address the problem described above have led to the development of the MUltiple SIgnal Classification- (MUSIC-) type algorithm to find unknown inhomogeneities and the algorithm has been applied to various problems, for example, detection of small inhomogeneities in homogeneous space [6–9], location identification of small inhomogeneities embedded in a half-space or multilayered medium [10–12], reconstructing perfectly conducting cracks [13, 14], imaging of internal corrosion [15], shape recognition of crack-like thin inhomogeneities [16–18] and volumetric extended targets [19–21], and application to the biomedical imaging [22]. We also refer to [23, 24] for a detailed and concise description of MUSIC. Several research efforts have contributed to confirming that MUSIC is a fast and stable algorithm that can easily be extended to multiple inhomogeneities and that does not require specific regularization terms that are highly dependent on the problem at hand. However, its feasibility is only confirmed when the background medium is homogeneous; that is, the imaging performance of MUSIC when unknown inhomogeneities are surrounded by random scatterers remains unknown. In several works [25–28], an inverse scattering problem in random media has been concerned. Specially, mathematical theory of MUSIC for detecting point-like scatterers embedded in an inhomogeneous medium has been concerned in [29]. Motivated by these remarkable works, a more careful investigation of the mathematical theory is still required.Motivated by the above, MUSIC algorithm has been applied for detecting the locations of small electromagnetic inhomogeneities when they are surrounded by electromagnetic random scatterers and confirmed that it can be applied satisfactorily. However, this only relied on the results of numerical simulations, that is, a heuristic approach to some extent, which is the motivation for the current work. In this contribution, we carefully analyze the mathematical structure of MUSIC-type imaging function and discover some properties. This work is based on the relationship between the singular vectors associated with nonzero singular values of a multistatic response (MSR) matrix and asymptotic expansion formula due to the existence of small inhomogeneities; refer to [23].This paper is organized as follows. Section2 introduces the two-dimensional direct scattering problem and an asymptotic expansion formula in the presence of small inhomogeneities. In Section 3, MUSIC-type imaging function is introduced. In Section 4, we analyze the mathematical structure of the MUSIC-type imaging function and discuss its properties. In Section 5, we present the results of numerical simulations to support the analyzed structure of MUSIC and Section 6 presents a short conclusion.


## 2. Two-Dimensional Direct Scattering Problem

In this section, we survey a two-dimensional direct scattering problem and introduce an asymptotic expansion formula. For a more detailed description we recommend [18, 23, 30]. Let Σ
m, m
=
1,2
,
…
,
M, be an electromagnetic inhomogeneity with a small diameter r
m in two-dimensional space R
2. Throughout this paper, we assume that every Σ
m is expressed as(1)
Σ
m
=
z
m
+
r
m
B
m
,where z
m denotes the location of Σ
m and B
m is a simple connected smooth domain containing the origin. For the sake of simplicity, we let Σ be the collection of Σ
m. Throughout this paper, we assume that inhomogeneities are well separated from each other such that(2)
ω
z
m
-
z
m
′
≫
1
-
1
4
=
0.75
,for all m
,
m
′
=
1,2
,
…
,
M and m
≠
m
′.Let us denoteΔ
s, s
=
1,2
,
…
,
S, as the random scatterer with small radius r
s
<
r and let Δ be the collection of Δ
s. Similarly, we assume that Δ
s is of the form(3)
Δ
s
=
y
s
+
r
s
B
m
.As before, suppose that Δ
s
∩
Δ
s
′
=
∅ for all s
,
s
′
=
1,2
,
…
,
S and s
≠
s
′ and the positions of y
s are random but they are fixed for all frequencies discussed later.In this work, we assume that every inhomogeneity is characterized by its dielectric permittivity and magnetic permeability at a given positive angular frequencyω
=
2
π
/
λ, where λ denotes the wavelength. Let ɛ
m, ɛ
s, and ɛ
0 be the electric permittivities of Σ
m, Δ
s, and R
2, respectively. Then, we can introduce the piecewise-constant electric permittivity ɛ
(
x
) and magnetic permeability μ
(
x
) such that(4)
ɛ
x
=
ɛ
m
for
x
∈
Σ
m
ɛ
s
for
x
∈
Δ
s
ɛ
0
for
x
∈
R
2
∖
Σ
¯
∪
Δ
¯
,
μ
x
=
μ
m
for
x
∈
Σ
m
μ
s
for
x
∈
Δ
s
μ
0
for
x
∈
R
2
∖
Σ
¯
∪
Δ
¯
,respectively. For the sake of simplicity, we let ɛ
0
=
μ
0
=
1, ɛ
m
>
ɛ
s, and μ
m
>
μ
s for all m and s. Hence, we can set the wavenumber k
=
ω
ɛ
0
μ
0
=
ω.For a given fixed frequencyω, we denote(5)
u
i
n
c
x
,
θ
=
e
i
ω
θ
·
xto be a plane-wave incident field with the incident direction θ
∈
S
1, where S
1 denotes the two-dimensional unit circle. Let u
(
x
,
θ
) denote the time-harmonic total field that satisfies the following Helmholtz equation(6)
∇
·
1
μ
x
∇
u
x
,
θ
+
ω
2
ɛ
x
u
x
,
θ
=
0with transmission conditions on the boundaries of Σ
m and Δ
s. This configuration is associated with a scalar scattering problem for E-polarized (Transverse Magnetic (TM) polarization, corresponding to dielectric contrasts) field; the H-polarized (Transverse Electric (TE) polarization, corresponding to magnetic contrasts) case could be dealt with per duality. It is well known that u
(
x
,
θ
) can be decomposed as(7)
u
x
,
θ
=
u
i
n
c
x
,
θ
+
u
s
c
a
t
x
,
θ
,where u
s
c
a
t
(
x
,
θ
) denotes the unknown scattered field that satisfies the Sommerfeld radiation condition(8)
lim
x
→
0
⁡
x
∂
u
s
c
a
t
x
,
θ
∂
x
-
i
ω
u
s
c
a
t
x
,
θ
=
0uniformly in all directions ϑ
=
x
/
x
∈
S
1. The far-field pattern u
f
a
r
(
ϑ
,
θ
) of the scattered field u
s
c
a
t
(
x
,
θ
) is defined on S
1. It can be expressed as(9)
u
s
c
a
t
x
,
θ
=
e
i
ω
x
x
u
f
a
r
ϑ
,
θ
+
o
1
x
,
x
⟶
+
∞
.Then by virtue of [31], the far-field pattern u
f
a
r
(
ϑ
,
θ
) can be written as the following asymptotic expansion formula, which plays a key role in the MUSIC-type algorithm that will be designed in the next section:(10)
u
f
a
r
ϑ
,
θ
=
ω
2
1
+
i
4
ω
π
∑
m
=
1
M
r
m
2
ɛ
m
-
ɛ
0
B
m
-
μ
0
μ
m
+
μ
0
2
ϑ
·
2
θ
e
i
ω
θ
-
ϑ
·
z
m
+
∑
s
=
1
S
r
s
2
ɛ
s
-
ɛ
0
ɛ
0
μ
0
B
s
-
μ
0
μ
s
+
μ
0
2
ϑ
·
2
θ
e
i
ω
θ
-
ϑ
·
y
s
.


## 3. MUSIC-Type Imaging Algorithm

In this section, we introduce the MUSIC-type algorithm for detecting the locations of small inhomogeneities. For the sake of simplicity, we exclude the constant termω
2
(
1
+
i
)
/
4
ω
π from (10). For this, let us consider the eigenvalue structure of the MSR matrix(11)
K
=
u
f
a
r
ϑ
1
,
θ
1
u
f
a
r
ϑ
1
,
θ
2
⋯
u
f
a
r
ϑ
1
,
θ
N
u
f
a
r
ϑ
2
,
θ
1
u
f
a
r
ϑ
2
,
θ
2
⋯
u
f
a
r
ϑ
2
,
θ
N
⋮
⋮
⋱
⋮
u
f
a
r
ϑ
N
,
θ
1
u
f
a
r
ϑ
N
,
θ
2
⋯
u
f
a
r
ϑ
N
,
θ
N
.Suppose that ϑ
j
=
-
θ
j for all j; then K is a complex symmetric matrix but not a Hermitian. Thus, instead of eigenvalue decomposition, we perform singular value decomposition (SVD) of K (see [24], for instance)(12)
K
≈
∑
m
=
1
3
M
σ
m
U
m
V
m
∗
+
∑
s
=
3
M
+
1
3
M
+
3
S
σ
s
U
s
V
s
∗
,where superscript ∗ is the mark of a Hermitian. Then, U
1
,
U
2
,
…
,
U
3
M
+
3
S is the orthogonal basis for the signal space of K. Therefore, one can define the projection operator onto the null (or noise) subspace, P
n
o
i
s
e
:
C
N
×
1
→
C
N
×
1. This projection is given explicitly by(13)
P
n
o
i
s
e
≔
I
N
-
∑
m
=
1
3
M
+
3
S
U
m
U
m
∗
,where I
N denotes the N
×
N identity matrix. For any point x
∈
R
2 and suitable vectors c
n
∈
R
3
∖
0, n
=
1,2
,
…
,
N, define a test vector f
(
x
)
∈
C
N
×
1 as(14)
f
x
=
c
1
·
1
,
θ
1
T
e
i
ω
θ
1
·
x
,
c
2
·
1
,
θ
2
T
e
i
ω
θ
2
·
x
,
…
,
c
N
·
1
,
θ
N
T
e
i
ω
θ
N
·
x
T
.Then, by virtue of [23], there exists N
0
∈
N such that, for any N
≥
N
0, the following statement holds:(15)
f
x
∈
Range
K
K
¯
iff
x
∈
z
m
,
y
sfor m
=
1,2
,
…
,
M and s
=
1,2
,
…
,
S. This means that if x
∈
Σ
m or x
∈
Δ
s then P
n
o
i
s
e
(
f
(
x
)
)
=
0. Thus, the locations of Σ
m and Δ
s follow from computing the MUSIC-type imaging function(16)
F
x
=
1
P
n
o
i
s
e
f
x
.The resulting plot of F
(
x
) will have peaks of large magnitudes at z
m
∈
Σ
m and y
s
∈
Δ
s.Remark 1.
Based on several works [17, 18, 20], selection of c
n in (14) is highly depending on the shape of Σ
m. Unfortunately, the shape of  Σ
m is unknown; it is impossible to find proper vectors c
n. Due to this fact, following [20], we assume that c
n
·
[
1
,
θ
n
]
T
=
1 for all n; that is, we consider the following test vector instead of (14):(17)
f
x
=
1
N
e
i
ω
θ
1
·
x
,
e
i
ω
θ
2
·
x
,
…
,
e
i
ω
θ
N
·
x
Tand we analyze the mathematical structure of F
(
x
).


## 4. Structure of Imaging Function

Henceforth, we analyze the mathematical structure ofF
(
x
) and examine certain of its properties. Before starting, we recall a useful result derived in [32].Lemma 2.
Assume thatθ
n
:
n
=
1,2
,
…
,
N spans S
1. Then, for sufficiently large N, ξ
∈
S
1, and x
∈
R
2, the following relation holds:(18)
1
N
∑
n
=
1
N
e
i
ω
θ
n
·
x
=
1
2
π
∫
S
1
e
i
ω
θ
·
x
d
S
θ
=
J
0
ω
x
,
1
N
∑
n
=
1
N
ξ
·
θ
n
e
i
ω
θ
n
·
x
=
1
2
π
∫
S
1
ξ
·
θ
e
i
ω
θ
·
x
d
S
θ
=
i
x
x
·
ξ
J
1
ω
x
,where J
n denotes Bessel function of order n of the first kind.Now, we introduce the main result.Theorem 3.
For sufficiently largeN
>
3
M
+
3
S and ω, F
(
x
) can be represented as follows: for e
1
=
[
1,0
]
T and e
2
=
[
0,1
]
T, (19)
F
x
≈
1
-
∑
m
=
1
M
J
0
ω
x
-
z
m
2
-
∑
m
=
1
M
∑
h
=
1
2
x
-
z
m
·
e
h
x
-
z
m
2
J
1
ω
x
-
z
m
2
-
∑
s
=
1
S
J
0
ω
x
-
y
s
2
-
∑
s
=
1
S
∑
h
=
1
2
x
-
y
s
·
e
h
x
-
y
s
2
J
1
ω
x
-
y
s
2
-
1
/
2
.Proof.
Based on the asymptotic expansion formula (10) and results in [13], P
n
o
i
s
e can be represented as(20)
P
n
o
i
s
e
=
I
N
-
∑
m
=
1
3
M
U
m
U
m
∗
-
∑
s
=
1
3
S
U
s
U
s
∗
≈
I
N
-
∑
m
=
1
M
∑
h
=
1
3
W
h
z
m
W
h
z
m
∗
-
∑
s
=
1
S
∑
h
=
1
3
W
h
y
s
W
h
y
s
∗
,where(21)
W
1
x
=
1
N
e
i
ω
θ
1
·
x
,
e
i
ω
θ
2
·
x
,
…
,
e
i
ω
θ
N
·
x
T
,
W
2
x
=
2
N
e
1
·
θ
1
e
i
ω
θ
1
·
x
,
e
1
·
θ
2
e
i
ω
θ
2
·
x
,
…
,
e
1
·
θ
N
e
i
ω
θ
N
·
x
T
,
W
3
x
=
2
N
e
2
·
θ
1
e
i
ω
θ
1
·
x
,
e
2
·
θ
2
e
i
ω
θ
2
·
x
,
…
,
e
2
·
θ
N
e
i
ω
θ
N
·
x
T
.With this, applying (18) and performing a tedious calculation, we arrive at(22)
P
n
o
i
s
e
f
x
=
f
x
-
1
N
N
∑
m
=
1
M
A
z
m
+
B
1
z
m
+
B
2
z
m
-
∑
s
=
1
S
A
y
s
+
B
1
y
s
+
B
2
y
s
,where(23)
A
ξ
≔
e
i
ω
θ
1
·
ξ
J
0
ω
x
-
ξ
e
i
ω
θ
2
·
ξ
J
0
ω
x
-
ξ
⋮
e
i
ω
θ
N
·
ξ
J
0
ω
x
-
ξ
,
B
h
ξ
≔
i
e
h
·
θ
1
x
-
ξ
·
e
h
x
-
ξ
e
i
ω
θ
1
·
ξ
J
1
ω
x
-
ξ
i
e
h
·
θ
2
x
-
ξ
·
e
h
x
-
ξ
e
i
ω
θ
2
·
ξ
J
1
ω
x
-
ξ
⋮
i
e
h
·
θ
N
x
-
ξ
·
e
h
x
-
ξ
e
i
ω
θ
N
·
ξ
J
1
ω
x
-
ξfor ξ
∈
R
2 and h
=
1,2. By implementing elementary calculus, we can show that(24)
P
n
o
i
s
e
f
x
2
=
P
n
o
i
s
e
f
x
P
n
o
i
s
e
f
x
¯
=
1
N
∑
n
=
1
N
1
-
∑
h
=
1
8
Φ
h
,where(25)
Φ
1
=
∑
m
=
1
M
e
i
ω
θ
n
·
x
-
z
m
+
e
-
i
ω
θ
n
·
x
-
z
m
J
0
ω
x
-
z
m
,
Φ
2
=
-
∑
m
=
1
M
e
i
ω
θ
n
·
z
m
J
0
ω
x
-
z
m
∑
m
′
=
1
M
e
-
i
ω
θ
n
·
z
m
′
J
0
ω
x
-
z
m
′
,
Φ
3
=
-
i
∑
m
=
1
M
∑
h
=
1
2
e
h
·
θ
n
x
-
z
m
·
e
h
x
-
z
m
e
i
ω
θ
n
·
x
-
z
m
-
e
-
i
ω
θ
n
·
x
-
z
m
J
1
ω
x
-
z
m
,
Φ
4
=
-
∑
m
=
1
M
∑
h
=
1
2
e
h
·
θ
n
x
-
z
m
·
e
h
x
-
z
m
e
i
ω
θ
n
·
z
m
J
1
ω
x
-
z
m
∑
m
′
=
1
M
∑
h
=
1
2
e
h
·
θ
n
x
-
z
m
′
·
e
h
x
-
z
m
′
e
i
ω
θ
n
·
z
m
′
J
1
ω
x
-
z
m
′
,
Φ
5
=
∑
s
=
1
S
e
i
ω
θ
n
·
x
-
y
s
+
e
-
i
ω
θ
n
·
x
-
y
s
J
0
ω
x
-
y
s
,
Φ
6
=
-
∑
s
=
1
S
e
i
ω
θ
n
·
y
s
J
0
ω
x
-
y
s
∑
s
′
=
1
S
e
-
i
ω
θ
n
·
y
s
′
J
0
ω
x
-
y
s
′
,
Φ
7
=
-
i
∑
s
=
1
M
∑
h
=
1
2
e
h
·
θ
n
x
-
y
s
·
e
h
x
-
y
s
e
i
ω
θ
n
·
x
-
y
s
+
e
-
i
ω
θ
n
·
x
-
y
s
J
1
ω
x
-
y
s
,
Φ
8
=
-
∑
s
=
1
M
∑
h
=
1
2
e
h
·
θ
n
x
-
y
s
·
e
h
x
-
y
s
e
i
ω
θ
n
·
y
s
J
1
ω
x
-
y
s
∑
s
′
=
1
M
∑
h
=
1
2
e
h
·
θ
n
x
-
z
s
′
·
e
h
x
-
z
s
′
e
i
ω
θ
n
·
z
s
′
J
1
ω
x
-
z
s
′
.First, applying (18), we can obtain(26)
1
N
∑
n
=
1
N
∑
m
=
1
M
e
i
ω
θ
n
·
x
-
z
m
J
0
ω
x
-
z
m
=
∑
m
=
1
M
J
0
ω
x
-
z
m
2
.This leads us to(27)
1
N
∑
n
=
1
N
Φ
1
=
1
N
∑
n
=
1
N
∑
m
=
1
M
e
i
ω
θ
n
·
x
-
z
m
+
e
-
i
ω
θ
n
·
x
-
z
m
J
0
ω
x
-
z
m
=
2
∑
m
=
1
M
J
0
ω
x
-
z
m
2and similarly to(28)
1
N
∑
n
=
1
N
Φ
5
=
1
N
∑
n
=
1
N
∑
s
=
1
S
e
i
ω
θ
n
·
x
-
y
s
+
e
-
i
ω
θ
n
·
x
-
y
s
J
0
ω
x
-
y
s
=
2
∑
s
=
1
S
J
0
ω
x
-
y
s
2
.
Next, based on the orthonormal property of singular vectors, relations (2) and (18), and the following asymptotic form(29)
J
0
ω
z
-
z
m
′
≈
2
ω
z
-
z
m
′
π
cos
⁡
ω
z
-
z
m
′
-
π
4
,we can derive(30)
1
N
∑
n
=
1
N
Φ
2
=
-
1
N
∑
n
=
1
N
∑
m
=
1
M
e
i
ω
θ
n
·
z
m
J
0
ω
x
-
z
m
∑
m
′
=
1
M
e
-
i
ω
θ
n
·
z
m
′
J
0
ω
x
-
z
m
′
=
-
∑
m
=
1
M
∑
m
′
=
1
M
1
N
∑
n
=
1
N
e
i
ω
θ
n
·
z
-
z
m
′
J
0
ω
x
-
z
m
J
0
ω
x
-
z
m
′
=
-
∑
m
=
1
M
∑
m
′
=
1
M
J
0
ω
z
-
z
m
′
J
0
ω
x
-
z
m
J
0
ω
x
-
z
m
′
=
-
∑
m
=
1
M
J
0
ω
x
-
z
m
2and similarly(31)
1
N
∑
n
=
1
N
Φ
6
=
-
∑
s
=
1
S
J
0
ω
x
-
y
s
2
.
For evaluatingΦ
3, let us perform an elementary calculus (32)
1
N
∑
n
=
1
N
i
∑
m
=
1
M
∑
h
=
1
2
e
h
·
θ
n
e
i
ω
θ
n
·
x
-
z
m
x
-
z
m
·
e
h
x
-
z
m
J
1
ω
x
-
z
m
=
∑
m
=
1
M
∑
h
=
1
2
i
1
N
∑
n
=
1
N
e
h
·
θ
n
e
i
ω
θ
n
·
x
-
z
m
x
-
z
m
·
e
h
x
-
z
m
J
1
ω
x
-
z
m
=
-
∑
m
=
1
M
∑
h
=
1
2
x
-
z
m
·
e
h
x
-
z
m
2
J
1
ω
x
-
z
m
2
.Then, we can conclude that(33)
1
N
∑
n
=
1
N
Φ
3
=
2
∑
m
=
1
M
∑
h
=
1
2
x
-
z
m
·
e
h
x
-
z
m
2
J
1
ω
x
-
z
m
2
,
1
N
∑
n
=
1
N
Φ
7
=
2
∑
s
=
1
S
∑
h
=
1
2
x
-
y
s
·
e
h
x
-
y
s
2
J
1
ω
x
-
y
s
2
.
Finally, forΦ
4, by applying following integral, for θ
n
,
θ
,
ξ
∈
S
1,(34)
1
N
∑
n
=
1
N
θ
n
·
ξ
2
≈
1
2
π
∫
S
1
θ
·
ξ
2
d
θ
=
1
2
,we can derive the following:(35)
1
N
∑
n
=
1
N
Φ
4
=
-
1
N
∑
n
=
1
N
∑
m
=
1
M
∑
h
=
1
2
e
h
·
θ
n
x
-
z
m
·
e
h
x
-
z
m
e
i
ω
θ
n
·
z
m
J
1
ω
x
-
z
m
∑
m
′
=
1
M
∑
h
′
′
=
1
2
e
h
′
′
·
θ
n
x
-
z
m
′
·
e
h
′
′
x
-
z
m
′
e
i
ω
θ
n
·
z
m
′
J
1
ω
x
-
z
m
′
=
-
∑
m
=
1
M
1
N
∑
n
=
1
N
∑
h
=
1
2
e
s
·
θ
n
2
2
∑
h
=
1
2
x
-
z
m
·
e
s
x
-
z
m
J
1
ω
x
-
z
m
2
=
-
∑
m
=
1
M
∑
h
=
1
2
x
-
z
m
·
e
h
x
-
z
m
J
1
ω
x
-
z
m
2
.Correspondingly,(36)
1
N
∑
n
=
1
N
Φ
8
=
-
∑
s
=
1
S
∑
h
=
1
2
x
-
y
s
·
e
h
x
-
y
s
J
1
ω
x
-
y
s
2
.Hence, by combining (27)–(36), we can obtain the following mathematical structure: (37)
P
n
o
i
s
e
f
r
2
=
1
-
∑
m
=
1
M
J
0
ω
x
-
z
m
2
-
∑
m
=
1
M
∑
h
=
1
2
x
-
z
m
·
e
h
x
-
z
m
2
J
1
ω
x
-
z
m
2
-
∑
s
=
1
S
J
0
ω
x
-
y
s
2
-
∑
s
=
1
S
∑
h
=
1
2
x
-
y
s
·
e
h
x
-
y
s
2
J
1
ω
x
-
y
s
2
.This enables us to obtain the desired result. This completes the proof.Remark 4 (applicability of MUSIC).
SinceJ
0
(
0
)
=
1, the value of F
(
x
) will be sufficiently large when x
=
z
m or y
s for all m and s. Hence, based on the result in Theorem 3, the locations of Σ
m and Δ
s can be identified via the map of F
(
z
). This is the reason why it is possible to detect the locations of small inhomogeneities as well as random scatterers. Note that, for a successful detection, based on the hypothesis in Theorem 3, the value of N (at least, greater than 3
M
+
3
S) and ω must be sufficiently large enough. If applied frequency is low or total number of N is small, poor result would appear in the map of F
(
x
).Remark 5 (discrimination of singular values).
Theoretically, if the size, permittivity, and permeability of the random scatterers are smaller than those of the inhomogeneities, thenσ
s
<
σ
m for all m and s. This means that if it were possible to discriminate singular values associated with small inhomogeneities, then the structure of F
(
x
) would become(38)
F
x
≈
1
-
∑
m
=
1
M
J
0
ω
x
-
z
m
-
∑
m
=
1
M
∑
h
=
1
2
x
-
z
m
·
e
h
x
-
z
m
2
J
1
ω
x
-
z
m
2
-
1
/
2
.Hence, it is expected that more good results can be obtained. Our approach presents an improvement. However, if the relation σ
s
<
σ
m were no longer valid, the locations of random scatterers would have to be identified via MUSIC such that poor results would appear in the map of F
(
x
).


## 5. Results of Numerical Simulations

Selected results of numerical simulations are presented here to support the identified structure of the MUSIC-type imaging function. In this section, we only consider the dielectric permittivity contrast case; that is, we setɛ
m
=
3, ɛ
0
=
1, and μ
m
=
μ
s
=
μ
0 for all m and s. The radius of all Σ
m and Δ
s is set to 0.1 and 0.05, respectively. The applied angular frequency is ω
=
2
π
/
λ and a total of N number of incident directions is applied such that(39)
θ
j
=
-
cos
⁡
2
π
j
-
1
N
,
2
π
j
-
1
N
T
,
j
=
1,2
,
…
,
N
.M
=
3 small inhomogeneities are selected with locations z
1
=
[
0.25,0
]
T, z
2
=
[
-
0.4,0.5
]
T, and z
3
=
[
-
0.3
,
-
0.7
]
T. We set S
=
100 number of small scatterers as being randomly distributed in Ω
=
[
-
1,1
]
×
[
-
1,1
]
⊂
R
2 such that(40)
y
s
=
η
1
-
1,1
,
η
2
-
1,1
Tfor all s and also select the permittivities randomly as(41)
ɛ
s
=
η
3
1,2
,where η
p
(
a
,
b
), p
=
1,2, and 3, is an arbitrary real value within [
a
,
b
]. Refer to Figure 1 for a sketch of the distribution of the three inhomogeneities and random scatterers.Figure 1
Distribution of inhomogeneities (red-colored dots) and random scatterers (blue-colored “×” mark).The far-field elements of MSR matrixK are generated by means of the Foldy-Lax framework to avoid aninverse crime. After the generation, a singular value decomposition of K is performed via the MATLAB command svd. The nonzero singular values of K are discriminated as follows: first, a 0.1-threshold scheme (by first choosing the j singular values σ
j such that σ
j
/
σ
1
≥
0.1) is applied based on [18], and second, the first 3-singular values are selected.Figure2 exhibits the distribution of the normalized singular values of K and maps of F
(
x
) with the 0.1-threshold scheme and with selection of the first 3-singular values when λ
=
0.3 and N
=
32. Note that due to the huge number of artifacts it is very hard to identify the locations of Σ
m with the 0.1-threshold scheme but, fortunately in this example, one can discriminate three nonzero singular values such that, based on Remark 5, the locations of Σ
m can be identified more clearly. This result supports the derived mathematical structure in Theorem 3.Figure 2
Distribution of normalized singular values (a, c) and maps ofF
(
x
) with first 3-singular values (b) and with 0.1-threshold scheme (d).
(a)
(b)
(c)
(d)Now, let us examine the effect of total number of directionsN in the extreme cases. Figure 3 exhibits normalized singular values and map of F
(
x
) with small number of N
=
5 when λ
=
0.4. Based on Remark 4, the value of N must be sufficiently large so, as we expected, locations of Σ
m cannot be identified via the map of F
(
x
) with small N.Figure 3
Distribution of normalized singular values (a) and map ofF
(
x
) with first 2-singular values (b).
(a)
(b)Opposite to the previous result, Figure4 displays normalized singular values and maps of F
(
x
) with large number of N
=
256 when λ
=
0.4. Similar to the results in Figure 2, locations of Σ
m can be examined clearly via the selection of first 3-singular values. Applying 0.1-threshold, it is very hard to identify locations of Σ
m but, opposite to the result in Figure 2, their locations can be recognized even though some artifacts still exist.Figure 4
Distribution of normalized singular values (a, c) and maps ofF
(
x
) with first 3-singular values (b) and with 0.1-threshold scheme (d).
(a)
(b)
(c)
(d)On the basis of recent works [13, 20], it has been confirmed that MUSIC is robust with respect to the random noise. In order to examine the robustness, assume that 10
d
B Gaussian random noise is added to the unperturbed data u
f
a
r
(
ϑ
j
,
θ
l
). Throughout results in Figure 5 when N
=
32 and λ
=
0.3, although some blurring appears in the map of F
(
x
), we can easily find proper singular values and obtain an accurate image. It is interesting to observe that, opposite to the results in Figure 2, locations of Σ
m can be detected despite existence of some artifacts.Figure 5
Distribution of normalized singular values (a, c) and maps ofF
(
x
) with first 3-singular values (b) and with 0.1-threshold scheme (d) when N
=
32 and λ
=
0.3, and collected far-field data is perturbed by a white Gaussian random noise.
(a)
(b)
(c)
(d)From the above results, we can examine that, by having small perturbations of random scatterersΔ
s, their effects to the scattered fields are quite small so that Σ
m can be discriminated very accurately. Opposite to the this examination, let us consider the effect of Δ
s when their size and permittivities satisfy r
s
=
0.1 and ɛ
s
=
η
(
2.5,3
), respectively (remember that r
m
=
0.1 and ɛ
m
≡
3 for all m). In this example, it is very hard to discriminate nonzero singular values associated with Σ
m so that it is impossible to detect their exact locations; refer to Figure 6 when N
=
32 and λ
=
0.4.Figure 6
Distribution of normalized singular values (a) and map ofF
(
x
) with 0.1-threshold scheme (b).
(a)
(b)It is well-known that using multifrequency improves the imaging performance; refer to [13, 32–34]. At this moment, we consider multifrequency MUSIC-type imaging in order to compare the imaging performance against the traditional single-frequency one. For given F-different frequencies 0
<
ω
1
<
ω
2
<
⋯
<
ω
F, SVD of MSR matrix K
(
ω
f
) is(42)
K
ω
f
≈
∑
m
=
1
3
M
σ
m
ω
f
U
m
ω
f
V
m
∗
ω
f
+
∑
s
=
3
M
+
1
3
M
+
3
S
σ
s
ω
f
U
s
ω
f
V
s
∗
ω
f
.Then, by choosing test vector(43)
f
x
;
ω
f
=
1
N
e
i
ω
f
θ
1
·
x
,
e
i
ω
f
θ
2
·
x
,
…
,
e
i
ω
f
θ
N
·
x
T
,we can survey the projection operator onto the null (or noise) subspace such that(44)
P
n
o
i
s
e
f
x
;
ω
f
≔
I
N
-
∑
m
=
1
3
M
+
3
S
U
m
ω
f
U
m
∗
ω
f
f
x
;
ω
f
,and correspondingly multifrequency MUSIC-type imaging function Q
(
x
;
F
) can be introduced as(45)
Q
x
;
F
=
1
F
∑
f
=
1
F
P
n
o
i
s
e
f
x
;
ω
f
-
1
.Figure7 shows maps of Q
(
x
;
10
), where ω
f
=
2
π
/
λ
s. Here, N
=
32 directions are applied and λ
f are equidistributed in the interval [
λ
F
,
λ
1
] with λ
1
=
0.7 and λ
F
=
0.3. By comparing results in Figure 2, we can observe that unexpected artifacts have been eliminated so that applying multiple frequencies yields a more accurate result and then single frequency.Figure 7
Maps ofQ
(
x
;
10
) with first 3-singular values (a) and with 0.1-threshold scheme (b).
(a)
(b)


## 6. Concluding Remarks

The mathematical structure of MUSIC-type imaging function is carefully identified by establishing a relationship with integer ordered Bessel functions. This is based on the fact that the elements of the MSR matrix can be expressed by an asymptotic expansion formula. The identified structure explains some unexplained phenomena and provides a method for improvements.Based on recent work [7], the electric field E in the existence of small inhomogeneity with radius r can be expressed as follows:(46)
E
x
=
E
0
x
+
r
3
B
m
∑
m
=
1
M
k
2
3
ɛ
m
-
ɛ
0
ɛ
m
+
2
ɛ
0
G
x
,
z
m
·
E
0
z
m
-
i
ω
μ
0
3
μ
m
-
μ
0
μ
m
+
2
μ
0
∇
×
G
x
,
z
m
·
H
0
z
m
+
O
r
4
,where electromagnetic fields (
E
0
,
H
0
) are the solutions of the Maxwell equations(47)
∇
×
E
0
=
i
ω
μ
0
H
0
in
R
3
,
∇
×
H
0
=
-
i
ω
ɛ
0
E
0
+
J
0
in
R
3
,
lim
r
→
∞
⁡
r
∇
×
E
0
-
i
k
r
r
×
E
0
=
0
,
lim
r
→
∞
⁡
r
∇
×
H
0
-
i
k
r
r
×
H
0
=
0and G is Green’s function(48)
G
x
,
z
m
≔
1
0
0
0
1
0
0
0
1
+
∇
∇
k
2
e
i
k
x
-
z
m
4
π
x
-
z
m
.Thus, by applying above asymptotic expansion formula and through the similar process in Theorem 3, the result in this paper can be extended to the three-dimensional problem so that MUSIC will be applicable for detecting three-dimensional inhomogeneities surrounded by random scatterers.In comparison with the MUSIC, other closely related reconstruction algorithms such as linear sampling method [35–37], subspace migration [32, 33, 38], and direct sampling method [39–41] will be applicable for detecting inhomogeneities in random medium. Analysis of imaging functions and exploring their certain properties will be the forthcoming work.



---

*Source: 7872548-2016-02-02.xml*