- Email: [email protected]
Inverse electromagnetic scattering for a locally perturbed perfectly conducting plate
Wave Motion 60 (2016) 95–107
Contents lists available at ScienceDirect
Wave Motion journal homepage: www.elsevier.com/locate/wavemoti
Inverse electromagnetic scattering for a locally perturbed perfectly conducting plate Heping Dong a,∗ , Fuming Ma a , Xiaokai Yuan a , Yukun Guo b a
Institute of Mathematics, Jilin University, Changchun, 130012, China
b
Department of Mathematics, Harbin Institute of Technology, Harbin, 150001, China
highlights • The linear sampling method is applied to reconstruct the shape of a local perturbation of a perfectly conducting plate from far-field measurements of time harmonic electromagnetic fields.
• The linear sampling method with near-field data is also discussed. • Numerical examples in case of far-field data, limited aperture data and near-field data are included to illustrate the feasibility of this method.
article
info
Article history: Received 9 January 2015 Received in revised form 17 July 2015 Accepted 14 September 2015 Available online 25 September 2015 Keywords: Locally rough surface Inverse electromagnetic scattering problem Maxwell’s equations Linear sampling method
abstract We consider a simple but fully three-dimensional inverse problem to determine the shape of a local perturbation of a perfectly conducting plate from far-field measurements of time harmonic electromagnetic fields. For this purpose we reformulate the model problem as an exterior Maxwell problem for a symmetric domain, and prove an equivalence between the model problem and its reformulation. Then, linear sampling method is applied to solve the reformulated problem. We illustrate the feasibility of this method by some numerical examples. © 2015 Elsevier B.V. All rights reserved.
1. Introduction In recent years, there has been increased interest in inverse electromagnetic scattering from unbounded obstacles. It can be classified roughly into four classes: locally rough surface problem, rough surface problem, cavity problem and grating problem. These problems arise from many applications such as radar, remote sensing, medical imaging and nondestructive testing. In this paper, our attention is restricted to the discussion of locally rough surface problem. In [1], for the local perturbation that is both star-like and above the infinite plane case, the scattering problem is reformulated as an exterior Dirichlet problem for a symmetric planar domain with corners, and then a Newton method is proposed to reconstruct the local rough surface from the far-field pattern. In [2], the scattering problem is reformulated as an equivalent boundary value problem in a bounded domain with a DtN map, and a continuation method is applied to reconstruct the local disturbance with multiple frequency data. Both of the papers [3,4] are concerned with the near-field imaging of locally rough surface.
∗
Corresponding author. E-mail address: [email protected] (H. Dong).
http://dx.doi.org/10.1016/j.wavemoti.2015.09.004 0165-2125/© 2015 Elsevier B.V. All rights reserved.
96
H. Dong et al. / Wave Motion 60 (2016) 95–107
Based on the plane wave decomposition of the scattered wave and boundary integral equations, an optimization method is introduced in [3] to identify the profile of a reflective local rough surface by the scattered wave measured in the near field. In [4], the surface displacement is reconstructed by solving a minimization problem from the near-field data under the condition that the local surface displacement is above an infinite ground plane. In [5], an integral equation formulation defined on a bounded curve with two corners is designed to solve the direct scattering problem, and a Newton method with multifrequency far field data is applied to reconstruct the local perturbed surface. In contrast to [2], the integral equation formulation does not involve any infinite integral or a DtN map, and therefore it can be solved by using the Nyström method efficiently. Compared with [1], the authors in [6] make use of linear sampling method [7,8] to solve the inverse problem under the condition that the symmetric domain is with no corners. Then, [9] gives the uniqueness of the inverse scattering problem by using the generalization of Green’s theorem to nonsmooth domain. The well-posedness of direct scattering problems is established in [2,10–13]. For the rough surface problems we refer to [14–16]. For the inverse cavity problems, the reader is referred to [17,18]. For the grating problems, we refer to [19] and references therein. Motivated by [1,6] for the two dimensional case, we are currently trying to extend the technique to a simple but fully three dimensional model for the electromagnetic exploration under the condition that the local perturbation is above an infinite perfectly conducting plate. We reformulate the scattering problem as an exterior Maxwell problem for a symmetric domain, and then use linear sampling method to reconstruct the shape of local perturbation. For the application of linear sampling method in inverse electromagnetic scattering problems, we refer to [20]. This paper is organized as follows. In Section 2 we describe the model problem in detail, and transform the model problem for unbounded obstacle into a new scattering problem for bounded obstacle by symmetric continuation. Then, we prove an equivalence between the direct scattering problem and its reformulation. Section 3 formulates the inverse scattering problem, and describes the implementation details of linear sampling method. Finally, we show some numerical examples in this section. 2. Presentation of the forward and inverse scattering problem In this section, a direct scattering problem is considered for Maxwell’s equations within the upper half-space R3+ = {x ∈ R : x3 > 0}. We decompose the space R3 = R3+ ∪ R30 ∪ R3− in a hyperplane R30 = {x ∈ R3 : x3 = 0} corresponding to the surface where the local perturbation above on, and R3− = {x ∈ R3 : x3 < 0} corresponding to the lower half-space. Let D+ ⊂ R3+,0 := {x ∈ R3 : x3 ≥ 0} be a bounded domain such that ∂ D+ ⊂ R3+ is a smooth surface. Suppose ΣD,0 = D+ ∩ R30 is also bounded such that the boundary ∂ ΣD,0 is a smooth closed curve. Denote by (e1 , e2 , e3 ) the usual Cartesian basis of R3 , by x = (x1 , x2 , x3 )⊤ a point in R3 , and by ν+ the unit outward normal to the boundary. Throughout let x · y and x × y be the scalar product and the vector product of x, y ∈ R3 , respectively, and let |x| denote the Euclidean norm of x. For x ∈ R3 we denote by 3
T1 x =
1 0 0
0 1 0
0 0 x = (x1 , x2 , −x3 )⊤ −1
its reflection about the plane R30 , and define T2 x = −T1 x = (−x1 , −x2 , x3 )⊤ . For convenience and also for later use in this section, we give that (T1 x) · y = x · (T1 y ), (T1 x) × (T1 y ) = T2 (x × y ) and (T2 x) × (T1 y ) = T1 (x × y ). Given as incident fields i
∇ × ∇ × pexp(ikx · d ) = ik(d × p) × dexp(ikx · d ), k H i (x) = ∇ × pexp(ikx · d ) = ikd × pexp(ikx · d ), E i (x) =
where the wave number k is positive, d is a unit vector describing the direction of propagation, and p is the polarization vector. The solution of scattering problem for the half-space, i.e. absence of D+ , is given by the reflected wave E r (x) = T2 E i (T1 x), H r (x) = T1 H i (T1 x). For convenience, we denote E i,r = E i + E r , H i,r = H i + H r and f+ := −ν+ × E i,r = 0, for x ∈ R30 . Then the direct scattering problem for the scattering of a time harmonic electromagnetic plane wave by a locally s s perturbed perfectly conducting surface ∂ D+ is to find an electric field E+ = E i,r + E+ and a magnetic field H+ = H i,r + H+ such that
∇ × E+ − ikH+ = 0,
in R3+ \ D+ ,
(1)
∇ × H+ + ikE+ = 0,
in R+ \ D+ ,
(2)
ν+ × E+ = 0,
\ ΣD,0 ) ∪ ∂ D+ .
(3)
on (
R30
3
s s The unknown scattered fields E+ and H+ are required to satisfy the Silver–Müller radiation condition s lim r (H+ × xˆ − E+s ) = 0,
r →∞
uniformly in xˆ = x/|x| ∈ ∂ B := {x ∈ R3 : |x| = 1}, r = |x|. The geometry of this model problem is shown in Fig. 2.1.
(4)
H. Dong et al. / Wave Motion 60 (2016) 95–107
97
Fig. 2.1. Geometry of the model problem.
The corresponding inverse scattering problem can be stated as: for given wave number k > 0 and incident field (E i , H i ), ∞ reconstruct the unknown boundary ∂ D+ , by use of the electric far field pattern E+ (ˆx, d , p), ∀ˆx ∈ ∂ B+ , ∀d ∈ ∂ B− , ∀p ∈ R3 . In order to reformulate this direct scattering problem in half-space as an exterior Maxwell problem, we denote by ∂ D− = {x ∈ R3 : T1 x ∈ ∂ D+ } the reflection of ∂ D+ about R30 . Then ∂ D− ∪ ∂ D+ is the boundary ∂ D of a bounded domain D. Let E s , H s be the radiating solution to the Maxwell equations satisfying
∇ × E s − ikH s = 0,
¯, in R3 \ D
(5)
∇ × H + ikE = 0, in R3 \ D¯ , ν × E s = f , on ∂ D, lim r (H s × xˆ − E s ) = 0,
(6)
s
s
(7) (8)
r →∞
where f+ (x),
x3 ≥ 0,
T1 f+ (T1 x),
x3 < 0,
f (x) =
and
ν(x) =
ν+ (x), T1 ν+ (T1 x),
x3 ≥ 0 , x3 < 0.
Then we can prove that Theorem 1. Problem (1)–(4) is equivalent to problem (5)–(8). s Proof. Assume that (E+ , H+s ) is a radiating solution to the problem (1)–(4), and extend (E+s , H+s ) by s E+ (x),
x3 ≥ 0,
s T2 E+ (T1 x),
x3 < 0,
E (x) = s
s H+ (x),
x3 ≥ 0 ,
s T1 H + (T1 x),
x3 < 0.
H (x) = s
and
Straightforward calculations show that (E , H ) is the solution of the problem (5)–(8). On the other hand, if E s (x), H s (x) solves the problem (5)–(8), we then obtain that T2 E s (T1 x), T1 H s (T1 x) also solves the problem (5)–(8). From the uniqueness for solutions to exterior Maxwell boundary value problems [21] we conclude that s E s (x) = T2 E s (T1 x), which implies that the first two components of vector E s vanish at x ∈ R30 . Consequently, ν+ × E+ = 0, 3 s s for x ∈ R0 . Thus it can be verified that (E+ , H+ ) is the solution of the problem (1)–(4). s
s
Theorem 2. The electric far field pattern E∞ (ˆx, d, p) for the scattering problem (5)–(8) satisfies the following relations (1) E∞ (T1 xˆ , d, p) = T2 E∞ (ˆx, d, p); (2) E∞ (ˆx, T1 d, p) = −E∞ (ˆx, d, T1 p); (3) E∞ (T1 xˆ , T1 d, p) = T1 E∞ (ˆx, d, T1 p). Proof. (1) Since by the Huygens’ principle ([21, Th 6.24]), the far field pattern is given by E (ˆx) = ∞
ik 4π
xˆ × ∂D
ν(y) × H (y) × xˆ exp(−ikxˆ · y)ds(y).
Define R(ˆx) =
∂D
Q (y, xˆ )exp(−ikxˆ · y)ds(y),
where Q (y, xˆ ) = P (y) × xˆ , P (y) = ν(y) × H (y). Now the far field pattern is written as E ∞ (ˆx) =
ik 4π
xˆ × R(ˆx).
98
H. Dong et al. / Wave Motion 60 (2016) 95–107
Noting that H i,r (T1 x) = T1 H i,r (x) and H s (T1 x) = T1 H s (x), we have H (T1 x) = T1 H (x). Using the symmetry of scatterer D, we obtain that ν(T1 x) = T1 ν(x). Hence, E ∞ (T1 xˆ ) =
= = = = = =
ik 4π ik 4π ik 4π ik 4π ik 4π ik 4π ik 4π
(T1 xˆ ) × (T1 xˆ ) ×
ν(y) × H (y) × (T1 xˆ )exp −ik(T1 xˆ ) · y ds(y) ∂D
T1 ν(T1 y) × T1 H (T1 y) × (T1 xˆ )exp −ikxˆ · (T1 y) ds(y)
∂D
× (T1 xˆ )exp −ikxˆ · (T1 y) ds(y) T2 ν(T1 y) × H (T1 y) (T1 xˆ ) × ∂D
(T1 xˆ ) ×
(T1 xˆ ) ×
∂D
∂D
(T1 xˆ ) × T1 R(ˆx)
T1 Q (T1 y, xˆ )exp −ikxˆ · (T1 y) ds(y)
T2 xˆ × R(ˆx)
T2 P (T1 y) × (T1 xˆ )exp −ikxˆ · (T1 y) ds(y)
= T2 E ∞ (ˆx). (2) In terms of the incident field E i and the reflected wave E r , we have E i (x, T1 d , T1 p) = T1 E i (T1 x, d , p),
E r (x, T1 d , T1 p) = −E i (x, d , p),
and T1 E i,r (T1 x, d , p) = E i,r (x, T1 d , T1 p). Since E i,r (T1 x, d , p) = T2 E i,r (x, d , p), we get E i,r (x, T1 d , T1 p) = −E i,r (x, d , p). According to the boundary condition (7), one can derive that
ν × E s (x, T1 d , T1 p) = −ν × E i,r (x, T1 d , T1 p) = −ν × −E i,r (x, d , p) = ν × E i,r (x, d , p). Thus,
ν × −E s (x, T1 d , T1 p) = −ν × E i,r (x, d , p) = ν × E s (x, d , p). From the uniqueness for solutions to exterior Maxwell boundary value problems [21] we get
−E s (x, T1 d , T1 p) = E s (x, d , p). Consequently, the corresponding far field pattern satisfies E ∞ (ˆx, d , p) = −E ∞ (ˆx, T1 d , T1 p). (3) Combining the previous results (1) and (2), we finally obtain the following results
T1 E ∞ (T1 xˆ , T1 d , T1 p) = T1 −E ∞ (T1 xˆ , d , p)
= T2 E ∞ (T1 xˆ , d , p) = E ∞ (ˆx, d , p). The proof of the theorem is now completed.
∞ This theorem tells us that for given the far field pattern E+ (ˆx, d , p), for xˆ ∈ ∂ B+ , d ∈ ∂ B− , the far field pattern
E ∞ (ˆx, d , p),
∀ˆx ∈ ∂ B, ∀d ∈ ∂ B,
corresponding the scattered field E s of problem (5)–(8) can be obtained. Hence, we are able to deal with this problem by using linear sampling method.
H. Dong et al. / Wave Motion 60 (2016) 95–107
99
3. Linear sampling method and numerical implementation To analyze the inverse scattering problem we consider the following Hilbert spaces: H (curl, D) := {a ∈ (L2 (D))3 : curla ∈ (L2 (D))3 }, − 21
Ht
1 3 (∂ D) := a ∈ H − 2 (∂ D) : ν · a = 0 ,
L2t (∂ B) := {a ∈ (L2 (∂ B))3 : a · xˆ = 0, xˆ ∈ ∂ B}, H
− 12
− 21
(Div, ∂ D) := a ∈ Ht
(∂ D) : Diva ∈ H
− 12
(∂ D)
with norm
∥a∥2H −1/2 (Div,∂ D) := ∥a∥2(H −1/2 (∂ D))3 + ∥Diva∥2H −1/2 (∂ D) . ¯ ) denotes the space of all fields a : R3 \ D¯ → C3 such that a ∈ H (curl, (R3 \ D¯ ) ∪ BR ) for Furthermore, Hloc (curl, R3 \ D all open balls BR containing the closure of D. 3.1. Far-field data We are now in a position to describe the numerical algorithm. To this end, we introduce the far field operator F : L2t (∂ B) → L2t (∂ B) such that
(F g )(ˆx) :=
∂B
E ∞ (ˆx, d , g (d ))ds(d ),
xˆ ∈ ∂ B.
(9)
The linear sampling method is based on seeking a function g (·, z , p) ∈ L2t (∂ B), such that for fixed p and z,
F g (·, z , p) (ˆx) = Ee∞ (ˆx, z , p)
=
ik 4π
(ˆx × p) × xˆ exp(−ikxˆ · z ) −
ik 4π
(ˆx × T1 p) × xˆ exp(−ikxˆ · T1 z ),
(10)
where Ee∞ is the electric far field pattern of the electric dipole Ee . In contrast to [22], the electric dipole Ee has to be modified to satisfy ν × Ee = 0 on R30 . In this case, we have Ee (x, z , p) =
i k
i ∇x × ∇x × pΦ (x, z ) − ∇x × ∇x × T1 pΦ (x, T1 z ) , k
where
Φ (x, z ) =
1 exp(ik|x − z |) 4π
|x − z |
,
x, z ∈ R3 , x ̸= z . 1
Moreover, we introduce a bounded linear operator A : H − 2 (Div, ∂ D) → L2t (∂ B) which maps the boundary data f ∈H
− 21
(Div, ∂ D) to the far field pattern E ∞ of the radiating solution E s ∈ Hloc (curl, R3 \ D¯ ) of (5)–(8). Then we have that
F g = −A(ν × Eg ),
(11)
where Eg (x) = ik
∂B
eikx·d g (d )ds(d ) + ik
∂B
eik(T1 x)·d T2 g (d )ds(d ).
We can now prove the main result of this paper. Theorem 3. Assume that k is not a Maxwell eigenvalue for D and let F be the far field operator (9) for the scattering problem (5)–(8). Then the following hold: (1) For z ∈ D and a given ε > 0 there exists a gεz ∈ L2t (∂ B) such that
∥F gεz − E∞ e (·, z, p)∥L2 (∂ B) < ε t
(12)
and the corresponding Herglotz function Egεz converges to the solution of interior Maxwell problem
∇ × E − ikH = 0, ν × (E + Ee ) = 0, in H (curl, D) as ε → 0.
∇ × H + ikE = 0 in D, on ∂ D,
(13) (14)
100
H. Dong et al. / Wave Motion 60 (2016) 95–107
¯ and a given ε > 0, every gεz ∈ L2t (∂ B) that satisfies (12) is such that (2) For z ∈ R3 \ D lim ∥gεz ∥L2 (∂ B) = ∞.
ε→0
t
Proof. For convenience, we denote
∇x × ∇x × pΦ (x, z ) , k Eg1 (x) = ik eikx·d g (d )ds(d ), i
Ee1 (x, z , p) =
i Ee2 (x, z , p) = − ∇x × ∇x × k Eg2 (x) = ik
∂B
∂B
T1 pΦ (x, T1 z ) ,
eik(T1 x)·d T2 g (d )ds(d ).
Then Ee = Ee1 + Ee2 , Eg = Eg1 + Eg2 . Under the assumption on k we have well-posedness of the interior Maxwell problem
in H (curl, D). Given ε > 0, from the denseness of Herglotz operator(Corollary 7.19 in [21]) we can choose gzε ∈ L2t (∂ B) such that
ν(x) × Eg1zε (x, z , p) + Ee1 (x, z , p)
H −1/2 (Div,∂ D)
<
ε . 2∥A∥
Define W (x, z , p) = ν(x) × Eg1ε (x, z , p) + Ee1 (x, z , p) , the above equation implies
z
∥W (·, z , p)∥H −1/2 (Div,∂ D) <
ε 2∥A∥
.
(15)
Straightforward calculations show that
ν(x) × Eg2zε (x, z , p) + Ee2 (x, z , p) = T1 ν(T1 x) × T2 Eg1zε (T1 x, z , p) + T2 Ee1 (T1 x, z , p) = T2 ν(T1 x) × Eg1zε (T1 x, z , p) + Ee1 (T1 x, z , p) = T2 W (T1 x, z , p). With the aid of (15) and the symmetry of scatterer D, it can be verified that
∥T2 W (T1 ·, z , p)∥H −1/2 (Div,∂ D) = ∥W (·, z , p)∥H −1/2 (Div,∂ D) <
ε 2∥A∥
.
Thus, we get
ν(x) × Egzε (x, z , p) + Ee (x, z , p)
H −1/2 (Div,∂ D)
≤ ∥W (·, z , p)∥H −1/2 (Div,∂ D) + ∥T2 W (T1 ·, z , p)∥H −1/2 (Div,∂ D) <
ε . ∥A∥
Then (12) follows by observing Eq. (11). Now if z ∈ D, from the well-posedness of the interior Maxwell problem,
lim ν(x) × Egzε (x, z , p) + Ee (x, z , p)
ε→0
H −1/2 (Div,∂ D)
=0
implies lim Egzε − E H (curl,D) = 0,
ε→0
where E is the solution to (13)–(14). By arguments similar to those in Theorem 7.21 of [21], the second statements can be established by noting (15).
In the following, we briefly describe the discretization of the far field equation and the details may be found in [20,22,23]. Let (ˆx, e1 (ˆx), e2 (ˆx)) denote an orthonormal basis of R3 with xˆ ∈ ∂ B. Note that E ∞ (ˆx, d , p) · xˆ = 0, the far field equation (10) is then equivalent to
∞ ∞ e1 (ˆx) · E (ˆx, d , g (d , z , p))ds(d ) = e1 (ˆx) · Ee (ˆx, z , p), ∂ B e2 (ˆx) · E ∞ (ˆx, d , g (d , z , p))ds(d ) = e2 (ˆx) · Ee∞ (ˆx, z , p). ∂B
By using the reciprocity relation p · E ∞ (ˆx, d , q) = q · E ∞ (−d , −ˆx, p),
∀ˆx, d ∈ ∂ B, ∀p, q ∈ C3 ,
(16)
H. Dong et al. / Wave Motion 60 (2016) 95–107
101
and the numerical quadrature formula
∂B
f (d )ds(d ) ≈
N
ωj f (dj ),
j =1
where ωj are the weights associated with the mesh triangular meshing of the unit sphere ∂ B, dj are the vertices of mesh triangles. The system (16) can be transformed into the following discrete version
N ωj E ∞ (−dj , −ˆx, e1 (ˆx)) · g (dj , z , p) = e1 (ˆx) · Ee∞ (ˆx, z , p), j =1
N ωj E ∞ (−dj , −ˆx, e2 (ˆx)) · g (dj , z , p) = e2 (ˆx) · Ee∞ (ˆx, z , p). j =1
In view of the fact that the unknown g (dj , z , p) is a tangential vector field, we can set g (dj , z , p) = g1 (dj , z , p)e1 (dj ) + g2 (dj , z , p)e2 (dj ) with g1 (dj , z , p) and g2 (dj , z , p) ∈ C. Thus, we obtain the fully discrete linear equations
N 2 ωj gl (dj , z , p)E ∞ (−dj , −ˆx, e1 (ˆx)) · el (dj ) = e1 (ˆx) · Ee∞ (ˆx, z , p), j =1 l =1
(17)
N 2 ωj gl (dj , z , p)E ∞ (−dj , −ˆx, e2 (ˆx)) · el (dj ) = e2 (ˆx) · Ee∞ (ˆx, z , p). j =1 l =1
In our numerical experiments, the far field pattern E ∞ (−dj , −ˆx, el (ˆx)), l = 1, 2 for N measurement points {ˆxj }Nj=1 is computed by boundary integral equation method [24]. These measurement points are roughly uniformly distributed on the unit sphere ∂ B, and they will also serve as incident directions {dj }Nj=1 . Thus, the quadrature weights ωj can be seen as a
constant. Let qˆ be some unit vector in R3 satisfying qˆ × xˆ j ̸= 0 for any j. In our examples, we take qˆ = (1, 1, 1)⊤ . Then we set the two orthogonal polarizations e1 (ˆxj ) = qˆ × xˆ j /|ˆq × xˆ j |
and e2 (ˆxj ) = qˆ × (ˆxj × qˆ )/|ˆq × (ˆxj × qˆ )|.
Now, the system (17) can be rewritten as A∞ g⃗ (z , p) = ⃗ f (z , p)
(18)
where A∞ is a 2N × 2N complex matrix corresponding to the given far field data, vector g⃗ = (g1 , g2 ) unknowns and ⃗ f is a vector given by the right hand side of (17). In order to avoid inverse crimes, we corrupt the matrix A∞ with random noise
⊤
represents the
A∞ := A∞ (I + ϵ R1 + iϵ R2 ) where R1 and R2 are two matrices of real random numbers in [−1, 1] and ϵ is a parameter. In the following examples, we have chosen ϵ = 0.01 and ϵ = 0.1. Since the ill-posedness of this problem, instead of Eq. (18), we solve the regularized equation
γ g⃗ + A∗∞ A∞ g⃗ = A∗∞ ⃗f
(19)
is the adjoint of A∞ and γ is the regularization parameter. Let p1 = (1, 0, 0) , p2 = (0, 1, 0) , p3 = (0, 0, 1)⊤ . For each z in the sampling grid, we compute the indicator function
where A∗∞
G(z ) =
⊤
1
∥⃗g (z , p1 )∥−2 + ∥⃗g (z , p2 )∥−2 + ∥⃗g (z , p3 )∥−2
21
⊤
.
3 The reconstruction corresponds to the isosurfaces of G(z ) close to zero. Consequently, we will choose a constant C and plot surfaces where G(z ) = C . In this section, we present and discuss three numerical examples. In each example, we use 50 and 100 incoming waves respectively to reconstruct the shape of local perturbation, and these incoming waves are generated by Gmsh which is a 3D finite element grid generator. Because of the symmetry of this problem, we only need to put sampling points in upper space. Example 1. In the first example, we consider a specific case in which the shape of local perturbation ∂ D+ is a half sphere with radius R = 1. The wave number is k = 3. We take (41)3 sampling points within a cube prescribed as Ω = [−2, 2] × [−2, 2] × [0.01, 1.5]. As can be seen from Fig. 3.1, the local perturbation can be reconstructed when illuminations and measurements N = 50 and N = 100 respectively, provided a good isosurface of G(z ) and regularization parameter are chosen.
102
H. Dong et al. / Wave Motion 60 (2016) 95–107
Fig. 3.1. Reconstruction of the half ball with radius R = 1 and the exact surface.
Example 2. For the second example, we wish to reconstruct a half cushion-shaped surface with the parametrization z (θ , ϕ) =
0.27 + 0.065(cos 2ϕ − 1)(cos 4θ − 1)r (θ , ϕ),
where spherical coordinates r (θ , ϕ) = (sin θ cos ϕ, sin θ sin ϕ, cos θ ),
θ ∈ [0, π /2], ϕ ∈ [0, 2π ].
We choose wave number k = 4 and (41) sampling points within Ω = [−1, 1] × [−1, 1] × [0.01, 0.8]. Fig. 3.2 shows the results. 3
Example 3. Our third example is intended to reconstruct a half bean-like surface with the parametrization
0.5 1 − 0.1 cos(π cos θ ) sin θ cos ϕ , z (θ , ϕ) = 0.5 sin θ sin ϕ + 0.3 cos(π cos θ ) 0.7 cos θ
θ ∈ [0, π /2], ϕ ∈ [0, 2π ].
We assume wave number k = 7 and the sampling domain is Ω = [−1, 1] × [−1, 1] × [0.01, 1] with (41)3 sampling points. Results are shown in Fig. 3.3. Moreover, one can observe how the increase of N helps a better resolution of the reconstruction. The above numerical results illustrate that the linear sampling method gives a feasible reconstruction of the locally perturbed perfectly conducting plate even in the presence of 10% noise in measurements.
H. Dong et al. / Wave Motion 60 (2016) 95–107
103
Fig. 3.2. Reconstruction of the half cushion and the exact surface.
3.2. Limited aperture data In this section, we restrict the set of observation points to only half of the upper half sphere Γ ⊂ ∂ B. This means we know the half of measurement data. In this case the data are E ∞ (ˆxi , dj , p) for some (ˆxi , dj , p) ∈ Γ × (−Γ ) × R3 . Figs. 3.4 and 3.5 display reconstructions that correspond to the three examples. We observe that the quality of the reconstruction deteriorates as the aperture decreases. However, one notices in all cases that we have a rather good approximation of the projection of the scatterer on the plane located in the illuminated region. 3.3. Near-field data In this section we consider the linear sampling method for the near field measurements. Let G+ ⊂ R3+,0 be a bounded
domain enclosing the local perturbation such that ∂ G+ ⊂ R3+ is a smooth surface. The incident field is considered to be the electric field of an electric dipole located at y ∈ ∂ G+ with polarization p ∈ R3 , and is given by Eei (x, y , p) =
i k
∇x × ∇x × pΦ (x, y ) = Ee1 (x, y , p).
The inverse scattering problem we are interested in is to determine the unknown boundary ∂ D+ from knowledge of the s tangential components ν+ ×E+ of the scattered electric field measured on ∂ G+ for all incident point sources Eei (x, y , p), ∀y ∈ 3 ∂ G+ , ∀p ∈ R .
104
H. Dong et al. / Wave Motion 60 (2016) 95–107
Fig. 3.3. Reconstruction of the half bean and the exact surface.
Note that the reflected wave is Eer (x, y , p) = T2 Eei (T1 x, y , p) = Ee2 (x, y , p). Analogously to the case of far-field data, we ¯ with E i,r replaced can reformulate the scattering problem as an exterior Maxwell problem (5)–(8) for a symmetric domain D by Ee . Then, we have the following result. Theorem 4. For x, y ∈ ∂ G, the tangential component ν × Es of the scattered electric field satisfies the following relations (1) (ν × Es )(T1 x, y, p) = T2 (ν × Es )(x, y, p); (2) (ν × Es )(x, T1 y, p) = −(ν × Es )(x, y, T1 p); (3) (ν × Es )(T1 x, T1 y, p) = T1 (ν × Es )(x, y, T1 p). Now we reformulate the linear sampling method in terms of the near field data. To this end we consider the near field ¯ is contained in G, and the near field operator F : L2t (∂ G) → L2t (∂ G) where G = G+ ∪ G− is a bounded domain such that D equation
(F ϕz )(x) :=
∂G
ν(x) × E s (x, y , ϕz ,p (y ))ds(y ) = ν(x) × Ee (x, z , p),
x ∈ ∂ G.
(20)
For convenience, ∂ G is chosen to be a sphere with radius r = 2 in Example 1 and r = 1 in Examples 2 and 3. By using the symmetry relation (Theorem 6.32 in [21]) q · E s (x, y , p) = p · E s (y , x, q),
¯, x, y ∈ R3 \ D
H. Dong et al. / Wave Motion 60 (2016) 95–107
Fig. 3.4. Reconstruction for limited aperture case.
Fig. 3.5. Reconstruction for limited aperture case.
105
106
H. Dong et al. / Wave Motion 60 (2016) 95–107
Fig. 3.6. Reconstruction for near-field data.
the near field equation (20) can be transformed into the fully discrete version
N 2 ωj ϕzl ,p (yj )E s (yj , xi , e2 (ˆxi )) · el (ˆyj ) = e2 (ˆxi ) · Ee (xi , z , p), j =1 l =1
N 2 ωj ϕzl ,p (yj )E s (yj , xi , e1 (ˆxi )) · el (ˆyj ) = e1 (ˆxi ) · Ee (xi , z , p), j =1 l =1
where we have set ϕz ,p (yj ) = ϕz1,p (yj )e1 (ˆyj ) + ϕz2,p (yj )e2 (ˆyj ) and the measurement points {yj }Nj=1 also serve as the locations {xi }Ni=1 of the point sources. Similarly as far-field data case, we plot the isosurfaces of G(z ) for ϕz ,p . In the numerical experiments, we keep the settings of wave number, sampling domain and the number of sampling points in each example. Fig. 3.6 shows the result of reconstruction.
4. Conclusions In this paper we have presented a numerical method for solving a three-dimensional inverse problem to determine the shape of a local perturbation of a perfectly conducting plate. We have applied the linear sampling method to the far-field data case, limited aperture case and near-field data case. The numerical experiments show that satisfactory reconstructions can be obtained in the three cases. Acknowledgments This work was supported by the National Science Foundation of China (Grant Nos. 91123031-3, 11371172, 41474102, 91230119) and the Scientific Research Foundation of Jilin Province (Grant No. 20140520055JH).
References [1] R. Kress, T. Tran, Inverse scattering for a locally perturbed half-plane, Inverse Problems 16 (2000) 1541–1559. [2] G. Bao, J. Lin, Imaging of local surface displacement on an infinite ground plane: the multiple frequency case, SIAM J. Appl. Math. 71 (2011) 1733–1752.
H. Dong et al. / Wave Motion 60 (2016) 95–107
107
[3] G. Bao, J. Lin, Imaging of reflective surfaces by near-field optics, Opt. Lett. 37 (2012) 5027–5029. [4] G. Bao, J. Lin, Near-field imaging of the surface displacement on an infinite ground plane, Inverse Probl. Imaging 7 (2013) 377–396. [5] H.W. Zhang, B. Zhang, A novel integral equation for scattering by locally rough surfaces and application to the inverse problem, SIAM J. Appl. Math. 73 (2013) 1811–1829. [6] L.X. Feng, F.M. Ma, A sampling method for solving inverse scattering problems with a locally perturbed half plane, Northeast. Math. J. 19 (2003) 1–4. [7] D. Colton, A. Kirsch, A simple method for solving inverse scattering problems in the resonance region, Inverse Problems 12 (1996) 383–393. [8] D. Colton, H. Haddar, M. Piana, The linear sampling method in inverse electromagnetic scattering theory, Inverse Problems 19 (2003) S105–S137. [9] G.Z. Yan, Uniqueness of the inverse scattering problem of a locally perturbed half-plane, Comput. Math. Appl. 48 (2004) 411–418. [10] M. Durán, I. Muga, J.C. Nédélec, The Helmholtz equation in a locally perturbed half-plane with passive boundary, IMA J. Appl. Math. 71 (2006) 853–876. [11] M. Durán, I. Muga, J.C. Nédélec, The Helmholtz equation in a locally perturbed half-space with non-absorbing boundary, Arch. Ration. Mech. Anal. 191 (2009) 143–172. [12] S.N. Chandler-Wilde, A.T. Peplow, A boundary integral equation formulation for the Helmholtz equation in a locally perturbed half-plane, ZAMM Z. Angew. Math. Mech. 85 (2005) 79–88. [13] H. Ammari, C. Latiri-Grouz, An integral equation method for the electromagnetic scattering from a scatterer on an absorbing plane, Integral Equations Operator Theory 39 (2001) 159–181. [14] Y. He, P.J. Li, J. Shen, A new spectral method for numerical solution of the unbounded rough surface scattering problem, J. Comput. Phys. 275 (2014) 608–625. [15] S.N. Chandler-Wilde, C. Lines, A time domain point source method for inverse scattering by rough surfaces, Computing 75 (2005) 157–180. [16] C. Burkard, R. Potthast, A multi-section approach for rough surface reconstruction via the Kirsch–Kress scheme, Inverse Problems 26 (2010) 045007. [17] G. Bao, J. Gao, P. Li, Analysis of direct and inverse cavity scattering problems, Numer. Math. Theory Methods Appl. 4 (2011) 419–422. [18] P.J. Li, An inverse cavity problem for Maxwell’s equations, J. Differential Equations 252 (2012) 3209–3225. [19] G. Bao, P.J. Li, Convergence analysis in near-field imaging, Inverse Problems 30 (2014) 085008. [20] F. Cakoni, D. Colton, P. Monk, The Linear Sampling Method in Inverse Electromagnetic Scattering, SIAM, Philadelphia, 2011. [21] D. Colton, R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, third ed., Springer, New York, 2013. [22] P. Monk, Finite Element Methods for Maxwell’s Equations, Oxford, New York, 2003. [23] D. Colton, H. Haddar, P. Monk, The linear sampling method for solving the electromagnetic inverse scattering problem, SIAM J. Sci. Comput. 24 (2002) 719–731. [24] M. Ganesh, S.C. Hawkins, A high-order tangential basis algorithm for electromagnetic scattering by curved surfaces, J. Comput. Phys. 227 (2008) 4543–4562.
Contents lists available at ScienceDirect
Wave Motion journal homepage: www.elsevier.com/locate/wavemoti
Inverse electromagnetic scattering for a locally perturbed perfectly conducting plate Heping Dong a,∗ , Fuming Ma a , Xiaokai Yuan a , Yukun Guo b a
Institute of Mathematics, Jilin University, Changchun, 130012, China
b
Department of Mathematics, Harbin Institute of Technology, Harbin, 150001, China
highlights • The linear sampling method is applied to reconstruct the shape of a local perturbation of a perfectly conducting plate from far-field measurements of time harmonic electromagnetic fields.
• The linear sampling method with near-field data is also discussed. • Numerical examples in case of far-field data, limited aperture data and near-field data are included to illustrate the feasibility of this method.
article
info
Article history: Received 9 January 2015 Received in revised form 17 July 2015 Accepted 14 September 2015 Available online 25 September 2015 Keywords: Locally rough surface Inverse electromagnetic scattering problem Maxwell’s equations Linear sampling method
abstract We consider a simple but fully three-dimensional inverse problem to determine the shape of a local perturbation of a perfectly conducting plate from far-field measurements of time harmonic electromagnetic fields. For this purpose we reformulate the model problem as an exterior Maxwell problem for a symmetric domain, and prove an equivalence between the model problem and its reformulation. Then, linear sampling method is applied to solve the reformulated problem. We illustrate the feasibility of this method by some numerical examples. © 2015 Elsevier B.V. All rights reserved.
1. Introduction In recent years, there has been increased interest in inverse electromagnetic scattering from unbounded obstacles. It can be classified roughly into four classes: locally rough surface problem, rough surface problem, cavity problem and grating problem. These problems arise from many applications such as radar, remote sensing, medical imaging and nondestructive testing. In this paper, our attention is restricted to the discussion of locally rough surface problem. In [1], for the local perturbation that is both star-like and above the infinite plane case, the scattering problem is reformulated as an exterior Dirichlet problem for a symmetric planar domain with corners, and then a Newton method is proposed to reconstruct the local rough surface from the far-field pattern. In [2], the scattering problem is reformulated as an equivalent boundary value problem in a bounded domain with a DtN map, and a continuation method is applied to reconstruct the local disturbance with multiple frequency data. Both of the papers [3,4] are concerned with the near-field imaging of locally rough surface.
∗
Corresponding author. E-mail address: [email protected] (H. Dong).
http://dx.doi.org/10.1016/j.wavemoti.2015.09.004 0165-2125/© 2015 Elsevier B.V. All rights reserved.
96
H. Dong et al. / Wave Motion 60 (2016) 95–107
Based on the plane wave decomposition of the scattered wave and boundary integral equations, an optimization method is introduced in [3] to identify the profile of a reflective local rough surface by the scattered wave measured in the near field. In [4], the surface displacement is reconstructed by solving a minimization problem from the near-field data under the condition that the local surface displacement is above an infinite ground plane. In [5], an integral equation formulation defined on a bounded curve with two corners is designed to solve the direct scattering problem, and a Newton method with multifrequency far field data is applied to reconstruct the local perturbed surface. In contrast to [2], the integral equation formulation does not involve any infinite integral or a DtN map, and therefore it can be solved by using the Nyström method efficiently. Compared with [1], the authors in [6] make use of linear sampling method [7,8] to solve the inverse problem under the condition that the symmetric domain is with no corners. Then, [9] gives the uniqueness of the inverse scattering problem by using the generalization of Green’s theorem to nonsmooth domain. The well-posedness of direct scattering problems is established in [2,10–13]. For the rough surface problems we refer to [14–16]. For the inverse cavity problems, the reader is referred to [17,18]. For the grating problems, we refer to [19] and references therein. Motivated by [1,6] for the two dimensional case, we are currently trying to extend the technique to a simple but fully three dimensional model for the electromagnetic exploration under the condition that the local perturbation is above an infinite perfectly conducting plate. We reformulate the scattering problem as an exterior Maxwell problem for a symmetric domain, and then use linear sampling method to reconstruct the shape of local perturbation. For the application of linear sampling method in inverse electromagnetic scattering problems, we refer to [20]. This paper is organized as follows. In Section 2 we describe the model problem in detail, and transform the model problem for unbounded obstacle into a new scattering problem for bounded obstacle by symmetric continuation. Then, we prove an equivalence between the direct scattering problem and its reformulation. Section 3 formulates the inverse scattering problem, and describes the implementation details of linear sampling method. Finally, we show some numerical examples in this section. 2. Presentation of the forward and inverse scattering problem In this section, a direct scattering problem is considered for Maxwell’s equations within the upper half-space R3+ = {x ∈ R : x3 > 0}. We decompose the space R3 = R3+ ∪ R30 ∪ R3− in a hyperplane R30 = {x ∈ R3 : x3 = 0} corresponding to the surface where the local perturbation above on, and R3− = {x ∈ R3 : x3 < 0} corresponding to the lower half-space. Let D+ ⊂ R3+,0 := {x ∈ R3 : x3 ≥ 0} be a bounded domain such that ∂ D+ ⊂ R3+ is a smooth surface. Suppose ΣD,0 = D+ ∩ R30 is also bounded such that the boundary ∂ ΣD,0 is a smooth closed curve. Denote by (e1 , e2 , e3 ) the usual Cartesian basis of R3 , by x = (x1 , x2 , x3 )⊤ a point in R3 , and by ν+ the unit outward normal to the boundary. Throughout let x · y and x × y be the scalar product and the vector product of x, y ∈ R3 , respectively, and let |x| denote the Euclidean norm of x. For x ∈ R3 we denote by 3
T1 x =
1 0 0
0 1 0
0 0 x = (x1 , x2 , −x3 )⊤ −1
its reflection about the plane R30 , and define T2 x = −T1 x = (−x1 , −x2 , x3 )⊤ . For convenience and also for later use in this section, we give that (T1 x) · y = x · (T1 y ), (T1 x) × (T1 y ) = T2 (x × y ) and (T2 x) × (T1 y ) = T1 (x × y ). Given as incident fields i
∇ × ∇ × pexp(ikx · d ) = ik(d × p) × dexp(ikx · d ), k H i (x) = ∇ × pexp(ikx · d ) = ikd × pexp(ikx · d ), E i (x) =
where the wave number k is positive, d is a unit vector describing the direction of propagation, and p is the polarization vector. The solution of scattering problem for the half-space, i.e. absence of D+ , is given by the reflected wave E r (x) = T2 E i (T1 x), H r (x) = T1 H i (T1 x). For convenience, we denote E i,r = E i + E r , H i,r = H i + H r and f+ := −ν+ × E i,r = 0, for x ∈ R30 . Then the direct scattering problem for the scattering of a time harmonic electromagnetic plane wave by a locally s s perturbed perfectly conducting surface ∂ D+ is to find an electric field E+ = E i,r + E+ and a magnetic field H+ = H i,r + H+ such that
∇ × E+ − ikH+ = 0,
in R3+ \ D+ ,
(1)
∇ × H+ + ikE+ = 0,
in R+ \ D+ ,
(2)
ν+ × E+ = 0,
\ ΣD,0 ) ∪ ∂ D+ .
(3)
on (
R30
3
s s The unknown scattered fields E+ and H+ are required to satisfy the Silver–Müller radiation condition s lim r (H+ × xˆ − E+s ) = 0,
r →∞
uniformly in xˆ = x/|x| ∈ ∂ B := {x ∈ R3 : |x| = 1}, r = |x|. The geometry of this model problem is shown in Fig. 2.1.
(4)
H. Dong et al. / Wave Motion 60 (2016) 95–107
97
Fig. 2.1. Geometry of the model problem.
The corresponding inverse scattering problem can be stated as: for given wave number k > 0 and incident field (E i , H i ), ∞ reconstruct the unknown boundary ∂ D+ , by use of the electric far field pattern E+ (ˆx, d , p), ∀ˆx ∈ ∂ B+ , ∀d ∈ ∂ B− , ∀p ∈ R3 . In order to reformulate this direct scattering problem in half-space as an exterior Maxwell problem, we denote by ∂ D− = {x ∈ R3 : T1 x ∈ ∂ D+ } the reflection of ∂ D+ about R30 . Then ∂ D− ∪ ∂ D+ is the boundary ∂ D of a bounded domain D. Let E s , H s be the radiating solution to the Maxwell equations satisfying
∇ × E s − ikH s = 0,
¯, in R3 \ D
(5)
∇ × H + ikE = 0, in R3 \ D¯ , ν × E s = f , on ∂ D, lim r (H s × xˆ − E s ) = 0,
(6)
s
s
(7) (8)
r →∞
where f+ (x),
x3 ≥ 0,
T1 f+ (T1 x),
x3 < 0,
f (x) =
and
ν(x) =
ν+ (x), T1 ν+ (T1 x),
x3 ≥ 0 , x3 < 0.
Then we can prove that Theorem 1. Problem (1)–(4) is equivalent to problem (5)–(8). s Proof. Assume that (E+ , H+s ) is a radiating solution to the problem (1)–(4), and extend (E+s , H+s ) by s E+ (x),
x3 ≥ 0,
s T2 E+ (T1 x),
x3 < 0,
E (x) = s
s H+ (x),
x3 ≥ 0 ,
s T1 H + (T1 x),
x3 < 0.
H (x) = s
and
Straightforward calculations show that (E , H ) is the solution of the problem (5)–(8). On the other hand, if E s (x), H s (x) solves the problem (5)–(8), we then obtain that T2 E s (T1 x), T1 H s (T1 x) also solves the problem (5)–(8). From the uniqueness for solutions to exterior Maxwell boundary value problems [21] we conclude that s E s (x) = T2 E s (T1 x), which implies that the first two components of vector E s vanish at x ∈ R30 . Consequently, ν+ × E+ = 0, 3 s s for x ∈ R0 . Thus it can be verified that (E+ , H+ ) is the solution of the problem (1)–(4). s
s
Theorem 2. The electric far field pattern E∞ (ˆx, d, p) for the scattering problem (5)–(8) satisfies the following relations (1) E∞ (T1 xˆ , d, p) = T2 E∞ (ˆx, d, p); (2) E∞ (ˆx, T1 d, p) = −E∞ (ˆx, d, T1 p); (3) E∞ (T1 xˆ , T1 d, p) = T1 E∞ (ˆx, d, T1 p). Proof. (1) Since by the Huygens’ principle ([21, Th 6.24]), the far field pattern is given by E (ˆx) = ∞
ik 4π
xˆ × ∂D
ν(y) × H (y) × xˆ exp(−ikxˆ · y)ds(y).
Define R(ˆx) =
∂D
Q (y, xˆ )exp(−ikxˆ · y)ds(y),
where Q (y, xˆ ) = P (y) × xˆ , P (y) = ν(y) × H (y). Now the far field pattern is written as E ∞ (ˆx) =
ik 4π
xˆ × R(ˆx).
98
H. Dong et al. / Wave Motion 60 (2016) 95–107
Noting that H i,r (T1 x) = T1 H i,r (x) and H s (T1 x) = T1 H s (x), we have H (T1 x) = T1 H (x). Using the symmetry of scatterer D, we obtain that ν(T1 x) = T1 ν(x). Hence, E ∞ (T1 xˆ ) =
= = = = = =
ik 4π ik 4π ik 4π ik 4π ik 4π ik 4π ik 4π
(T1 xˆ ) × (T1 xˆ ) ×
ν(y) × H (y) × (T1 xˆ )exp −ik(T1 xˆ ) · y ds(y) ∂D
T1 ν(T1 y) × T1 H (T1 y) × (T1 xˆ )exp −ikxˆ · (T1 y) ds(y)
∂D
× (T1 xˆ )exp −ikxˆ · (T1 y) ds(y) T2 ν(T1 y) × H (T1 y) (T1 xˆ ) × ∂D
(T1 xˆ ) ×
(T1 xˆ ) ×
∂D
∂D
(T1 xˆ ) × T1 R(ˆx)
T1 Q (T1 y, xˆ )exp −ikxˆ · (T1 y) ds(y)
T2 xˆ × R(ˆx)
T2 P (T1 y) × (T1 xˆ )exp −ikxˆ · (T1 y) ds(y)
= T2 E ∞ (ˆx). (2) In terms of the incident field E i and the reflected wave E r , we have E i (x, T1 d , T1 p) = T1 E i (T1 x, d , p),
E r (x, T1 d , T1 p) = −E i (x, d , p),
and T1 E i,r (T1 x, d , p) = E i,r (x, T1 d , T1 p). Since E i,r (T1 x, d , p) = T2 E i,r (x, d , p), we get E i,r (x, T1 d , T1 p) = −E i,r (x, d , p). According to the boundary condition (7), one can derive that
ν × E s (x, T1 d , T1 p) = −ν × E i,r (x, T1 d , T1 p) = −ν × −E i,r (x, d , p) = ν × E i,r (x, d , p). Thus,
ν × −E s (x, T1 d , T1 p) = −ν × E i,r (x, d , p) = ν × E s (x, d , p). From the uniqueness for solutions to exterior Maxwell boundary value problems [21] we get
−E s (x, T1 d , T1 p) = E s (x, d , p). Consequently, the corresponding far field pattern satisfies E ∞ (ˆx, d , p) = −E ∞ (ˆx, T1 d , T1 p). (3) Combining the previous results (1) and (2), we finally obtain the following results
T1 E ∞ (T1 xˆ , T1 d , T1 p) = T1 −E ∞ (T1 xˆ , d , p)
= T2 E ∞ (T1 xˆ , d , p) = E ∞ (ˆx, d , p). The proof of the theorem is now completed.
∞ This theorem tells us that for given the far field pattern E+ (ˆx, d , p), for xˆ ∈ ∂ B+ , d ∈ ∂ B− , the far field pattern
E ∞ (ˆx, d , p),
∀ˆx ∈ ∂ B, ∀d ∈ ∂ B,
corresponding the scattered field E s of problem (5)–(8) can be obtained. Hence, we are able to deal with this problem by using linear sampling method.
H. Dong et al. / Wave Motion 60 (2016) 95–107
99
3. Linear sampling method and numerical implementation To analyze the inverse scattering problem we consider the following Hilbert spaces: H (curl, D) := {a ∈ (L2 (D))3 : curla ∈ (L2 (D))3 }, − 21
Ht
1 3 (∂ D) := a ∈ H − 2 (∂ D) : ν · a = 0 ,
L2t (∂ B) := {a ∈ (L2 (∂ B))3 : a · xˆ = 0, xˆ ∈ ∂ B}, H
− 12
− 21
(Div, ∂ D) := a ∈ Ht
(∂ D) : Diva ∈ H
− 12
(∂ D)
with norm
∥a∥2H −1/2 (Div,∂ D) := ∥a∥2(H −1/2 (∂ D))3 + ∥Diva∥2H −1/2 (∂ D) . ¯ ) denotes the space of all fields a : R3 \ D¯ → C3 such that a ∈ H (curl, (R3 \ D¯ ) ∪ BR ) for Furthermore, Hloc (curl, R3 \ D all open balls BR containing the closure of D. 3.1. Far-field data We are now in a position to describe the numerical algorithm. To this end, we introduce the far field operator F : L2t (∂ B) → L2t (∂ B) such that
(F g )(ˆx) :=
∂B
E ∞ (ˆx, d , g (d ))ds(d ),
xˆ ∈ ∂ B.
(9)
The linear sampling method is based on seeking a function g (·, z , p) ∈ L2t (∂ B), such that for fixed p and z,
F g (·, z , p) (ˆx) = Ee∞ (ˆx, z , p)
=
ik 4π
(ˆx × p) × xˆ exp(−ikxˆ · z ) −
ik 4π
(ˆx × T1 p) × xˆ exp(−ikxˆ · T1 z ),
(10)
where Ee∞ is the electric far field pattern of the electric dipole Ee . In contrast to [22], the electric dipole Ee has to be modified to satisfy ν × Ee = 0 on R30 . In this case, we have Ee (x, z , p) =
i k
i ∇x × ∇x × pΦ (x, z ) − ∇x × ∇x × T1 pΦ (x, T1 z ) , k
where
Φ (x, z ) =
1 exp(ik|x − z |) 4π
|x − z |
,
x, z ∈ R3 , x ̸= z . 1
Moreover, we introduce a bounded linear operator A : H − 2 (Div, ∂ D) → L2t (∂ B) which maps the boundary data f ∈H
− 21
(Div, ∂ D) to the far field pattern E ∞ of the radiating solution E s ∈ Hloc (curl, R3 \ D¯ ) of (5)–(8). Then we have that
F g = −A(ν × Eg ),
(11)
where Eg (x) = ik
∂B
eikx·d g (d )ds(d ) + ik
∂B
eik(T1 x)·d T2 g (d )ds(d ).
We can now prove the main result of this paper. Theorem 3. Assume that k is not a Maxwell eigenvalue for D and let F be the far field operator (9) for the scattering problem (5)–(8). Then the following hold: (1) For z ∈ D and a given ε > 0 there exists a gεz ∈ L2t (∂ B) such that
∥F gεz − E∞ e (·, z, p)∥L2 (∂ B) < ε t
(12)
and the corresponding Herglotz function Egεz converges to the solution of interior Maxwell problem
∇ × E − ikH = 0, ν × (E + Ee ) = 0, in H (curl, D) as ε → 0.
∇ × H + ikE = 0 in D, on ∂ D,
(13) (14)
100
H. Dong et al. / Wave Motion 60 (2016) 95–107
¯ and a given ε > 0, every gεz ∈ L2t (∂ B) that satisfies (12) is such that (2) For z ∈ R3 \ D lim ∥gεz ∥L2 (∂ B) = ∞.
ε→0
t
Proof. For convenience, we denote
∇x × ∇x × pΦ (x, z ) , k Eg1 (x) = ik eikx·d g (d )ds(d ), i
Ee1 (x, z , p) =
i Ee2 (x, z , p) = − ∇x × ∇x × k Eg2 (x) = ik
∂B
∂B
T1 pΦ (x, T1 z ) ,
eik(T1 x)·d T2 g (d )ds(d ).
Then Ee = Ee1 + Ee2 , Eg = Eg1 + Eg2 . Under the assumption on k we have well-posedness of the interior Maxwell problem
in H (curl, D). Given ε > 0, from the denseness of Herglotz operator(Corollary 7.19 in [21]) we can choose gzε ∈ L2t (∂ B) such that
ν(x) × Eg1zε (x, z , p) + Ee1 (x, z , p)
H −1/2 (Div,∂ D)
<
ε . 2∥A∥
Define W (x, z , p) = ν(x) × Eg1ε (x, z , p) + Ee1 (x, z , p) , the above equation implies
z
∥W (·, z , p)∥H −1/2 (Div,∂ D) <
ε 2∥A∥
.
(15)
Straightforward calculations show that
ν(x) × Eg2zε (x, z , p) + Ee2 (x, z , p) = T1 ν(T1 x) × T2 Eg1zε (T1 x, z , p) + T2 Ee1 (T1 x, z , p) = T2 ν(T1 x) × Eg1zε (T1 x, z , p) + Ee1 (T1 x, z , p) = T2 W (T1 x, z , p). With the aid of (15) and the symmetry of scatterer D, it can be verified that
∥T2 W (T1 ·, z , p)∥H −1/2 (Div,∂ D) = ∥W (·, z , p)∥H −1/2 (Div,∂ D) <
ε 2∥A∥
.
Thus, we get
ν(x) × Egzε (x, z , p) + Ee (x, z , p)
H −1/2 (Div,∂ D)
≤ ∥W (·, z , p)∥H −1/2 (Div,∂ D) + ∥T2 W (T1 ·, z , p)∥H −1/2 (Div,∂ D) <
ε . ∥A∥
Then (12) follows by observing Eq. (11). Now if z ∈ D, from the well-posedness of the interior Maxwell problem,
lim ν(x) × Egzε (x, z , p) + Ee (x, z , p)
ε→0
H −1/2 (Div,∂ D)
=0
implies lim Egzε − E H (curl,D) = 0,
ε→0
where E is the solution to (13)–(14). By arguments similar to those in Theorem 7.21 of [21], the second statements can be established by noting (15).
In the following, we briefly describe the discretization of the far field equation and the details may be found in [20,22,23]. Let (ˆx, e1 (ˆx), e2 (ˆx)) denote an orthonormal basis of R3 with xˆ ∈ ∂ B. Note that E ∞ (ˆx, d , p) · xˆ = 0, the far field equation (10) is then equivalent to
∞ ∞ e1 (ˆx) · E (ˆx, d , g (d , z , p))ds(d ) = e1 (ˆx) · Ee (ˆx, z , p), ∂ B e2 (ˆx) · E ∞ (ˆx, d , g (d , z , p))ds(d ) = e2 (ˆx) · Ee∞ (ˆx, z , p). ∂B
By using the reciprocity relation p · E ∞ (ˆx, d , q) = q · E ∞ (−d , −ˆx, p),
∀ˆx, d ∈ ∂ B, ∀p, q ∈ C3 ,
(16)
H. Dong et al. / Wave Motion 60 (2016) 95–107
101
and the numerical quadrature formula
∂B
f (d )ds(d ) ≈
N
ωj f (dj ),
j =1
where ωj are the weights associated with the mesh triangular meshing of the unit sphere ∂ B, dj are the vertices of mesh triangles. The system (16) can be transformed into the following discrete version
N ωj E ∞ (−dj , −ˆx, e1 (ˆx)) · g (dj , z , p) = e1 (ˆx) · Ee∞ (ˆx, z , p), j =1
N ωj E ∞ (−dj , −ˆx, e2 (ˆx)) · g (dj , z , p) = e2 (ˆx) · Ee∞ (ˆx, z , p). j =1
In view of the fact that the unknown g (dj , z , p) is a tangential vector field, we can set g (dj , z , p) = g1 (dj , z , p)e1 (dj ) + g2 (dj , z , p)e2 (dj ) with g1 (dj , z , p) and g2 (dj , z , p) ∈ C. Thus, we obtain the fully discrete linear equations
N 2 ωj gl (dj , z , p)E ∞ (−dj , −ˆx, e1 (ˆx)) · el (dj ) = e1 (ˆx) · Ee∞ (ˆx, z , p), j =1 l =1
(17)
N 2 ωj gl (dj , z , p)E ∞ (−dj , −ˆx, e2 (ˆx)) · el (dj ) = e2 (ˆx) · Ee∞ (ˆx, z , p). j =1 l =1
In our numerical experiments, the far field pattern E ∞ (−dj , −ˆx, el (ˆx)), l = 1, 2 for N measurement points {ˆxj }Nj=1 is computed by boundary integral equation method [24]. These measurement points are roughly uniformly distributed on the unit sphere ∂ B, and they will also serve as incident directions {dj }Nj=1 . Thus, the quadrature weights ωj can be seen as a
constant. Let qˆ be some unit vector in R3 satisfying qˆ × xˆ j ̸= 0 for any j. In our examples, we take qˆ = (1, 1, 1)⊤ . Then we set the two orthogonal polarizations e1 (ˆxj ) = qˆ × xˆ j /|ˆq × xˆ j |
and e2 (ˆxj ) = qˆ × (ˆxj × qˆ )/|ˆq × (ˆxj × qˆ )|.
Now, the system (17) can be rewritten as A∞ g⃗ (z , p) = ⃗ f (z , p)
(18)
where A∞ is a 2N × 2N complex matrix corresponding to the given far field data, vector g⃗ = (g1 , g2 ) unknowns and ⃗ f is a vector given by the right hand side of (17). In order to avoid inverse crimes, we corrupt the matrix A∞ with random noise
⊤
represents the
A∞ := A∞ (I + ϵ R1 + iϵ R2 ) where R1 and R2 are two matrices of real random numbers in [−1, 1] and ϵ is a parameter. In the following examples, we have chosen ϵ = 0.01 and ϵ = 0.1. Since the ill-posedness of this problem, instead of Eq. (18), we solve the regularized equation
γ g⃗ + A∗∞ A∞ g⃗ = A∗∞ ⃗f
(19)
is the adjoint of A∞ and γ is the regularization parameter. Let p1 = (1, 0, 0) , p2 = (0, 1, 0) , p3 = (0, 0, 1)⊤ . For each z in the sampling grid, we compute the indicator function
where A∗∞
G(z ) =
⊤
1
∥⃗g (z , p1 )∥−2 + ∥⃗g (z , p2 )∥−2 + ∥⃗g (z , p3 )∥−2
21
⊤
.
3 The reconstruction corresponds to the isosurfaces of G(z ) close to zero. Consequently, we will choose a constant C and plot surfaces where G(z ) = C . In this section, we present and discuss three numerical examples. In each example, we use 50 and 100 incoming waves respectively to reconstruct the shape of local perturbation, and these incoming waves are generated by Gmsh which is a 3D finite element grid generator. Because of the symmetry of this problem, we only need to put sampling points in upper space. Example 1. In the first example, we consider a specific case in which the shape of local perturbation ∂ D+ is a half sphere with radius R = 1. The wave number is k = 3. We take (41)3 sampling points within a cube prescribed as Ω = [−2, 2] × [−2, 2] × [0.01, 1.5]. As can be seen from Fig. 3.1, the local perturbation can be reconstructed when illuminations and measurements N = 50 and N = 100 respectively, provided a good isosurface of G(z ) and regularization parameter are chosen.
102
H. Dong et al. / Wave Motion 60 (2016) 95–107
Fig. 3.1. Reconstruction of the half ball with radius R = 1 and the exact surface.
Example 2. For the second example, we wish to reconstruct a half cushion-shaped surface with the parametrization z (θ , ϕ) =
0.27 + 0.065(cos 2ϕ − 1)(cos 4θ − 1)r (θ , ϕ),
where spherical coordinates r (θ , ϕ) = (sin θ cos ϕ, sin θ sin ϕ, cos θ ),
θ ∈ [0, π /2], ϕ ∈ [0, 2π ].
We choose wave number k = 4 and (41) sampling points within Ω = [−1, 1] × [−1, 1] × [0.01, 0.8]. Fig. 3.2 shows the results. 3
Example 3. Our third example is intended to reconstruct a half bean-like surface with the parametrization
0.5 1 − 0.1 cos(π cos θ ) sin θ cos ϕ , z (θ , ϕ) = 0.5 sin θ sin ϕ + 0.3 cos(π cos θ ) 0.7 cos θ
θ ∈ [0, π /2], ϕ ∈ [0, 2π ].
We assume wave number k = 7 and the sampling domain is Ω = [−1, 1] × [−1, 1] × [0.01, 1] with (41)3 sampling points. Results are shown in Fig. 3.3. Moreover, one can observe how the increase of N helps a better resolution of the reconstruction. The above numerical results illustrate that the linear sampling method gives a feasible reconstruction of the locally perturbed perfectly conducting plate even in the presence of 10% noise in measurements.
H. Dong et al. / Wave Motion 60 (2016) 95–107
103
Fig. 3.2. Reconstruction of the half cushion and the exact surface.
3.2. Limited aperture data In this section, we restrict the set of observation points to only half of the upper half sphere Γ ⊂ ∂ B. This means we know the half of measurement data. In this case the data are E ∞ (ˆxi , dj , p) for some (ˆxi , dj , p) ∈ Γ × (−Γ ) × R3 . Figs. 3.4 and 3.5 display reconstructions that correspond to the three examples. We observe that the quality of the reconstruction deteriorates as the aperture decreases. However, one notices in all cases that we have a rather good approximation of the projection of the scatterer on the plane located in the illuminated region. 3.3. Near-field data In this section we consider the linear sampling method for the near field measurements. Let G+ ⊂ R3+,0 be a bounded
domain enclosing the local perturbation such that ∂ G+ ⊂ R3+ is a smooth surface. The incident field is considered to be the electric field of an electric dipole located at y ∈ ∂ G+ with polarization p ∈ R3 , and is given by Eei (x, y , p) =
i k
∇x × ∇x × pΦ (x, y ) = Ee1 (x, y , p).
The inverse scattering problem we are interested in is to determine the unknown boundary ∂ D+ from knowledge of the s tangential components ν+ ×E+ of the scattered electric field measured on ∂ G+ for all incident point sources Eei (x, y , p), ∀y ∈ 3 ∂ G+ , ∀p ∈ R .
104
H. Dong et al. / Wave Motion 60 (2016) 95–107
Fig. 3.3. Reconstruction of the half bean and the exact surface.
Note that the reflected wave is Eer (x, y , p) = T2 Eei (T1 x, y , p) = Ee2 (x, y , p). Analogously to the case of far-field data, we ¯ with E i,r replaced can reformulate the scattering problem as an exterior Maxwell problem (5)–(8) for a symmetric domain D by Ee . Then, we have the following result. Theorem 4. For x, y ∈ ∂ G, the tangential component ν × Es of the scattered electric field satisfies the following relations (1) (ν × Es )(T1 x, y, p) = T2 (ν × Es )(x, y, p); (2) (ν × Es )(x, T1 y, p) = −(ν × Es )(x, y, T1 p); (3) (ν × Es )(T1 x, T1 y, p) = T1 (ν × Es )(x, y, T1 p). Now we reformulate the linear sampling method in terms of the near field data. To this end we consider the near field ¯ is contained in G, and the near field operator F : L2t (∂ G) → L2t (∂ G) where G = G+ ∪ G− is a bounded domain such that D equation
(F ϕz )(x) :=
∂G
ν(x) × E s (x, y , ϕz ,p (y ))ds(y ) = ν(x) × Ee (x, z , p),
x ∈ ∂ G.
(20)
For convenience, ∂ G is chosen to be a sphere with radius r = 2 in Example 1 and r = 1 in Examples 2 and 3. By using the symmetry relation (Theorem 6.32 in [21]) q · E s (x, y , p) = p · E s (y , x, q),
¯, x, y ∈ R3 \ D
H. Dong et al. / Wave Motion 60 (2016) 95–107
Fig. 3.4. Reconstruction for limited aperture case.
Fig. 3.5. Reconstruction for limited aperture case.
105
106
H. Dong et al. / Wave Motion 60 (2016) 95–107
Fig. 3.6. Reconstruction for near-field data.
the near field equation (20) can be transformed into the fully discrete version
N 2 ωj ϕzl ,p (yj )E s (yj , xi , e2 (ˆxi )) · el (ˆyj ) = e2 (ˆxi ) · Ee (xi , z , p), j =1 l =1
N 2 ωj ϕzl ,p (yj )E s (yj , xi , e1 (ˆxi )) · el (ˆyj ) = e1 (ˆxi ) · Ee (xi , z , p), j =1 l =1
where we have set ϕz ,p (yj ) = ϕz1,p (yj )e1 (ˆyj ) + ϕz2,p (yj )e2 (ˆyj ) and the measurement points {yj }Nj=1 also serve as the locations {xi }Ni=1 of the point sources. Similarly as far-field data case, we plot the isosurfaces of G(z ) for ϕz ,p . In the numerical experiments, we keep the settings of wave number, sampling domain and the number of sampling points in each example. Fig. 3.6 shows the result of reconstruction.
4. Conclusions In this paper we have presented a numerical method for solving a three-dimensional inverse problem to determine the shape of a local perturbation of a perfectly conducting plate. We have applied the linear sampling method to the far-field data case, limited aperture case and near-field data case. The numerical experiments show that satisfactory reconstructions can be obtained in the three cases. Acknowledgments This work was supported by the National Science Foundation of China (Grant Nos. 91123031-3, 11371172, 41474102, 91230119) and the Scientific Research Foundation of Jilin Province (Grant No. 20140520055JH).
References [1] R. Kress, T. Tran, Inverse scattering for a locally perturbed half-plane, Inverse Problems 16 (2000) 1541–1559. [2] G. Bao, J. Lin, Imaging of local surface displacement on an infinite ground plane: the multiple frequency case, SIAM J. Appl. Math. 71 (2011) 1733–1752.
H. Dong et al. / Wave Motion 60 (2016) 95–107
107
[3] G. Bao, J. Lin, Imaging of reflective surfaces by near-field optics, Opt. Lett. 37 (2012) 5027–5029. [4] G. Bao, J. Lin, Near-field imaging of the surface displacement on an infinite ground plane, Inverse Probl. Imaging 7 (2013) 377–396. [5] H.W. Zhang, B. Zhang, A novel integral equation for scattering by locally rough surfaces and application to the inverse problem, SIAM J. Appl. Math. 73 (2013) 1811–1829. [6] L.X. Feng, F.M. Ma, A sampling method for solving inverse scattering problems with a locally perturbed half plane, Northeast. Math. J. 19 (2003) 1–4. [7] D. Colton, A. Kirsch, A simple method for solving inverse scattering problems in the resonance region, Inverse Problems 12 (1996) 383–393. [8] D. Colton, H. Haddar, M. Piana, The linear sampling method in inverse electromagnetic scattering theory, Inverse Problems 19 (2003) S105–S137. [9] G.Z. Yan, Uniqueness of the inverse scattering problem of a locally perturbed half-plane, Comput. Math. Appl. 48 (2004) 411–418. [10] M. Durán, I. Muga, J.C. Nédélec, The Helmholtz equation in a locally perturbed half-plane with passive boundary, IMA J. Appl. Math. 71 (2006) 853–876. [11] M. Durán, I. Muga, J.C. Nédélec, The Helmholtz equation in a locally perturbed half-space with non-absorbing boundary, Arch. Ration. Mech. Anal. 191 (2009) 143–172. [12] S.N. Chandler-Wilde, A.T. Peplow, A boundary integral equation formulation for the Helmholtz equation in a locally perturbed half-plane, ZAMM Z. Angew. Math. Mech. 85 (2005) 79–88. [13] H. Ammari, C. Latiri-Grouz, An integral equation method for the electromagnetic scattering from a scatterer on an absorbing plane, Integral Equations Operator Theory 39 (2001) 159–181. [14] Y. He, P.J. Li, J. Shen, A new spectral method for numerical solution of the unbounded rough surface scattering problem, J. Comput. Phys. 275 (2014) 608–625. [15] S.N. Chandler-Wilde, C. Lines, A time domain point source method for inverse scattering by rough surfaces, Computing 75 (2005) 157–180. [16] C. Burkard, R. Potthast, A multi-section approach for rough surface reconstruction via the Kirsch–Kress scheme, Inverse Problems 26 (2010) 045007. [17] G. Bao, J. Gao, P. Li, Analysis of direct and inverse cavity scattering problems, Numer. Math. Theory Methods Appl. 4 (2011) 419–422. [18] P.J. Li, An inverse cavity problem for Maxwell’s equations, J. Differential Equations 252 (2012) 3209–3225. [19] G. Bao, P.J. Li, Convergence analysis in near-field imaging, Inverse Problems 30 (2014) 085008. [20] F. Cakoni, D. Colton, P. Monk, The Linear Sampling Method in Inverse Electromagnetic Scattering, SIAM, Philadelphia, 2011. [21] D. Colton, R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, third ed., Springer, New York, 2013. [22] P. Monk, Finite Element Methods for Maxwell’s Equations, Oxford, New York, 2003. [23] D. Colton, H. Haddar, P. Monk, The linear sampling method for solving the electromagnetic inverse scattering problem, SIAM J. Sci. Comput. 24 (2002) 719–731. [24] M. Ganesh, S.C. Hawkins, A high-order tangential basis algorithm for electromagnetic scattering by curved surfaces, J. Comput. Phys. 227 (2008) 4543–4562.