Spatial Sampling in Monostatic Subsurface Radar Imaging

—This paper deals with microwave subsurface imag- ing obtained by a migration-like inversion scheme, for a 2 D monostatic scalar conﬁguration and a two-layered background medium. The focus is on the determination of a data sampling strategy which allows to reduce the number of required mea- surements and at the same time keep the same performance in the reconstructions. To this end, the measurement points are determined in order to approximate the point-spread function corresponding to the ideal continuous case (i.e., the case in which the data space is not sampled at all). Basically, thanks to suitable variable transformations the point-spread functions is recast as a Fourier-like operator and this provides insight to devise the sampling scheme. It is shown that resulting measurement spatial positions are non-uniformly arranged across the measurement domain and their number can be much lower than the one provided by some literature standard sampling criteria. The study also contains a comparison with the free-space case so as to highlight the role played by the half-space that schematized the subsurface scattering scenario on the number and the locations of the measurement points. Numerical examples are also included to check the theoretical arguments.


I. INTRODUCTION
In this paper we consider microwave imaging achieved by inverting the linearized scattering integral equation for subsurface prospecting. Accordingly, a reflection mode configuration is addressed. Moreover, data are assumed to be collected under a monostatic strategy [1].
Regardless of the inversion method one may choose to employ, the scattering problem must be first translated into its finite dimensional counterpart through a suitable discretization process. This step is crucial and must be designed in order to trade-off different aspects of the problem. On the one hand, to keep the number of measurements (especially in terms of the spatial locations) as low as possible which favourably affects the data acquisition time and the computational and the storage burden. On the other hand, the number and the way data are sampled impacts on the eigenspectrum of the matrix model arising from discretization and hence on the achievable performance and the resilience against noise and uncertainties [2].
The discretization step can be seen as a sensors' selection problem [3] which can be roughly stated as follows: given an N × M matrix model to be inverted, where N is the number of measurements and M the number of unknowns, select L < N so that the achievable performance is optimized. The selection problem presents a combinatorial complexity M. Maisto, R. Pierri and R. Solimene are with the Diparitmento di Ingegneria, Università degli Studi della Campania -Luigi Vanvitelli, 81031 Aversa, Italy, e-mail: mariaantonia.maisto@unicampania.it. and hence cannot be in practice addressed by an exhaustive exploration in the data space. To overcome this drawback a number of methods basically based on convex optimization, greedy methods and heuristics have been developed in order to optmize the selection against some error metrics, such as the confidence ellipsoid, the mean square error [4], the frame potential [5], etc. [6], [7], which are all linked to the singular values of the problem and try to shape their behavior [8]. All these methods still require iterative procedures and have the merit to have a clear connection to the error metric. However, they usually employ a fine sampling grid (i.e., large N ) and an a priori fixed size, L, of the data to be retained. The latter matter is of course connected to the Number of Degrees of Freedom of the scattered field which can be used to estimate L [9], [10] and can also provide insight for sampling the field data [11].
A more classical way to determine the data set is to employ the properties of the kernel function of the scattering integral equation and related bandwidth arguments to devise samplinglike schemes. In this case, no iterative procedures are required and the number of data to be employed in the reconstructions arise as the sampling points following within the measurement aperture and the adopted frequency band. For example, planewave expansion of the Green function, once evanescent waves are neglected, leads to a λ/4 spatial sampling step, λ being the wavelength. Stationary phase arguments are employed in [12]- [15]; the resulting spatial sampling step depends on the extent of the spatial region to be imaged and is lower than λ/4 so as to save many data points.
In view of the mentioned advantages provided by the sampling approach, this is the strategy adopted in this paper to find the sensors' locations. In particular, we progress with respect to the results reported in [12], [13] by refining the theory upon which sampling is based, which is also able to take into account the spatially varying filtering introduced by the propagator in near-field configurations, the latter being the ones relevant for subsurface prospecting. In particular, it is shown that the sampling points must be non-uniformly arranged across the measurement domain according to a law that also depends on the dielectric permittivity of the halfspace which models the subsurface scattering scenario. What is more, the number of required points is lower than the one predicted by previous approaches in [12], [13].
The inverse problem related to subsurface imaging is illposed, or equivalently said, the singular values of the scattering operators tend to zero as their index increases [2]. Therefore, during the reconstruction some regularization scheme must be employed. Considering a finite number of measurements implicitly regularizes the inverse problem; indeed, the sensors' selection methods choose the measurements in order to shape the singular values in order to optimize some error metric. However, when L is large, it may happen that low singular values are retained and regularization is still necessary. This question is less important if the inverse operator is approximated by the adjoint one (i.e., adjoint based inversion). This is the case for migration inversion schemes which basically obtain the reconstructions by back-propagating the scattered field data according to the background Green function [16], [17]. Indeed, by the adjoint inversion, even a large L does not represent a problem since stability against the noise is implicitly obtained by the filtering imposed by the singular values themselves. However, the resulting reconstructions do not optmize some error metrics, for example they do not minimize the mean square error [17]. Accordingly, error metrics cannot be used to check the goodness of the sensors' selection by the sampling approach.
A common way to estimate the achievable performance, even by adjoint methods, is in terms of the point-spread function, which of course is related to the resolution that can be obtained in the reconstructions. Accordingly, in this manuscript, we derive the data sampling scheme in order to approximate the point-spread function corresponding to the continuous case, that is the ideal case in which data are not assumed to be sampled.
The rest of the paper is organized as follows. In Section II, the mathematical formulation of the problem is introduced. Section III is devoted to deriving the sampling scheme which is then checked by numerical examples in Section ??. In this section we consider the effect of the half-space dielectric permittivity. In particular, the case of free-space is also included and used as benchmark to better highlight the particular features due to the subsurface scenario. Finally, conclusions are briefly reported.

II. PROBLEM DESCRIPTION
The study is developed for the 2D scalar scattering problem sketched in Fig. 1 with invariance assumed along the y-axis. The background medium consists of two homogeneous nonmagnetic (i.e., the magnetic permeability is everywhere the one of the free-space µ 0 ) half-spaces separated by a planar interface at z = 0. In particular, the upper half-space is assumed as the free space and its dielectric permittivity is denoted by u = 0 ; the lower one schematizes the subsurface region and is electromagnetically denser with l > 0 .
The unknown targets are embedded in the lower halfspace and assumed to reside within the rectangular region SD = [z min , z max ] × [−X s , X s ] (i.e., the investigation domain). The incident field is radiated by an y-polarized line source of unitary amplitude at different frequencies within the wavenumber band k 0 ∈ Ω = [k 0min , k 0max ], k 0 being the wavenumber if free-space. A monostatic configuration is considered, so that the scattered field is collected at the same position as the source while the latter can assume different positions across the measurement aperture. In particular, the measurement domain is assumed to be the segment OD = [−X 0 , X 0 ] of the x-axis located at the height z o ≥ 0. The Under the Born approximation [18], the contrast function χ and the scattered field collected over Σ = OD × Ω are linked through the following linear scattering operator where L 2 (SD) and L 2 (Σ) represent the set of square integrable functions supported over SD and Σ, respectively. The χ function is defined as ( s − l )/ l , with s being the dielectric permittivity of the unknown scatterer. The explicit form of the operator A is given as where ω is the angular frequency, n = l / u the refractive index and k l = nk 0 the wavenumber in the lower half-space medium. The function G(·) is the Green function pertinent to the two-layered background medium and appears squared because of the considered monostatic configuration.
It is assumed that |z min | > λ l , λ l being the wavelength in the lower half-space. Accordingly, the Green function in (2) can be approximated as [19] G where h(x o , r, k 0 ) takes into account the relevant amplitude factors, and φ(x o , r) = (R u + nR l ) is the phase term that takes into account the propagation path. In particular, r = (x, y) and r o = (x o , z o ) are the target and the field points, and the paths travelled by the waves in the upper and lower half-spaces, respectively, Finally, x m (x o , r) is the refraction point at the half-spaces' interface which is given by the Snell's law as When the observation point is just located at the separation interface then x m = x o and z o = 0 so that the phase term simplifies as φ( Moreover, for the case of a homogeneous free-space background medium, the functions in (3) In order to perform the reconstruction, equation (2) should be inverted for the χ function. In this paper we choose to achieve inversion through the adjoint of the scattering operator, that isχ with This is a very common inversion strategy found in literature where it is addressed as migration, backpropagation, timereversal and so on [20]- [22].
The point-spread function obtained by such an inversion scheme is the so-called model resolution kernel [23] which links the actual unknown to the reconstructed one, that iŝ where The previous formulation applies to the ideal cases where data are continuously collected over the frequency band and the measurement domain. Of course, practical scenarios requires to discretize the problem by properly sampling the spatial and the frequency variables x o and k o . This aim is pursued in the next section where a sampling scheme is devised in order to approximate the point-spread function.

III. SAMPLING SCHEME
In order to derive the sampling scheme, we mainly focus on the spatial integral appearing in the point-spread function expression (8). Hence, let us rewrite this term as with k 0 considered as a parameter. Note that in (9) the dependence on z has been omitted since it has considered z = z min . The latter being the most critical range coordinate for finding the discrete approximation of the pointspread function. Indeed, the scattered field arising from targets located at such a range has the largest spatial bandwidth and hence requires more samples for a good point-spread function discrete approximation. Now, psf k is rewritten in a more convenient way so that Fourier based arguments can be used to gain insights on the spatial sampling step. In particular, the developments and tools presented in [24], [25] are adapted to the problem at hand.
More in detail, as a first step, we introduce a transformation η(x) that maps the interval [−X s , X s ] into [η(−X s ), η(X s )]. At this juncture, such a function does not need to be defined yet and its usefulness will be clear later on.
Then, introduce the following spectral variable dν (10) with η = η(x) and η = η(x ). Accordingly, the phase argument in (9) can be rewritten as The function w(η, η , x o ) is continuous and monotonic with respect to x o (hence invertible) ∀η, η . This allows to replace the integration in x o with the integration in w. In particular, on denoting and w avg (η, η , X 0 ) = w(η, η , −X 0 ) + w(η, η , X 0 ) 2 (13) and by setting w =w + w avg , expression (9) can be rewritten as where K(w, η, η , k 0 ) = −A(x o (w), η, η , k 0 ) dxo dw . Note that, for the sake of notation simplicity, we omitted to explicitly report the variables which Ωw and w avg depend on. However, this dependence is indicated in (12) and (13) and highlights the spatially varying behaviour of the point-spread function which, as opposed to far-field configurations, is a characteristic feature of near-field configurations.
At this point, the transformation η(x) can be chosen in order to simplify the matter a little bit. In particular, it has been shown that by setting [24] η the point-spread function bandwidth Ωw can be made constant and equal to 1 so that the spatially varying behaviour mentioned above results now encoded into the non-linear link between the variables x and η, which clearly implies a spatially varying resolution. Also, by noting that the leading order contribution in (14) occurs for η − η = 0 [26], the amplitude factor is approximated as K(w, η, η , k 0 ) ≈ K(w, η , η , k 0 ) = K(w, η , k 0 ). Accordingly, (14) becomes We assume to collect data over a set of spatial positions so thatw results uniformly sampled and say ∆w = ∆w the corresponding sampling step. Then, Now, standard Fourier arguments, as the one invoked for avoiding greeting lobes in array antenna theory, suggest that the sampling step forw should be where the outer right hand side term arises because η(−X s ) = −η(X s ). By making explicit the relationship betweenw n and the corresponding sampling points in x o , i.e. x on , (18) rewrites as Moreover, since (19) must be satisfied for each η, η ∈ [−η(X s ), η(X s )], it must be imposed that In particular, the maximum on the left side can be easily evaluated from (10) and when inserted in (20) yields so that, eventually, the equation for finding the sampling points In particular, eq. (22) entails the non-uniform sampling of the observation variable x o .

A. Discussion
A number of considerations can be added to complement the previous derivation.
First, as shown in (22), the sampling points are frequency dependent. Therefore, a conservative choice is to set the sampling points in correspondence to the highest adopted frequency; hence, (22) simply becomes Second, the previous derivation only dealt with the sampling of x o . The sampling of the wavenumber can be achieved according to standard arguments based on the range extent of the area to be imaged, that is as ∆k o = π/n(z max − z min ). Basically, the sampling has been tackled by a sort of factorised approach. This is indeed a common way to go through the sampling problem because of its simplicity. Of course, a better (in the sense that a lower number of samples is expected) and more rigours way should be to simultaneously address the sampling in the x o − k 0 domain. With some complications the previous approach could be generalized to address this case as well. However, this would lead to less practical measurement configurations where spatial sampling in general changes with the frequency. As a result, that number of spatial samples could be even larger than the one predicted by previous derivation. Since, we believe that reducing the spatial measurements positions is more important than reducing frequencies, the full 2D sampling approach has been ruled out in this study.
Third, the derived procedure is indeed general and not dependent on the background medium. Actually, the medium enters only through the specif expression of the phase function φ(·). Of course, the sampling step and hence the number of samples vary accordingly to the medium features.
Moreover, it is interesting to compare the number of required spatial samples obtained by the previous derivation to literature results. To this end, we consider the estimation reported in [12], [13]. In particular, since those studies refer to a free-space configuration, here we just pursue the comparison for such a scenario. In particular, it is assumed z o = 0. Also, in those papers two sampling criteria were suggested. One, which besides the different configuration parameters, also depends on the measurement aperture; and another one from which, after signal compression, such a dependence is removed. The second one requires fewer spatial samples.
According to (20), the number of required spatial samples is given as The estimations reported in [13] instead return and after signal compression Of course, N xo > N c xo . Now, using easy math, the following conclusions can be drawn: • For fixed X 0 and X s , N w ≤ N c xo with equality reached when z min increases. More in detail, the two estimations coincide for Fresnel zone configuration; • For fixed z min and increasing X 0 , N w grows more slowly than N c xo and what is more it is always bounded. Eventually, the proposed spatial sampling scheme proves to be more effective, especially for near-zone configurations.
As a final remark, we point out that, while using (23), in some cases a slight oversampling could be required in order to ensure the for each target positions, the point-spread function side-lobes do not start to increase. This is especially true for a small relative bandwidth. For large bandwidths, having set the sampling step in correspondence of k max , guarantees a certain degree of oversampling for most frequencies.

IV. NUMERICAL ANALYSIS
In this section we presents a few numerical examples in order to check the sampling scheme presented above. In particular, the examples focus on the comparison between the 'exact' point-spread function and the one obtained by employing the sampling scheme in (23). More in detail, by exact point-spread function we mean (8) numerically implemented by employing a very fine and uniform grid of points within the data space Σ; the point-spread function that takes into account the proposed sampling strategy is still obtained by implementing (8) but the spatial points are non-uniformly deployed according to (23). In particular, according to the formulation introduced in Section III, in this case the point-spread function is computed as It must be noted that, while in Section III η was highlighted as a function of only x because z was set at z min (i.e., the most critical case for sampling), here, η turns out to be dependent on x and z.
As a final remark before switching to address the numerical examples, we advise the reader that in the following, the approximate Green function introduced in Section II is not exploited. Indeed, the Green function is computed by using its plane-wave expansion (Weyl expansion) [18].
The numerical examples deals basically with three scattering scenarios: the free-space case (scenario 1), the half-space case with measurements taken away from the separation interface (scenario 2) and the half-space case with measurements just collected over the separation interface (scenario 3). In particular, the free-space case is addressed in order to have a reference for the subsurface configuration. Also, it allows an easy vis a vis comparison for the estimation of the spatial samples' number reported in [13], which has been indeed worked out for such a case. Of course, N xo and N c xo can be   generalized to the half-space case. In particular, for the case the scattered field is observed just over the separation interface, the number of samples is computed by using the wavenumber of the lower half-space as the scattering phenomenon was taking place in a homogeneous medium. However, for the case of the stand off configuration we did not pursue the generalization.
We start by showing in Fig. 2 the number and the spatial measurement deployment across the measurement line [−X s , X s ] = [−2, 2]m for the considered scattering scenarios. In particular, panels b) and c) both address the case of scenario 2 but with a different stand off distance from the interface. From such a figure the expected non-uniform sensor arrangement can be appreciated. Moreover, N w results much lower than N c xo (see panels a) and d) for which the comparison has been pursued). Also, it can be noted that the required spatial samples increase while moving from case a) to case d). This indeed could be expected since this is somehow like switching from the free space to a denser medium corresponding to the lower half-space.
In Fig. 3 the normalized point-spread functions correspond-ing to the cases depicted in Fig. 2 are shown. As can be seen, while moving towards the case relative to the scenario 3, the main beam of point-spread function tends to become narrower. This entails an improvement of the achievable resolution. This is a very well known fact and it is consistent with the higher number of samples that are required. A detailed estimation of how the half-space configuration parameters affect resolution has been reported in [27]. Here, what matters is that the exact point-spread function and the one returned by (27) practically overlap. This means that the estimated number of samples works very well. What is more, this holds true even though N w is lower than N c xo .

V. CONCLUSION
In this contribution, we have addressed the problem of devising a sampling scheme in order to collect the data to be used in a migration-like reconstruction scheme for subsurface prospecting.
To this end, thanks to suitable variable transformation, the inversion method point-spread function has been recast as a Fourier type integral operator and its explicit form used to obtain insight concerning the spatial sampling scheme. Instead, frequencies have been sampled according to a standard criterion.
It has been shown that the resulting spatial sampling points are to be arranged non-uniformly across the measurement aperture. Also, the required number of spatial samples is usually lower than the one dictated by commonly used literature results; this being a remarkable advantage since it allows to use fewer sensors for real aperture radar systems or to reduce the time for data collection in synthetic aperture radar systems. A few numerical examples confirmed the goodness of the proposed sampling scheme. Also, the role of the dielectric permittivity of the lower half-space schematizing the subsurface scattering scenario has been highlighted. As expected, compared to the homogeneous free-space case, more spatial points are needed and they must be placed in general on different positions, depending on the medium dielectric features.
With reference to the sampling there is, of course, the question of the achievable resolution. Indeed, configurations that allow for better resolution (i.e., point-spread function having narrower main beam) require a denser sampling. This has been clearly verified in the reported numerical examples. Nonetheless, we did not thoroughly addressed this point in the paper since the focus was on the sampling scheme. A detailed study of the achievable resolution and related closedform estimation for the half-space configuration considered herein has been reported in [27].