International Journal of Antennas and Propagation
Volume 2018 (2018), Article ID 3020417, 15 pages
https://doi.org/10.1155/2018/3020417

Research Article

Embedding Approach to Modeling Electromagnetic Fields in a Complex Two-Dimensional Environment

Anton Tijhuis,¹ Ann Franchois,² and Jean-Michel Geffrin³

¹Department of Electrical Engineering, Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, Netherlands
²Department of Information Technology, Ghent University, iGent-Technologiepark-Zwijnaarde 15, 9052 Gent, Belgium
³Aix Marseille University, CNRS, Centrale Marseille, Institut Fresnel, Marseille, France

Correspondence should be addressed to Ann Franchois; ann.franchois@ugent.be

Received 31 October 2017; Revised 19 February 2018; Accepted 6 March 2018; Published 10 May 2018

Academic Editor: Paolo Burghignoli

Copyright © 2018 Anton Tijhuis et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

An approach is presented to combine the response of a two-dimensionally inhomogeneous dielectric object in a homogeneous environment with that of an empty inhomogeneous environment. This allows an efficient computation of the scattering behavior of the dielectric cylinder with the aid of the CGFFT method and a dedicated extrapolation procedure. Since a circular observation contour is adopted, an angular spectral representation can be employed for the embedding. Implementation details are discussed for the case of a closed 434 MHz microwave scanner, and the accuracy and efficiency of all steps in the numerical procedure are investigated. Guidelines are proposed for choosing computational parameters such as truncation limits and tolerances. We show that the embedding approach does not increase the CPU time with respect to the forward problem solution in a homogeneous environment, if only the fields on the observation contour are computed, and that it leads to a relatively small increase when the fields on the mesh are computed as well.

1. Introduction

In almost any computational approach to solving nonlinear inverse-scattering problems, a discretized configuration is introduced that depends on a fixed number of parameters. Subsequently, a cost functional is defined in terms of simulated and known scattered fields. Here, two different strategies can be distinguished. Conventionally, the corresponding forward problem is treated as an auxiliary problem, which is solved exactly for successive approximate configurations [1–6]. For multidimensional problems, this requires a number of field computations for a varying physical parameter such as frequency or source position. The cost function then refers to the known measured field information and preferably includes a regularizing function of the configuration parameters [7–9]. In the so-called modified gradient method and subsequent generalizations, the configuration and the unknown fields are determined simultaneously [10, 11]. The conventional approach has the advantage that the formulation of the inverse problem directly relates the parameterized configuration to the known field data. From a practical point of view, however, it is sometimes considered as less feasible because of the computational effort required in the repeated field computations. The argument is that it is not needed to compute the field with full accuracy in a configuration that still deviates considerably from the actual one.

For the case of an inhomogeneous, lossy dielectric cylinder in a homogeneous surrounding medium, however, it was demonstrated that a highly efficient implementation is obtained when the fields are computed by solving a contrast-source integral equation with a combination of the conjugate-gradient FFT (CGFFT) method and a special extrapolation procedure [12]. The extrapolation can be performed for almost any physical parameter [13], such as frequency or source position. Thus, the forward scattering problem can be solved for each new value of the physical parameter in a few iterations of the CGFFT procedure. This technique has been demonstrated successfully in the context of Newton-type inverse scattering [5, 6]. It is the authors’ experience that these schemes, with only the parameters in the profile parameterization as fundamental unknowns, are generally more efficient than schemes where the field and the profile are determined simultaneously.

A special feature of our implementation is that its efficiency is based on the circumstance that the dielectric cylinder is embedded in a homogeneous surrounding medium. This means that Green’s function in the integral equation above exhibits convolution symmetry. Preserving that symmetry in the relevant space discretization allows the application of FFT operations in evaluating the operator products in the conjugate gradient method [14–16]. In practical experiments, however, the surrounding medium may be inhomogeneous and the symmetry is broken. In that case, FFT operations are no longer applicable. The same problem also arises in the modified gradient method, where FFT operations are used as well to compute the field updates.

To circumvent this problem, we use the feature that the scattering operator characterizes the complete electromagnetic response of the region inside a closed observation contour. Hence, it must be possible to determine the scattered-field data from a cylinder in an arbitrary environment from the scattering operator for the same object in a homogeneous environment. That data, in turn, can then be obtained with the existing implementation. We introduced this so-called embedding approach for a cylinder inside a circular observation contour in [17, 18]. The choice of this particular configuration was inspired by the experimental research with a 434 MHz scanner of the third author [19]. Several authors have shown interest in quantitative imaging with a circular scanner with metallic enclosure, including [20–22]. It is well-known that employing the 2D Green’s function of the empty casing is computationally expensive [19, 20].

In the present paper, we formulate the embedding approach in the angular spectral representation for a general surrounding 2D medium and subsequently specialize to the case where the dielectric cylinder is surrounded by a perfectly conducting circular container. Besides a more comprehensive theoretical formulation than in [17, 18], we provide details on the numerical implementation and on its performance in accuracy and speed as a function of various parameters. We show that, with well-chosen values for these parameters, the embedding approach does not increase the CPU time as compared to the forward problem solution in a homogeneous environment, if only the fields on the observation contour are computed, and that it leads to a relatively small increase in CPU time, when the fields in the object are needed as well, for example, to compute the Jacobian matrix in a Newton-type inversion scheme.

The embedding approach relies on the identification of the scattering and reflection operators for the dielectric cylinder and the empty microwave scanner, respectively. The idea of using such an operator to characterize scattering properties has a long tradition in the electromagnetic literature [23–26]. The use of a numerically computed scattering operator, however, is new and originates from the availability of the “march in source position” method [13]. A generalization for multiple interacting domains of arbitrary shapes is given in [27, 28].

Finally, it should be remarked that in [29], a procedure is proposed based on reciprocity that is capable of converting the field in the complete configuration into the field in a homogeneous environment, that is, the reverse procedure from what is proposed in the present paper. The suggestion is to perform the profile inversion on the thus corrected data, using inverse-profiling algorithms for objects in a homogeneous background. However, this idea has two possible drawbacks. First, in order to carry out the conversion from one environment to another, complete data on a contour surrounding the scatterer must be available. In an actual experiment, such data may not always be available, while theoretical results for an estimated configuration can always be computed. Second, the conversion renormalizes the experimental data including the measurements, while the present procedure allows a comparison with the actual data. This makes it easier to account for the accuracy of these data, for example, by including appropriate weighting coefficients in the cost functional.

The paper is organized as follows. In Section 2, we describe the scanner configuration and its mathematical idealization. Section 3 summarizes the field computation for an object in a homogeneous environment. The scattering operators are introduced in Section 4 and used to formulate the embedding approach in Section 5. Section 6 presents the computational details for a homogeneous environment, the empty scanner and the complete configuration. In Section 7, the computational complexity of the algorithm is analyzed, and a procedure is given for tuning the computational parameters. The conclusions are drawn in Section 8.

2. Formulation of the Problem

In the present paper, we describe and investigate an efficient procedure to calculate the electric field inside a cylindrical scanner for a given permittivity profile. This field may then be used in inverse-profiling algorithms. For our numerical experiments, we adopted the configuration of the scanner described in [19] and shown in Figure 1, which was developed to conduct biomedical imaging experiments. This scanner comprises a circular array of 64 transmitting/receiving conical dipole antennas, which operate at 434 MHz in a multi-incidence mode, that is, one antenna at a time is transmitting and the others are receiving. The array has a radius of 27.6 cm and is placed inside a water-filled metal casing with a slightly larger radius of 29.0 cm. Measurements of the relative permittivity of the water typically yielded , which corresponds to a wavelength  cm in the water. The diameter of the T/R circle thus is about , the antennas are spaced apart about , and they are at a distance of about from the casing.

Figure 1: Synthetic dielectric object in the scanner.

Our modeling assumptions are that we may assume the fields to be two-dimensional and that mutual coupling between the antenna elements can be neglected. We therefore consider the configuration shown in Figure 2. An inhomogeneous, lossy dielectric cylinder in an observation domain is excited by a time-harmonic electric line source on a circular contour with radius outside the cylinder. “Complete” scattering data need to be determined, that is, the scattered electric field must be calculated on a circle for a line source anywhere on that contour. The “environment” in is two-dimensional and linearly reacting, that is, the reflected field due to sources in satisfies the superposition principle. The observation contour is located in a homogeneous region with relative permittivity . The radii and bound the cylinder and the environment from the outside and inside, respectively. For now, we leave the environment unspecified, but in the calculations, it will be either a homogeneous space or a conducting cylinder.

Figure 2: Two-dimensional model of object in scanner.

The electrically polarized field caused by a line source at source point can be identified as Green’s function, that is, the solution of the second-order differential equation

for that satisfies the proper conditions in . In (1), is a two-dimensional position vector and is a complex frequency with . In the present paper, we describe an approach for solving (1) for the case where the scattering configuration is located inside a general environment. This solution will be denoted as . In our approach, we first determine Green’s function for the same configuration in a homogeneous environment, which will be denoted as . The differential equation (1) applies to both problems. For the object in a homogeneous environment, satisfies the radiation condition as . For the object in a general environment, the boundary condition depends on the properties of the exterior medium. For the special case where the environment is a metal wall with inner radius , we have at that wall. For the general case, a reflection operator will be introduced further on that can be employed to formulate a boundary condition on the observation contour .

3. Homogeneous Environment

The feasibility of our approach depends strongly on the availability of a fast procedure for determining the field for the scattering configuration in a homogeneous environment. Such a scheme is available from [12]. We briefly repeat its main features. Starting point is the contrast-source integral equation

where , is a finite domain in which , and is Green’s function of the surrounding homogeneous dielectric medium:

In (3), denotes the modified Bessel function of the second kind of order zero. In this integral equation, the incident field and Green’s function are available in closed form.

The square region and , in which is enclosed, is subdivided into subregions with mesh size . The grid points of the square mesh are located at , with for and for . Solving (2) now amounts to determining an approximation of at the grid points .

The space discretization of the integral in the right-hand side of (2) has two special aspects. First, the logarithmically singular behavior of as is subtracted by breaking up the integral over into

Second, the discretization of the integrals in (4) is based on approximating suitable parts of both integrands by piecewise-linear approximations and integrating analytically over polygons determined by the boundary of and the grid. This results in a discretized integral equation of the form

where is a sampled, filtered version of . In (5), the convolution-type structure of the continuous equation (2) has been preserved. This makes this equation suitable for the application of the CGFFT method. In addition, it is second-order accurate in the mesh size . More information on the discretization and the corresponding error estimate can be found in [12].

The initial estimate for the CGFFT procedure is obtained by taking a linear combination of previous “final” solutions and determining the coefficients by minimizing the squared error for the problem at hand. In [13], a more detailed explanation is given, as well as several examples of physical parameters for which the effectiveness of this extrapolation has been demonstrated. In the present context, we extrapolate in source position [6]. Since , the physical parameter that is varied here is the angle , which explains the denomination “marching on in angle.”

4. Scattering Operators

Before we consider the complete problem, we first introduce the scattering operators for the individual building blocks of the configuration. Since both the observation contour and the casing have circular symmetry, it is convenient to carry out a Fourier transformation with respect to the angle and to carry out the analysis in the spectral domain. This leads to a matrix formulation in terms of the angular spectral coefficients.

4.1. Scattering by a Dielectric Cylinder

We start by considering the dielectric cylinder. For a general excitation outside , we can write the field in in spectral form as

The coefficients represent the incident field and the coefficients the field scattered by the dielectric object. Since the cylinder is linearly reacting, these coefficients must be related by a linear scattering operator:

The value of the elements in (7) can be obtained from Green’s function for the cylinder in a homogeneous environment. In the special case described in Section 3, the incident field is the field caused by a line source in a completely homogeneous background and can be expressed in spectral form as

where and . This expression is valid for all , that is, for .

Now, for , the incident fields specified in (6) and (8) must be identical. Since the expansion functions are linearly independent, comparing both expressions directly leads to the identification

Using the definition (7) then leads to the following expression for the scattered field for a receiver position for a cylinder in a homogeneous environment

This means that the individual elements that define the scattering operator can be obtained from by applying Fourier transformations with respect to and . The modified Bessel function of the second kind has no zeros in the right half-plane . Therefore, no problems with division by such functions will be encountered.

4.2. Reflection by a Cylindrical Casing

Next, we consider the environment. Writing the field in in the spectral form (6), where the coefficients correspond to a field radiated by sources in and where the coefficients represent a source-free field reflected by the environment, leads to the introduction of the reflection operator:

This operator can be determined from the field that is excited by a line source in a configuration where the actual environment in surrounds an “empty” observation domain that contains a homogeneous dielectric medium with .

In the computations, we will consider the special case where the environment is a perfectly conducting circular container with inner radius . In that case, it follows from the boundary condition that is a diagonal operator:

where is the Kronecker symbol. The function can have zeros on the imaginary -axis. Such a zero corresponds to the situation where the empty casing has an interior mode. However, interior modes occur only when the medium inside the casing is lossless. In the experimental setup, the water in which the dielectric object is immersed is lossy. In the present study, we have therefore not considered special measures to handle interior modes.

5. The Complete Configuration

Now that the individual building blocks of the configuration have been analyzed, we can combine the results to obtain the field in the complete configuration.

5.1. Field on the Observation Contour

We first determine the field on the observation contour . To that end, we divide the domain , in which the spectral analysis is carried out, in two subdomains that are separated by .

In the complete region , we have the spectral representation

where the coefficients and depend on the angle .

For the “interior subdomain” , we may envisage the situation as scattering by the dielectric cylinder. We treat the field due to the line source and the field reflected from the casing (given by the coefficients ) as primary fields that generate a secondary field (given by the coefficients ) that propagates in the direction of increasing . With the aid of the definition (7), we then directly arrive at

For the “exterior subdomain” , we similarly treat the line source and the field scattered by the dielectric cylinder as primary fields that generate a secondary field reflected from the metal casing. This leads to

Combining both results finally gives a linear equation for the coefficients :

This equation can be solved by truncating the summation over and and inverting the resulting matrix equation. This must be carried out for varying , that is, for multiple right-hand sides. Once the coefficients are found, we can use expression (14) to determine . Substitution of the values for these coefficients in (13) then gives the total field in .

5.2. Operator Formulation

The equation found in (16) can also be written in operator form. From a spectral point of view, the line source generates a source-free incident field in and a radiating field in , with spectral amplitudes

respectively. Identifying these amplitudes as the elements of two excitation vectors and then results in the operator equation

where the unknown vector contains the spectral amplitudes for the secondary field that originates from the casing. From energy considerations, it follows that the norm of the operator product is sufficiently small that the solution of (18) may formally be written as a geometrical series:

This confirms that repeated scattering and diffraction effects are accounted for in solving the system of equations given by (16) and (18). Moreover, (19) constitutes the basis for a Neumann-type iterative scheme that can be used for the solution of these equations. In the same notation, we have from (14)

which completes the formulation in operator form.

5.3. Field in the Observation Domain

The previous analysis allows us to compute the field for and on . However, in inverse-scattering algorithms, we must also determine for and for [5, 18]. This field is needed to determine the “profile update” in Newton-type optimization.

To this end, we use the equivalence principle. In , the total field may be envisioned as a response to the current at the line source and the induced surface current on the casing. Both these currents radiate a field that is incident on the dielectric cylinder. Alternatively, we may treat the second constituent as originating from an equivalent surface current on . This means that we replace the right-hand side of (1) by

and compute the field generated by this source for a dielectric cylinder in a homogeneous environment. In , we will then find the correct field . This conclusion is in fact merely a special formulation of Huygens’ principle. Obviously, the field in domain of this equivalent configuration will deviate from the actual field.

Determining the incident field by separation of variables and comparing the result for with the terms involving in (13) lead to the definition

Now, the uniqueness of the solution of the homogeneous wave equation implies that, for a given incident field in , the corresponding total field in is fixed. Therefore, we may identify the total field in as a superposition of fields generated by line sources on in the homogeneous embedding

which of course holds only for . In (23), is a point on , characterized by the angle . With (23), we can now compute directly from with .

For a general environment, an analogous procedure is available for the field in . We then express the field in the environment in terms of the line-source response of the empty casing. For a circular metal casing, such a procedure is not needed, since the spectral representation (13) remains valid up to .

5.4. Direct Reflection from the Casing

In the practical configuration shown in Figure 1, the line source is located close to the interface at . This means that the direct response, corresponding to the term in (19), must compensate the logarithmically singular behavior of the incident field as . In fact, for and , the field directly reflected from the casing may be approximated by that of an image source with opposite sign at and . Therefore, the convergence of the angular series slows down as . Now, it is observed from (12) that the computation of the reflection coefficients is much easier than the computation of the elements of the scattering matrix . Therefore, it makes sense to extract the direct reflection from the casing, which gives rise to the term , out of the total casing response determined by the coefficients . To this end, we rewrite the operator equation as

which leads to the power-series solution

The solution (25) separates the field that originates from the casing into a field that would be present in an empty casing and an additional field scattered by the dielectric. Both the field excited by the line source and the field that results from a direct reflection by the casing are considered as fields that are incident on the dielectric cylinder. This combined incident field is scattered at least once by the dielectric cylinder and reflected at least once at the casing before it contributes to the regularized spectral coefficient . Since the distance between this cylinder and the observation contour is considerably larger than the distance between the casing and the observation contour, the resulting diffracted field is much smoother than the fields that originate from the line source and from direct reflection at the casing.

To facilitate the discussion of the numerical aspects in the upcoming sections, we express the operator form of the decomposition (25) of the field in the region in the spectral form used in Section 5.1. Further, we restrict ourselves to the diagonal operator specified in (12). Combining (17) and (24) leads to the form

The second sum in (26) represents the field that is directly reflected from the casing, with given by (12) and given by (17). The field can also be written as the sum

where is the field in the empty casing—the first two sums in (26)—and is the difference field due to the insertion of the object—the last two sums in (26). In an experiment, is obtained by subtracting the measured fields with and without the object in place.

In a similar fashion, the equivalent surface current in (21) is decomposed as

where

corresponds to the direct reflection from the casing and where is obtained by replacing by in (22). generates the source-free constituent of the difference field, that is, the third sum in (26).

6. Numerical Study

We discuss the numerical implementation of the embedding approach, and we examine its performance in accuracy and speed as a function of various parameters for the configuration of the 434 MHz scanner [19]. All simulations were performed on a SUN ULTRA HPC 4000 workstation. In our code, we have used the public domain software packages LAPACK [30], for the linear system solutions; AMOS [31], for the Bessel function computations; and NMS [32], for the 2D FFTs. In particular, the accuracy of the various building blocks in the embedding approach is evaluated with tests on a number of homogeneous circular dielectric cylinders, specified in Table 1, for which accurate field solutions are available from spectral representations. This accuracy is quantified by means of the normalized root-mean-square error (NRMSE), defined as

Table 1: Diameter, relative permittivity and mesh size of the test objects: homogeneous dielectric circular cylinders Muscle1, Muscle4, Air1 and Air2, and an inhomogeneous dielectric cylinder Leg2. The permittivity for muscle is from [19], that for bone and fat is based on [33].

(a)

Cylinder
Muscle1
Muscle4
Air1
Air2

Diameter (cm)
8.8
35.2
8.8
17.6

54.2 + i38.4
54.2 + i38.4
1
1

2ℓ (cm)
9
36
9
18

(b)

Cylinder
Leg2: 3 layers

Fat
Muscle
Bone

Diameter (cm)
17.6
15.6
4.4

6 + i2.5
54.2 + i38.4
6 + i2.5

2ℓ (cm)

18

Besides the homogeneous cylinders of muscle and air, with approximate diameters and , there is also an inhomogeneous leg object with approximate diameter , which consists of a circular cylinder of muscle, covered with a 1 cm thick layer of fat and containing a bone with approximate diameter . The bone is decentered 1.5 cm in the -direction. The homogeneous cylinders are centered at the origin (i.e., the center of the casing), while the leg is decentered by −0.5 cm in the -direction. On the one hand, the discretization cell size and hence the number of grid points determine the accuracy of on the grid, of on and of all quantities derived from these in the embedding approach. On the other hand, a sufficient number of line sources—assumed to be equally spaced on —are needed to obtain accurate representations for the scattering matrix (10), say , and for the equivalent surface current on (21), say . The number of forward problem solutions may differ from the actual number of sources used in the cost function for an imaging experiment. In this section, we illustrate the influence of and on the computational performance.

6.1. Homogeneous Environment

6.1.1. Discretization Error—Choice of

We compared the total field in the grid points, computed by solving (5) with the CGFFT method for the source position , with the exact solution as a function of the discretization cell size . For the nonsymmetrical Leg2 object, the exact solution was chosen as the discretized solution for a much smaller cell size. We also compared the corresponding scattered field in the 64 receiver points with the exact solution. Figure 3(a) shows the discretization error on the total field, with ranging from 0.0352 cm to 2.25 cm . Depending on the diameter of the object, we used a mesh with side 9 cm, 18 cm, or 36 cm, and we varied between 4 and 512 (or between 25 and 263,169). Figure 3(b) shows the resulting discretization error on the scattered field. It is clearly illustrated that the discretization errors are of as , and it furthermore appears that for fixed , the variation in the errors for the different cylinders considered is limited. Let us add here that for highly contrasting objects, such as air cylinders, extremely small cells may be needed in areas where the field rapidly varies—this occurs, for example, when the air interface is in the neighborhood of the source. For objects with characteristics within the ranges considered in Table 1, Figure 3(b) is helpful for choosing the cell size corresponding to a given acceptable error on the scattered field, and Figure 3(a) then gives an indication for the CGFFT stop criterion. Proceeding with the iterations when the CGFFT NRMSE (with respect to the incident field on the mesh) is smaller than a certain fraction of the total field discretization error is a waste of computational effort.

(a)

(b)

Figure 3: The NRMSE for (a) on the mesh and for (b) on as a function of the cell size for the objects of Table 1. In presence of the casing and for proper choices of and , the values in (a) are multiplied with a factor between 1 and 2 for the NRMSE of , and (b) remains valid for the NRMSE of .

6.1.2. CPU Time

The CPU times for the different steps in the forward problem are illustrated in Figure 4 for Leg2 as a function of . Most of the time is taken up by the CGFFT iterations to solve (5) for the different angles of incidence; only a few percent of this time is needed to compute, respectively, the initial estimates for the total field by means of the “marching on in angle” procedure, the Green functions in (5), and the scattered fields in the receivers. The second most time-consuming step is the computation of the incident fields on the grid; the values presented here correspond to a worst-case scenario, since we did not exploit the symmetries in the grid and T/R points. Note that the Green functions can be stored for a given configuration of the mesh, transmitters/receivers, and exterior medium permittivity.

Figure 4: CPU times for the different steps in the forward problem for Leg2 with sources: (1) solution of equation (5) for with CGFFT and “marching on in angle”; (2) computation of the incident field ; (3) computation of the scattered field ; (4) one CGFFT iteration; (5) the summation (42) to obtain the total field . The curve (6) gives the CPU time for step (1) with the incident field as initial guess.

6.1.3. Marching on in Angle

We tested the efficiency of “marching on in angle,” where we used the total field solutions from the three previous excitations to compute the initial estimate, by comparing the total number of CGFFT iterations to those in a conventional CGFFT forward problem solution, where the initial estimate is chosen equal to the incident field. Such a study has not yet been reported for the case of a relatively large value of the exterior medium permittivity . When the contrast is small, only a limited number of CGFFT iterations is needed to obtain convergence.

In all our tests, “marching on in angle” was faster, except for the case where the spacing between the sources was larger than . With proper choices for the number of sources , for and for the CGFFT stop criterium, “marching on in angle” leads to a significant reduction in the computational effort. For example, with corresponding to a scattered field discretization error of a few percent (see Figure 3(b)), the relative reduction typically is 20 for Muscle1 , 52 for Leg2 —see also Figure 4 where the squares are below the dotted curve—and 32 for Muscle4 , see Table 2, which gives the relative reduction in the total number of CGFFT iterations as a function of and . The gain in efficiency is not as spectacular as with some examples in [13], for which a much larger number of CGFFT iterations is needed to obtain convergence, which leaves more room for improving the computational efficiency.

Table 2: Comparison of the total number of CGFFT iterations as a function of N and K for the conventional (i.e., incident field initial guess) and “marching on in angle” approaches. The relative reduction in the number of iterations, the CGFFT stop criterium (NMRS CGFFT), and the resulting NRMSE for on the grid also are indicated.

N
16
32
64
128
256
512

Muscle1

K = 16
No marching
92
108
112
136

Marching
96
111
127
139

Reduction (%)
−4
—
—
—

K = 32
No marching
188
220
224
280
312
344

Marching
151
163
216
256
267
320

Reduction (%)
20
26
4
9
14
7

K = 64
No marching
380
444
448
568
632
696

Marching
176
200
325
390
460
518

Reduction (%)
54
55
27
31
27
26

NRMSE CGFFT
2 × 10⁻³
10⁻³
2 × 10⁻⁴
5 × 10⁻⁵
10⁻⁵
2 × 10⁻⁶

NRMSE
5.9 × 10⁻³
1.7 × 10⁻³
4.4 × 10⁻⁴
1.1 × 10⁻⁴
2.8 × 10⁻⁵
7.1 × 10⁻⁶

Leg2

K = 32
No marching
1043
1272
1571

Marching
780
1015
1406

Reduction (%)
25
20
11

K = 64
No marching
2095
2555
3170
3489
3755
3981

Marching
878
1216
1896
2673
3239
3607

Reduction (%)
58
52
40
23
14
9

K = 128
No marching
4199
5112
6352

7533

Marching
1856
2466
3135

4596

Reduction (%)
56
52
49

39

NRMSE CGFFT
10⁻²
5 × 10⁻³
10⁻³
2 × 10⁻⁴
5 × 10⁻⁵
10⁻⁵

NRMSE
1.1 × 10⁻¹
3.2 × 10⁻²
7.9 × 10⁻³
1.9 × 10⁻³
5.1 × 10⁻⁴
1.2 × 10⁻⁴

Muscle4

K = 64
No marching
700
832
1128
1376
1648
1752

Marching
619
761
1032
1302
1563
1678

Reduction (%)
12
9
9
10
5
4

K = 128
No marching
1404
1664
2256
2760
3312

Marching
791
911
1531
2000
2586

Reduction (%)
44
45
32
28
22

K = 256
No marching
2804
3328
4512
5536

Marching
1034
1103
1679
2565

Reduction (%)
63
67
63
54

NRMSE CGFFT
2 × 10⁻²
10⁻²
2 × 10⁻³
5 × 10⁻⁴
10⁻⁴
5 × 10⁻⁵

NRMSE
6.4 × 10⁻²
2.6 × 10⁻²
6.7 × 10⁻³
1.7 × 10⁻³
4.5 × 10⁻⁴
1.3 × 10⁻⁴

It appears that the relative reduction decreases with increasing values of , that is, finer meshing. This can be explained as follows. The dimension of the vector space in which the solution needs to be determined increases as , while the discretization error, hence the recommended CGFFT stop criterion, decreases as . The error on the “marching on in angle” initial estimate does not depend on , while the number of iterations, which is needed to further reduce the error, increases and becomes less dependent on the choice of the initial guess.

It can furthermore be seen that the relative reduction increases with increasing numbers of sources . For example, when is doubled in the aforementioned examples, the relative reduction typically is 54 for Muscle1, 52 for Leg2, and 63 for Muscle4. The sources are closer, such that the error in the “marching on in angle” initial estimate is reduced. This is also illustrated in Figure 5, which shows the number of CGFFT iterations, with and without “marching on in angle” as a function of for Muscle1.

Figure 5: The total number of CGFFT iterations as a function of the number of sources for the “marching on in angle” and conventional approaches, for Muscle1 with .

Let us conclude by stressing the importance of choosing the largest possible value for the CGFFT stop criterion (NRMSE CGFFT), in order to get the most benefit from “marching on in angle.” When, for example, for Leg2 , the CGFFT stop criterion is reduced by a factor of 10, the number of CGFFT iterations is almost doubled, from 1216 to 2148, while the resulting reduction in the NRMSE on the total field, from 3.2 to 2.8, is not significant.

6.1.4. Scattering Matrix—Choice of

For the computation of the scattering matrix, (10) is expressed in transmitter/receiver positions on , and the summations are truncated to terms,

with , and where should be large enough for the aliasing effects to be negligible or for . The elements

are computed from (31) by means of a 2D FFT. The computational effort of this step is negligible: 10 ms for the example of Figure 4. For a centered circular homogeneous cylinder, the scattering matrix is diagonal with given by

with and the radius of and the wave speed in the cylinder, respectively. For high orders , the elements (33) decrease with increasing order at a rate which primarily depends on , for a given . This is illustrated in Table 3, which shows that the truncation errors in (31) are negligible when for the smallest cylinders, which yields double precision, and for the largest cylinder, which is largely sufficient to obtain single precision. The accuracy of the scattering matrix then is determined by the discretization error on .

Table 3: The orders and for which the elements of the scattering matrix have dropped to fractions 10⁻⁷ (single precision) and 10⁻¹⁴ (double precision), respectively, of their maximum absolute value, for the homogeneous cylinders of Table 1.

Cylinder
Muscle1
Air1
Air2
Muscle4

9
9
14
24

14
14
20
36

6.2. Direct Reflection by the Casing

The direct reflection from the casing, which is treated as a separate constituent in the expressions for the field on and for the equivalent surface current, does not depend on the object. Hence, these contributions can be computed once beforehand for a given water-filled scanner geometry. Let us therefore consider the case of a line source on in the casing without object. The spectral representation of the field in the complete region is then given by the first two sums in (26).

6.2.1. Exact Solution

The convergence of the first sum in (26) is extremely slow for observation points on . This is visible in Figure 6, in Table 4, and from applying the asymptotic approximation ([34], 9.3.1) for high orders ,

Figure 6: Absolute values of the terms in (26) for : (1) , (2) , computed with equation (40) for (3) Muscle1, (4) Leg2, and (5) Muscle4. The terms are indicated with +.

Table 4: The orders and for which the terms in the series in (8) and (35) have dropped to fractions 10⁻⁷ (single precision) and 10⁻¹⁴ (double precision), respectively, of their maximum absolute values, for different .

(cm)
27.6
25
15
5

(8):

1
0.91
0.54
0.18

10⁷
130
30
15

10¹⁴
283
52
22

(35):

0.91
0.82
0.49
0.16

97
60
27
15

211
124
44
22

Consequently, it is more efficient to compute with (3) or

where the reflected-field constituent is truncated at order . The convergence of the spectral representation of the direct casing reflection in (35) is better than that of (8) in the region ; see Figure 6 and Table 4. In this case, the asymptotic approximation of the term , for , yields

which is identical to that of a line source located on a contour with radius in a homogeneous background— thus is the geometric mean of and . We observed that the direct casing reflection reaches stable and accurate DP values in the 64 receivers on when choosing .

6.2.2. Expression from Equivalent Line Sources

Alternatively, in the region , the field can be regarded as if it originated from an equivalent surface current (21), which we replace with a discrete set of equivalent line sources with spacing on in a homogeneous background, denoting by the complex amplitude of the equivalent source at due to an excitation at . Imposing the identity of the reflected field constituent in (35) and the field generated by these sources, for ,

where we used the spectral representation (8) for , we obtain the following expression for the complex amplitude

According to (23) and since in this section , the field in the region then is expressed as the linear combination

The truncation to terms in the RHS of (37) may lead to significant errors on the field (39) in grid points that are near to . In Table 5, we compared the field on the mesh (39), generated by different numbers of equivalent line sources, with the exact solution (35) for meshes with sides 9 cm, 18 cm, and 36 cm. It follows that for the mesh with side 9 cm, already yields a very high precision; for the mesh with side 36 cm, yields a moderate precision, and is needed for high precision computations. Figure 7 shows an image of the amplitude and phase of the field (35) in the water-filled casing.

Table 5: The NRMSE of the field (39) in the empty casing as a function of the number of equivalent sources for different sizes of the mesh.

(cm)

9
5 × 
2 × 
2 × 
2 × 

18
6 × 
7 × 
4 × 
4 × 

36

2 × 
3 × 
5 × 

(a)

(b)

Figure 7: Amplitude (a) and phase (b) of the field (35) in the empty water-filled casing. An upper limit of 0.2 was chosen in the image of the amplitude.

6.3. Complete Configuration

Finally, we look into the computation of the fields for the complete configuration of an object in the casing. It is shown that the embedding approach maintains the accuracy of the forward problem solution in homogeneous space, if (order of the scattering matrix) and (number of equivalent line sources) are properly chosen. The additional computational effort for the embedding as such is also examined.

6.3.1. Scattered Field on the Receivers

Instead of solving (16) for , we take into account the separation of the direct casing reflection (26) and solve the resulting set of equations, truncated to terms, for

where and with . The elements then are computed with

The absolute values of and rapidly decrease as a function of , as is shown in Figure 6. As a consequence, the convergence of the spectral representation (26) of the field is much better than that of . Provided that is chosen according to Table 3, the error on only depends on the discretization error in the scattered field . We compared the values of with the exact solution as a function of for Muscle1 with , Leg2 with , and Muscle4 with , and we observed that the resulting NRMSE is almost identical to that for the homogeneous solutions; see Figure 3(b).

6.3.2. Total Field on the Grid

The total field on the mesh is computed as a linear combination of homogeneous solutions for equivalent line sources

where is given by (38) and where

The accuracy of in (42) depends on the number of equivalent sources , as was shown for the case of the empty casing in Section 6.2.2, on the discretization error of and on . We compared the values of with the exact solution as a function of : for Muscle1 with , the NRMSE is almost identical to that in Figure 3(a); for Leg2 with and for Muscle4 with , the NRMSE is approximately 1.5 times as high as the values in Figure 3(a). When the number of equivalent sources for Muscle4 is reduced to , the NRMSE for increases from 1.6 × 10⁻⁴ to 3.3 × 10⁻⁴, as could be expected from Table 5. Figure 8 shows the exact solution for the total field in the casing for Leg2.

(a)

(b)

(c)

(d)

Figure 8: Amplitude (a, c) and phase (b, d) of the field for Leg2 in water without casing (a, b) and with casing (c, d). The field is displayed over a square subregion of width 36 cm.

6.3.3. CPU Time

All summations containing complex exponentials, such as (35), (38), (43) and the right hand side in (40), are computed by means of FFTs; hence, the computational effort involved in these steps is negligible. The effort for the computation of the field on the T/R circle is primarily determined by the solution of (40), for which we used LU factorization: we observed CPU times of 5 ms for Muscle1 with , 30 ms for Leg2 with , and 0.2 s for Muscle4 with . The effort needed for the summations in (42) to compute the field on the mesh is much more important but remains lower than that needed for the CGFFT solution, as can be seen from the circles in Figure 4. It can be concluded that the embedding approach does not increase the CPU time with respect to the forward problem solution in a homogeneous environment, when only the fields on the observation contour are computed, and that it leads to a relatively small increase when the fields on the mesh are computed as well. The relative increase is less than 10 for the previously specified examples Muscle1 and Leg2 and less than 50 for Muscle 4.

7. Computational Procedure

We conclude this paper with a summary of the computational details. In analyzing the computational complexity, it should be kept in mind that the approach described in this paper was devised for use as forward scheme in inverse profiling, where an unknown configuration is reconstructed by matching the corresponding scattered field to a known measured field by linear or nonlinear optimization. For each new estimate of the configuration, the field caused by sources on the observation contour must be determined. When the optimization converges, the successive estimates gradually approach the desired optimum.

7.1. Homogeneous Environment

To demonstrate the efficiency of the scheme, we compare it with a straightforward implementation of the method of moments. For the object in a homogeneous environment, the first advantage is the second-order accuracy of the space discretization. In Figure 3, the conclusion from [12] that the error in the computed fields is of for decreasing or equivalently of for increasing , was confirmed for our test objects. Second, as shown in (5), the convolution structure of the continuous equation (2) was preserved. A straightforward evaluation of a matrix-vector product requires an effort of in each CG iteration step for unknown field values. Replacing these multiplications by two-dimensional FFT operations reduces the computational complexity to per step. Third, marching on in angle reduces the number of iterations. From Table 2, an acceleration by about 50 is observed. In all cases, the computational procedure is considerably more efficient than a straightforward matrix inversion, which requires a computational effort of , followed by matrix-vector computations at an effort of . In fact, the motivation for treating only the profile parameters as independent variables during the optimization in [5, 6, 18] was the efficiency of this forward scheme.

7.2. Complete Configuration

For the object inside the scanner, the conventional approach requires evaluating the field due to line sources in an empty scanner for at mesh nodes, each for observation points at mesh nodes. For each pair of nodes, the modified Bessel functions have to be computed for orders, which leads to a matrix fill time of . In the algorithm described in this paper, the Bessel functions are computed once for the argument , and we invert the truncated version (40) of (16). Moreover, for each line source, we synthesize the actual field (42) at points from fields in a homogeneous environment, which amounts to a total effort of . As mentioned towards the end of Section 6.3.3, only the last step leads to a relatively small increase in computation time. The proposed embedding approach thus is significantly more efficient than computing the fields in the complete configuration by means of Green’s functions of the empty casing [19, 20].

7.3. Guidelines

Last but not least, we enumerate the various steps of our algorithm, giving some guidelines for efficient application. The goal is to compute with accuracy the difference field on the T/R circle and the total field on the mesh for an object with a maximum size of . Based on the results given in this paper, we recommend the following procedures:

(1)Choose the smallest possible value for , which yields an accuracy , with the aid of Figure 3(b).(2)Choose the smallest possible value for , which yields a (much) better accuracy than on the mesh, with the aid of Table 5.(3)Choose the smallest possible value for , which yields a (much) better accuracy than for the given object size, with the aid of Table 3.(4)In general, ; hence, choose the number of forward problem solutions . For convenience, we have chosen in all our examples.(5)Compute ; ; and, if also the total field on the T/R circle is needed, .(6)Compute Green’s functions in a homogeneous environment.(7)Choose the CGFFT stop criterion, with the aid of Figure 3(a), and solve the forward problem in a homogeneous environment with “marching on in angle.”(8)Compute the elements of the scattering matrix and the elements and .(9)Compute the field on the T/R circle.(10)Compute and the field on the mesh.

8. Conclusion

In this paper, we have described a procedure to decompose the computation of electromagnetic fields in a relatively complicated configuration. The procedure allows (re)computing the field in part of the configuration, while the remaining part of the configuration and its electromagnetic response is left unchanged. By considering scattering in a homogeneous environment, the efficiency of the CGFFT procedure is exploited. Generating the initial estimate by “marching on in angle” accelerates the convergence of this procedure significantly.

The procedure has been applied to a standard scanner configuration for 2D inverse profiling. In the model, we have neglected the influence of the finite length of the antennas, the mutual coupling, and the variation in the properties of the individual antennas. Previous expertise [35] has shown that such effects can be handled in the calibration of the results, which is needed anyway. Reconstruction results for the idealized configuration have already been described in [18].

In the present paper, we have addressed the efficiency and accuracy of the forward algorithm and described the influence of the different tuning parameters in the algorithm. Results have been presented and discussed for canonical objects with representative values for the permittivity and the object dimension. A systematic procedure has been proposed for choosing computational parameters such as truncation limits and tolerances.

Appendix

A.1. Direct Reflection from Casing

Additional confirmation of the interpretation given in Section 5.4 can be obtained by considering the coefficients in (14) and (15) as the fundamental unknowns. The choice for in Section 5.1 was motivated by the circumstance that these coefficients determine the secondary contribution to the field incident on the dielectric cylinder in the equivalence principle discussed in Section 5.3. From a computational point of view, however, there is no preference. Repeating the analysis of Sections 5.1 and 5.2 for the coefficients now results in the operator equation

which leads to the power-series solution

which has a similar interpretation as the solution (25). For completeness, it should finally be mentioned that the coefficient vector can now be obtained from

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors would like to thank Dr. Kamal Belkebir and Dr. Amélie C. S. Litman, presently at the Institut Fresnel of the Université Paul Cézanne Aix Marseille III, France, and Dr. Bastiaan P. de Hon of the Eindhoven University of Technology for contributions to and discussions about the work presented in this paper. In particular, Dr. Belkebir was the first to suggest developing a generalization of existing “free-space” inverse-scattering algorithms to handle experimental data from the 434 MHz scanner developed by the third author.

References

W. C. Chew and Y. M. Wang, “Reconstruction of two-dimensional permittivity distribution using the distorted Born iterative method,” IEEE Transactions on Medical Imaging, vol. 9, no. 2, pp. 218–225, 1990.

N. Joachimowicz, C. Pichot, and J.-P. Hugonin, “Inverse scattering: an iterative numerical method for electromagnetic imaging,” IEEE Transactions on Antennas and Propagation, vol. 39, no. 12, pp. 1742–1753, 1991.

P. M. Meaney, K. D. Paulsen, and T. P. Ryan, “Two-dimensional hybrid element image reconstruction for TM illumination,” IEEE Transactions on Antennas and Propagation, vol. 43, no. 3, pp. 239–247, 1995.

A. Franchois and C. Pichot, “Microwave imaging – complex permittivity reconstruction with a Levenberg-Marquardt method,” IEEE Transactions on Antennas and Propagation, vol. 45, no. 2, pp. 203–215, 1997.

A. G. Tijhuis, K. Belkebir, A. C. S. Litman, and B. P. de Hon, “Theoretical and computational aspects of 2-D inverse profiling,” IEEE Transactions on Geoscience and Remote Sensing, vol. 39, no. 6, pp. 1316–1330, 2001.

J. De Zaeytijd, A. Franchois, C. Eyraud, and J.-M. Geffrin, “Full-wave three-dimensional microwave imaging with a regularized Gauss–Newton method— theory and experiment,” IEEE Transactions on Antennas and Propagation, vol. 55, no. 11, pp. 3279–3292, 2007.

J. De Zaeytijd and A. Franchois, “Three-dimensional quantitative microwave imaging from measured data with multiplicative smoothing and value picking regularization,” Inverse Problems, vol. 25, no. 2, article 024004, 2009.

A. H. Golnabi, S. D. Geimer, P. M. Meaney, and K. D. Paulsen, “Comparison of no-prior and soft-prior regularization in biomedical microwave imaging,” Journal of Medical Physics, vol. 36, no. 3, pp. 159–170, 2011.

F. Bai, A. Franchois, and A. Pižurica, “3D microwave tomography with Huber regularization applied to realistic numerical breast phantoms,” Progress In Electromagnetics Research, vol. 155, pp. 75–91, 2016.

R. E. Kleinman and P. M. van den Berg, “A modified gradient method for two- dimensional problems in tomography,” Journal of Computational and Applied Mathematics, vol. 42, no. 1, pp. 17–35, 1992.

A. Abubakar, P. M. van den Berg, and J. J. Mallorqui, “Imaging of biomedical data using a multiplicative regularized contrast source inversion method,” IEEE Transactions on Microwave Theory and Techniques, vol. 50, no. 7, pp. 1761–1771, 2002.

Z. Q. Peng and A. G. Tijhuis, “Transient scattering by a lossy dielectric cylinder: marching-on-in-frequency approach,” Journal of Electromagnetic Waves and Applications, vol. 7, no. 5, pp. 739–763, 2012.

A. G. Tijhuis, M. C. van Beurden, and A. P. M. Zwamborn, “Iterative solution of field problems with a varying physical parameter,” Elektrik, Turkish Journal of Electrical Engineering & Computer Sciences, vol. 10, pp. 163–183, 2002.

T. Sarkar, E. Arvas, and S. Rao, “Application of FFT and the conjugate gradient method for the solution of electromagnetic radiation from electrically large and small conducting bodies,” IEEE Transactions on Antennas and Propagation, vol. 34, no. 5, pp. 635–640, 1986.

P. Zwamborn and P. M. van den Berg, “The three dimensional weak form of the conjugate gradient FFT method for solving scattering problems,” IEEE Transactions on Microwave Theory and Techniques, vol. 40, no. 9, pp. 1757–1766, 1992.

Z. Q. Zhang and Q. H. Liu, “Three-dimensional weak-form conjugate- and biconjugate-gradient FFT methods for volume integral equations,” IEEE Transactions on Biomedical Engineering, vol. 29, no. 5, pp. 350–356, 2001.

A. G. Tijhuis, A. Franchois, W. H. A. B. Janssen, and A. P. M. Zwamborn, “Two-dimenional inverse profiling in a complex environment,” in Proceedings of the Microwave Imaging Methods and Techniques Workshop, European Microwave Week, Paris, October, 2000.

A. Franchois and A. G. Tijhuis, “A quasi-Newton reconstruction algorithm for a complex microwave imaging scanner environment,” Radio Science, vol. 38, no. 2, pp. VIC 12-1–VIC 12-13, 2003.

J. M. Geffrin, Imagerie Microonde: Etude Dun Scanner à 434 MHz Pour Applications Biomédicales, Ph.D, Thesis of the University of Paris XI, Orsay, France, 1995.

C. Gilmore and J. LoVetri, “Enhancement of microwave tomography through the use of electrically conducting enclosures,” Inverse Problems, vol. 24, no. 3, article 035008, 2008.

L. Crocco and A. Litman, “On embedded microwave imaging systems: retrievable information and design guidelines,” Inverse Problems, vol. 25, no. 6, article 065001, 2009.

R. Lencrerot, A. Litman, H. Tortel, and J. -M. Geffrin, “Measurement strategies for a confined microwave circular scanner,” Inverse Problems in Science and Engineering, vol. 17, no. 6, pp. 787–802, 2009.

P. C. Waterman, “Matrix formulation of electromagnetic scattering,” Proceedings of the IEEE, vol. 53, no. 8, pp. 805–812, 1965.

G. Kristensson and P. C. Waterman, “The T matrix for acoustic and electromagnetic scattering by circular disks,” The Journal of the Acoustical Society of America, vol. 72, no. 5, pp. 1612–1625, 1982.

B. Peterson and S. Ström, “T Matrix for electromagnetic scattering from an arbitrary number of scatterers and representations of E(3),” Physical Review D, vol. 8, no. 10, pp. 3661–3678, 1973.

W. C. Chew, L. Gürel, Y.-M. Wang, G. Otto, R. L. Wagner, and Q. H. Liu, “A generalized recursive algorithm for wave-scattering solutions in two dimensions,” IEEE Transactions on Microwave Theory and Techniques, vol. 40, no. 4, pp. 716–723, 1992.

A. M. van de Water, B. P. de Hon, M. C. van Beurden, A. G. Tijhuis, and P. de Maagt, “Linear embedding via Green’s operators: a modeling technique for finite electromagnetic band-gap structures,” Physical Review E, vol. 72, no. 5, 2005.

S. Mokhlespour, V. Lancellotti, and A. G. Tijhuis, “An extension of the linear embedding via Green’s operators method for the analysis of disconnected finite antenna arrays,” Progress In Electromagnetics Research M, vol. 49, pp. 141–151, 2016.

P. M. van den Berg and J. T. Fokkema, “Removal of undesired wavefields related to the casing of a microwave scanner,” IEEE Transactions on Microwave Theory and Techniques, vol. 51, no. 1, pp. 187–192, 2003.

E. Anderson, Z. Bai, C. Bischof et al., LAPACK Users’ Guide, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, Third Edition edition, 1999.

D. E. Amos, “A remark on Algorithm 644: a portable package for Bessel functions of a complex argument and nonnegative order,” ACM Transactions on Mathematical Software, vol. 21, no. 4, pp. 388–393, 1995.

D. Kahaner, C. Moler, and S. Nash, Numerical Methods and Software, Prentice-Hall, Inc., Upper Saddle River, NJ, USA, 1988.

M. A. Stuchly and S. S. Stuchly, “Dielectric properties of biological substances - tabulated,” Journal of Microwave Power, vol. 15, no. 1, pp. 19–25, 2016.

M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, Dover Publications, New York, NY, USA, 1970.

A. G. Tijhuis, K. Belkebir, A. C. S. Litman, and B. P. de Hon, “Multiple-frequency distorted-wave Born approach to 2D inverse profiling,” Inverse Problems, vol. 17, no. 6, pp. 1635–1644, 2001.