Three-dimensional reconstruction of non-homogeneous dielectric objects by the coupled-dipole method

REGULAR ARTICLES

Three-dimensional reconstruction of non-homogeneous dielectric objects by the coupled-dipole method

Amin Bassrei^I; Thierry J. Lemaire^II

^IUniversidade Federal da Bahia, Instituto de Física & Centro de Pesquisa em Geofísica e Geologia Campus Universitário de Ondina 40170-290 Salvador BA, Brazil

^IIUniversidade Estadual de Feira de Santana Departamento de Física Av. Universitária, s/n, Campus Universitário, km 03, BR 116 44031-460 Feira de Santana BA, Brazil

ABSTRACT

In this paper we present an inversion procedure for electromagnetic scattering, based on the coupled-dipole method (CDM) combined with an inversion algorithm making use of the singular value decomposition procedure associated with a regularization factor. This method permits to obtain images of non-homogeneous dielectric objects whose dimensions are comparable to the incident wavelength. The feasibility of this method is showed in two synthetic examples using the CDM with 257 and 515 dipoles, for spherical objects with a spherical inclusion. The method also works if the scattered electric field, which is the input data, is corrupted with Gaussian noise.

Keywords: Electromagnetic scattering; Coupled-dipole method; Inverse problems; Singular value decomposition; Regularization

I. INTRODUCTION

Electromagnetic scattering phenomena have been shown to be of great importance because of their various applications in many areas [1-4]. The ability of electromagnetic waves to penetrate into objects without making alterations of their structure is an important property and permits to obtain information about the studied scattering object, i.e., it is possible to make a non-destructive imaging analysis.

Several methods were developed to solve inverse electromagnetic problems. Most of the methods dealing with microwave imaging are based on the discretization of an integral equation for the electric field or the current density. This leads to a linear system of equations which is ill-posed and consequently the solution is not unique.

Recently, a new numerical technique [5] was introduced in order to obtain the complex refractive index of a homogeneous scattering object in the scope of active microwave imaging. The modeling of the scattering phenomenon is based on the well-known coupled-dipole method (CDM) [6,7] and simulations showed the feasibility of this technique even with corrupted data. In this paper are presented new results for non-homogeneous scattering objects with artificial data (noise free and also corrupted by noise), using the truncated singular value decomposition (SVD) method to make the data inversion. The method applied to non-homogeneous spheres shows a reasonable resolution to detect the inclusion.

This non invasive and non destructive technique is very useful to determine the dielectric properties of an object [4]. The input physical quantity of interest is the electric field and the searched characteristics of the object are the permittivity (generally a function of space coordinates) and geometry. Our approach is based on the inversion of the CDM which permits to describe the scattered electric field. This inversion technique [5] using the CDM is shown to be efficient when applied to a homogeneous and isotropic scatterer of size parameter less than 10.0.

This inversion scheme is tested with synthetic data for simple geometries of the scatterer, and noisy data are used to illustrate its stability. The a priori information used in this model is the external shape of the studied object. This paper is organized as follows. In Section II we briefly describe the CDM. In Section III, the inversion procedure is introduced and numerical results are given in Section IV. Section V is the conclusion.

II. COUPLED-DIPOLE METHOD

The theoretical description of the scattering phenomenon is done with the well known CDM. The polarizability of each induced dipole located at the sites of a cubic lattice is given by the Clausius-Mossotti (or Lorentz-Lorenz) prescription [1]. We remind that the basic equation to be inverted is the relation between the measured scattered electric field E_scat(r_i) at r_i, and the N electric dipolar moments p_j,(j = 1,...,N) of the dipolar units, located at r_j,

where the 3 ×3 complex matrix P(r) is given by [8] (for a time factor exp(-iwt)),

"Ä" is the Kronecker product, and I is the unit 3×3 matrix. The data to be inverted are the scattered electric fields, and the inversion output is the dipolar moment of each dipolar unit. This last quantity is used to obtain the polarizability of each dipolar unit and, consequently, the refractive index of each site of the cubic lattice. However, the Freedholm equation of first kind (1) is ill-conditioned making the inversion process not trivial.

III. DATA INVERSION

The inversion procedure is a 3-step algorithm. The first one leads to the inversion of a system of linear equations (1) obtained by writing the scattered field at each point where it is measured. The second step leads to the estimation of the polarizability at each site of the lattice by use of the equation [5],

where E₀ is the incident field. The last step leads to the determination of the complex refractive index at each site from the Clausius-Mossotti formula. Let us notice that the main computational time is spent at the first step of the algorithm because of the large size of the matrix to be inverted.

As mentioned above, we will perform the inverse procedure by the SVD technique. Our basic equation here is

which is a linear transformation on p. If the vector E_scat describes the observed actual output of the system, the problem is to "choose" the vector of model parameters p in order to minimize, in some sense, the difference between the observed E_scat and the prescribed output of the system Pp. In order to solve the system of linear equations, we have first to invert the matrix P, obtaining the pseudo-inverse P⁺, then multiply the pseudo-inverse by the scattered field, obtaining finally the recovered vector of moments p^est which permits the computation of the refractive index.

The generalized inverse matrix is frequently used in inverse problems, and its solution has the minimum Euclidean norm. In this case the objective function to be minimized is

where t is the vector of Lagrange multipliers. The minimization yields

The above equation is equivalent to the so-called "pseudo-inverse" for underdetermined systems developed by [9] and later by [10]. One versatile way to calculate the pseudo-inverse is through SVD [11] where the square kernel matrix is expressed as

where U is the matrix which contains the orthonormalized eigenvectors of PP^T, V the orthonormalized eigenvectors of P^TP and the matrix S is formed by the singular values of P. The pseudo-inverse P⁺ will be given by

where S^-1 = (1/s_ii), i = 1; ...;k; and 1/s_ll = 0, for l > k.

IV. NUMERICAL SIMULATIONS

In order to show the feasibility of this method, we present two synthetic examples where a non-homogeneous sphere with a spherical inclusion is described by 257 and by 515 dipoles. In both examples the refractive index of the homogeneous background is 1.1 and the index of the inclusion is 1.4. The choice of this simple shape is not a limitation and the extension to other scatterer geometry is straightforward. In order to obtain the synthetic data to be inverted, which is the scattered electric field, we calculate the dipolar moments using the CDM approach. Each dipole is associated with a vector of six components: three components of the real dipolar moment, and three components of the imaginary one. Thus, in the inversion procedure, we have 6 × 257 = 1542 model parameters for the first example and 6 × 515 = 3090 for the second one.

The acquisition geometry is performed in two planes around the scatterer. Using standard notation for spherical coordinates, for one plane we have f = 0 and for the other f = p/2. The angle q varies from 0 to p. For the first example one plane has 128 observation points and the other has 129, while for the second example one plane has 257 observation points and the other has 258. This means that for the first situation there are 257 observation points and for the second one there are 515 observation points.

Since the scattered electric field used as the input for the inverse procedure is also a complex vector, we have 6 × 257 = 1542 data parameters for the first example and 6 × 515 = 3090 for the second one. The number of data parameters is then equal to the number of model parameters, making the system determined although ill-posed. The matrix P is then square, but this is not a limitation, that is, we could have more data parameters than model parameters, making the problem overdetermined, or the opposite situation, with more model parameters than data parameters, making the problem in this case underdetermined.

As explained before for the data inversion we used the SVD technique. The sphere has a size parameter equal to 2.0, with a spherical inclusion of size parameter equal to 1.2, centered at (x₁, y₁, z₁) with (k x₁ = 0.4, k y₁ = 0.4, k z₁ = 0.4). The singular values can be seen in Fig. 1 for 257 dipoles example. The behavior of singular value decay was quite similar for the simulations with 515 dipoles.

The ill-posedness of the problem suggests us the need of some kind of regularization. We have employed the damped least squares approach, as proposed by [12] and [13]. Here the matrix P^TP is damped with a positive constant, in order to stabilize the inversion. Other kinds of regularization can be applied like filtering the matrix P^TP with a second or fourth order derivative matrix as proposed by [14-18]. Let us notice that the choice of the regularization factor l is a problem by itself. Several techniques have been suggested in the literature. Since the purpose of this work is not to explore this search for the optimum regularization factor, we have adopted the following procedure: for each example, we used several noise levels, and for each noise level we performed 35 inversions, each one with a different l, computing the root mean square (rms) error for the data parameters and for the model parameters. The range for l was rather wide, from l_min = 1 ×10^-15 to l_max = 1 ×10⁺¹⁹.

The results of the simulations are summarized in Tables I and ^II. Table I is for 257 dipoles and ^{Table II} for 515 dipoles. In each table the first column gives the value of noise factor b which is added to the scattered electric field (), so that

where ran_j is a sequence of random numbers within the range [0, 1].

The second column gives the noise estimate e_rms, which is the relative rms error between E_scat and :

As mentioned above for each noise level we run 35 inversions, where each simulation had a different regularization factor l, used to compute the estimated moment vector:

The best l, that is, the l which provides the minimum rms error between the vector of true moment and the vector of estimated moment, is given in the third column. This relative rms error is denoted by (see the fourth column in Tables I and ^II) and expressed by

Note that the value for the best l is kept constant with different noise levels, not showing fluctuation.

We can compute once more the scattered field associated with the estimated moment vectors by applying the forward modeling procedure,

In the tables, the fifth column gives the relative rms error between the (observed) scattered electric field and the calculated scattered electric field, denoted by :

The inversion output, p^est, is used to calculate the polarizability a^est, which is then used to calculate the estimated refraction index n^est, with the Clausius-Mossotti formula. Finally the last column in the tables shows the relative rms error between the true refraction index and the estimated one,

Due to space limitations we show only some results from the two examples. In each figure we show, for a given slice or projection, the 2-D image with gray scale as well as the contour map. Just the noise free results are showed in the figures. A compilation of the results with all noise levels tested can be seen in the above tables. Also, only two projections are showed: kz = 0.0, which is at the equator, and kz = 0.75, which is located at 3/4 from the equator.

In the first example the sphere is formed by 257 dipoles. This means that at the equator the model has a diameter formed by 9 dipoles. The inclusion, with refraction index equal to 1.4, has a radius equal to 60 % of the "homogeneous" sphere in such a way that its maximum diameter at the equator is formed by 5 dipoles. The features above explained can be seen in Fig. 2, which we call the true model. Note the maximum diameter, at the equator, is a. This value is not specified, but is related to the wavenumber k by the relation k a = x. In our examples x = 2.0. The acquisition is done at circles of radius equal to k R = 4.0. The noise free inversion for k z = 0.0 is shown in Fig. 3. The next two figures repeat the set, now for k z = 0.75. Fig. 4 shows the true model, Fig. 5 presents the reconstruction. In both cases we can see a good agreement in relation to the location and intensity but a poor shape definition.

In the second example the sphere is formed by 515 dipoles, which implies that at the equator the model has a diameter formed by 11 dipoles. Again, the inclusion with refraction index equal to 1.4 is larger than the radius, with its maximum diameter at the equator formed by 7 dipoles. The true model for k z = 0.0 is presented in Fig. 6, and the noise free inversion in Fig. 7. For the k z = 0.75 projection, Fig. 8 shows the true model and Fig. 9 the noise free inversion. Again, in both cases we can see a good agreement in relation to the location and intensity but a poor shape definition.

In a previous work [5] each inversion with 171 dipoles demanded around 80 min of CPU time in a RISC workstation. Now using a Pentium IV with a 2.4 GHz clock and 2 Gb of RAM memory each inversion demanded around 5 min for 257 dipoles and around 100 min for 515 dipoles.

The reconstruction is indeed three-dimensional. In both examples there are only two acquisition planes, in such a way that the data acquisition could essentially be considered two-dimensional. The reconstruction could be better, that is, with smaller , with the inclusion of more acquisition planes. In this case the data acquisition would be properly considered three-dimensional. Although the matrix P is square there is a shortage of data. This can be clearly seen in Fig. 1, where the singular values are rapidly attenuated. In all the images displaying the reconstructed model the location of the scatterer and the value of the refraction index are well given, but the shape of the scatterer is not well defined.

V. CONCLUSIONS

In this paper, we show the feasibility of the inversion method for non-homogeneous scattering objects. The inversion can be done in practice for objects having a real part of the complex refractive index close to 1 and with weak absorption. From the simulations with synthetic examples with ill-conditioned kernel matrices and data corrupted by noise, we showed that the proposed algorithm is feasible for the inversion of electromagnetic data. In general the estimated models indicate the presence of a inclusion inside the scatterer and give a reasonable estimate of its refractive index, although the shape is not well defined.

Received on 9 February, 2006; revised version received on 5 November, 2006

• [1] H. C. Van de Hulst, Light scattering by small particles (Dover, New York, 1981).
• [2] B. T. Draine, Astrophysical Journal 333, 848 (1988).
• [3] M. Baribaud, Journal of Physics D 23, 269 (1990).
• [4] S. Caorsi, G. L. Gragnani, and M. Pasorino, IEEE Transactions on Antennas and Propagation 42, 581 (1994).
• [5] T. J. Lemaire and A. Bassrei, Applied Optics 39, 1272 (2000).
• [6] E. M. Purcell and C. R. Pennypacker, Astrophysical Journal 186, 705 (1973).
• [7] B. T. Draine and P. J. Flatau, Journal of the Optical Society of America A 11, 1491 (1994).
• [8] P. Chiappetta, Journal of Physics A 13, 2101 (1980).
• [9] E. H. Moore, Bulletin of the American Mathematical Society 26, 394 (1920).
• [10] R. Penrose, Proceedings of the Cambridge Philosophical Society 51, 406 (1955).
• [11] C. Lanczos, Linear Differential Operators (Van Nostrand, London, 1961).
• [12] K. Levenberg, Quartelly of Applied Mathematics 2, 164 (1944).
• [13] D. W. Marquardt, Journal of the Society for Applied and Industrial Mathematics 11, 431 (1963).
• [14] D. L. Phillips, Journal of the Association Computers Manufacturers 9, 84 (1962).
• [15] B. C. Cook, Nuclear Instrumentation and Methods 24, 256 (1963).
• [16] A. N. Tikhonov, Soviet Mathematics Doklady 4, 1035 (1963).
• [17] A. N. Tikhonov, Soviet Mathematics Doklady 4, 1624 (1963).
• [18] S. Twomey, Journal of the Association of Computers Manufacturers 10, 97 (1963).

Publication Dates

• Publication in this collection
  17 July 2007
• Date of issue
  June 2007

History

• Received
  09 Feb 2006
• Accepted
  11 May 2006