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.
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.
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
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)
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.
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.
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.
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 .
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 ,
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.
(a)
(b)
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−4 to 3.3 × 10−4, 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)
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.