Fast quantitative microwave imaging with resolvent kernel extracted from measurements

A quantitative imaging method is proposed based on microwave measurements where a direct inversion in real space is employed. The electrical properties of penetrable objects are reconstructed using a resolvent kernel in the forward model, which is extracted from calibration measurements. These measurements are performed on two known objects: the reference object (RO) representing the scatterer-free measurement and the calibration object representing a small scatterer embedded in the RO. Since the method does not need analytical or numerical approximations of the forward model, it is particularly valuable in short-range imaging, where analytical models of the incident field do not exist while the fidelity of the simulation models is often inadequate. The experimentally determined resolvent kernel inherently includes the particulars of the measurement setup, including all transmitting and receiving antennas. The inversion is fast, allowing for quasi-real-time image reconstruction. The proposed technique is demonstrated and validated through examples using simulated and experimental data. Its performance with noisy data is also examined. The concept of experimentally determined resolvent kernel in the forward model may be valuable in other imaging modalities such as ultrasound, photonic imaging, electrical-impedance tomography, etc.


Introduction
Recently, active microwave systems have been deployed in a variety of short-range imaging applications such as concealed weapon detection, through-the-wall imaging, nondestructive testing, medical diagnostics, etc; see, e.g., [1][2][3][4][5][6][7][8]. In these applications, the inspected object is close to the microwave source and the sensing antennas, i.e., at a distance comparable or smaller than at least one of the following three measures: (i) the target's size, (ii) the size of the antenna or the sensor array, (iii) the wavelength. The goal is to reconstruct the dielectric properties (the real and imaginary parts of the complex permittivity) as functions of position. In quantitative imaging, these functions yield values that estimate the actual property distributions in the imaged object. In contrast, qualitative imaging only aims at localizing the abnormalities, where the properties differ from those in the normal state of the object.
Among the quantitative methods, the Born iterative methods [9,10] and the model-based optimization methods [11] are the most common in microwave near-field imaging. All of these employ iterative strategies where numerical forward models are updated to match the measured data. Thus, the reconstruction is often time-consuming and its convergence is critically dependent on the fidelity of the forward model. In the near-field imaging of complex objects such as tissues, luggage or structural components, achieving adequate fidelity of the forward model is challenging [12]. The major obstacle is not the numerical accuracy of the simulation model as this can be ensured by proper discretization. The model fidelity is mostly corrupted by the incomplete representation of the acquisition setup and factors such as fabrication tolerances, electromagnetic interference, positioning errors, environmental uncertainties, etc. Such errors, commonly referred to as modeling errors, have a detrimental impact on the inverse-problem solution similar to that of the noise in the measured data.
The poor performance of the direct inversion methods in near-field imaging applications is due to the assumption that the incident field and/or the Green function are in the form of a plane, cylindrical or spherical waves. This is rarely valid in the near zone of the transmitting and receiving antennas. In [15][16][17] and [21,22], this limitation is overcome by making use of simulation models of the incident field and the Green function distributions. Unfortunately, the fidelity of the simulation models, although better than the analytical approximations, may also be questionable for the same reasons as the simulations used in the optimization-based quantitative reconstruction methods.
Here, we propose a direct-inversion imaging method, which yields quantitative images while offering real-time performance similar to that of the qualitative approaches. The quantitative reconstruction becomes possible due to a calibration measurement of a known electrically very small (point) scatterer, which is placed in the known environment formed by the background medium and the acquisition system. Here, this environment is referred to as the reference object (RO), the measurements of which provide the baseline (or RO) data. When the point scatterer is placed in the RO, the resulting object is referred to as the calibration object (CO), the measurements of which yield the calibration (or CO) data.
The CO response is in effect the system point-spread function (PSF). Mathematically, it is the product of the RO Green function, the total field at the point scatterer and the known contrast of the point scatterer. The total field in turn is proportional to the incident field (the field in the RO) according to the localized nonlinear (LN) approximation [23]. Thus, the point scatterer in the CO effectively serves as a scattering probe inside the RO that provides sampling of the resolvent kernel of the scattering problem specific to the acquisition setup. The measured PSF and its respective resolvent kernel characterize the imaging system quantitatively. Moreover, they provide higher fidelity than any analytical or simulation model.
As a first-order approximation, the proposed method views the scattering from the object under test (OUT) as a convolution of its contrast function with the resolvent kernel. Unlike existing deconvolution techniques, e.g., [28,29], here the inversion is performed in real space, i.e., Fourier transforms are not employed. Most importantly, the images are quantitative, which, to our knowledge, is not achieved by any of the existing direct-inversion methods.
The utility of the proposed method is two-fold: (i) fast quantitative solution of weakscattering problems where multiple scattering is negligible, and (ii) providing a good starting point for an iterative reconstruction of complex targets.
The paper is organized as follows. The forward model based on the measured system PSF is formulated in section 2. The PSF acquisition through the CO measurements is presented in section 3. Section 4 describes the raster-scanning measurement setup used in all validation examples. The image formation is summarized in section 5. Section 6 presents the validation examples and the noise analysis. Discussion and conclusions are given in section 7.

Theory
2.1. Forward model based on the measured system PSF The microwave imaging system employs a network of N r receivers (Rxs) and N t transmitters (Txs), all of which are equipped with antennas for signal acquisition in the desired frequency band and over the desired surfaces (typically of canonical shapes such as planes, cylinders or spheres). In frequency-sweep acquisitions, the data are collected at N f sample frequencies f m (m = 1, …, N f ), one frequency at a time. The scattering parameters (or S-parameters) S ij m ( ) (i = 1, …, N r and j = 1, …, N t ) are the responses of interest. They represent the complex ratios of the received-to-transmitted signals, where the subscripts i and j indicate the Rx and the Tx antennas, respectively. At each frequency, the N t Txs are switched on, one at a time, while the N r Rxs record the scattered signals. Thus, at each (mth) frequency, = × N N N s t r responses are collected denoted as S n m ( ) (n = 1, 2, …, N s ). Note that each n corresponds to a specific (i, j) pair.
When the jth Tx antenna is excited by a voltage V , where r denotes position in the observation domain, R i m ( ) is the linear operator, which describes how the ith Rx antenna converts the field distribution in its vicinity into a voltage signal, and Ē r ( ) To simplify the notations, in the following, the normalized field distribution will be simply denoted as A set of responses S n m RO, ( ) is first obtained by measurements of the known RO (target-free host environment). These serve as the baseline data. Let us describe the dielectric-property distribution of an imaged object in volume V through the relative complex permittivity  Here, ̲ I is the identity tensor and δ is the Dirac delta function. With reference to figure 1, the imaged volume V is discretized uniformly into N v voxels, the volume Ω v of which is electrically small and is less than the expected spatial resolution limit 1 . Let the center of the pth voxel be at ′ ∈ V r , p where p can assume any integer value Figure 1. Illustration of the discretization of the imaged volume. The acronyms Rx and Tx stand for receiver and transmitter, respectively. 1 The spatial resolution can be analytically estimated by making use of the far-zone assumption. This estimate, however, may be too conservative as the spatial resolution is expected to be better in near-field measurements where evanescent field exists. Thus, the size of the small scatterer in the CO should be as small as possible, yet, it has to be sufficient to ensure detectability. of the pth CO can be obtained by substituting (7) into (4):  From (8) and (11), it follows that This is the PSF of the nth response of the imaging system for a point scatterer at the pth voxel. If the system is invariant upon translation or rotation, this PSF can be used to characterize the nth system response to point scatterers at other locations by making use of coordinate translations.
Further, because the voxel-size scatterer is small, the internal total field ′ E r ( )  which is in accordance with the well-known LN [23] and quasi-analytic (QA) [30][31][32] approximations. In (13), the tensor ̲ Γ m CO ( ) depends on the contrast δε m CO ( ) and the shape of the small target. It tends to the identity tensor ̲ I asymptotically as δε m CO ( ) decreases, leading to the linearized Born model of scattering. Also, in the static limit of the LN approximation and in the QA approximation, ) is a scalar. Using (13), (12) can be written as Next, the forward model of the scattering from an OUT can be written using (1) and (4) as: Neglecting the coupling among the scattering voxels and employing the linear relationship (13) between the total field in the OUT and that in the RO, (15) can be expressed in terms of the resolvent kernel (14) as  ( ) is to unity. The forward model in (17) inherently incorporates the Green function of the particular acquisition setup. As such, it is more reliable than any analytical or simulation model.

Power maps
and discretizing the integral into a sum over the contributions of all voxels, (17) can be written as    The power map bears similarity to the 'sensitivity map' in sensitivity-based imaging [21,22]-it can be viewed as the sensitivity map scaled by the permittivity contrast of the small scatterer in the CO. In [21,22], however, the sensitivity maps are built using incidentfield information acquired from simulations or semi-analytical models, which are not likely to possess sufficient fidelity for near-field measurements.
The power map also bears mathematical similarity to back-propagation or time-reversal schemes, although here the Tx and Rx arrays need not be co-located. This becomes clear if (20) is viewed as the system of data and A is the × N N s v system matrix built from the measured PSFs of the acquisition system. In the context of time reversal, A can be viewed as a discrete representation of the measured background coherent PSF [24]. Therefore, the quantity ⋅ A b * , equivalent to the power map definition in (21), can be viewed as an estimate of the projection of the data onto the space defined by the system singular vectors.
To this end, the power maps yield only qualitative imaging results, i.e., there is no property estimation involved. With wideband data, the OUT power maps at all frequencies are usually combined (often after proper normalization [21,22]) to obtain a single qualitative image of the OUT.
To gain a deeper insight into the power map, we substitute (20) into (21) and obtain is not necessarily zero but it is expected to be small in magnitude. Increasing the distance between the pth and qth voxels leads to decrease in the magnitude of On the other hand, when p = q, (24) becomes First, the linear system (23) is written in a matrix form, consisting of the CO power maps, m m ( ) is an N v × 1 vector containing the OUT power map, and τ m ( ) is an N v × 1 vector of the unknown property contrast ratios defined by (19). To build M , m ( ) all necessary CO measurements must be performed and the respective CO power maps to be built using (24). To build m , m ( ) the OUT measurements are needed from which the OUT power map is built using (21). Since the small scatterer in the CO has known contrast δε , of the OUT can be computed using (19). Finally, the actual property values of the OUT can be calculated from Δε

⎧ ⎨ ⎩
Such power map would imply infinitesimal spatial resolution (e.g., the system employs a very short wavelength) and noise-free measurements. In this ideal case, M m ( ) would be a diagonal matrix, which is well-conditioned and easily invertible. In the extreme opposite case, where the whole imaged volume is below the spatial resolution limit (e.g., the system employs a very low frequency), all elements of become similar. Then, M m ( ) is illconditioned and cannot be inverted.
Also, due to the large matrix size, the solution of (26) through direct inversion or LU decomposition is not practical. Here, (26) is solved as a linear least-square (LS) problem, which is cast as Here, || ⋅ || is the l 2 norm and τ m Once the solution τ m * ( ) of (31) is available, the estimated property distribution of the OUT at the mth frequency is calculated as In order to evaluate the accuracy of the property estimation in all subsequent examples, we define a relative root-mean-square error (RRMSE) as is the true property distribution of the OUT.
We note that, in principle, the τ m ( ) coefficients could also be found by solving (20) directly, i.e., without constructing the power maps first. Note that (20) can be cast as a rectangular system τ = b A . However, if the data b are not entirely in the range of the forward operator (represented by A), then such direct solution will not produce a valid result. On the other hand, under the same conditions, the system τ = A A A b * * , which is in effect (26), may have a solution and this solution will minimize the least-square error of (20); see, for example [24]. This also explains why it is advantageous to obtain A experimentally as this provides the physically correct functional space spanning all possible response functions acquired with the particular setup.
In addition, constructing the power maps brings the following three advantages. (i) The images provided by the absolute values of the power maps are generated practically instantaneously, which enables real-time qualitative imaging. (ii) A good initial guess of the OUT property distribution can be obtained once the power maps are generated; see (32). This is critical in dealing with nonuniqueness and in accelerating the iterative solution to (31). (iii) Multi-frequency data combination is efficiently done using the power maps (see section 5.3 below) thus enabling fast processing of wideband data where we may be dealing with thousands of frequency points. Such number of data points would render the direct solution of (20) impractical.

The limitations of the proposed method
The reconstruction method described above rests on the forward model in (18) obtained through two assumptions: (i) the mutual coupling between the scattering centers and the multiple scattering effects in the OUT are negligible, which reduces (15) to (16), and (ii) the ratio Γ ( ) is close to unity, which reduces (17) to (18). The first limitation is typical for all direct-inversion methods. If multiple-scattering effects are prominent in the OUT, they would lead to image artifacts. Such artifacts could be suppressed if the property reconstruction obtained with the proposed method is further refined using an iterative technique such as the iterative Born method [9], the distorted Born iterative method [10], or a model-based optimization method [11].
We next examine the limitations stemming from the assumption that the ratio ( ) is close to unity. As shown in [23], at the low-frequency (or dc) limit of the LN approximation, the coefficient Γ ,   These conditions can serve as guidelines when choosing the RO and CO properties for the system calibration. If the conditions in (37) and (38) are not observed, the ratio in (36) should be used in conjunction with the forward model in (17).
Assuming that the estimation inaccuracy stems only from the approximation of and not from multiple scattering effects), an acceptable estimation range for the OUT complex permittivity distribution can be obtained within a specified error level α, Finally, we emphasize that size of the small scatterer in the CO limits the resolution of the images, i.e., resolution better than the CO size cannot be achieved. This is why it is critical to use a CO of the smallest possible size (yet sufficiently large for acceptable SNR of the calibration data) so that the reconstruction algorithm can achieve its best performance in terms of resolution. We note that, in near-field imaging, the resolution is strongly influenced by the variability of the near-zone incident field rather than the wavelength. From that point of view, since the proposed approach utilizes experimentally acquired resolvent kernel, it does allow for the retrieval of near-zone scattering information as long as the CO size is sufficiently small.

CO measurements
The crux of the proposed method is in the measurement-based characterization of the imaging system through the acquisition of its PSF for the given RO. For a heterogeneous RO, the PSF is a function of the position of the voxel-size scatterer, which would require a measurement On the other hand, if the background is assumed homogeneous in all three dimensions or invariant in at least one of the coordinate variables, the number of the required calibration measurements is drastically reduced. This is because the PSF corresponding to a scatterer position along any coordinate variable, with respect to which the background properties are constant, can be obtained from a single PSF measurement via coordinate translation. The assumption of a homogeneous or layered background is common in imaging and it is also exploited here for the case of a planar raster-scanning acquisition system. As explained next, our system is characterized only by N z calibration measurements, where N z is the number of the imaged layers along the range (or depth) direction, here assumed to be along z.
The number of image samples along range must be such that the sampling interval Δz is equal or smaller than the expected range resolution δz. The latter can be analytically estimated in the case of far-zone measurements and it is given by δz = c/2B [16,33], where c is the speed of light in the respective medium and B is the bandwidth of the acquired signals. As the field variation with range is faster in the near zone than in the far zone, it is recommended to choose Δz < δz. Insufficient sampling along range mostly impacts the range resolution in the final images. Note that the range sampling positions determine the range locations of the small scatterer used in the CO measurements. Figure 2 shows an example configuration of the planar raster-scanning acquisition system. The space between the scanning planes (separated by a distance d) is the volume where the imaged object, i.e., the RO, the CO or the OUT, is placed. Two antennas (e.g., half-wavelength dipoles) face each other along boresight (the z-direction) and form a two-port network. In general, the antennas can be switched antenna arrays or dual-polarization antennas, in which case the system is multi-port. Each antenna is able to work either in Rx or in Tx mode. In a two-antenna system as the one shown in

System arrangement
, 1 where u = 1, …, N x , v = 1, …, N y , are the indices of the samples along x and y, respectively while N x and N y are the total numbers of sampling points along the x and y directions.
The number of acquired signals at each sampling location with our two-antenna system is This number can be further increased if more Tx and Rx antennas are used in a switched array. At each sample frequency, the scanning system provides a total of

Calibration data acquisition with planar raster scanning
When the RO is homogeneous or layered, the responses measured by the antennas are invariant to translations along x and y. On the other hand, the CO consists of the RO and a small scatterer embedded in it. Thus, in order to obtain the system PSF for all points in a ′ = z const plane, only one scan is necessary. Typically, the small scatterer is embedded in the

Summary of the image-formation procedure
The image-formation procedure is summarized here for the case of a planar raster-scanning acquisition system. It can be easily adjusted for other scanning schemes. There are two main procedural stages: the system calibration and the OUT imaging.
The system calibration consists of: (i) acquisition of the RO responses, (ii) acquisition of the CO responses, and (iii) the generation and storing of the CO power maps. The CO measurements are carried out as described in section 4.2, with only one scan per range location.
The OUT imaging consists of: (i) acquiring the responses of the OUT, (ii) generation of the OUT power map, (iii) generation of the OUT permittivity map. The detailed steps of the procedure are summarized below.

System calibration
Step 1 Acquire the baseline signals S n m RO, ( ) (RO measurement).
Step 2 Place the small scatter at the center ′ r 0 of the desired ′ = z const planes respectively and acquire the signals S n m CO, 0, ( ) (CO measurement).
Step 10 Solve (26)  Step 11 Obtain ε ′ r ( ) for the linear systems in (26). Once built, these matrices are repeatedly used to process the data acquired with various OUTs.

Multi-frequency combination (used with wideband data)
The multi-frequency strategy is recommended when thousands of frequency samples may need to be processed. First, the power maps of the OUT and the CO at all N f frequencies are combined via The normalized CO and OUT maps obtain through (47) are next used to fill the matrix M and the vector m in (26), respectively. Following step 9 to step 11 from the procedure at a single frequency, (26) We note that the multi-frequency combination can be done via simple averaging of the permittivity distributions obtained at all single-frequency solutions. However, the method described above produces better results. This is because the imaging resolution of the directaveraging approach is badly influenced by the low-frequency solutions. This is not the case in the proposed multi-frequency scheme, which reduces the impact of the low-resolution images at an earlier stage by combining the power maps. Also, notice that the LS solver is only called once with the multi-frequency combination scheme. The planar raster-scanning acquisition system described in section 3 is emulated with a full-wave electromagnetic numerical solver based on the method of moments [34]. We attempt to image an F-shaped dielectric bar (i.e., the OUT) in the z = 0 plane; see The RO is set to be air, i.e., ε ε ε Since the background medium is homogeneous, in the absence of numerical errors, the S-parameters should be the same at all the sampled positions. Thus, at each frequency, only one simulation of the RO measurement is performed.
In the simulation of the CO measurement, a small dielectric cube of 1 cm 3 is placed at the origin (the center of the z = 0 plane). The small scatterer has a property contrast of 0.1 in the real permittivity only, i.e., δε = 0.1, m CO ( ) which is set to be frequency independent. For another z′ = const plane, the dielectric cube is placed in the center of this plane and the CO data acquisition is repeated.
In the simulation of the OUT measurement, the F-shaped dielectric bar consists of 16 dielectric cubes. Each cube has the same size of 1 cm 3 and has a property contrast of 0.2 in the real permittivity only, i.e., ε m The OUT contrast is set to be frequency independent as well.
The scanning aperture in the OUT acquisition is 20 × 20 cm 2 with N x = 21 and N y = 21, respectively. In order to perform the coordinate translation given in section 4.2, the scanning aperture is enlarged to 40 × 40 cm 2 in the CO acquisition providing (2N x − 1) and (2N y − 1) sample points along x and y, respectively. The sampling step size along x and y is 1 cm.
The volume between the scanning planes is divided evenly into five layers along the z axis. Each layer is 1 cm thick. The middle three layers, i.e., the z = −1, 0, 1 planes, are being reconstructed.
6.1.2. Power maps of COs. Using the S-parameters acquired with the CO simulations, the database of CO power maps is built using (24); see also step 5 in the system calibration stage Figure 3. Configuration of the simulation setup for the OUT measurement. The OUT is an F-shaped dielectric bar that is placed at the middle layer. described in section 5. These complex-valued maps represent the system PSF and they are crucial for the subsequent property estimation obtained by solving (26). Figure 4 shows the normalized magnitude plots (in dB) of the power maps of the CO at 3, 7, 11 and 15 GHz. The small scatterer is placed at the center (0, 0, 0) cm, where the peak values are clearly observed. If the scatterer is repositioned, the map peak would move in correspondence with its location. It is observed that the map approaches the ideal delta function distribution as the frequency increases, which indicates the spatial resolution is better at higher frequencies. Due to the source polarization (the dipoles are x-polarized), steeper slopes can be observed along the y direction. The improved spatial resolution in the y direction is consistent with the fact that the dipoles' patterns in this direction are broader, which is another prerequisite for better spatial resolution. On the other hand, there are more artifacts along the y direction. This is due to the finite size of the scanned surface which supports an angle (with respect to the center of the imaged volume) that is narrower than the antenna beam in the = x 0 plane. It is also observed that the cross-range resolution (along x or y) is better than the range resolution (along z), which is consistent with the fact that the bandwidth of the acquired data (which determines the range resolution) is smaller than the maximum frequency (which determines the cross-range resolution). 6.1.3. Power maps of the F-shaped dielectric bar. Figure 5 shows the normalized magnitude plots (in dB) of the power maps of the OUT at 3, 7, 11, and 15 GHz. These images are only qualitative. It is observed that the cross-range resolution (along x and y) improves as the frequency increases, which is expected. Due to the different beamwidths in the xz and yz antenna radiation patterns, the resolution performance along x and y is somewhat different; i.e., the branches of the F-shape along x have somewhat sharper contrast compared with that along y. It is also observed that the dependence of the range resolution on frequency is not significant. This is expected as in planar scanning it is the bandwidth which impacts the range resolution.  Figure 6 shows the initial guess for the real permittivity (ε′) distribution at 3, 7, 11 and 15 GHz. Figure 7 shows the initial guess for the imaginary permittivity (ε″) at the same frequencies. As mentioned, the F-shaped bar in the OUT has a complex relative permittivity of + 1.2 i0 with a property contrast in real permittivity only with respect to the RO (air). This difference between the real and the imaginary permittivity distributions in the OUT is already well captured in the initial guess shown. However, the initial permittivity estimate has substantial errors and can only serve as a good starting point for the LS solution of (26). Figure 8 shows the LS solutions for the real permittivity distribution of the OUT at 3, 7, 11 and 15 GHz. Figure 9 shows the respective imaginary permittivity distribution. The LS solutions use the LS solver lsqlin of MATLAB [35] with physical constraints on the real and imaginary permittivities: ε′ ⩾ 1 and ε″ ⩾ 0. Compared with the initial guess (see figures 6 and 7), the estimation accuracy has been improved dramatically. Figure 10 shows the results obtained with the multi-frequency combination in section 5.3. It is worth noting that these The RRMSE defined in (34) for the LS solution is given in figure 11 (see the red curve with the scale on the right). No additional noise is added at this stage. It is worth mentioning that the raw simulation data are not purely noiseless; it contains numerical errors at a very low level. The average mesh convergence error for the S-parameters in the simulations is around 0.02 (added uncertainty of ±0.01 for all reflection and transmission coefficients) over the entire frequency band. It is evident that the RRMSE of the LS solutions over the whole frequency band is well below 0.015 and orders of magnitude better than the initial guess. Generally, the accuracy improves as the frequency increases.
The choice of the solver and the initial point are critical in dealing with the nonuniqueness of the solution. As an example, figure 12 shows the real permittivity solutions of the direct (or linear) Born inversion using the Moore-Penrose pseudoinverse of the equations defined by (20). Note that the pseudoinverse solver (the MATLAB function pinv) allows neither for an initial point setting nor for constraints. The F-shaped target is detected well along the cross-range but not along the range. The pseudoinverse solution can be applied to the equations defined by (26) as well. The result is comparable to the one shown in figure 12 (not shown here for brevity) and much worse than the iterative LS solution where we employed the initial guess (32) along with physically based constraints on the complex permittivity values.
Finally, we note that the direct solution of either (20) or (26) with multi-frequency wideband data is impractical when the number of frequency samples is on the order of hundreds or thousands. In this case, the intermediate step of constructing multi-frequency power maps is critical as it provides computational speed and an adequate initial point for the linear LS solution. with the MATLAB function awgn [35]. The signal-to-noise ratio (SNR) is set to 60, 40, 20, 10, 3 and 0 dB, respectively. Figure 11 shows the RRMSE of the LS solutions obtained using simulated data with the added GWN. As the SNR decreases, the RRMSE increases significantly. Figure 13 gives an ( ) approaches infinity. In this situation, the addition of noise improves the conditioning by introducing differences in the matrix element; however, this does not improve the image quality since the original data have been corrupted by noise. When the SNR is poor, the image quality might be improved by using regularization techniques, e.g., Tikhonov regularization [36]. Figure 15 shows the LS solutions of the real permittivity obtained using Tikhonov regularization with the regularization parameter of 0.01 when SNR = 20 dB. Compared with figure 13(c), the F-shaped target is better captured in figure 15(c), although the accuracy is still poor.  6.2. Reconstruction of a lossy RO with two dielectric cylinders using experimental data 6.2.1. Experimental setup. The experiment is performed with a planar raster-scanning system shown in figure 16, which is far more complex than the synthetic simulated example. The Rx is equipped with an electronically switched antenna array of nine bow-tie elements. The array has five bow-tie elements polarized along x and four polarized along y. The Tx antenna is a TEM (transverse electromagnetic) horn [37] polarized long x. The distance between the top and the bottom scanning planes is about 6 cm. Both the Rx and Tx antennas are at a distance of less than 2 mm from the objects surface. Together with the imaged object, the Tx and Rx antennas form a 10-port microwave network allowing for the measurement of nine transmission S-parameters at each sampling location. Measurements are also performed with the imaged object flipped over, which is equivalent to switching the Tx and Rx locations. Thus, the total number of acquired signals at each sample point is N s = 18. The frequency sweep employs N f = 86 equally spaced sampling points from 3 to 9 GHz. The RO is made of five absorber sheets 2 , each of which has the dimensions of 20 × 20 × 1 cm, as shown in figure 17(a).
The CO is the same as the RO except that a dielectric cylinder (see footnote 2) is inserted in the center of the middle layer of the RO, as shown in figure 17(b). The dielectric cylinder has a diameter of 1 cm and a height of 1 cm. The band-average complex relative permittivities of the absorber sheets and the dielectric cylinders are − 10 i5 and − 15 i0.003, respectively. The dispersion of these materials is negligible.
The OUT is the same as the RO except that two dielectric cylinders are inserted in the middle layer of the RO, as shown in figure 17(c). The cylindrical inclusions have the same size and properties as those of the scatterer in the CO. They are placed next to each other at an edge-to-edge distance of 1 cm, which is roughly equal to the anticipated spatial resolution. The imaged area is a square of 5.5 × 5.5 cm; see figure 17(c).
The experimental signals contain noise caused by many factors, e.g., mechanical vibrations and uncertainties in the positioning, electronic noise of the instruments (the vector- network analyzer), electromagnetic interference, other wireless signals, etc. Figure 18 shows the averaged SNR for the 9 channels of the scanning system. The SNR is calculated from the RO measurement data in resemblance with the way the MATLAB function awgn adds noise (see the appendix). The averaged system SNR varies roughly from 2 to 13 dB over the frequency band, which implies very low signal quality. This SNR is typical for the imaging in a lossy environment-a scenario which arises in applications such as tissue imaging, nondestructive testing and underground surveillance.
In this example, 3D images cannot be generated due to limitations of the hardware, which is incapable of acquiring the reflected signals. With planar raster scanning, the lack of reflected signals leads to loss of range resolution; see [15,16,33,38,39]. Here, only the 2D results for the middle layer of the OUT are presented.
6.2.2. Power maps of the lossy RO with two dielectric cylinders. Figure 19 shows the normalized magnitude plots (in dB) of the power maps for the middle layer of the OUT at different frequencies. The two dielectric cylinders embedded in the middle layer are clearly detected. As the frequency increases, the cross-range resolution improves significantly.
It has been attempted to reconstruct the OUT by using a simulated resolvent kernel, computed from the incident field distributions [21,22] associated with each of the antennas in the simulation model of the experiment. This attempt failed, i.e., the obtained power maps did not show the targets. This failure is due to the inability of the model to represent adequately the complexities of the experiment despite the fact that the computations are extremely timedemanding. 6.2.3. Dielectric-property estimation of the lossy RO with two dielectric cylinders. The LS problem is solved using Tikhonov regularization with the regularization parameter of 0.001. Figure 20 shows the LS solutions for the real permittivity distribution in the middle layer of the OUT. Figure 21 shows the plots of the imaginary permittivity distributions. It is worth mentioning that, without regularization, the LS solution fails to localize and quantify the dielectric cylinders at low frequencies, e.g., 3.56 and 5.03 GHz. With regularization, acceptable estimates are obtained at more sample frequencies. At some frequencies, the two   cylinders are not only well localized but also their estimated permittivity is closer to the true value of − 15 i0.003; e.g., see figure 21(c) as well as figure 22(c). In both cases, the estimated permittivities of the dielectric cylinders and the background medium are about − 12.5 i2 and − 10 i5; this is impressive given the poor SNR level.
It is instructive to compare the imaging results of this experiment (with SNR levels roughly between 2 and 13 dB) with those in the simulation example in section 6.1 where the simulated data have been corrupted by similar levels of noise (see figures 11 and 13). As is evident from figure 13(d), at SNR = 10 dB, the reconstruction in the simulation example fails completely-neither the location nor the permittivity of the object are recovered. Yet, in the experiment, which has similar or lower SNR, the targets are not only localized but also a fairly good estimate of the complex permittivity is obtained without regularization or constraints employed in the LS solution. This is mainly due to the fact that in the experiment we use nine Rx sensors with one Tx antenna allowing for the acquisition of 18 transmission coefficients at each frequency and sampling point (on both sides of the RO). In contrast, the simulation setup employs two dipoles, each acting as both a Tx and an Rx antenna. This implies only three independent responses (taking into account reciprocity) at each sampling position and frequency. Increasing the number of independent responses is a well-known method to counteract the detrimental effects of noise; see e.g., [21,22]. Also, the high losses in the RO help reduce the modeling errors associated with multiple scattering. Figure 22 shows the images generated by further applying the proposed multi-frequency combination. Given the true properties of the dielectric targets and the background medium, the overall-frequency results are obviously better than that of some single-frequencies; see figures 22(a) and 20(b). Finally, figure 23 shows the condition number of the LS problem with the experimental data, which appears to vary with frequency in accordance with the averaged SNR curve in figure 18. This illustrates one more time the impact of noise on the uniqueness and the well-posedness of the LS problem.

Discussion and conclusions
A fast imaging method is proposed, which achieves quantitative reconstruction through the direct inversion of the scattering problem in real space (as opposed to Fourier or k space). The method does not require analytical or numerical approximations of the forward model since the acquisition system is fully characterized through the calibration measurements of two known objects: the RO and the CO.
First, the baseline measurement is that of the RO, i.e., the assumed known background medium together with the acquisition apparatus. It is essential in estimating the incident field and, therefore, in extracting the scattered signal component from the total signal. We emphasize that the main goal of the RO measurement is to describe properly all the features of the particular acquisition system (e.g., transmitting and receiving antennas, the fixtures and shielding) rather than possible complexities of the OUT. The OUT features are the subject of the reconstruction, not the system calibration. The construction of a proper physical RO depends on the expected averaged OUT permittivity distribution. The closer the RO properties are to those of the OUT, the more successful the reconstruction would be as discussed in section 2.4. If no additional prior knowledge is available about the expected OUTs, the RO is fabricated as a homogeneous object, which is also done here. Second, the CO measurement takes the above standard calibration procedure only one step further by employing a small known object in the already available RO. As in the RO measurement, the CO measurement aims at describing properly the PSFs of the particular sensors rather than attempting to mimic the OUT. The acquired CO signals are critically important for the proposed quantitative inversion as they comprise the functional basis for quantifying the OUT responses.
The direct inversion without the need for Fourier transforms enables real-time reconstruction. Acquiring the system's resolvent kernel experimentally, instead of using analytical  or simulation models, improves drastically the image quality, especially in the case of nearfield imaging where these models become inadequate. The proposed method is not subject to the typical limitation associated with the linear Born or the Rytov approximations. However, its accuracy is shown to improve if the dielectric properties of the small known target in the CO are similar to those in the OUT. In addition, the method does not account for multiple scattering. These limitations cannot be ignored in the imaging of heterogeneous objects with wide contrast variations. In such cases, the method can provide a good initial estimate of the object's permittivity distribution toward subsequent refinement through a nonlinear iterative solution.

Appendix. Calculation of the SNR of the measurement data
The SNR of the measured data are calculated from the data acquired with the RO. Assume that at the mth sample frequency there is a total of N acquired signals S ,