Orthogonality Sampling Method for the Electromagnetic Inverse Scattering Problem

TL;DR

This paper develops a theoretical foundation and resolution analysis for the orthogonality sampling method in electromagnetic inverse scattering, demonstrating its effectiveness for locating and shaping anisotropic scatterers from far field data.

Contribution

It provides a new theoretical basis and resolution analysis for the orthogonality sampling method, linking it to the direct sampling and factorization methods.

Findings

The method is fast, robust, and effective for 3D anisotropic scatterers.

Theoretical resolution limits are established.

Numerical examples validate the method's performance.

Abstract

This paper is concerned with the electromagnetic inverse scattering problem that aims to determine the location and shape of anisotropic scatterers from far field data (at a fixed frequency). We study the orthogonality sampling method which is a simple, fast and robust imaging method for solving the electromagnetic inverse shape problem. We first provide a theoretical foundation for the sampling method and a resolution analysis of its imaging functional. We then establish an equivalent relation between the orthogonality sampling method and direct sampling method as well as resolution analysis for the latter. The analysis used to justify the Factorization method for the far field operator plays an important role in the justifications. Finally, we present some numerical examples to validate the performance of the sampling methods for anisotropic scatterers in three dimensions.

Equations

$$H(\mathrm{curl},\mathcal{O})$$

$$H_{\mathrm{loc}}(\mathrm{curl},{\mathbb{R}}^3)$$

$$(\cdot ,\cdot {)}_{H(\mathrm{curl},\mathcal{O})}=(\mathrm{curl}\cdot ,\mathrm{curl}\cdot )+(\cdot ,\cdot ).$$

$$\mathrm{curl}\mathbf{E}-i\omega {\mu }_0\mathbf{H}=0\text{ and }\mathrm{curl}\mathbf{H}+i\omega \epsilon \mathbf{E}=0\text{in }{\mathbb{R}}^3.$$

$${\epsilon }_{\mathrm{r}}=\epsilon /{\epsilon }_0,k=\omega \sqrt{{\epsilon }_0{\mu }_0}.$$

$$\mathrm{curl}\mathrm{curl}\mathbf{E}-k^2{\epsilon }_{\mathrm{r}}\mathbf{E}=0,\text{in }{\mathbb{R}}^3.$$

$$\nu \times {\mathbf{E}}_{+}=\nu \times {\mathbf{E}}_{-}\text{ and}\nu \times \mathrm{curl}{\mathbf{E}}_{+}=\nu \times \mathrm{curl}{\mathbf{E}}_{-},\text{on }\partial \Omega .$$

$$\mathrm{curl}{\mathbf{H}}_{\mathrm{in}}+i\omega {\epsilon }_0{\mathbf{E}}_{\mathrm{in}}=0\text{ and}\mathrm{curl}{\mathbf{E}}_{\mathrm{in}}-i\omega {\mu }_0{\mathbf{H}}_{\mathrm{in}}=0,\text{in }{\mathbb{R}}^3.$$

$$\mathrm{curl}\mathrm{curl}{\mathbf{E}}_{\mathrm{in}}-k^2{\mathbf{E}}_{\mathrm{in}}=0,\text{in }{\mathbb{R}}^3$$

$$\mathrm{curl}\mathrm{curl}\mathbf{u}-k^2{\epsilon }_{\mathrm{r}}\mathbf{u}=k^{2}P{\mathbf{E}}_{\mathrm{in}}\text{in }{\mathbb{R}}^3,$$

$$P:={\epsilon }_{\mathrm{r}}-I.$$

$$\mathrm{curl}\mathbf{u}\times \frac{\mathbf{x}}{|\mathbf{x}|}-ik\mathbf{u}=\mathcal{O}(|\mathbf{x}{|}^{-2})\text{as }|\mathbf{x}|\to \infty ,$$

$${\mathbb{S}}^2=\{\mathbf{x}\in {\mathbb{R}}^3:|\mathbf{x}|=1\}\text{ and }{L}_t^2({\mathbb{S}}^2)=\{\mathbf{v}\in L^2({\mathbb{S}}^2{)}^3:\hat{\mathbf{x}}\cdot \mathbf{v}(\hat{\mathbf{x}})=0,\hat{\mathbf{x}}\in {\mathbb{S}}^2\}.$$

$$\mathbf{u}(\mathbf{x},\mathbf{d},\mathbf{q})=\frac{e^{ik|\mathbf{x}|}}{|\mathbf{x}|}\left({\mathbf{u}}^{\infty }(\hat{\mathbf{x}},\mathbf{d},\mathbf{q})+O\left(\frac{1}{|\mathbf{x}{|}^2}\right)\right)\text{ as }|\mathbf{x}|\to \infty ,$$

$${\mathbf{u}}^{\infty }(\hat{\mathbf{x}},\mathbf{d},\mathbf{q})={\mathbf{u}}^{\infty }(\hat{\mathbf{x}},\mathbf{d})\mathbf{q}$$

$${\mathbf{u}}^{\infty }(\hat{\mathbf{x}},\mathbf{d})=[{\mathbf{u}}^{\infty }(-\mathbf{d},-\hat{\mathbf{x}}){]}^{\top }.$$

$$(\mathcal{F}\mathbf{g})(\hat{\mathbf{x}})={\int }_{{\mathbb{S}}^2}{\mathbf{u}}^{\infty }(\hat{\mathbf{x}},\mathbf{d})\mathbf{g}(\mathbf{d})\mathrm{d}s(\mathbf{d}).$$

$$H:L_t^2({\mathbb{S}}^2)\to L^2(\Omega {)}^3\text{ is given by }(H\mathbf{g})(\mathbf{x})={\int }_{{\mathbb{S}}^2}\mathbf{g}(\mathbf{d})e^{ik\mathbf{x}\cdot \mathbf{d}}\mathrm{d}s(\mathbf{d}),\mathbf{x}\in \Omega$$

$$G:L^2(\Omega {)}^3\to L_t^2({\mathbb{S}}^2)\text{ is given by }G\mathbf{f}={\mathbf{v}}^{\infty }$$

$$H^{\ast }\mathbf{f}(\mathbf{d})=\mathbf{d}\times \left({\int }_{\Omega }\mathbf{f}(\mathbf{x})e^{-ik\mathbf{x}\cdot \mathbf{d}}\mathrm{d}\mathbf{x}\right)\times \mathbf{d},\mathbf{d}\in {\mathbb{S}}^2.$$

$$T:L^2(\Omega {)}^3\to L^2(\Omega {)}^3\text{ is given by }(T\mathbf{f})(\mathbf{x})=k^{2}P(\mathbf{f}+\mathbf{v})$$

$$\mathrm{curl}\mathrm{curl}\mathbf{w}-k^2{\epsilon }_r\mathbf{w}=0\text{ and }\mathrm{curl}\mathrm{curl}\mathbf{\phi }-k^2\mathbf{\phi }=0\text{ in}\Omega$$

$$\nu \times \mathbf{w}=\nu \times \mathbf{\phi }\text{ and }\nu \times \mathrm{curl}\mathbf{w}=\nu \times \mathrm{curl}\mathbf{\phi }\text{on }\partial \Omega$$

$${\mathcal{I}}_{OSM}({\mathbf{y}}_s):={\int }_{{\mathbb{S}}^2}{\left|{\int }_{{\mathbb{S}}^2}{\mathbf{u}}^{\infty }(\hat{\mathbf{x}},\mathbf{d})(\mathbf{d}\times \mathbf{p})\times \mathbf{d}\cdot (\hat{\mathbf{x}}\times \mathbf{p})\times \hat{\mathbf{x}}{e}^{ik\hat{\mathbf{x}}\cdot {\mathbf{y}}_s}\mathrm{d}s(\hat{\mathbf{x}})\right|}^2\mathrm{d}s(\mathbf{d}).$$

$${\mathcal{I}}_{OSM}({\mathbf{y}}_s)=\Vert \mathbf{p}\cdot \mathcal{F}{\varphi }_{{\mathbf{y}}_s}{\Vert }_{L^2({\mathbb{S}}^2)}^2,$$

$${\varphi }_{{\mathbf{y}}_s}(\mathbf{d})=(\mathbf{d}\times \mathbf{p})\times \mathbf{d}{e}^{-ik\mathbf{d}\cdot {\mathbf{y}}_s}.$$

$${\mathcal{I}}_{OSM}({\mathbf{y}}_s)$$

$$={\int }_{{\mathbb{S}}^2}{\left|{\int }_{{\mathbb{S}}^2}[{\mathbf{u}}^{\infty }(-\mathbf{d},-\hat{\mathbf{x}}){]}^{\top }(\mathbf{d}\times \mathbf{p})\times \mathbf{d}\cdot (\hat{\mathbf{x}}\times \mathbf{p})\times \hat{\mathbf{x}}{e}^{ik\hat{\mathbf{x}}\cdot {\mathbf{y}}_s}\mathrm{d}s(\hat{\mathbf{x}})\right|}^2\mathrm{d}s(\mathbf{d})$$

$$={\int }_{{\mathbb{S}}^2}{\left|{\int }_{{\mathbb{S}}^2}[{\mathbf{u}}^{\infty }(\mathbf{d},\hat{\mathbf{x}}){]}^{\top }(\mathbf{d}\times \mathbf{p})\times \mathbf{d}\cdot (\hat{\mathbf{x}}\times \mathbf{p})\times \hat{\mathbf{x}}{e}^{-ik\hat{\mathbf{x}}\cdot {\mathbf{y}}_s}\mathrm{d}s(\hat{\mathbf{x}})\right|}^2\mathrm{d}s(\mathbf{d})$$

$$={\int }_{{\mathbb{S}}^2}{\left|{\int }_{{\mathbb{S}}^2}{\mathbf{u}}^{\infty }(\mathbf{d},\hat{\mathbf{x}})((\hat{\mathbf{x}}\times \mathbf{p})\times \hat{\mathbf{x}}{e}^{-ik\hat{\mathbf{x}}\cdot {\mathbf{y}}_s})\cdot (\mathbf{d}\times \mathbf{p})\times \mathbf{d}\mathrm{d}s(\hat{\mathbf{x}})\right|}^2\mathrm{d}s(\mathbf{d})$$

$$={\int }_{{\mathbb{S}}^2}{\left|(\hat{\mathbf{x}}\times \mathbf{p})\times \hat{\mathbf{x}}\cdot {\int }_{{\mathbb{S}}^2}{\mathbf{u}}^{\infty }(\hat{\mathbf{x}},\mathbf{d})(\mathbf{d}\times \mathbf{p})\times \mathbf{d}{e}^{-ik\mathbf{d}\cdot {\mathbf{y}}_s}\mathrm{d}s(\mathbf{d})\right|}^2\mathrm{d}s(\hat{\mathbf{x}})$$

$$={\Vert (\hat{\mathbf{x}}\times \mathbf{p})\times \hat{\mathbf{x}}\cdot \mathcal{F}{\varphi }_{{\mathbf{y}}_s}\Vert }_{L^2({\mathbb{S}}^2)}^2.$$

$$(\hat{\mathbf{x}}\times \mathbf{p})\times \hat{\mathbf{x}}=\mathbf{p}(\hat{\mathbf{x}}\cdot \hat{\mathbf{x}})-\hat{\mathbf{x}}(\hat{\mathbf{x}}\cdot \mathbf{p})=\mathbf{p}-\hat{\mathbf{x}}(\hat{\mathbf{x}}\cdot \mathbf{p}),$$

$${\mathcal{I}}_{OSM}({\mathbf{y}}_s)=O\left(\frac{1}{\mathrm{dist}({\mathbf{y}}_s,\Omega {)}^2}\right),\text{as }\mathrm{dist}({\mathbf{y}}_s,\Omega )\to \infty .$$

Full text

Keywords. orthogonality sampling method, inverse electromagnetic scattering, direct sampling method, Maxwell’s equations, anisotropic media

**AMS subject classification. ** 35R30, 35R09, 65R20

1 Introduction

In this paper, we consider the inverse shape problem that is derived from the time-harmonic electromagnetic scattering of an inhomogeneous anisotropic medium. In many physical applications such as non-destructive testing and medical imaging one wishes to infer the shape and/or material properties of the scatterer from measured electromagnetic data. We assume that the far field measurements are known and we wish to analyze two sampling methods for recovering the scatterer. Sampling methods generally fall under the category of qualitative (otherwise known as non-iterative or direct) reconstruction techniques. These methods are advantageous to use since they require little a-prior information to implement and are computationally simple. Qualitative methods have been used to solve multiple inverse shape problems in electromagnetic scattering (see for e.g. [3, 15, 10] and the references therein). These methods have also been extended to inverse problems in the time domain. In [6, 4] the linear sampling and factorization methods are extended to inverse scattering problem problems in the time domain and in [5] the MUSIC algorithm is studied for recovering small volume scatterers for the time-dependent acoustic scattering problem. Here we rigorously analyze both the orthogonality sampling method (OSM) and direct sampling method (DSM) for recovering a penetrable inhomogeneous anisotropic medium from electromagnetic far field data.

In general, sampling methods allow one to recover the scatterer by connecting the scatterer to the solution of a linear ill-posed equation involving the far field operator. Roughly speaking, the linear sampling method gives that the so-called far field equation $\mathcal{F}g=\phi_{\mathbf{y}_{s}}$ is only solvable (via a regularization strategy) provided the sampling point is contained in the scatterer. Here the righthand side $\phi_{\mathbf{y}_{s}}$ is known and depends on the sampling point $\mathbf{y}_{s}\in\mathbb{R}^{3}$. This allows one to define an imaging functional that is the reciprocal of the norm of the solution to the far field equation which should only be non-zero as the regularization tends to zero for sampling points in the scatterer. See [3] for the analysis of the linear sampling method for the electromagnetic scattering problem. The factorization method gives that $\phi_{\mathbf{y}_{s}}$ is in the range of a positive self-adjoint compact operator defined by the far field operator if and only if the sampling point is in the scatterer. By appealing to Picard’s criteria one can derive an imaging functional using the spectral decomposition of the far field operator see [15].

The OSM was first introduced in [20] for the inverse acoustic scattering from sound soft scatterers. Comparing with classical sampling methods the OSM is simpler to implement, can image (small) scatterers with only one incident field, and its stability can be easily justified. However, its mathematical foundation was only partly known. For the Helmholtz equation case, the method was rigorously justified in [20] for small scatterers and in [17] for scatterers with arbitrary shape using multi-static data. It was also first proved in [17] that the OSM is equivalent to the DSM studied in this cited paper via a remarkable connection to the analysis of the factorization method. Recently, in [11, 16] this DSM was studied in connection to the spectral decomposition of the far field operator. We also refer to [9] for the analysis of a multifrequency OSM. The OSM for Maxwell’s equations and its analysis for the case of small scatterers have been recently established in [19]. Motivated by these recent works we study in this paper the OSM for anisotropic Maxwell’s equations and provide a theoretical foundation of the method for scatterers with arbitrary shape as well as a resolution analysis for its imaging functional. Furthermore, we establish an equivalent relation between the OSM and the DSM as well as resolution analysis for the latter. The factorization analysis for the far field operator as well as the Funk-Hecke formula play an important role in our theory. We also provide numerical results for three-dimensional anisotropic scatterers to validate the efficiency of the sampling methods.

We want to mention that the DSM in [17] was extended to the electromagnetic case in [1]. This extension relies on the Factorization method analysis in [15] for an isotropic medium and the decay rate of the imaging functional was not established. Another DSM which is related to the OSM was studied in [13] for small electromagnetic scatterers. We also refer to [12, 2, 8] and references therein for results on coefficient reconstruction for the isotropic inverse electromagnetic scattering problem.

We now introduce some basic notations for the paper. Let $\mathcal{O}\subset\mathbb{R}^{3}$ be a domain (connected and open) with Lipschitz boundary, we indistinctly denote by $(\cdot,\cdot)$ the inner product of $L^{2}(\mathcal{O})$ or $L^{2}(\mathcal{O})^{3}$ and by $\|\cdot\|$ the associated norms. We further denote

[TABLE]

where $H(\mathrm{curl}\,,\mathcal{O})$ is equipped by usual inner product

[TABLE]

For the following sections we will first rigorously formulate the direct and inverse electromagnetic scattering problems under consideration in section 2. Section 3 is dedicated to an analysis of the far field operator and the factorization method which is necessary for the study of the OSM. We establish the main theoretical results of the paper in section 4. More precisely, we define the imaging functional for the OSM, prove its resolution and stability, and an equivalent relation to the imaging functional of the DSM. Lastly, numerical examples are given where we reconstruct bounded anisotropic scatterers in $\mathbb{R}^{3}$ using the electromagnetic far field data. We see that the sampling methods are robust reconstruction methods that can recover scatterers of many different shapes and sizes.

2 Direct and inverse problem formulation

We consider the scattering of time-harmonic electromagnetic waves at positive frequency $\omega$ from a non-magnetic inhomogeneous medium. Suppose that there is no free charge and current density. Then, the $H_{\mathrm{loc}}(\mathrm{curl}\,,\mathbb{R}^{3})$ electric field ${\mathbf{E}}$ and the magnetic field ${\mathbf{H}}$ satisfy the Maxwell’s equations

[TABLE]

Here we assume that $\varepsilon$ is the electric permittivity and $\mu_{0}$ (positive constant) is the magnetic permeability of the medium. The permittivity $\varepsilon$ is assumed to be a bounded matrix-valued function. Let $\Omega$ be a bounded domain occupied by the non-magnetic inhomogeneous medium. The medium outside of $\Omega$ is assumed to be homogeneous. This means that there is a positive constant $\varepsilon_{0}$ such that $\varepsilon=\varepsilon_{0}I$ outside of $\Omega$, where $I$ is the $3\times 3$ identity matrix. We define the relative material parameter and the wave number as

[TABLE]

Eliminating magnetic field ${\mathbf{H}}$ from (1) we obtain

[TABLE]

The transmission conditions across the boundary of $\Omega$ are given by

[TABLE]

We denote by $\mathbf{F}_{+}$ and $\mathbf{F}_{-}$ the traces on $\partial\Omega$ from the exterior and interior of the domain $\Omega$ for a vector-valued function $\mathbf{F}$ respectively, and $\nu$ is the unit outward normal vector on $\partial\Omega$. Assume that we illuminate the inhomogeneous anisotropic medium with the electric and magnetic incident fields ${\mathbf{E}}_{\mathrm{in}}$ and ${\mathbf{H}}_{\mathrm{in}}$, respectively, satisfying

[TABLE]

Then there arises the scattered electric field $\mathbf{u}$, defined by $\mathbf{u}:={\mathbf{E}}-{\mathbf{E}_{\mathrm{in}}}$. Since the incident field ${\mathbf{E}_{\mathrm{in}}}$ satisfies the homogeneous Maxwell equation with wave number $k$ given by

[TABLE]

subtracting this equation from (2) we can conclude that the scattered field $\mathbf{u}\in H_{\mathrm{loc}}(\mathrm{curl}\,,\mathbb{R}^{3})$ is the solution to

[TABLE]

where the contrast $P$ is defined by

[TABLE]

Therefore, by definition we have that the support of the contrast $P$ is given by $\Omega$. Note that we also have the corresponding transmission conditions for the scattered field $\mathbf{u}$ following from (3). We complete the scattering problem by the Silver-Müller radiation condition for the scattered field $\mathbf{u}$ given by

[TABLE]

which is assumed to hold uniformly with respect to $\mathbf{x}/|\mathbf{x}|$.

It is known (see for e.g. [18]) that (4)–(5) is well-posed provided that the contrast $P$ is bounded with non-negative real and imaginary parts with support $\Omega$ provided that the only solution to the homogeneous problem (i.e. ${\mathbf{E}_{\mathrm{in}}}=0$) is trivial. We will assume that the homogeneous scattering problem only admits the trivial solution. This gives that the mapping ${\mathbf{E}_{\mathrm{in}}}\mapsto\mathbf{u}$ is linear and bounded from $L^{2}(\Omega)^{3}$ into $H_{\mathrm{loc}}(\mathrm{curl}\,,\mathbb{R}^{3})$.

Now we can talk about the inverse problem. To this end, we define $\widehat{\mathbf{x}}=\mathbf{x}/|\mathbf{x}|$,

[TABLE]

We consider the incident plane wave ${\mathbf{E}_{\mathrm{in}}}(\mathbf{x},\mathbf{d},\mathbf{q})=\mathbf{q}e^{ik\mathbf{x}\cdot\mathbf{d}}$, where the vector $\mathbf{d}\in\mathbb{S}^{2}$ indicates the direction of the incident propagation and $\mathbf{q}\in\mathbb{R}^{3}$ is the polarization vector such that $\mathbf{q}\cdot\mathbf{d}=0$. It’s well-known that we can express the corresponding scattered wave in terms of the asymptotic expansion

[TABLE]

uniformly in all observation directions $\widehat{\mathbf{x}}\in\mathbb{S}^{2}$. The function $\mathbf{u}^{\infty}(\widehat{\mathbf{x}},\mathbf{d},\mathbf{q})$ belonging to $L_{t}^{2}(\mathbb{S}^{2})$ for each incident and observation direction is called the far field pattern.

Inverse problem: Determine the shape and location of the scatterer $\Omega$ given the far field pattern $\mathbf{u}^{\infty}(\widehat{\mathbf{x}},\mathbf{d},\mathbf{q})$ for all $\widehat{\mathbf{x}},\mathbf{d}\in\mathbb{S}^{2}$ for a single wave number.

3 The far field operator and its factorization analysis

In this section, we will define and analyze the far field operator corresponding to (4)–(5). The analysis in this section will be used to derive sampling methods to solve the inverse shape problem of recovering the scatterer $\Omega$ from the far field data. It is well-known that the far field pattern $\mathbf{u}^{\infty}(\widehat{\mathbf{x}},\mathbf{d},\mathbf{q})$ is linear in $\mathbf{q}$ and can be written as

[TABLE]

where $\mathbf{u}^{\infty}(\widehat{\mathbf{x}},\mathbf{d})$ is a $3\times 3$ matrix and $\mathbf{q}\cdot\mathbf{d}=0$, see [7]. The following reciprocity relation is important in our analysis and its proof can be found in [7, Theorem 6.30].

Theorem 1**.**

For all $\widehat{\mathbf{x}},\mathbf{d}\in\mathbb{S}^{2}$, the far field pattern $\mathbf{u}^{\infty}(\widehat{\mathbf{x}},\mathbf{d})$ satisfies a reciprocity relation

[TABLE]

We now define the far field operator $\mathcal{F}:L^{2}_{t}(\mathbb{S}^{2})\to L^{2}_{t}(\mathbb{S}^{2})$ as

[TABLE]

In order to derive our sampling methods we will need to factorize the far field operator $\mathcal{F}$. To this end, it has been shown in [14] that the far field operator has the following factorization $\mathcal{F}=GH$. Here the operator

[TABLE]

and is the superposition of incident plane waves. It is easy to see that $H$ is compact and injective. The data to far field pattern operator

[TABLE]

where $\mathbf{v}\in H_{\mathrm{loc}}(\mathrm{curl}\,,\mathbb{R}^{3})$ is the solution to (4)–(5) with ${\mathbf{E}_{\mathrm{in}}}=\mathbf{f}$. The adjoint operator $H^{*}$ of $H$ is given by

[TABLE]

From the form of the data to far field pattern operator a factorization of $G$ in [14] gives that $G=H^{*}T$ with

[TABLE]

where again $\mathbf{v}\in H_{\mathrm{loc}}(\mathrm{curl}\,,\mathbb{R}^{3})$ is the solution to (4)–(5) with ${\mathbf{E}_{\mathrm{in}}}=\mathbf{f}$. Due to the well-posedness of the direct problem and the boundedness of $P$ it is clear that $T$ is a bounded linear operator on $L^{2}(\Omega)^{3}$, and hence $G=H^{*}T$ is also a compact operator. This gives that $\mathcal{F}=H^{*}TH$. This factorization is used to derive the Factorization method for reconstruction $\Omega$ in [14]. In order to continue we first make the following assumptions on the domain and coefficients.

Assumption 2**.**

We will assume that $\Omega$ is a Lipschitz bounded domain in $\mathbb{R}^{3}$. The contrast is assumed that $P\in L^{\infty}(\Omega)^{3\times 3}$. We lastly assume that the wave number $k$ is not a transmission eigenvalue i.e. the only solution to the homogeneous system in $L^{2}(\Omega)^{3}\times L^{2}(\Omega)^{3}$

[TABLE]

is trivial.

It is known that the set of real transmission eigenvalues is at most discrete for a real-valued permittivity and is empty for a complex-valued permittivity. See [3] and the references therein for the analysis of transmission eigenvalue problems. The following results can be found in [14]. These results will be critical in the later section where we study the orthogonality and direct sampling methods.

Theorem 3**.**

Let the operators $H$ and $T$ be as defined in (6)–(8) and let $\mathbf{p},\,\mathbf{y}_{s}\in\mathbb{R}^{3}$.

1. a.

The operator $H$ is compact and injective. 2. b.

The far field operator $\mathcal{F}$ is compact and injective. 3. c.

The far field operator has the factorization $\mathcal{F}=GH$ where $G=H^{*}T$. 4. d.

Let $\phi_{\mathbf{y}_{s}}(\mathbf{d})=(\mathbf{d}\times\mathbf{p})\times\mathbf{d}\,e^{-ik\mathbf{d}\cdot\mathbf{y}_{s}}$ then $\phi_{\mathbf{y}_{s}}\in\mathrm{Range}(H^{*})\iff\mathbf{y}_{s}\in\Omega$.

The following coercivity result from [10] is important for our theoretical analysis of the sampling methods.

Theorem 4**.**

Let the operator $T$ be as defined in (8). If the matrix-valued contrast $P\in L^{\infty}(\Omega)^{3\times 3}$ satisfies either

The imaginary part $\mathrm{Im}\,(P)$ is uniformly positive definite or 2. 2.

There is a constant $\alpha\geq 0$ such that $\mathrm{Re}\,(P)+\alpha\mathrm{Im}\,(P)$ is uniformly positive definite and $\mathrm{Im}\,(P)$ positive semidefinite.

Then provided that $k$ is not a transmission eigenvalue the operator $T$ is coercive on $\overline{\mathrm{Range}(H)}$.

Proof.

See the proof of Theorem 10 in [10] for details. ∎

4 Orthogonality sampling method

This section is dedicated to studying the OSM for solving the inverse problem. We first define the imaging functional for the OSM and then prove a resolution analysis as well as stability for the imaging functional which is the main theoretical result of this section. After that we show an equivalence between the imaging functionals of the OSM and the DSM.

The imaging functional. Let $\mathbf{y}_{s}$ be the sampling points in the imaging process and $\mathbf{p}\in\mathbb{R}^{3}$ is a fixed vector. We are interested in imaging of the scatterer $\Omega$ given the far field pattern $\mathbf{u}^{\infty}(\widehat{\mathbf{x}},\mathbf{d},\mathbf{q})$, for all $\widehat{\mathbf{x}},\mathbf{d}\in\mathbb{S}^{2}$ and polarization vector $\mathbf{q}=(\mathbf{d}\times\mathbf{p})\times\mathbf{d}$. We define the imaging functional $\mathcal{I}_{OSM}$ for the OSM as

[TABLE]

The use of $(\mathbf{d}\times\mathbf{p})\times\mathbf{d}$ in the imaging functional is useful for the analysis in this section and to have the polarization vector belonging to $L^{2}_{t}(\mathbb{S}^{2})$ for any choice of $\mathbf{p}\in\mathbb{R}^{3}$. We now wish to write the imaging functional in terms of the far field operator. This will allow us to study the the imaging functional using the factorization analysis in the previous section.

Lemma 5**.**

The imaging functional $\mathcal{I}_{OSM}$ satisfies

[TABLE]

where $\phi_{\mathbf{y}_{s}}\in L^{2}_{t}(\mathbb{S}^{2})$ is given by

[TABLE]

Proof.

We use the reciprocity relation and interchange the roles of $\widehat{\mathbf{x}}$ and $\mathbf{d}$

[TABLE]

Now using the identity

[TABLE]

and the fact that $\mathcal{F}\phi_{\mathbf{y}_{s}}(\widehat{\mathbf{x}})\cdot\widehat{\mathbf{x}}=0$, proves the claim. ∎

We have shown that the imaging functional can be represented in terms of the far field operator. We now turn our attention to studying its properties. Just as in the previous section we will assume that the contrast satisfies the assumptions of Theorem 4. This is to insure the coercivity property of the middle operator $T$.

Theorem 6**.**

For every $\mathbf{y}_{s}\in\Omega$ the imaging functional $\mathcal{I}_{OSM}(\mathbf{y}_{s})$ is bounded from below by a positive constant. Moreover, for $\mathbf{y}_{s}\notin\Omega$ the imaging functional satisfies

[TABLE]

Proof.

We first observe that using again $\mathcal{F}\phi_{\mathbf{y}_{s}}(\widehat{\mathbf{x}})\cdot\widehat{\mathbf{x}}=0$ and the identity $(\widehat{\mathbf{x}}\times\mathbf{p})\times\widehat{\mathbf{x}}=\mathbf{p}-\widehat{\mathbf{x}}(\widehat{\mathbf{x}}\cdot\mathbf{p})$ we have

[TABLE]

Therefore, Lemma 5 and Cauchy-Schwarz inequality deduce that

[TABLE]

Here $|\mathbb{S}^{2}|$ is the surface area of $\mathbb{S}^{2}$. From the factorization $\mathcal{F}=H^{*}TH$ and the coercive property of $T$ we have that

[TABLE]

where $C$ is some positive constant. Now let $\mathbf{y}_{s}\in\Omega$ then we have that $\phi_{\mathbf{y}_{s}}\in\mathrm{Range}(H^{*})$ by Theorem 3. Therefore, we have $\phi_{\mathbf{y}_{s}}=H^{*}\varphi_{\mathbf{y}_{s}}$ for some $\varphi_{\mathbf{y}_{s}}\neq 0$, and

[TABLE]

proving the first statement of the theorem.

We now show that the imaging functional decays as $\text{dist}(\mathbf{y}_{s},\Omega)\to\infty$. To do so, we first observe from Lemma 5 and the factorization $\mathcal{F}=GH$ that

[TABLE]

Therefore, we can show that \big{\|}H\phi_{\mathbf{y}_{s}}\big{\|}^{2} satisfies the decay property and use the upper bound on the imaging functional given above. Here we let $Y^{m}_{\ell}$ denote the spherical harmonics which form a complete orthonormal system on $L^{2}(\mathbb{S}^{2})$. Now, recall the Funk-Hecke formula (see for e.g. [7])

[TABLE]

where $j_{\ell}$ is the first kind spherical Bessel function of order $\ell$. In particular,

[TABLE]

Just as in recent works [11, 19] we will use the Funk-Hecke formula to show that \big{\|}H\phi_{\mathbf{y}_{s}}\big{\|}^{2} decays as $\text{dist}(\mathbf{y}_{s},\Omega)\to\infty$. Indeed, using the Funk-Hecke formula for $m=\ell=0$ with the formula $\mathrm{curl}\,_{\mathbf{z}}\mathrm{curl}\,_{\mathbf{z}}(\mathbf{p}e^{-ik\mathbf{d}\cdot\mathbf{z}})=-k^{2}(\mathbf{d}\times\mathbf{p})\times\mathbf{d}e^{-ik\mathbf{d}\cdot\mathbf{z}}$ straightforward calculations gives

[TABLE]

where, for $j=1,2,3$,

[TABLE]

It is obvious that $v_{j}=O(1/|\mathbf{z}|)$ as $|\mathbf{z}|\to\infty$. Therefore

[TABLE]

which leads to

[TABLE]

proving the claim. ∎

This resolution analysis implies the for any sampling points $\mathbf{y}_{s}\in\Omega$ the imaging functional is strictly positive and will decay as $\mathbf{y}_{s}$ moves away from the scatterer $\Omega$. Therefore, one can plot $\mathcal{I}_{OSM}(\mathbf{y}_{s})$ in order to recover the scatterer. Using the imaging functional $\mathcal{I}_{OSM}(\mathbf{y}_{s})$ has the advantage that one does not have to solve an ill-posed equation at each sampling point. Other sampling methods such as the Linear sampling method [3] and Factorization method [15] requires one to solve an ill-posed equation involving the far field operator at each sampling point. This requires one to compute the singular-value decomposition of the far field operator as well as applying a regularization scheme. Also stability of these methods with respect to noise added to the data is not justified where as the orthogonality sampling method only requires one to compute an inner-product and by the following theorem we see that the imaging functional is stable with respect to noise added to the far field data.

In many physical applications one only knows the far field data up to some ‘small’ perturbation. Now, assume that we only know the far field operator up to a perturbation, that means we have access to its noisy version $\mathcal{F}_{\delta}$, which satisfies

[TABLE]

for some $\delta>0$. Here $\mathcal{F}_{\delta}$ represents the measured far field operator physical experiments. We now give a stablity estimate for the imaging functional $\mathcal{I}_{OSM}(\mathbf{y}_{s})$.

Theorem 7** (stability estimate).**

Denote by $\mathcal{I}_{OSM,\delta}$ the imaging functional corresponding to noisy far field operator $\mathcal{F}_{\delta}$. Then

[TABLE]

where $|\mathbb{S}^{2}|$ is again the surface area of $\mathbb{S}^{2}$, and $\mathbf{p}\in\mathbb{R}^{3}$ is the polarization vector in the definition of the imaging functional (11).

Proof.

Using Lemma 5, the triangle inequality and the Cauchy-Schwarz inequality we have

[TABLE]

proving the theorem. ∎

Motivated by the work in [17] we prove that the imaging functionals for the OSM and DSM are equivalent. The DSM was considered in [1] for the far field operator corresponding to the magnetic field. In their analysis they assume that the contrast is scalar-valued and use the factorization established in [15]. The analysis presented here is valid for a matrix-valued contrast. Moreover, we give an explicit decay rate for the imaging function. To this end, we again let the sampling points $\mathbf{y}_{s}$ and fixed vector $\mathbf{p}\in\mathbb{R}^{3}$, the imaging functional for the DSM is defined as

[TABLE]

where

[TABLE]

Now, similarly to Theorem 6 we will represent the functional in terms of the far field operator. This will be used to prove the equivalence with OSM.

Lemma 8**.**

The imaging functional $\mathcal{I}_{DSM}$ for the DSM satisfies

[TABLE]

Proof.

We use the reciprocity relation and interchange the roles of $\widehat{\mathbf{x}}$ and $\mathbf{d}$

[TABLE]

proving the claim. ∎

From Theorem 6 and 8 we have all we need to derive the equivalence of the two imaging functionals studied in this section. Again we assume that the contrast satisfies the assumptions of Theorem 4. The equivalence of the two sampling methods is proven in the following result.

Theorem 9**.**

There exists positive constants $c_{1}$ and $c_{2}$ such that

[TABLE]

Proof.

From Lemma 8, the factorization $\mathcal{F}=H^{*}TH$ and the coercive property of $T$ we obtain

[TABLE]

where $C$ is some positive constant. The first inequality of the theorem hence follows from the estimate

[TABLE]

The following estimate is from the beginning of the proof of Theorem 6

[TABLE]

This estimate and Lemma 8 implies the second inequality of the theorem. ∎

Following the method for proving Theorem 6 we can obtain similar resolution analysis for the direct sampling method.

Theorem 10**.**

For every $\mathbf{y}_{s}\in\Omega$ the imaging functional $\mathcal{I}_{DSM}(\mathbf{y}_{s})$ is bounded from below by a positive constant. Moreover, for $\mathbf{y}_{s}\notin\Omega$ the imaging functional satisfies

[TABLE]

We also have the stability of the DSM that can be proved using the triangle and the Cauchy-Schwarz inequalities. The analysis in this section allows one to solve the inverse shape problem for electromagnetic scattering by plotting either $\mathcal{I}_{DSM}$ or $\mathcal{I}_{OSM}$. This amounts to a fast yet stable reconstruction method that only relies on the knowledge of the far-field data. Unlike some traditional sampling methods one does not need to minimize a non-linear functional at each sampling point.

5 Numerical examples

We present in this section several numerical examples to validate the performance of the OSM as well as the DSM. The simulations were carried on a Quad Core 3.6 GHz machine with 32GB RAM and the implementation was done using the computing software Matlab. The synthetic data are generated by numerically solving the direct problem with the spectral solver studied in [19]. We solve the direct problem (4)–(5) with incident field

[TABLE]

where $\mathbf{p}=(1/\sqrt{3},-1/\sqrt{3},1/\sqrt{3})^{\top}$ and $N_{\mathbf{d}}$ is the number of directions $\mathbf{d}_{j}$ that is specified below for each numerical example. Likewise, we denote by $N_{\widehat{\mathbf{x}}}$ the number of points ${\widehat{\mathbf{x}}}\in\mathbb{S}^{2}$ where the far field pattern data are collected. The points $\widehat{\mathbf{x}}$ and $\mathbf{d}$ that are chosen for generating the scattering data are almost uniformly distributed on $\mathbb{S}^{2}$. Let

[TABLE]

be our synthetic far field data that has three components $\mathcal{D}_{n},n=1,2,3,$ and $\mathcal{D}_{n}$ can be considered as an $N_{\widehat{\mathbf{x}}}\times N_{\mathbf{d}}$ matrix. To consider noisy data, we add artificial noise to our synthetic data. More precisely, a complex-valued noise matrix $\mathcal{N}$ containing random numbers that are uniformly distributed in the complex square

[TABLE]

is added to the data matrix $\mathcal{D}_{n}$. Denoting by $\delta$ the noise level, the noisy data matrix $\mathcal{D}_{n,\delta}$ is then given by

[TABLE]

where $\|\cdot\|_{2}$ is the matrix 2-norm. To define the anisotropic contrasts for the numerical examples, we need the following diagonal matrix

[TABLE]

5.1 OSM and DSM (Figure 1).

We present in Figure 1 the reconstruction results of the OSM and DSM for an anisotropic scatterer including three different balls. The scatterer is characterized by smoothly varying contrast $P(\mathbf{x})$ defined by

[TABLE]

where $A$ is given by (16).

Here we choose wave number $k=12$ (the corresponding wavelength is about 0.52) and $N_{\widehat{\mathbf{x}}}\times N_{\mathbf{d}}=325\times 325$. The far field data are perturbed by 30$\%$ noise. In Figure 1 we also present 3D-visualizations of the exact geometry and its reconstructions using isosurface in Matlab. The isovalue for the isosurface plotting is chosen to be $1/3$ of the maximal value of the computed imaging functionals $\mathcal{I}_{OSM}$ and $\mathcal{I}_{DSM}$. It can be seen in Figure 1 that both the OSM and the DSM are robust with noise added to the data and provide reasonable reconstructions of the scatterer. We also observe that the reconstruction from $\mathcal{I}_{OSM}$ is more similar to that from $(\mathcal{I}_{DSM})^{2}$ rather than the one from $\mathcal{I}_{DSM}$. This might reflect the estimate $\mathcal{I}_{DSM}\leq c_{2}\sqrt{\mathcal{I}_{OSM}}$ from the equivalent relation in Theorem 9.

5.2 OSM for highly noisy data (Figure 2).

The second example is presented in Figure 2 where we focus on performance of the sampling method on data perturbed by high amounts of noise.

Here consider the scatterer as in the previous example (17), $k=12$ and $N_{\widehat{\mathbf{x}}}\times N_{\mathbf{d}}=325\times 325$. From Figure 2 we can see that, even there are high levels of noise $\delta=0.6$ and $0.9$ in the far field data, that the OSM is still able to provide reasonable reconstructions. The computed images are not very different for 60$\%$ and 90$\%$ noise in the data. The solid performance of the method on noisy data can be justified by the stability of the method that is discussed in Theorem 7.

5.3 OSM with less data (Figure 3).

This is the focus of the third example that is to examine the performance of the OSM on a smaller wave number and a smaller set of scattering data. We consider the same scatterer as in the first and second examples (17). The data in this example are perturbed by 30$\%$ noise. We can see in Figure 3(b) that for $N_{\widehat{\mathbf{x}}}\times N_{\mathbf{d}}=325\times 325$ and wave number $k=6$ (wavelength is about 1.04) the reconstruction is not as good as that of the case $k=12$ (Figure 3(c)). Although we can see two components of the scatterers in the reconstruction the shape and locations are not very accurate. The result is even worse if we have less data, see Figure 3(d). More precisely, for $N_{\widehat{\mathbf{x}}}\times N_{\mathbf{d}}=91\times 91$ and $k=12$, the reconstruction is no longer reasonable.

5.4 OSM for other types of scatterers (Figure 4).

In this example we present in Figure 4 the reconstruction results for three different types of scatterers. We consider scatterers with more complicated shapes and non-smooth geometries. For the scatterer in Figure 4(a) the anisotropic contrast $P(\mathbf{x})$ is again a smoothly varying function defined by

[TABLE]

where matrix $A$ is given by (16). The contrast $P(\mathbf{x})$ for the scatterers in Figures 4(d) and 4(g) are respectively equal to $A/2$ and $A/4$ in $\Omega$ and 0 outside of $\Omega$. Again, $N_{\widehat{\mathbf{x}}}\times N_{\mathbf{d}}=325\times 325$, $k=12$ and the data are perturbed by 30$\%$ noise. The pictures show that with the right amount of data the sampling method is able to provide good reconstruction results for scatterers with more complicated shapes. We also observe that for scatterers with non-smooth geometries like in Figures 4(d) and 4(g), the isovalue in the isosurface plotting should be 1/2 of the maximal value of the computed imaging functional $\mathcal{I}_{OSM}$ to give a better three-dimensional image.

6 Summary

We study the OSM for solving the electromagnetic inverse scattering problem for anisotropic media with far field data. We propose an imaging functional for the OSM which is able to compute the location and shape of electromagnetic scatterers in a fast and robust way. Using tools of the Factorization method analysis and the Funk-Hecke formula we are able to establish a rigorous justification and resolution analysis for the proposed imaging functional. We also prove that this functional is equivalent to that of the DSM and that our resolution analysis for the OSM can be directly applied to the DSM. Numerical results for three-dimensional anisotropic scatterers are presented. Together with our recent work in [19], where the OSM is justified for the electromagnetic inverse scattering with one incident plane wave, we have provided a versatile approach for solving the electromagnetic inverse scattering problem.

Acknowledgement. The work of DLN is partially supported by NSF grant DMS-1812693.