Accurate Range Migration for Fast Quantitative Fourier-Based Image Reconstruction With Monostatic Radar

Range migration (or range focusing) techniques are widely used in optical, acoustic, and microwave real-time image reconstruction methods. They have been successfully applied to far-field 3-D imaging where they rely on plane-wave assumptions, which ignore the data amplitude variation over the acquisition aperture. The accuracy of the plane-wave assumption, however, quickly degrades in close-range imaging, where amplitude variations are significant and where the range to the target is on the order of the range sampling step. Here, we present a range-focusing method of improved accuracy, which is applicable to both far-zone and close-range monostatic radar. It refocuses the system point-spread function (PSF) to any range location, taking into account both magnitude and phase changes. The approach can be applied with any Fourier-based imaging algorithm utilizing the Lippmann–Schwinger equation as the underlying scattering model. Here, it is validated through examples based on simulated and measured data where the images are reconstructed with quantitative microwave holography (QMH). QMH employs measured PSFs to achieve quantitative imaging in real-time. The proposed range-migration method is applicable with measured PSFs, too, leading to reduced system-calibration effort and the ability to focus an image at any desired range.


I. INTRODUCTION
F OURIER-BASED image reconstruction is used in numerous applications of wave scattering, including groundpenetrating radar, nondestructive testing, biomedical imaging, concealed weapon detection, target localization, and tracking [1]- [11]. By acquiring the wavefront's magnitude and phase across a finite surface (the acquisition aperture) and across a range of frequencies, target localization, and image generation is possible in a 3-D volume [12]. The benefit of performing the reconstruction in Fourier (or wavenumber) space is the remarkable computational speed. Since the inverse problem is cast as a deconvolution, solving it in Fourier space leads to drastic reduction of the computational complexity compared to the real-space solution [13]. These inversion methods solve a Lippmann-Schwinger type of equation, which is an integral equation of acoustic or electromagnetic (EM) scattering employing a linearizing approximation of its resolvent kernel [14]. A common challenge among the Fourier-based imaging algorithms is how to resolve objects along the range (or depth) dimension (see Fig. 1). Range migration (or range focusing) techniques have been previously developed for varying scenarios [15]- [19]. One common method in synthetic aperture radar (SAR) exploits the Stolt mapping [20]- [22], which maps the data frequency (ω) dependence into k z dependence of the form ∼ exp (−ik z z) (range migration), where k z is the Fourier variable corresponding to the range (z) realspace variable and i = √ −1. This allows for target reflectivity recovery via 3-D inverse Fourier transform (FT). Range stacking [23]- [26] avoids the somewhat computationally expensive Stolt mapping but it, too, needs the range migration factor to represent analytically the z-dependence of the scattered field.
Range migration approaches depend on the resolvent kernel of the chosen forward model of scattering. This kernel is proportional to the product of the total internal field (due to the transmitting (Tx) antenna) and the background Green's function [13]. Since range migration is common in fast real-time image reconstruction, it is usually applied with linearized scattering models where the total internal field is replaced by the incident field due to the Tx antenna as per the zeroth-order Born approximation. On the other hand, Green's function is equivalent to the incident field due to the receiving (Rx) This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ antenna if this antenna were to transmit [13], [27]. In far-field imaging, both of these fields are approximated as plane waves, leading to the simple range-migration factor of ∼ exp (−ik z z) in Fourier space; see [28].
An improvement in accuracy is achieved by an assumption that, while the transmitted field is a plane wave, the scattering is a superposition of spherical waves emanating from the scattering centers building the object under test (OUT); see [17], [26]. By reciprocity, this also means that Green's function is in the form of a spherical wave. The improvement is due to the ability of the spherical-wave model to account better for the amplitude variations of the measured data over the acquisition aperture, especially in close-range applications, e.g., the whole-body imaging for concealed weapon detection.
The data amplitude dependence on the range to target is the most pronounced in near-field imaging where the extent of the acquisition aperture is much larger than the range to the target. In order to make Fourier-based real-time 3-D imaging applicable to these scenarios, it is imperative to devise a range-migration strategy that takes into account the spherical spread of both the transmitted field and Green's function. Here, we propose an accurate range-migration formula, which is applicable directly in the Fourier-transformed data space. The result is analytical and thus it does not increase the running time.
The proposed range-migration approach is applicable not only with the conventional qualitative Fourier-based image reconstruction. Recently, Fourier-based quantitative imagereconstruction methods have been proposed [9], [29]- [31]. The quantitative estimation of the dielectric profile of an OUT becomes possible with the use of measured system point-spread functions (PSFs) in place of analytical/simulated resolvent kernels. The PSF is the measured response to an electrically small scattering probe (SP) and its acquisition is part of the system calibration. The importance of the measured PSFs as a means of calibrating the image reconstruction cannot be overstated. The PSFs represent the measured impulse response function of the imaging system. As such, they avoid the significant inaccuracies in analytical approximations or the modeling errors in simulations of the imaging setup. They capture the characteristics of the employed antennas along with the influence of the entire measurement environment. They also capture quantitatively the system sensitivity to the contrast of a point scatterer. However, the use of a measured PSF is not without its drawbacks. One important drawback is the significant number of PSF measurements required for 3-D imaging with dense sampling along range. Since each range slice in the image requires a PSF, the calibration effort may render the approach impractical. Thus, a method to analytically refocus a measured PSF (with the probe at one range) to any other range is desired.
The contribution of this work is the derivation of an accurate analytical representation of the range dependence of the resolvent kernel of the monostatic-radar scattering model in Fourier space. This representation can be employed directly in the qualitative Fourier-based image-reconstruction algorithms for improved accuracy without increasing the computational requirements. In quantitative Fourier-based imaging employing measured PSFs, the representation leads to a magnitude-and-phase scaling factor that is applied to the Fourier-transformed PSF in order to refocus it to any desired range-slice location. This range-migration technique reduces the calibration effort to a single PSF measurement, which is critically important for practical 3-D imaging applications. A comparison with the conventional far-field range-migration approach confirms the improvement in the re-focused PSF accuracy and the final image structural and quantitative accuracy. We also highlight the limitations of the method related to the initial PSF measurement position and give recommendations for Fourier-domain filtering based on the maximum viewing angle at the refocused SP range location. The method is validated with both simulated and experimental data used in image reconstruction with quantitative microwave holography (QMH) [9], [29], [32].

II. RANGE MIGRATION FOR MONOSTATIC RADAR
In the following derivations, we employ the stationary phase approximation (SPA), which is a method for approximating integrals in the form [33] ∛ Note that the approximation becomes more accurate as p − → ∛, but notations are simplified by setting p = 1 [33].

A. Range Migration for Monostatic Scattering
Consider an SP at the center of the imaged volume (x , y , z ) = (0, 0, 0) and a monostatic measurement of the respective PSF, where the Tx and Rx antennas are co-located and scanning together on an acquisition plane atz = const (z > 0) along x and y. In order to employ the well-known Green's function of the scalar Helmholtz equation, we make the following assumptions about the measured PSF data: 1) it captures correctly the signal strength at the rangez as determined by the system transmitter; 2) it also captures the antennas' polarization and directivity properties, which are not expected to vary significantly when the range position of the point scatterer changes within the image z extent; and 3) the evanescent field components are negligible. The remaining amplitude and phase dependence S(x, y) of the PSF on the distance between the Tx/Rx antennas and the SP is then where r = (x 2 + y 2 +z 2 ) 1/2 , k = 2π/λ is the background wavenumber k ∈ R, and λ is the respective wavelength. The 2-D FT of S(x, y) is Taking into account (1) and (2), the phase function for the given k is determined to be and its derivatives are From (8) and (9), the stationary phase positions are obtained where r 0 = (x 2 0 + y 2 0 +z 2 ) 1/2 . The phase φ(x 0 , y 0 ) at the stationary position is now obtained as We next introduce the spectral variable k z as Using (10) leads to a relation between z and k z , which is analogous to those for x and y in (10), namely z Note that (10) along with (13) imply that x z Now the stationary phase is cast as a function of k z andz, which, in the case ofz > 0, is It is now clear that if k z is imaginary, i.e., k 2 x + k 2 y > (2k) 2 , the phase function becomes imaginary and violates the SPA requirement that the phase function is real-valued. An imaginary k z corresponds to modes evanescent along z and, indeed, such modes have been dismissed in the assumptions leading to the use of the monostatic resolvent kernel in (5). Hereafter, only real-valued k z are considered, so that (15) is written as where To determine the value of the constant K in (4), we need the determinant of (x 0 , y 0 ) in (3), which is found to be Furthermore, the trace tr is found as It then follows from (4) that K = −i, and the FT of S(x, y, k) for a probe at (0, 0, 0) is found as Here, the last argument in (k x , k y , k; 0) emphasizes the probe's range position. The expression in (20) is the FT of the resolvent kernel in (5). As such it can be used directly in qualitative Fourier-based reconstruction algorithms (e.g., microwave holography and range stacking) bearing in mind thatz represents the range from the probe to the acquisition aperture. Note that, unlike previous range migration formulas, where S ∼ e −ik zz or S ∼ e −ik zz /k z , the expression in (20) also containsz in the denominator. This is a consequence of the kernel's amplitude dependence ∼ 1/r 2 in real space. Consider now the case where the PSF has been measured with the probe at (x , y , z ) = (0, 0, 0) but the PSF for a probe at (0, 0, z ) is desired. Now, the range between the acquisition plane and the SP isz − z . Similar to the case of z = 0, we setz − z > 0, i.e., the probe remains on the same side of the acquisition plane. The FT of the corresponding kernel is Dividing (21) with (20) provides the range-migration scaling factor that can be applied to refocus any measured PSF to a new range location This is the factor used to refocus measured PSFs, such as those utilized in QMH. The QMH algorithm is summarized in Appendix. Its PSFs, denoted as where z 0 is the aperture-to-probe distance at which the PSF is measured and z is the increase (z > 0) or decrease (z < 0) of this distance when the PSF is refocused. Compare this to the conventional plane-wave range migration technique derived from the angular spectrum representation [33] H It is obvious that if the shift along range is insignificant relative to z 0 , (24) is acceptable. Yet, in close-range imaging, where the imaged volume is sampled at range intervals of a similar order of magnitude to z 0 , the magnitude scaling factor is not negligible. This becomes especially important in quantitative imaging, where the accuracy of the approximation directly impacts the accuracy of the permittivity estimates.

B. Analytical Example With Finite Apertures
This example validates the result in (23) and illustrates the inability of the method to account for the case of imaginary k z . It is also used to study the impact of the aperture size. The FT in (6) is first computed via the discrete Fourier transform (DFT) with the frequency set to 3 GHz (λ ≈ 0.1 m). When computing the real-space function (5), the x and y spatial sampling steps are set to 0.025 m. The aperture extent is 12 m along x and y. The probe is positioned 0.2 m away from the aperture. This DFT serves as a reference. The SPA solution (20) is also computed (with all imaginary values of k z set to 0) and compared with the DFT result in Fig. 2.
It is clear that the SPA and the DFT results match well. A circle circumscribes the 4k 2 = k 2 x +k 2 y region (of radius 2k), separating the propagating and evanescent regions in k-space. At this circle (for the DFT case), minor Gibb's effects due to the FT of a finite aperture appears to create rippling artifacts at k x = ±125, where the magnitude plot transitions rapidly [6]. It is especially clear in the DFT case that the contributions in the propagating region are much more dominant in magnitude than the evanescent region. By setting the evanescent region of k-space to zero, a close agreement between the SPA and DFT results is achieved since the attenuation (in the DFT case) is extremely rapid.
In practice, large apertures such as the one considered here (of extent 120λ in both x and y) are rarely implemented in near-field imaging. Consider instead a situation where the aperture is limited to a maximum viewing angle of 45 • (aperture size of 0.4 m × 0.4 m). The results are shown in Fig. 3. The Fourier domain sampling step is tied to the aperture size, and thus the lower k-space sampling step becomes apparent in Fig. 3. 1 Note also that the accuracy of the approximation is reduced, as large Gibb's artifacts are now present and distorting the DFT result.
It is well known that the k-space region of the measured scattered-field propagating modes is also limited by the maximum target viewing angle α [13]. The viewing angle in monostatic measurements is the angle of arrival of the scattered wave relative to the aperture's unit normal. A reception at grazing angles (α ≈ ±90 • ) ensures the availability of propagating modes with wavenumbers k x and k y as high as to fulfill k 2 x + k 2 y ≈ (2k) 2 . However, often such reception is not achievable due to the finite size of the aperture or the limited beamwidth of the antennas. In this case, α < 90 • , and the values of k x and k y are limited by the projection of the wave vector k onto the acquisition aperture, i.e., k α = k sin α. In monostatic radar, this sets the k-space region to a circle of radius 2k α The two circles drawn on Fig. 3(a) and (b) highlight the 2k and 2k α radius circles which contain the propagating wave modes. This example highlights an important aspect of the PSF range migration when the fixed-size aperture does not ensure large viewing angles. In the case where the probe is close to the aperture, the maximum viewing angle approaches 90 • , and thus 2k α ≈ 2k. However, as the probe is shifted further from the aperture, the maximum viewing angle decreases, and thus the region containing the propagating modes shrinks. This observation is important in close-range imaging with measured PSFs. When a measured PSF is refocused from a further to a closer position (relative to the aperture), due to an increasing viewing angle, information beyond the 2k α boundary is required but is not available. This suggests that the initial probe location, at which the system PSF is measured, should be selected to maximize the viewing angle (the closest range location) and thus maximize the propagation-mode information within the 2k-radius region.
To summarize, the SPA approximation agrees well with the DFT result as long as filters are applied to restrict the SPA result to a region which corresponds to 2k α . Note that Fig. 3 appears to show disagreement, especially in the magnitude response [see Fig. 3(a) and (e)], but these rippling artifacts are due to the FT itself. These Gibb's effects (ripple artifacts) are always present with limited apertures since the signal does not decay to zero at the aperture edges. Since QMH operates in the Fourier domain, these artifacts must be suppressed, which is achieved with apodization filtering. Apodization acts as a smooth spatial window, reducing the signal at the aperture edges to zero [32].

A. Simulation
An imaging experiment is simulated using the FEKO full-wave EM simulator [34]. Five half-lambda dipole antennas are positioned 10 mm apart in an "X" configuration [see Fig. 4(a)], and perform a raster scan across a 60 mm × 60 mm aperture in 2 mm increments. A frequency sweep is performed from 25 to 40 GHz in 1 GHz increments, capturing five sets of monostatic measurements. The OUT comprises three structures, which can be seen in Fig. 4. Their permittivities are provided in the caption of Fig. 4. The background medium is vacuum.
To perform the reconstruction, the range migration factor in (23) is applied to the closest PSF measurement (probe at z = 5 mm) in order to refocus it to the slices at z = 15 mm, After the range migration of the PSF, an eighth-order Butterworth low-pass filter is applied in the Fourier domain to restrict the spectrum to the region accessible within the viewing angle of the aperture (25). The system of equations is then inverted to extract the Fourier-domain reflectivity function at each desired range slice. Before casting this result back to real space via 2-D inverse FT, a low-pass 16th order Butterworth filter is applied at each range slice. The filter has a 3 dB value at approximately 1.6k α whereas its 20 dB cutoff is at the 2k α limit for the highest frequency in the measurement as per (25). The suppression beyond 2k α is critically important. The fast FT of the data and the PSFs acquired on a finite aperture produces high-spatial-frequency artifacts. As a result, the Fourier-domain inversion produces erroneous reflectivity values at k > 2k α , where signal levels are low relative to the artifacts. Moreover, these erroneous reflectivity values may be large due to ill-conditioned system matrices. Thus, suppression by more than 20 dB is desirable. The choice of filters is discussed in more detail in [32].
The resultant reconstructions can be seen in Fig. 5. The quantitative estimates of both the large central object and the small cubicles are reasonably accurate when using the four measured PSFs for the four range slices The two closely spaced probes ( r = 1.5) in the z = 5 mm plane are barely distinguishable from each other. The farfield monostatic-radar resolution limit with α max = 80.54 • is calculated as [13] ξ = λ min 4 sin α max ≈ 2 mm, ξ ≡ x, y which is exactly the size of these probes and the sampling step in the acquisition. Their estimated real permittivity is 1.345 in the range migration approach and 1.33 in the measured PSF approach (see Fig. 6). The larger object at z = 15 mm is structurally identified in both approaches, and has an average real permittivity of 1.361 in the range migration approach and 1.327 in the measured PSF approach. The furthest components at z = 35 mm highlight the impact of the aperture limit on resolution. At this range location, the two cubes in corner-tocorner contact appear as one larger cube of higher permittivity, and the isolated cube appears as a lower permittivity, larger cube structure.

B. Experiment
To further validate the approach, a planar raster scanning chamber is configured with a waveguide rectangular (WR)-28 horn antenna [35] with a beamwidth of approximately 54 • . The antenna is connected to an Agilent E8363B vector network analyzer capturing S 11 data from 26 to 40 GHz in 100 MHz increments. While the antenna remains fixed, the platform carrying the imaged object shifts in 2 mm increments across a 30 cm by 30 cm aperture. The antenna is positioned at approximately 4.5 mm above the object. Overall, there are 151 (x) by 151 (y) by 141 (ω) data points. The number of spatial and frequency samples is determined from the extent of the acquisition aperture (determined by the lateral extent of the target), the necessary frequency bandwidth (determined by the desired range resolution), the frequency sampling step (determined by the maximum range to a target in the imaged volume), and the spatial sampling step (determined by the anticipated lateral resolution as dependent on the shortest wavelength) [13], [32].
The imaged object is a stack of four 12.7 mm (0.5 in) thick and 500 mm wide Styrofoam sheets with an average estimated permittivity of r ≈ 1.175 − i0 across the frequency of interest (see Fig. 5). The size of the Styrofoam sheet was selected to extend beyond the beamwidth of the antenna which minimizes reflections caused by the air interface at the edge of the sheets. Several items are placed at different locations in the image volume. Two crosses, one made of carbon rubber and the other ceramic, are positioned in Layers 1 and 3 of the stack-up (see Fig. 7). Their permittivity values can be found in Table I. Note that the QMH algorithm generates an averaged permittivity estimate across the frequency band. Accordingly, Table I lists permittivity values that are averaged over all frequencies. These values are determined either via measurement with a dielectric probe or taken from [36]. Note that the measured permittivity of the components that have been starred (*) is only an estimate, likely lower than the actual  value. The slim-form dielectric probe [37] provides accurate measurements only in the case of mechanically penetrable materials such as liquids and gels where immersion of at least 5 mm is achievable. To achieve immersion in the ceramic or carbon-rubber material, a hole is drilled where the probe is inserted. The approach, however, is prone to errors due to the presence of air gaps. Also, the manufacturers of these materials do not provide permittivity data at the particular frequencies of interest. In Layers 2 and 4, a series of small targets (2.4 mm diameter Nylon balls and 3 mm × 2.6 × mm 2 mm carbon rubber prisms) are positioned laterally in pairs, each separated (edgeto-edge) by either 8, 4, or 2 mm. By doing so, we can evaluate the image resolution of the system relative to the shortest wavelength (7.5 mm at 40 GHz).
A series of PSF measurements are also performed on COs, which consist of a small SP embedded in a Styrofoam sheet. One of the carbon rubber prisms serves as a SP (3 mm × 2.6 mm × 2 mm, r ≈ 7.85 − i3.01). The probe is centered * These properties are underestimated due to the inability to properly measure them with the available dielectric measurement probe [37]. Documentation on the dielectric properties of these materials is not available above 20 GHz. laterally and is positioned at the top-most portion of each layer formed by a Styrofoam sheet. Note that this probe position is offset along range with respect to the cross-shaped components, which are vertically centered within the Styrofoam sheet.
To validate the proposed PSF range migration approach, the PSF measured at the plane closest to the antenna is refocused to all other desired range slices for image reconstruction. For comparison, this is repeated for the standard plane-wave range translation. The results are summarized in Figs. 9 and 10.
In the results using only measured PSFs (see Fig. 9), the images have good structural accuracy, i.e., the item's position and shape are reconstructed well. "Bleeding" artifacts along range occur throughout the image, due to both the large size and permittivity contrast of the cross-shaped components, as well as the misalignment with the PSF probe positioning in the CO measurements. The quantitative accuracy is low, especially in the regions occupied by the two cross-shaped components, which are not only high-contrast but also are not electrically small. The top-most carbon-rubber cross is reconstructed with some nonphysical (positive) imaginary permittivity values whereas the reconstructed ceramic cross in Layer 3 has a negative real part of the permittivity. Such errors are expected because of the linearizing approximation of the forward model [13]. Another potential contributing factor is the misalignment between the measured PSF range position and the range position of the cross-shaped objects. The measured PSFs are aligned with the top of each sheet, while the crosses occupy the entirety of the sheet thickness. Thus, repositioning analytically the PSFs closer to the center of the sheet may improve the permittivity estimate. The smaller probes (Nylon balls in Layer 2 and thin rubber prisms in Layer 4) are reconstructed at accurate range locations with reasonable permittivity estimates. Now consider the results generated using the range migration algorithm with a single measured PSF (z = 4.5 mm), shown in Fig. 10. Since the range migration allows for the PSF to be migrated to any position, slight adjustments are made to the PSF for Layer 3 to align it with the ceramic cross item. A 1.6 mm increase is applied to the range of the migrated PSF position. This increase was empirically determined to reduce the nonphysical permittivity estimates. Indeed, this analytical range tuning focuses the image, resulting in improved permittivity and shape representation of the actual object.
While the bleeding artifacts of the crosses are still present, the cross property estimation has improved and the nonphysical permittivity values have decreased. Interestingly, this result appears to suggest that the ceramic cross has significant loss, which has not been found in the material permittivity estimated through the dielectric probe measurement. Note also that the permittivity of both cross-shaped components appears to change across the surface of the object, which may be due to the different leveling with respect to the aperture plane.
The comparison of the images in Figs. 9 and 10 suggests that the spatial resolution is not impacted by the replacement of the measured PSFs with migrated ones. The smallest probe spacing is 2 mm in Layers 2 and 4, which was chosen to be close to λ min /4 ≈ 1.874 mm, where λ min is the shortest wavelength (at 40 GHz). Using (26), the viewing-angle limited resolution can also be determined. Here, α is limited by the beamwidth of the horn antenna which is ≈ 54 • [35], which would imply a resolution of 2.47 mm. Thus, the difficulty in distinguishing the close proximity probes in Layers 2 and 4 is expected. The real-permittivity slice image in Fig. 11 shows an enlarged section where the small probes reside in the layer z = 42.6 mm. It is observed that the image resolution is worse Fig. 10. QMH reconstructions of the simulated OUT shown in Fig. 4 using the combined Born/Rytov approach [9]. To reconstruct the images, a single PSF measurement is performed with the SP at z = 5 mm. The measured PSF is then migrated to z = 4.5, 16.0, 32.0, 42.6 mm. Each row shows plots of the real and imaginary parts of the OUT relative permittivity at a range slice: from z = 4.5 mm (top) to z = 42.6 mm (bottom). Fig. 11. Enlarged real-part of the permittivity of the z = 42.6 mm layer from Fig. 10, focusing on the array of carbon rubber probes. The probes that are the hardest to individually identify are the probes separated by ≈ λ min /4. The probes separated by ≈ λ min and ≈ λ min /2 are distinguishable. than λ min /4. However, the probes separated by distances of λ min /2 and λ min are distinguished very well.

IV. DISCUSSION AND CONCLUSION
A new accurate formulation of range migration is derived with the SPA for the case of monostatic measurements. The range migration is performed directly in Fourier domain. Compared to previous range-migration strategies, which adjust only the phase of the scattering kernel, this formulation has better accuracy, especially in near-field and close-range imaging scenarios, where the aperture size is much larger than the range to the target and the 3-D range resolution is in the same order of magnitude as the range to the target. The result is analytical and it does not increase the computational requirements of the Fourier-based image reconstruction algorithms.
The impact of the aperture size is also investigated with account for the maximum viewing angle of the imaging setup.
The results indicate that the best strategy for range migration of a measured PSF is to first acquire the PSF with a SP at the closest range position relative to the antenna(s) and then migrate away to all other desired range positions.
Both simulation and experimental examples using the QMH image-reconstruction algorithm demonstrate the image quality improvement when compared to the plane-wave migration approach. Most importantly, the accurate range-migration strategy allows for taking full advantage of measured PSFs with significantly reduced calibration effort. This is due to the reduction of the number of PSF measurements from N z (the number of imaged slices in a 3-D image) to just one.
It is worth emphasizing the importance of using a measured PSF whenever possible. It captures the properties of the used antennas and the lateral field distribution at the range where the SP is. Thus, the reconstruction algorithm (e.g., QMH) does not assume that the antennas act as point sources or pointlike receivers. Importantly, a measured PSF also captures the system response sensitivity to the contrast of a point scatterer thus enabling quantitative reconstruction. The analytical approximation of the scattering kernel behavior through the term ∼ exp(−i2kr )/r 2 comes only in the assumption that it captures the range dependence of the PSF since the measured PSF at a given range already captures the system-specific lateral distribution. In the far zone of the antenna, this approximation is accurate and with the antenna polarization and gain pattern already captured by the measured PSF, the range migrated PSFs are highly accurate. Most close-range imaging systems (e.g., whole-body imagers) operate on targets, which are in the far zone of the antennas but are close to the acquisition aperture (the extent of the aperture is greater than the range distance to the target) so that the spherical spread of the incident and scattered waves is not negligible. Nonetheless, as demonstrated in the simulated-data and the measured-data examples, near-field imaging benefits, too. This is because, in the radiative near zone (Fresnel zone) of the antenna, the field behavior represented by the term ∼ exp(−ikr )/r dominates. Range migration within the reactive near zone of the antennas, however, remains a challenge.
Finally, the limitations with regard to the imaging of electrically large high-contrast targets must be kept in mind. Such targets violate the Born and Rytov approximations that are the basis of the scattered-field data extraction strategies employed by many direct reconstruction algorithms [9], [39]. Additionally, the linearized model of scattering employed by such algorithms does not account for higher-order scattering effects such as multiple scattering and mutual coupling, which can be significant in the presence of large high-contrast objects. This leads to increased image artifacts and low fidelity of the quantitative reconstruction. The replacement of a measured PSF with a migrated PSF is not expected to have any impact, negative or positive, on these limitations.
Future work aims at extending this study to bistatic radar. Also, the efficacy of the proposed method is to be further validated with other quantitative Fourier-based imaging algorithms, e.g., [30]. APPENDIX QUANTITATIVE MICROWAVE HOLOGRAPHY QMH has been extensively derived in previous publications and is only briefly described here for contextual purposes [9], [10], [32].

A. Measurement Setup
The following derivation of QMH applies to planar apertures (along x, y, range along z), though formulations of QMH for cylindrical apertures are also available [40], [41]. QMH relies on the measurements of three separate objects for image reconstruction, the reference object (RO), the CO, and the OUT. The RO is a measurement of the uniform background medium with no embedded scatterers. It provides the incident-field response of the system and enables de-embedding of any system-specific artifacts that may be captured during the measurement. The CO is a measurement of an electrically small SP centered along the lateral scanning positions and at the desired reconstruction range location (0, 0, z ). The CO provides the data that produces the systemspecific PSF. Finally, the OUT is a measurement of the target to be reconstructed.

B. Forward Model of Scattering
The frequency-domain forward model of QMH has been derived with the use of scattering parameter (S-parameter) data instead of the conventional electric field data, enabling direct application to measurement data [27] S sc ik (r xy ,z; ω) ≈ V ρ(r )·H ik (r xy − r xy ,z; z ; ω)dx dy dz .
Here, S sc ik (r xy ,z; ω) is the S-parameter scattered response of the OUT measured by a Rx antenna i and a Tx antenna k at a given lateral position along the imaging aperture r xy = (x, y) and a fixed range positionz, and at a given angular frequency ω. ρ(r ) is the reflectivity function at a particular location within the imaged volume r = (x , y , z ), r ∈ V , where r xy = (x , y ), and H ik (r xy − r xy ,z; z ; ω) is the PSF response with a SP located at the center of the z plane r xy = (0, 0), and is assumed to be translationally invariant across the plane.
The reflectivity function is related to the permittivity of the imaged target via a measured SP [9], [42] ρ(r ) = r (r ) · ( r,SP SP ) −1 (28) where r (r ) = r − r,RO is the complex relative permittivity difference between the OUT and the background medium, r,SP = r,SP − r,RO is the complex relative permittivity difference between the SP permittivity and the background medium, and SP is the SP volume. The forward model in (27) relies on the external scattered field responses, whereas the measurements of the CO and OUT are of total external field responses. These can be extracted via either the subtraction approach of Born's approximation, or the logarithmic extraction process of Rytov's approximation, or a combination of the two [9], [13], [32], [43].
The forward model in (27) is a convolution along x and y in the real-domain, and through a FT, and discretization of the range integrand, becomes a multiplication in the frequency domain [44] S sc ik (κ κ κ,z; ω) ≈ z ρ(κ κ κ, z ) H ik (κ κ κ, z ; ω) v .
Here, tilde denotes the Fourier domain data, κ κ κ = (k x , k y ) is the spectral position, and v is the voxel volume. Equation (29) is a system of equations at each spectral position κ κ κ [44] A(κ κ κ) ρ ρ ρ(κ κ κ) = b(κ κ κ) where Here, N f is the number of frequency samples, and N z is the number of unique PSFs available for image reconstruction. These values are typically on the order of 10-100 and 3-20, respectively. However, the effort involved to increase the number of N z samples is substantially higher that N f , since each unique PSF requires a completely new measurement to be constructed, while N f can be increased by simply modifying the sampling settings on the vector network analyzer.
To reconstruct images, one can solve the inverse problem of (30) through a conventional inversion technique 2 Note that there are a unique system of equations for each κ κ κ, yet these equations are generally small. This makes QMH incredibly fast to solve, and parallelizable.