Inverse medium scattering for three-dimensional time harmonic Maxwell equations

Gang Bao^1,2 and Peijun Li¹

¹Department of Mathematics, Michigan State University, East Lansing, MI, 48824-1027, USA
²School of Mathematics, Jilin University, Changchun 130023, People's Republic of China

Email: bao@math.msu.edu and lipeijun@math.msu.edu

Received 16 December 2003
Published 22 January 2004

Abstract. A continuation method is developed for solving the inverse medium scattering problem of time harmonic Maxwell equations in . By using multi-frequency scattering data, our reconstruction algorithm first employs the Born approximation for an initial guess and proceeds via recursive linearization on the wavenumber k. At each linearization step, one forward and one adjoint state of the Maxwell equations are solved. Numerical examples are also presented and discussed.

1. Introduction

Consider the time harmonic Maxwell equations in three dimensions

where k is the wavenumber, and q(x)> - 1, which has a compact support, is the scatterer. The total electric field E^t consists of the incident field Eⁱ and the scattered field E:

Assume that the incident field is a plane wave

where is the propagation direction, and is the polarization vector satisfying . Evidently, such an incident wave satisfies the homogeneous equation

It follows from the equations (1.1) and (1.2) that the scattered field satisfies

In addition, the scattered field is required to satisfy the following Silver-Müller radiation condition:

where r  =  |x|. Let Ω be a bounded domain of with boundary Γ, which contains the compact support of the scatterer q(x). Denote by ν the outward unit normal to Γ. Computationally, it is convenient to reduce the problem to a bounded domain (cell) by imposing a suitable (artificial) boundary condition on Γ. For simplicity, we employ the first-order absorbing boundary condition [10]

Given the incident field Eⁱ, the direct problem is to determine the scattered field E for the known scatterer q(x), which has been well studied [4]. This work is devoted to the numerical solution of the inverse medium scattering problem, i.e., determining the scatterer q(x) from the measurements of near-field current densities, ν × E|_Γ. Although this is a classical problem in inverse scattering theory, little is known on reconstruction methods, due to the nonlinearity, ill-posedness, and large-scale computation associated with the inverse scattering problem. We refer the reader to [1, 5, 8, 9, 14] for related results on the inverse medium problem. See [4] for an account of the recent progress on the general inverse scattering problems.

Our goal of this work is to present a recursive linearization method that solves the inverse medium scattering problem of Maxwell's equations in three dimensions. The reader is referred to [2, 3] for recursive linearization approaches for solving the inverse medium scattering problems in two dimensions. The algorithm requires multi-frequency scattering data, and the recursive linearization is obtained by a continuation method on the wavenumber k. The algorithm first solves a linear equation (Born approximation) at the lowest wavenumber k. Updates are made by using the data at higher wavenumbers k sequentially. Following the idea of the Kaczmarz method [5, 12], we use partial data and solve an underdetermined minimal norm solution at each sweep. For each iteration, one forward and one adjoint state of the Maxwell equations are solved which may be implemented by using the symmetric second-order edge elements [11, 13].

2. Low frequency modes of the scatterer

Rewrite (1.3) as

where the incident wave is taken as . Consider a test function , where satisfying . Hence F satisfies (1.2).

Multiplying equation (2.1) by F, and integrating over Ω on both sides, we have

Integrating by parts and noting (1.2) for F, we deduce

Using the boundary condition (1.4) of the scattered field E, and the special form of the incident wave Eⁱ and F, we then get

A simple calculation yields

When the wavenumber k is small, the first part of the right-hand side of the above integral equation is dominant. Dropping the nonlinear (second) term, we obtain the linear integral equation

which is the Born approximation.

Since the scatterer q(x) we use a compact support, we use the notation

where is the Fourier transform of q(x) with . Choose

where θ_i,_i are the latitudinal and longitudinal angles respectively. It is obvious that the domain [0,π] × [0,2π] of (θ_i,_i),i  =  1,2, corresponds to the ball . Thus, the Fourier modes of in the ball can be determined. The scattering data with higher wavenumber k must be used in order to recover more modes of the true scatterer.

In practice, the Kaczmarz method [5, 12] is used to implement the linear integral equation (2.3) in order to reduce the computational cost and instability.

3. Recursive linearization

As discussed in the previous section, when the wavenumber k is small, the Born approximation allows a reconstruction of those Fourier modes less than or equal to 2k for the function q(x). We now describe a procedure that recursively determines q_k, an approximation of q(x) at k  =  k_j for j  =  1,2,..., with the increasing wavenumber. Suppose now that the scatterer has been recovered at some , and that k>0 is slightly larger than . We wish to determine q_k or, equivalently, to determine the perturbation

For the reconstructed scatterer , we can solve at the frequency k the forward scattering problem

For the scatterer q_k, we have

Subtracting (3.1), (3.2) from (3.3), (3.4), and omitting the second-order smallness in δq and in , we obtain

For the scatterer q_k, and the incident wave Eⁱ, we define the map S(q_k,Eⁱ) by

where E is the solution of (3.3), (3.4). Let γ be the trace operator to the boundary Γ of the bounded domain Ω. Define the scattering map

For simplicity, denote M(q_k,Eⁱ) by M(q_k). By the definition of trace operator, we have

We next examine the boundary data ν × E(x;θ₁,₁;k). Here, the variable x is the observation point which has two degrees of freedom because it is on the artificial boundary Γ. The terms θ₁,₁ are latitudinal and longitudinal angles of the incident wave Eⁱ, respectively. At each frequency, we have four degrees of freedom, and thus again data redundancy, which may be addressed by fixing one of the incident angles, say θ₁.

Use the notation (₁)_j  =  (j - 1)*2π/m,j  =  1,...,m, and the residual operator

where is the solution of (3.1), (3.2) with the incident longitudinal angle (₁)_j, and the scatterer . For each j, consider the minimal norm solution of the following problem:

which has the form

In practice, some regularization [6] would also be needed.

Lemma 3.1. Given residual , there exists a function F_j such that the adjoint Fréchet operator satisfies

where the bar denotes the complex conjugate, E_jⁱ is the incident wave with the longitudinal angle (₁)_j, and is the solution of (3.1), (3.2) with the incident wave E_jⁱ.

From the above lemma, it follows easily that

Therefore, for each incident wave with longitudinal angle (₁)_j, it is necessary to solve one forward and one adjoint problem for the Maxwell equations. Observe that the adjoint problem in fact takes a similar variational form to the forward problem. Essentially, we need to compute two forward problems at each sweep. Once δq_j is determined, is updated by . After the mth sweep is completed, we get the reconstructed scatterer q_k at the wavenumber k.

Figure 1. (a) Integrals with different wavenumbers k for fixed incident angle θ  =  π/3,  =  π/3. Solid curve: the exact integral value of the left-hand side of (2.2);  + : the numerical integral of the first term of the right-hand side of (2.2); *: the numerical integral of the second term of the right-hand side of (2.2); °: the numerical integral of the right-hand side of (2.2). (b) The numerical integral with different θ for fixed wavenumber k  =  2.0 and  =  π/3. Solid curve: the exact integral value of the left-hand side of (2.2); °: the numerical integral of the right-hand side of (2.2) (c) The numerical integral with different for fixed wavenumber k  =  2.0 and θ  =  π/3. Solid curve: the exact integral of the left-hand side of (2.2); °: the numerical integral of the right-hand side of (2.2).

4. Numerical experiments

In this section, we discuss the treatment of the forward scattering problem, and computational aspects of the recursive linearization algorithm.

For the forward solver, we adopt the symmetric second-order tetrahedral edge elements [11]. The reverse Cuthill-McKee ordering [7], the compressed row storage format, and the quasi-minimal residual algorithm with diagonal preconditioning are used in the assembly of unknowns, the storage of the coefficient matrix, and solving the linear system, respectively.

Consider a test problem with the exact scatterer

The compact support of this scatterer is an ellipsoid contained in the unit ball. For simplicity, we take , and in the test of the forward solver. Numerical results are shown in figure 1. In figure 1(a), for fixed incident latitudinal angle θ  =  π/3, and the longitudinal angle  =  π/3, the direct problem is solved at different wavenumbers k. In figures 1(b), and (c), the numerical results are shown with different latitudinal angles , and with different longitudinal angles , respectively. As is easily seen from figure 1(a), the first term of the right-hand side of the integral equation (2.2) is indeed dominant compared with the nonlinear second term when the wavenumber k is small.

For a simple stability analysis, some relative random noise is added to the data, e.g., the tangential component of the electric field is updated to

where rand gives normally distributed random numbers in [ - 1,1] and σ is an error parameter. In our numerical experiments, the latitudinal angle of the incident wave θ₁  =  π/2 is fixed, σ  =  0.02, the sweep number m  =  20, and some appropriately chosen regularization parameters β_k  =  0.8/k, where k is the wavenumber. Define the relative error by

where is the reconstructed scatterer, and q is the true scatterer. Figure 2 shows the slices of true scatterer, and figure 3 gives the reconstruction at the wavenumber k  =  5. The relative errors are shown in table 1 at different wavenumbers k.

This research was supported in part by the NSF grant DMS 01-04001 and the ONR grant N000140210365.

References

[1]

Ammari H and Bao G 2001 Analysis of the scattering map of a linearized inverse medium problem for electromagnetic wave Inverse Problems 17 219-34 

[2]

Bao G and Liu J 2003 Numerical solution of inverse scattering problems with multi-experimental limited aperture data SIAM J. Sci. Comput. 25 1102-17 

[3]

Chen Y 1997 Inverse scattering via Heisenberg uncertainty principle Inverse Problems 13 253-82 

[4]

Colton D and Kress R 1998 Inverse Acoustic and Electromagnetic Scattering Theory (New York: Springer) 

[5]

Dorn O, Bertete-Aguirre H, Berryman J G and Papanicolaou G C 1999 A nonlinear inversion method for 3D electromagnetic imaging using adjoint fields Inverse Problems 15 1523-58 

[6]

Engl H, Hanke M and Neubauer A 1996 Regularization of Inverse Problems (Dordrecht: Kluwer) 

[7]

Gibbs N E, Poole W GJr and Stockmeyer P K 1976 An algorithm for reducing bandwidth and profile of a sparse matrix SIAM J. Numer. Anal. 13 236-50 

[8]

Hohage T 2001 On the numerical solution of a three-dimensional inverse medium scattering problem Inverse Problems 17 1743-63 

[9]

Haddar H and Monk P 2002 The linear sampling method for solving the electromagnetic inverse medium problem Inverse Problems 18 891-906 

[10]

Jin J 2002 The Finite Element Methods in Electromagnetics (New York: Wiley) 

[11]

Kameari A 1999 Symmetric second order edge elements for triangle and tetrahedra IEEE Trans. Magn. 35 1394-7 

[12]

Natterer F and Wübbeling F 1995 A propagation-backpropagation method for ultrasound tomography Inverse Problems 11 1225-32 

[13]

Nédélec J C 1980 Mixed finite elements in  Numer. Math. 35 315-41 

[14]

Vögeler M 2003 Reconstruction of the three-dimensional refractive index in electromagnetic scattering by using a propagation-backpropagation method Inverse Problems 19 739-53