An Exponential Filtering Based Inversion Method for Microwave Imaging

In this paper, a new methodology based on the exponential filtering of singular values is adopted to solve the linear ill-posed problem of microwave imaging. This technique filters out the insignificant singular values and works as an efficient low pass filter to eliminate high-frequency noise from the estimated solution. Standard Tikhonov regularization has also proven to be a special case of this method. To show the effectiveness of this approach, various numerical examples of synthetic data and experimental data of Fresnel’s Institute are considered for the study. The reconstruction performance of this algorithm is quantified using the mean square error (MSE) and Pearson’s correlation coefficient (PCC). Further, the effect of noise on these metrics is presented. The results are compared with the standard Tikhonov regularization method, and it is observed that the proposed reconstruction algorithm provides accurate results compared to the standard Tikhonov regularization method.


Introduction
Microwave imaging (MI) is an important technology for estimating unknown objects by interrogating microwave signals [1]. It has wide applications in the fields of nondestructive testing and evaluation, biomedical engineering, subsurface sensing, hidden weapon detection, and remote sensing [1][2][3][4][5]. In general, these problems are non-linear and ill-conditioned in nature and require high computational time [6]. In scientific literature, several solution strategies have been proposed to solve such problems. Among them, the Born approximation is one of the most popular frameworks due to its speed and simplicity. This approach assumes that the relation between scattered field and permittivity profile is linear [7]. Therefore the imaging result is a solution to an ill-conditioned linear problem. To solve such problems, a regularization scheme needs to be introduced that stabilizes the solution by discarding the insignificant singular values [1], [8]. The solution of the regularized problem is well-behaved and close to the solution of the original illconditioned problem. Several techniques are available in the literature, such as spectral cutoff, Tikhonov regularization, and exponential filtering [9]. These schemes are called regularization by filtering because the only difference between them is the filter function. Since the main purpose of the filtering is to prevent the amplification of error in the problem due to insignificant singular values, the simplest approach to do this is to eliminate the contribution of insignificant singular values. The spectral cutoff is one of the simplest ways to do this. The main idea behind this method is to introduce a tolerance, for which all singular values below that tolerance are truncated [8]. However, finding the correct value of tolerance is a challenging task. The most commonly used regularization technique in MI is based on the Tikhonov method, which makes use of the L 2 -norm of the expected solution, thereby reconstructing the smooth image. The Tikhonov regularization method penalizes the residual norm and reduces the effect of small singular values. However, the exponential filtering system is considered to have very little bias towards the regularization parameter and is not subject to extra computational burden as it is implemented in the Tikhonov regularization [10]. The Tikhonov regularization method is an approximation of the exponential filtering technique [9]. Also, there is considerable literature on spectral cutoff and Tikhonov methods. None of the previous studies, however, investigated the exponential filtering method for the applications in MI. Therefore to further improve the imaging capabilities, this paper addresses this approach, which has hitherto been lacking in the scientific literature.
In this work, the reconstruction capability of exponential filtering is studied by analyzing filter factors. Further, the exponential filtering is applied to various numerical examples of synthetic data generated with dielectric and metallic objects. The merits of this approach are studied by comparing the results with standard Tikhonov regularization. Finally, the proposed algorithm is validated against an experimental database provided by the Institute of Fresnel, France [11]. In both cases, the results are encouraging and show improvement compared to the standard Tikhonov regularization.

Formulation
A schematic diagram of a two-dimensional (2D) microwave imaging configuration is shown in Fig. 1. The imaging region is represented by D and the observation region is represented by S. The unknown scatterer is illuminated using a set of incident fields generated by a single emitter (T x ) moving around the target. The scattered fields are collected with multiple receivers (R x ) located in the circular geometry that encloses the cross-section of the targets. The inverse scattering problem makes use of these measured scattered fields to determine the contrast function distribution.
Let us consider a 2D microwave imaging problem of transverse magnetic (TM) polarization. To determine the relationship between the contrast function and the scattered field, the scattering equations for non-magnetic materials can be represented as [1], where the spatial variables r, r ∈ (x, y) are the Cartesian coordinates of the measurement point and the source point, respectively. The terms e i and e t are the incident and the total electric fields in D, while e s is the scattered electric field in S. In equations (1) and (2), k 0 denotes the wavenumber of the background medium, G(r, r ) is the free space Green's function, and χ denotes the contrast function which carries information about the variation between the scatterer and the homogeneous background. It is defined as , where ε r and ε b denotes the relative permittivities of an object and the background, respectively. Equations (1) and (2) are solved by converting them to a discrete form with the help of method moments (MoM) [12]. The resulting matrix equations can be written as where G D and G S are internal and external radiation operators, respectively. The solution to find χ is non-linear [1]. Various methods are available in the literature to solve such problems. In this work, we consider the Born approximation (BA) model where the non-linear problem is approximated to a linear model. The BA is applicable when the object under test is a weak scatterer [5]. For strong scatterers, BA provides morphological features only. With this approach, the scattering equation (4) can be linearized by approximating the total field with the incident field, so that the only unknown in the equation is χ. Thus the resulting scattered field equation can be written as e s = G S e i χ.
This equation is linear and ill-posed, with the unknown parameter χ. In practice, scattered field data is contaminated by noise. Therefore, the direct inversion of this equation produces erroneous results. In this case, the regularization technique needs to be adopted to obtain a stable solution. Let us consider a summarized linear problem as where T is the ill-conditioned system matrix, and d is the scattered field vector. Here P represents the total number of observation points, and Q represents the number of cells in the test domain. In order to compute the optimal solution for this problem, the regularization strategy is needed that overcomes the instability of the problem. This work includes a description of two such techniques: Tikhonov regularization and exponential filtering. The characteristics of each method are discussed individually in the following section.

Tikhonov Regularization
Tikhonov regularization is a widely adopted technique for obtaining an approximate solution. In this regularization process, instead of directly solving the matrix equation (6) for the least-squares solution, we can solve the optimization problem that minimizes the objective function [13] where . 2 denotes the L 2 -norm, α is the regularization parameter, α χ 2 2 denotes the penalty/ regularization term used to obtain a stable solution that is more accurate and less sensitive to noise. The solution vector after solving this objective function can be written as where * denotes the conjugate transpose operator. By applying the singular value decomposition of system matrix as T = USV * , the regularized solution is given as where U and V are orthogonal matrices and are called as the left and right singular matrices of T, respectively, S is the diagonal singular value matrix, and S † is the diagonal inverse singular value matrix. We can rewrite the equation (9) as where γ denotes the system matrix rank, σ k is the k th singular value available along the diagonal of S, u k and v k denotes the k th left and right singular vectors respectively, and the term ψ k is called as the filter factor, which scales the solution component in the direction of v k such that the effect of error in the term u k , d is reduced. The filter factor for this method is given by The resulting diagonal elements of S † for the Tikhonov regularization approach is given by

Exponential Filtering
Even though the Tikhonov regularization is a widely used method, it is an approximation of the exponential filtering, often referred to as the asymptotic regularization or Showalter method [14]. The filter factor for this method is The eigenvector components corresponding to the insignificant eigenvalues are dampened by ψ EF k . The diagonal elements of the inverse singular value matrix for this case is given by The value of the filter factor is a function of α. Therefore, the selection of α is a key point for determining the quality of the reconstructed images. In this work, the value of α is evaluated using the state of the art regularized minimal residual method [15], [16], which solves the problem iteratively, and bypasses the calculation of an inverse matrix.
By performing the exponential series expansion, filter factors can be represented as By ignoring the higher-order terms for a case of σ 2 k α, we can rewrite the equation (15) as By applying the similar approximation to the Tikhonov filter factors of equation (11), we obtain This implies that the exponential filter factors reduce to Tikhonov filter factors when σ 2 α, which is the case for the last filter factors. Thus, Tikhonov regularization can be seen as a special case of exponential filtering.

Comparison of Filter Factors
The filter factors for Tikhonov regularization and exponential filtering are plotted as a function of singular values in Fig. 2. Filter factors vary depending on how quickly they reach saturation. Slow convergence indicates the significant smoothing in the obtained solution. A comparison plot shows that ψ EF k provides less smoothing than ψ TR k , and this implies that the solution obtained using exponential filtering is of higher quality than the Tikhonov regularized solution. Also, the effect of the regularization parameter is studied by considering different values of α (0.01, 0.1, and 0.5). It is quite evident from the figure that, for small values of α, the filter factors rate of attaining the saturation level (unity) is faster. This condition directly produces an inverse solution, but it results in inadequate filtering also. On the other hand, the value of α too high leads to a smooth curve, whose rate of reaching saturation is slower, and it yields in oversmoothing. The stability of the solution is ensured by the filter factor and the accuracy is controlled by the value of α. Thus the selection of proper values of ψ k and α is very crucial.
The filter factors are also plotted as a function of σ 2 /α in Fig. 3. Here, a regularization parameter of 0.01 is used. The figure clearly shows that the Tikhonov and exponential filtering become identical when σ 2 α (i.e. σ 2 /α 1). Both Tikhonov and Exponential filtering methods exhibit low pass filter characteristics. The limitation of Tikhonov regularization is that it filters out some of the significant singular values in the transition band. Hence, the obtained image is less accurate. In the transition band, the exponential filtering follows similar to the low pass filter. As a result, most significant singular values are passed, which contribute towards obtaining a more accurate solution.  By ignoring the higher-order terms for a case of we can rewrite the equation (15) as

Simulation Results
In this section, the performance of the proposed algorithm is validated through various numerical examples of synthetic and experimental data consisting of metallic as well as dielectric objects. In all the inversion examples, we will reconstruct the profiles for their location, shape as well as permittivity using the proposed reconstruction algorithm.

Synthetic Data Results
In this section, we show some reconstructions obtained using simulated data. The domain of interest is a 20 cm × 20 cm square, and it is divided into 2601 cells (51 × 51). In all the following examples, the measurement domain consists of 36 receivers that are uniformly placed in a circular geometry. The operating frequency can be chosen to be the frequency at which the length of the investigation domain is approximately one wavelength [17]. However, the resolution is improved when a higher frequency is used. Therefore, an operating frequency of 3 GHz is used.
The accuracy of the results can be quantified using the following figures of merit.

• Mean Square Error (MSE):
It is a measure of reconstruction quality based on comparison with the reference profile [18]. This error parameter is defined as where χ andχ denotes the actual and estimated contrast function values, respectively.

• Pearson Correlation Coefficient (PCC):
This coefficient represents the accurate detectability of the object by measuring the degree of correlation between the reference profile and the estimated profile [10]. It is given as where cov represents the covariance, and s represents the standard deviation. PCC value ranges from -1 to 1. Reconstructed images require a higher value of PCC.

Two Dielectric Objects
In this case, the value of ε r is selected in such a way that G D χ < 1 for the frequencies considered. Therefore, the target can be treated as a weak scatterer, and the BA can be applied for better reconstruction results [19]. As shown in Fig. 4(a), this example consists of two identical, circular cylinders of 2 cm radius and separated by a distance of 3 cm in free space (ε b = 1). The scatterers are dielectric with the permittivity value of ε r = 1.3 (χ = 0.3). Here, α is calculated using the regularized minimal residual method [15]. In this method, the residual is calculated at each iteration and the solution is updated accordingly. This process is repeated until the relative residual becomes less than or equal to τ = 10 −5 . At every iteration, the value of α is automatically chosen using Matlab's one-dimensional minimizer function fminbnd, which attempts to find the minimum of a singlevariable function within a fixed interval (here, the interval is between 0 to 1 for the considered examples).  Figure 4(e) shows the one dimensional cross-sectional plot, which gives the contrast function value retrieved along the x-axis (y = 0). Here, TR and EF symbolize Tikhonov regularization and exponential filtering, respectively. The results indicate that the accuracy of the exponential filtering is better than the standard Tikhonov method.
Although the results obtained using the proposed method are visually excellent, to objectively show the same, the figures of merit (MSE and PCC) are presented in Tab. 1. Here, the effect of noise on the quality of reconstruction is inspected by corrupting the simulated data with In this section, we show some reconstructions obtained using simulated data. The domain of interest is a 20 × 20 cm square and it is divided into 2601 cells (51 × 51). In all the following examples, the measurement domain consists of 36 receivers that are uniformly placed in a circular geometry. The operating frequency can be chosen to be the frequency at which the length of the investigation domain is approximately one wavelength [17]. However, the resolution is improved when a higher frequency is used. Therefore, an operating frequency of 3 GHz is used. × 100%. As can be seen, the PCC values for both methods are very close. This implies that the detectability of the proposed approach is almost identical to the standard Tikhonov method. However, the MSE values indicate that the proposed method is very much helpful in terms of improving the image accuracy. We observe that in all the cases, there is a reduction in error of approximately 83% for the proposed approach.

Austria Profile
In this example, the applicability of the proposed algorithm is tested on the Austria profile, which is the well-known target configuration in the inverse scattering community. Figure 5(a) shows the Austria profile (ε r = 1.3) with a homogeneous background (ε b = 1). It consists of two disks and one ring structure. The ring is centered at (0, -2) cm and the disks are located at (-3.5, 7) cm and (3.5, 7) cm. The radius of each disk of 2 cm, whereas the ring has an outer radius of 6 cm and an inner radius of 4 cm. With reference to Fig. 5(b)-(c), it is clear that the algorithm yields good reconstruction in terms of image accuracy. The MSE for the exponential filtering (0.8020) is less than the Tikhonov regularization (1.1422). Thus the proposed method achieves a statistically significant improvement compared to the standard Tikhonov method. Also, the PCC value is 0.7723 for the Tikhonov method and 0.7725 for exponential filtering.

Multilayer Dielectric Object
The multilayer dielectric cylinder configuration is often used for inversion procedures imaging applications. An example consists of two-layer circular non-concentric cylinders with radii of inner and outer cylinders, 2 cm and 6 cm, respectively. Furthermore, the relative permittivity values of the inner and outer cylinders are equal to 1.35 and 1.15, respectively, and are centered at (2, 0) cm and (0, 0) cm.
The actual profile of the non-concentric circular cylinder is given in Fig. 6(a), and the final reconstructed outputs (distribution of the contrast function) at 3 GHz are reported in Fig. 6(b)-(d). It is observed that the reconstructed image using the exponential filtering method is quite accurate in terms of target position, size, and dielectric property. The MSE value for the exponential filtering (0.2782) is less than the Tikhonov method (0.4329). Moreover, the PCC value is 0.9280 for the Tikhonov method and 0.9283 for the exponential filtering.

PEC Object
The proposed reconstruction algorithm is developed for dielectric scatterers, so it is not known how well this algorithm works for targets with perfect electric conductor (PEC) objects. In the case of PEC, the electromagnetic field does not enter into object and the scattered fields come from its boundary only [20], [21]. The imaging result is obtained by assuming the scatterer as a highly conducting dielectric cylinder so that the mathematical formulations derived for the dielectric object can be applied to the PEC objects [22]. Here, the target configuration consists of a circular PEC cylinder of radius 2 cm placed in a homogeneous free space at a distance of 3 cm from the center of the test domain, as shown in Fig. 7(a). The corresponding simulation results at a working frequency of 3 GHz are illustrated in Fig. 7(b)-(d). Since the PEC targets have infinite conductivity, only the boundary of these scatterers should be determined [22]. Therefore, the actual target permittivity value has not been defined in Fig. 7(d) and only the boundary of the scatterer has been indicated. It is quite evident from results that the exponential filtering provides better accuracy than the Tikhonov regularization. Thus the proposed reconstruction algorithm also works well for the PEC objects.

Experimental Data Results
In order to check the performance of the proposed algorithm under realistic conditions, experimental data acquired in the laboratory of Fresnel's institute, Marseille [11] is considered. The scattered field data is acquired in an anechoic chamber using the experimental setup shown in Fig. 8. In the setup, linearly polarized horn antennas are used. The transmitting antenna is fixed on the horizontal axis at a distance of 720 mm ± 3 mm, and the receiving antenna rotates around the target at a distance of 760 mm ± 3 mm. The target has an elongate shape along the vertical axis to form a two-dimensional structure. In Fig. 8, θ t represents the angle between the transmitting antenna and the target, whereas θ r denotes the angle between the transmitter and the receiver, whose angular range is from 60 • to 300 • . In this case, 36 transmission antenna positions are used. In the experimental case, four different imaging configurations (dielectric or metallic) are considered. The crosssectional dimensions of these scatterers are indicated in Fig. 9. The dielectric target consists of one ( Fig. 9(a)) or two ( Fig. 9(b)) filled dielectric cylinders, with a permittivity value of ε r = 3 ± 0.3, so it is considered a strong scatterer. Two metallic targets with rectangular and 'U-shaped' crosssection are also studied.

Single Dielectric
This example relates to the reconstruction of a single circular dielectric cylinder of 1.5 cm radius, located at a distance of 3 cm from the origin, as shown in Fig. 9(a). The reconstruction is carried out at an operating frequency of

Two Dielectrics
The dielectric profile of this example is shown in Fig. 9(b). It consists of two identical circular cylinders of a 1.5 cm radius and separated by 9 cm. The cylinders are separated along the horizontal axis by a distance of 9 cm. The imaging results for this target are illustrated in Fig. 11. Here again, better reconstruction is obtained using the proposed algorithm.

Metallic Rectangle
In this case, the target configuration consists of a metallic rectangular cylinder of dimensions (2.45 × 1.27) cm 2 , as shown in Fig. 9(c). The measurement is carried out at a frequency of 4 GHz, and the corresponding simulation results are reported in Fig. 12. As the object under inspection is small in size, the imaging domain is considered as 12 cm × 12 cm. It is observed that the profile of the target is clearly and quite accurately reconstructed by the exponential filtering approach.

Metallic 'U'
Finally, in this example, a more complicated object profile of 'U-shaped' metallic cylinder with a dimension of (8 × 5) cm 2 is inspected, and the corresponding results at an operating frequency of 4 GHz are presented in Fig. 13. Here, the test domain is taken as 16 cm × 16 cm. The results led to a similar conclusion where the exponential filtering results are more accurate than the Tikhonov regularization. x(cm)

Conclusion
This paper describes a new algorithm based on exponential filtering to solve the linear, discrete ill-posed problem of microwave imaging. The purpose of the filtering method is to suppress the contribution of lower singular values and keep significant singular values. As a result, the computed solution is more accurate. The linearized scattering model of Born approximation is considered. The required regularization parameter has been calculated using the regularized minimal residual method. The results are compared with the state of the art Tikhonov regularization for various numerical examples of synthetic data and experimental data (provided by Institut Fresnel). To show the effectiveness of the algorithm, results are also obtained with different signal to noise ratios of data. The results conclude that this method is able to reconstruct the dielectric scatterers with high accuracy. The same is validated for metallic objects. As future work, this model can be extended to the case of 3D.