1 Introduction
Abstract

In this paper, we consider two time-harmonic inverse scattering problems of reconstructing penetrable inhomogeneous obstacles from near field measurements. First we appeal to the Born approximation for reconstructing small isotropic scatterers via the MUSIC algorithm. Some numerical reconstructions using the MUSIC algorithm are provided for reconstructing the scatterer and piecewise constant refractive index using a Bayesian method. We then consider the reconstruction of an anisotropic extended scatterer by a modified linear sampling method and the factorization method applied to the near field operator. This provides a rigorous characterization of the support of the anisotropic obstacle in terms of a range test derived from the measured data. Under appropriate assumptions on the material parameters, the derived factorization can be used to determine the support of the medium without a-priori knowledge of the material properties.

Near field imaging of small isotropic and extended anisotropic scatterers

Isaac Harris111Department of Mathematics, Texas A M University, College Station, Texas 77843-3368, (Corresponding Author) E-mail: iharris@math.tamu.edu and Scott Rome222Email: romescott@gmail.com

Keywords: inverse scattering, factorization method, MUSIC algorithm.

AMS subject classifications: 35J05, 45Q05, 78A46

1 Introduction

Qualitative/Sampling methods have been used to solve multiple inverse problems such as parameter identification and shape reconstruction. These methods have been used to solve inverse boundary value problems for elliptic, parabolic and hyperbolic partial differential equation. See [21] and [20] for examples of the linear sampling methods applied to parabolic and hyperbolic equations, respectively. The vast available literature is a representative of the many directions that this research has taken (see [8], [15], [16], [25] and the references therein).

The factorization method is a qualitative method which develops a range test to determine the support of the unknown scatterer. In general, we have that

where is known and is a positive selfadjoint compact operator defined by the measurements. By appealing to the range test, we have that any two obstacles admitting the same data must be equal, proving uniqueness of the inverse problem. Since the support of the scatterer is connected to the range of a compact operator which is known, evaluating the indicator function derived from the factorization method amounts to applying Picard s criteria, which only requires the singular values and functions of a known operator. This makes the factorization method computationally cheap to implement and analytically rigorous. Using the characterization of the obstacle given by the factorization method we propose a modified linear sampling method as an alternative sampling algorithm for reconstructing the obstacle. The indicator function proposed for the modified linear sampling method can be shown to satisfy a similar analytically rigors characterization of the obstacle.

In this paper, we consider two problems associated with inverse scattering. The first problem we consider is the reconstruction of small isotropic obstacles and the refractive index from near field data. To do so, we derive the Multiple-Signal-Classification (MUSIC) algorithm which gives an analytically rigorous method of reconstructing small obstacles. Once we have reconstructed the obstacles, we will give a few example of how one can reconstruct constant refractive indices using a Bayesian method. Next we consider the reconstructing an anisotropic obstacle from near field data via the factorization method and modified linear sampling method. This problem has been considered in [9] and [26] for far field measurements using the linear sampling and factorization methods, respectively. The factorization method for an inhomogeneous isotropic media with near field measurements was studied in [22]. We analyze the factorization method applied to the near field operator for an anisotropic obstacle to derive a simple numerical algorithm to reconstruct the scatterer. For the case of an anisotropic media, the unique determination of the support is optimal because a matrix-valued coefficient is in general not uniquely determined by the scattering data (see for e.g. [18]). The analysis done in this paper is of a similar flavor to the work done in [12] and [29].

2 Recovering small volume isotropic scatterers and the refractive index

2.1 Scattering by small isotropic scatterers

In this section, we formulate the direct time-harmonic scattering problem in for small isotropic scatterers. The scattering obstacle may be made up of multiple simply components, but is not necessarily simply connected. We are interested in using near field measurements for ‘small’ obstacles embedded in a homogenous media. Let be the incident field that originates from the source point located on the curves . We will assume that the curve is class and that the incident field satisfies the Helmholtz equation in for or . Now let

where is a small parameter and is a domain with a piece-wise smooth boundary that is centered at the origin. Outside of the scattering obstacle(s) , the material parameters are homogeneous isotropic with refractive index scaled to one. Here we denote the refractive index of the homogenous background with the (possibly inhomogeneous) obstacle in by

where is the characteristic function. The refractive index of the obstacle(s) is represented by a continuous scalar function such that the real and imaginary parts satisfy

for almost every where we assume that or . Now the radiating scattered field given by the point source incident field is the unique solution to

(1)
(2)

where and the radiation condition (2) is satisfied uniformly in all directions . Recall the fundamental solution to Helmholtz equation

where is the first kind Hankel function of order zero. The corresponding total field for the scattering problem is given by satisfies

where the superscripts and for a generic function indicates the trace on the boundary taken from the exterior or interior of its surrounding domain, respectively. We have that the scattered field is given by the solution to the Lippmann-Schwinger Integral Equation (see for e.g. [15])

Assume that is sufficiently small such that

Therefore, we can conclude that the Born approximation, which is the first term in the Neumann series solution to the Lippmann-Schwinger Integral Equation, given by

(3)

is a ‘good’ approximation for the scattered field . Therefore, we assume that we have the measured near field data that is given by the Born approximation . Let be the bounded simply connected open set such that satisfies and . Now assume that we have the data set of near field measurements for all . The inverse problem we consider is to reconstruct the scatterer and the refractive index from a knowledge of the Born approximation of the near field measurements.

2.2 The MUSIC Algorithm

The MUSIC Algorithm can be considered as the discrete analogue of the factorization method (see for e.g. [14]). In particular, we connect the locations of the obstacles given by to the range of the so-called multi-static response matrix denoted by that is defined by the near field measurements for all . To arrive at the MUSIC algorithm, we will need to define our measurements operator. To this end, we let the incident field be give by . Notice that by (3) we obtain that

since the integrand is continuous with respect to the variable . Assume we have a finite number of incident and observation directions where , with being the number of small obstacles and being the number of incident and observation directions given by . Now define the multi-static response matrix

Notice that for sufficiently small we have that the multi-static response matrix is approximated by therefore we assume that is known. For convenience, we will assume that the sources and receivers are placed at the same points (i.e. ). We now factorize the multi-static response matrix, therefore we define the matrices and , where the matrix is given by

and the matrix where . Therefore, we have that and it follows that

(4)

Now define the vector for any point by

(5)

The goal of this section is to show that is in the range of the matrix if and only if , which is a discrete reformulation of the result of the factorization method stated in the Introduction. Since we are interested in finding obstacles , it is sufficient to prove the result only for values of . Recall that without loss of generality, we assume that the sources and receivers are placed at the same locations. We will follow the analytic framework in Section 4.1 of [25]. We now give a result that can be proven by using standard arguments from linear algebra (see proof of Theorem 3.1 in [17] for details).

Lemma 2.1.

Let the matrix have the following factorization where and with . Assume that

  1. the matrix has full rank

  2. the matrix is invertible

then RangeRange.

We now construct an indicator function derived from Lemma 2.1 to reconstruct the locations of the scattering obstacles. To this end, for each sampling point we will show that the vector is in the range of if and only if . We now prove an auxiliary result that connects the location of the scattering obstacles to the range of the matrix .

Theorem 2.1.

Assume that the set is dense in . Let then there is a number such that for all , we have that

  1. the matrix has full rank,

  2. Range if and only if .

Proof.

It is clear that is in the range of since is the -th column of To prove that for that is not in the range of for some sufficiently large, we proceed by way of contradiction. Let and assume that there does not exist such a , then we have that there exists , such that

(6)

and

We then conclude that (up to a subsequence) , , as . This gives that due to the density of and since is analytic with respect to for , we obtain that

Since the mapping

is a radiation exterior solution to the Helmholtz equation in we conclude that

by uniqueness to the exterior Dirichlet problem and unique continuation. Now letting gives that for , but by (6) and the convergence of the coefficients , as we have that

which gives a contradiction. Moreover, the fact that has full rank is a consequence of the given arguments. ∎

Now by combining this with Lemma 2.1 and Theorem 2.1 we have the following range test which can be used to reconstruct the set

Theorem 2.2.

Assume that the set is dense in . If then there is a number such that for all

Now we are ready to construct an indicator function which characterizes the locations of the small obstacles. To this end, let be the orthogonal projection onto . Therefore, by Theorem 2.2 we have that

Notice that since is a self-adjoint matrix we have that it is orthogonally diagonalizable. So we let be the -th orthonormal eigenvector of , we also let =Rank. This gives that for is an orthonormal basis for , and therefore we have that is given by

This leads to the following result:

Lemma 2.2.

Let be the -th orthonormal eigenvector of and let =Rank. Assume that the set is dense in . If then there is a number such that for all

(7)

2.3 Numerical validation of the MUSIC algorithm

We now present a few examples of reconstructing small obstacles using the MUSIC algorithm developed in this section in . To simulate the data, we use the Born approximation of the scattered field and a Gaussian quadrature to approximate the volume integral. We now let the collection curve be the boundary of a disk of radius one. In the following calculations we use different sources and receivers at the points

This leads to an approximated multi-static response matrix

where and are the weights and quadrature points for the numerical approximation of the volume integral by a 2 Gaussian quadrature method. We give examples with random noise added to the simulated data for . The random noise level is given by where the noise is added to the multi-static response matrix such that

The set can be visualized by plotting the corresponding indicator function given in (7)

where are the eigenvectors of the approximated matrix and =Rank is computed by the Matlab command rank. In the following examples, we fix the wave number .

Figure 1: Reconstruction of the location of the two scatterers, the circle of radius 0.2 centered at and the ellipse with and centered at . The refractive index in both scatterers. Contour plot of the indicator function where the red lines are the boundary of the scatterer.
Figure 2: Reconstruction of the location of one scatterer, the ellipse with and centered at . The refractive index in the scatterer. Contour plot of the indicator function where the red lines are the boundary of the scatterer.
Figure 3: Reconstruction of the location of the a square scatterer, . The refractive index in the scatterer. Contour plot of the indicator function where the red lines are the boundary of the scatterer. Remark: Even though both figures look identical, different data we used in the reconstructions.

2.4 Bayesian approach to approximating the Refractive Index

In this section, we consider the the second part of our inverse problem which is to reconstruct (i.e. approximate) the refractive index of the scatterer(s). In particular, the idea is to use a hierarchical Bayesian model for estimating the refractive index by a constant in the scatterer(s) which has been reconstructed, possibly by the MUSIC algorithm proposed in the previous section. Since we assume that is continuous in the scatterer the small size of the obstacle implies that there will be only small deviations from the value at the center of the obstacle.

Here we assume that we that an approximation of the scatterer is known using the reconstruction from the previous section or by another qualitative method. We use a Bayesian collocation method to approximate the value of the refractive index at the center of the obstacle. For this discussion, assume that the obstacle is centered at the point . Rather than solving for the value of the unknown exactly, we solve for a probability distribution of the possible values of . To do this, we must assume a prior distribution for possible values of and a data model for the scattered wave, which we will discuss below. Then, Markov chain Monte Carlo methods allow us to sample from the posterior distribution of given the scattered wave and the prior distribution. The mode of the resulting posterior distribution (called the Maximum a posteriori estimator or MAP). The MAP in certain settings can be shown to be equivalent to the standard linear regression estimator or the frequentist statistics maximum likelihood estimator. In this context, our approach provides more information on possible values of , and because the scatterer is small with a continuous refractive index.

To proceed, we define the necessary notation. Assume we have a reconstruction of via the MUSIC algorithm, denoted . Furthermore, we assume that the measured scattered field denoted is given by the Born approximation on , that is

where is the amount of error added for the reconstruction of and standard measurement error. We assume the error term is modeled by a normal distribution with zero mean and standard deviation , then we may write the data model, as a statistical distribution where values of are drawn from

(8)

Now, recall that we have readings of the scattered field at points with the location of the point source given by with being the collection curve discussed in the previous section. This forms the data set where indicates the number of sources and receivers.

We now motivate the subsequent Bayesian model. Since the goal to reconstruct the value of the refractive index at the center of the scatterer, using Taylor’s Theorem we have that where is small. Therefore, the constant term of the Taylor polynomial centered at the center of is approximately close to all values of on . Based on these assumptions, we define the following auxiliary function

(9)

Therefore, we describe the following Bayesian model to solve for the posterior distribution of given readings of the scatterer , i.e. we solve for the conditional pdf function

where and are the weights and quadrature points for the numerical approximation of the volume integral over . Here we assume that the error from the numerical integration is negligible since one can take a quadrature rule to approximate the integral to some preferred tolerance. By the definition of we have that

We may solve the model for posterior distribution for , which is given according to the proportion:

(10)

where is the likelihood function arising from the data model for and are the assumed priors given above. One may marginalize over obtain if desired. By using the Metropolis-Hastings algorithm, we can derive this posterior distribution, and thus values for . We note that this model is hierarchical with care taken to model the sources of error. This should lead to a more realistic solution.

For our calculation, we use the implementation in the Python package PyMC3 [32]. Therefore we will have which is the average value of which will be the approximation on . For our numerical example, it is a common practice to use a normal distribution with a large standard deviation () rather than a flat prior, . For our example we consider the square scatterer with refractive index where we add noise to the measurements. For the reconstruction we use the domain actual scatterer as well as which is an approximation taken from the MUSIC Algorithm’s reconstruction of the obstacle. Using the methods described in this section we are able to determine the distribution of possible values of which should approximately be .

Figure 4: Recovery from measurements in the case when . The distribution is the possible values of in the reconstructed scatterer .
Figure 5: Recovery from measurements in the case when . The distribution is the possible values of in the reconstructed scatterer .

Similar approaches utilizing Bayesian statistics and the Born approximation have been used when applying compressive sensing techniques to inverse scattering [33],[34],[35],[36]. In those context, there is an a-priori assumption on the sparsity of the scatterer. The approach has also been applied to continuous random media [37]. However, in this case, the scatterer is not sparse and it is not a continuous random medium. Moreover, we extend the approach by using the MUSIC algorithm to find the support of the scatterer, and use that information in our recovery of the refractive index.

3 Sampling methods for anisotropic scatterers

3.1 Scattering by extended anisotropic scatterers

In this section we formulate the direct and inverse time-harmonic scattering problems under consideration. To this end, we let the point source located at the point be given by the fundamental solution to Helmholtz equation where is assumed to be a closed curve (in ) or surface (in ) that is class or sufficiently smooth, such that is compactly embedded in . We consider the scattering of the (non-standard) incident field for by a penetrable anisotropic obstacle. Assume that the obstacle (for ) is a bounded simply connected open set having piece-wise smooth boundary with being the unit outward normal to the boundary.

We assume that the constitutive parameters of the scatterer are represented by a symmetric matrix and a scalar function such that

where as

for almost every and all . Outside the obstacle , the material parameters are homogeneous isotropic with refractive index scaled to one. We denote by and the constitutive parameters of the homogenous background with the anisotropic obstacle in by

where denotes the identity matrix and is the characteristic function on . Note that the support of the contrasts and is . The radiating scattered field given by the (non-standard) point source satisfies

(11)
(12)

where and (12) is satisfied uniformly in all directions . It can be shown that the scattering problem (11)-(12) has a unique solution by using a variational approach (see for e.g. [8]). The total field corresponding to the scattering problem (11) is defined as and satisfies

where the superscripts and for a generic function indicates the trace on the boundary taken from the exterior or interior of the domain, respectively. Now assume that we have the data set of near field measurements for all and define the Near Field operator

The inverse problem we consider is to reconstruct the scattering obstacle from a knowledge of for all . We will let be the bounded simply connected open set such that and .

In the coming sections, we will see that to obtain a symmetric factorization of the near field operator we require that the incident wave be the complex conjugate of a point source. However, it is known that these non-physical sources can be approximated by a linear combination of physical point sources (see [23]). In [22], the authors are able to avoid using non-physical sources but must construct the so-called outgoing-to-incoming operator. If is given by the boundary of a disk/sphere centered at the origin, then the outgoing-to-incoming operator is given by a boundary integral operator that does not depend on the type of scatterer.

3.2 Factorization of the near field operator

This section is dedicated to constructing a suitable factorization of the near field operator so that we are able to use the theoretical results in [24] and/or [27], in order to derive indicator functions for the support of in terms of the near field measurements. Here we consider a modified linear sampling method as well as the factorization method and just as in [4] we will us the range characterization from the factorization method to show that these two methods give equivalent indicator functions. To this end, motivated by the source problem given in equation (11)-(12) for the scattered field due to the (non-standard) point source, we consider the problem of finding for a given such that

(13)

It can be shown that for any given , the solution operator that maps is a bounded linear mapping from to . At this point, let us recall the exterior Dirichlet-to-Neumann mapping for Helmholtz equation

where satisfies

along with the radiation condition (12) where . With the help of the Dirichlet-to-Neumann operator we can write (13) in its equivalent variational form: find such that

(14)

where the boundary integral over is understood as the dual pairing between and . It is clear by well-posedness that if we take then we have that the scattered field for all . Now define the source to trace operator

as well as the bounded linear operator

By superposition, we have that near field operator is given by In order to use the theory developed in [25] and/or [27], we must construct a symmetric factorization of the near field operator. Therefore, we compute the adjoint of the operator which is given in the following result.

Theorem 3.1.

The operator is given by

where is the unique radiating solution to

(15)

for all .

Proof.

Let the radius of be sufficiently large so that . Given any we can construct a unique radiating field that satisfies (15). We now let be defined as

Therefore, we have that

Using integration by parts on the volume integrals and the fact that solves Helmholtz equation in we have that

Now using the jump relations for the normal derivative of the single layer potential for Helmholtz equation gives that

Notice that both and satisfy the radiation condition (12) therefore letting gives that

proving the claim. ∎

Just as in [12] notice that (13), implies that

(16)

where is the solution to (13) for a given . The variational form of (16) is given by

Now appealing to the Riesz representation theorem, we can define the bounded linear operator such that for all

(17)

Notice for the unique solution to (13) that by the definition of we have that . Now by the definition of and we conclude that . We conclude that for all . Now recalling that gives the desired factorization.

Lemma 3.1.

The near field operator for the scattering problem (11)-(12) has the factorization .

3.3 The solution to the inverse problem

The main goal of this section is to connect the support of to the range of an operator defined by the measured near field operator. We make this connection by analyzing the factorization of the near field operator developed in the previous section. We recall from the previous section that we have the following factorization . Under some appropriate assumptions on and , the factorization method gives that Range=. Here will be given by for a non-absorbing material or with for an absorbing material where

Furthermore for a generic self-adjoint compact operator on a Hilbert space , by appealing to the Hilbert-Schmidt theorem we can define in terms of the spectral decomposition such that

for all where is the orthonormal eigensystem of .

3.3.1 Analysis of the operator

Now we turn our attention to showing that the operator has the properties needed to apply Theorems A.1. To this end, let us define the interior transmission problem as finding a pair such that for a given pair satisfies

(18)
(19)
Definition 3.1.

The values of for which the homogeneous interior transmission problemi.e. (18)-(19) with , has nontrivial solutions are called transmission eigenvalues for and .

It is known that under appropriate assumptions on the coefficients that the interior transmission problem (18)-(19) is Fredholm of index zero, i.e if is not a transmission eigenvalue then (18)-(19) is well-posed (see for e.g. [7] and [11]). The following results are know for the transmission eigenvalue problem (see [7], [8] and [13])

  1. If the scatterer is absorbing such that and in then there are no real transmission eigenvalues.

  2. If the scatterer is non-absorbing (i.e. and ) then the set of real transmission eigenvalues is at most discrete provided that is positive (or negative) definite uniformly in .

  3. If the scatterer is non-absorbing then the set of real transmission eigenvalues is at most discrete provided that is positive (or negative) definite uniformly and is positive (or negative) in a neighborhood of the boundary .

Assumption 3.1.

The wave number is not a transmission eigenvalue for (18)-(19).

We are now ready to connect the support of to the range of . Since we assume that the obstacle is contained in the interior of we only need to consider sampling points .

Theorem 3.2.

Under Assumption 3.1, we have that for

Proof.

We start with the case where and assume that is such that . By the definition of there exists a satisfying (15) which implies that