Experimental Validation of a Microwave Imaging Method for Shallow Buried Target Detection by Under-Sampled Data and a Non-Coperative Source

In microwave imaging it is often of interest to inspect electrically large spatial regions. In these cases, data must be collected over a great deal of measurement points which entails long measurement time and/or costly, and often unfeasible, measurement configurations. In order to counteract such drawbacks, we have recently introduced a microwave imaging algorithm which looks for the scattering targets in terms of equivalent surface currents supported over a given reference plane. While this method is suited to detect shallowly buried targets, it allows to independently process each frequency data, hence the source and the receivers do not need to be synchronized. Moreover, spatial data can be reduced at large extent, without incurring in aliasing artefacts, by properly combining single-frequency reconstructions. In this paper, we validate such an approach by experimental measurements. In particular, the experimental test site consists of a sand box in open air where metallic plate targets are shallowly buried (few cm) under the air/soil interface. The investigated region is illuminated by a fixed transmitting horn antenna whereas the scattered field is collected over a planar measurement aperture at a fixed height from the air-sand interface. The transmitter and the receiver share only the working frequency information. Experimental results confirm the feasibility of the method.


Introduction
Microwave imaging, and more in general Radar imaging, is a mature research field that finds application in a number of different contexts, where pursuing non-destructive investigation is convenient or mandatory [1] - [11].
Ground Penetrating Radar (GPR) is a radar system that is properly conceived to address non-destructive imaging. Generally, GPRs work in contact with the interface between the air and the medium under investigation. However, there is a great interest for achieving target detection through non-contact measurement layouts, for example with GPRs mounted on a flying platform [12] [13]. Indeed, stand-off distance configurations allow for the investigation of regions which are not easily (or safely) accessible, as it happens, for instance, when one has to deal with mine or un-exploded device detection. Moreover, a flying GPR can allow for inspecting large areas quickly.
In this framework, however, the system cost and the achievable performance must be traded-off. This requires finding a compromise between the time needed to collected data, the number of sensors to be simultaneously deployed and the way transmitter and receivers "communicate". In this regard, a single-view/multistatic configuration seems convenient since only one sensor transmits and the others act like mere passive receivers, hence with reduced weight and cost. However, synchronization between TX and RXs is a critical issue when they are mounted on different platforms, since, differently from multi-monostatic arrangement, TX and RXs are no more co-located and hence do not share the same electronic system. Also, the number of measurement points is directly linked to the number of flying platforms and hence must be reduced as much as possible.
Recently, we have introduced a reconstruction scheme that allows to mitigate the previously mentioned drawbacks [14], [15]. This method relies on a certain scattering formulation where, according to the equivalence principle [16], the scattered field is modeled as being radiated by equivalent surface currents which are supported over the air/soil interface, or at some other reference plane whose depth is chosen. Basically, the main idea is that if the reference plane is close to the scattering target then the spatial support of the surface current gives information concerning the transverse location of the target. Here, the crucial point is the choice of the reference plane. If it is a priory known that the targets are in close proximity to the air/soil interface, like for mine [17] or un-exploded improvised device [18] detection (where the targets are often very shallowly buried just to hide from sight), then the reference plane can be set just at the air/soil interface. In this case the surface currents radiate in free-space and the related simple Green function can be considered as propagator. It is clear that by this model, we are renouncing to estimate the target depths. This is a minor drawback for targets close to the air/soil interface. On the other hand, the reconstruction is cast as a 2D inverse problems, since only the targets detection and their transverse locations are looked for, and hence single-frequency data can be employed. Accordingly, RXs do not need information about the TX, except the working frequency. Hence, the source can be considered as being non-cooperative (it does not share information with the receivers). However, it is not opportunistic as in most passive radar, since it is deliberately deployed in the scene. Multi-frequency data can be processed separately (i.e., incoherently) and then combined to counteract aliasing artefacts that can arise if the spatial measurement points are reduced (not properly sampled) [15]. Finally, depth can be explored by performing the reconstructions at different reference planes.
In previous contributions (see [14] and [15]) we have shown the feasibility of the method by employing synthetic numerical data and some measured data collected under lab condition for a free-space scattering scenario. In this contribution, we aim at checking the achievable performance for a more realistic scattering scenario where the background medium is actually in-homogeneous. Indeed, the test site mimics a realistic on-field situation, as it consists of an open air sand box. As to the RXs and the Txs, they are not mounted on flying platforms. However, they are at a stand-off distance from the air/soil interface and the TX signal is unknown to the RXs but the working frequency.
The rest of the paper is organized as follows. In section 2 the scattering model is briefly recalled and in section 3 the related reconstruction algorithm described. In section 4 the experimental set-up as well as some experimental results are presented and discussed. Conclusions end the paper.

Scattering model
In this section we briefly introduce the adopted scattering model and the necessary notation that is required in the following. A detailed derivation of the model can be found in [14]. Here, we adapt that derivation to the configuration used in the experimental set up (described in the sequel).
The background medium is two layered with the upper half-space being air whereas the lower one modeling the soil. The two half-spaces are separated by a planar interface (i.e., the air/soil interface) located at z = 0. The scattering targets are located in the lower half-space (i.e., for z < 0) and buried in close proximity to the separation interface. Moreover, the "transverse" investigation domain (i.e. the spatial region where the targets can belong to) is denoted as investigation domain is considered at a fixed depth z T . Generally, we will consider z T = 0 (that is just at the separation interface), but different depths z T < 0 can be considered as well.
The scattering scene is probed by a single source located in the upper half-space at some stand-off distance from the air/soil interface, h t , whereas the field scattered by the buried targets is collected over a set of sensors still located on the air side and all at the same height h O . Accordingly, the spatial measurement positions lay over a plane; say r n = (x n , y n , h O ), n = 1, 2, ..., N O , their positions, N O being 3 of 14 the number of spatial measurements. Figure 1 shows a pictorial view of the scattering configuration along with the adopted reference frame.
After the scene is illuminated by the incident field, the scattered field arises. Since all the targets are located in the half-space z < z T , by invoking the equivalence theorem [16] the scattered field can be considered as being radiated by equivalent surface currents supported over the plane at z = z T . In particular, by filling the half-space z < z T with a perfect electric conductor, only magnetic surface current survives. Say such a surface current. Note that, the magnetic equivalent current depends on the scattering scene and on the incident field as well. As such, it depends on the working frequency which is indeed highlighted in (1) in terms of the wavenumber k. In particular, if z T = 0 (i.e., just at the air/soil interface) the current in (1) radiates in free-space and hence the scattered field (in the upper half-space z > 0) can be written as where r = r t + z Tẑ is the source point, with r t = (x , y ), r is the field point and is the free-space 3D scalar Green function. G(r, r ; k) is the magnetic to electric dyadic Green function, whose expression is given as G xx (r, r ; k) G xy (r, r ; k) G xz (r, r ; k) G yx (r, r ; k) G yy (r, r ; k) G yz (r, r ; k) G zx (r, r ; k) G zy (r, r ; k) G zz (r, r ; k) (4) It must be remarked that, to be rigorous, surface integration in (2) should run over the entire plane z = 0. However, since the targets are very close to the air/soil interface (or in general to the reference plane z T ), it can reasonably be assumed that the current support is very similar to the target's cross section. Hence, D extent is chosen according to the size of the spatial region to be investigated. Again, in (2) we considered the free-space Green function. This is correct for z T = 0. When this is not the case, because the reconstruction at a different depth is required, we will still use the same Green function. Indeed, using the free-space Green function avoids to deal with the computation of the Green function pertaining to a layered background medium. Also, the latter requires the knowledge of the electromagnetic features (dielectric permittivity and conductivity) of the soil, which are not always available. Herein, such background medium electromagnetic parameters are assumed not known. By contrast, using the free-space Green function leads the targets appear more deeply located, because soil is electromagnetically denser than air. However, this is not a serious drawback because the targets of interest are in any case shallowly buried.
The magnetic surface current J eq has no component alongẑ. Accordingly, (2) particularizes as It is seen that E Sx is solely linked to J eqy and E Sy to J eqx . Therefore, if one collects separately such field components then the inverse problem in (5) splits in two identical scalar inverse problems from 4 of 14 which one can reconstruct the two source components independently. Then, these reconstructions can be combined as in [15]. However, in general this is not the case, even in view of the receiving antenna response. More precisely, what one can actually measure is the antenna output voltage. Therefore, in place of (5), the following equation should be considered where V is the voltage data and T a linear operator schematizing the antenna response. Eventually, this is the scattering model from which to start in order to perform the current reconstructions. Few details concerning the reconstruction algorithm along with some further simplifications to achieve such a task are provided in the next section.

Reconstruction Algorithm
According to (6), the magnetic current components cannot be in general separately reconstructed. Moreover, the antenna response should be known and put into the model. It would be useful to avoid both these issues. Indeed, looking for simultaneously both the source components means to deal with a doubled number of unknowns. Also, antenna response must be measured/known in advance.
As to the first question, if the receiving antenna is linearly polarized, for example, in the x − z plane, than its plane-wave spectrum vector belongs to the same plane. Accordingly, the main contribution to the voltage is due to J eqy , which is the one that contributes to E Sx . Similar considerations apply if the other source component is considered. Therefore, we make the problem scalar by approximating (6) as where r n are the measurement positions (antenna phase center) and H(r n , r ; k) is the kernel of the integral operator in (7), which is actually the approximation of the composition between the antenna response operator and the propagator, once the contribution due to J eqx has been neglected. However, the antenna response is still there. It is noted that both the data and the unknown depend on the frequency k. Hence, employing all the available multiple-frequency data to perform the reconstruction is not formally allowed. On the one hand, single-frequency reconstructions, as discussed above, do not allow to estimate target depths. This is however a minor drawbacks for shallowly buried targets and, as previously remarked, depth at which reconstruction is achieved can be changed. On the other hand, at single-frequency, the antenna response basically introduces a complex weight which is the same for each measurement position. This means that it does not affect source localization, which is instead related to the phase term that depends on r n − r . Hence, since we are mainly interested in the detection and the localization of the targets (i.e., we do not aim at quantitative reconstructions, even in view of the other approximations) antenna response is neglected in (7) while achieving single-frequency reconstructions. Note that this would be not the case for multi-frequency reconstructions since the complex weight in general changes with frequency.
In order to perform the reconstruction, the presented model (at a given frequency k l ) is discretized by representing the unknown current component J eqy as a truncated Fourier series, that is Accordingly, the unknowns of the problem become now the expansion coefficients I nm . The choice of N and M is linked to the so-called number of degrees of freedom of the problem, reflects the ill-posedness of the inverse problem at hand and depends on the operating frequency as well as the investigation and the observation domain extensions [14], [15]. In general, however, this discretization scheme allows to reduce at large extent the number of unknowns, as compared to a pixel based representation. The corresponding discrete model is then obtained as where V(k l ) ∈ C N O is the data vector at the l-th frequency, H ∈ C N O ×(2N+1)(2M+1) is the matrix version of the scattering operator and I ∈ C (2N+1)(2M+1) is the (vectorized) expansion coefficient matrix. Equation (9) is inverted for I by a standard Truncated Singular Value Decomposition (SVD) [19] of the relevant matrix operator. This allows to counteract the ill-posedness of the problem and to get a stable reconstruction.
Once the coefficients I nm have been recovered, the corresponding equivalent current J eqy (x, y; k l ) is computed by means of (8). Then, the support of such an equivalent surface current is provided by the image in the x, y investigation domain, I(x, y; k l ) = |J eqy (x, y; k l )| (10) In order to limit the system complexity the number of spatial data has to be reduced. This in general can lead to a reconstruction that is corrupted by aliasing artefacts that are difficult to distinguish from the actual current. To cope with this drawback, a simple strategy based on the combination of single-frequency reconstructions has been introduced in [15]. More in details, say N k the number of adopted frequencies, than the final reconstruction is obtained as The very basic idea behind (11) is that aliasing artefacts change positions with the working frequency whereas the actual source reconstruction does not. Therefore, (11) tends to mitigate all those peaks in the reconstruction that do not overlap (or overlap only partially) while the frequency changes. A criterion for the choice of the frequencies is provided in [15].
Summarizing, the algorithm presents the following steps: 1. Fix one frequency value; 2. Compute the scattering matrix model H; 3. Compute the SVD of H; 4. Fix a regularizing threshold for the normalized singular values of H (in the following experimental results, 20dB is used) and compute the unknown vector I via a TSVD inversion by retaining the data projection over the singular vectors corresponding to the singular values above the threshold; 5. Calculate I(x, y; k l ) by (8) and (10); 6. Repeat from point 1 by changing the frequency value; 7. Compute I(x, y) from (11).

Experimental Results
In this section we check the proposed algorithm against some experimental measurements collected in a semi-controlled scattering scenario. In particular, we first describe the test site and then show some reconstructions aiming at highlighting the role of the number of spatial measurements and of the employed frequencies.

Test site
The test site consists of a tank full of sand of about 3.5 m (length) 2.5 m (width) and 1.5 m (depth) in size. The tank is placed in the open air, so that the sand appears wet, apart the very surface layer which has been dried by sun. The electromagnetic features of the sand are unknown and are not estimated for detection purposes.
The transmitting antenna is a horn positioned at a h t = 1.5 m height from the sand floor and located at one of the end side of the tank. It is tilted so to pointing to the spatial region under investigation. The receiving antenna is still a horn and it is located on the other side of the tank. In particular, it is mounted on a wooden slide that allows to synthesize a planar measurement aperture at a fixed height from the air/sand interface (in the following examples h O = 80 cm or h O = 130 cm). Also, the receiving antenna is tilted toward the investigated spatial region and is linearly polarized. Fig.2 shows a schematic view of the measurement configuration along with some pictures of the test site. As a target a metallic rectangular plate 17.5 cm × 48 cm in size is considered.
A vector network analyzer is connected to the antennas by means of coaxial cables. Standard calibration at the end of each channel has been performed at the beginning of each measurement stage. Data have been acquired in the frequency band [2 − 9] GHz (201 equispaced frequencies) for each different position of the receiving antenna.
We performed measurements under two different conditions: flat and rough air/sand interface. Flatness was obtained by manually using a shovel for smoothing the sand floor. Of course, the obtained sand surface, though being smooth, it is far from being "ideally" flat. Roughness interface has been instead obtained by turning over the sand. For such scenarios, we took measurements with and without (background measurements) for comparison purposes.
It is worth remarking that the measurement scenario actually contains many features which are not accounted for in the scattering model used to developed for the detection algorithm. For example, though antennas are tilted towards the scattering scene they still present a direct link which implies direct coupling between the receiving and transmitting antenna. Also, because of the finite dimensions of the tank, the two-half space medium assumption clearly does not correspond to the actual background medium. This entails, that the received signal actually consists of different contributions besides the one expected from the targets. Also, note that, even though the medium were a precise two-layered medium with a flat separation interface, the air/sand interface reflection always superimposes to the target signals. Nonetheless, in order to mitigate such unwanted contributions we did not pre-processed data in the following reconstructions.

Detection results for flat air/soil interface
We start by considering the case of flat air/sand interface in the sense clarified above.   The investigation domain is a rectangle of 160 cm along x and 100 cm along y and it is located at z T = 0. Also, we introduce the two parameters o f f setx and o f f sety, which indicate the displacement, along x and y respectively, of the centre of the investigation domain with respect to the central point of the first measurement line (see Fig. 3a for the reference system). Basically, after acquiring the data, changing the investigation domain centre location (by varying the offset parameters) entails looking for the targets in different spatial regions. Accordingly, the same target will appear at different relative positions. This can be considered as a way to check reconstructions consistency and stability. A detection can be considered successful if the target localization point moves coherently with the change of the investigation domain centre position.
It must be remarked that the number of employed spatial data is already below the ones required if the field has to be properly spatially sampled (see [15]). This, of course, limits the highest frequency that can be used in the reconstructions. Indeed, we tested the inversion algorithm against different bands inside the available one of [2 − 9] GHz. As expected, when dealing with higher frequencies, even by our approach, detection is impaired because the reconstructions resulted crowed by a number of artefacts. Hence, we ended up by processing data collected only within the band [2 − 4] GHz. In particular, in such a band, N k = 11 frequencies was considered in order to be as close as possible to frequency selection criterion provided in [15].
Finally, the target is located as shown in Fig. 3b and roughly buried 2 cm below the air/sand interface.
In particular, Fig. 4 reports the reconstructions for a target located as in position 1 sketched in 3b, for different investigation domain centre offsets. As can be seen, a good detection is achieved with no significant artifacts corrupting the reconstructions. What is more, the reconstructed spot moves coherently with the investigation domain displacement. Some comments are here in order. Firstly, we remark that the inversion algorithm does not aim at providing the shape of the targets but rather it is intended for target detection. This is mainly due to the strategy we adopted to combine the different single frequency reconstructions. Indeed, the proposed multiplicative combination highlighted in (11) allows to mitigate aliasing artefacts but at the same time tends to enhance the strongest part of the reconstructions so that targets actually appear as hot spots.
Secondly, as we mentioned above, we did not processed data to counteract the direct link (form the transmitting antenna to the receiving one) or the air/sand reflection. Yet, reconstructions do not suffer from such spurious signals. This can be justified by observing that we have performed the reconstruction over a given spatial region: the investigation domain. Hence, the direct link should appear located outside such a region. As to the air/sand interface reflection, it actually enters in the investigation domain. However, it is not localized as the target contribution and tends to be spread over the whole investigation domain. The multiplicative combination strategy hence also enhances the target reconstruction against such a contribution. In Figs. 5 and 6 reconstructions corresponding to the target (still approximately buried at 2 cm) located as in position 2 (see 3b) are reported. In particular, while Fig. 5 has been obtained using the same source polarization as in the previous case, in Fig. 6 the transmitting antenna has been rotated 90 o so that the scene is illuminated by an orthogonal polarization. As can be seen, also in this case the target is clearly detected, regardless of the transmitting antenna features (in this case polarization).

Detection results for rough air/sand surface
We now turn to consider the case the air/sand interface has not smoothed. In this case we consider data collected over a grid of 8 × 7 positions with the same spatial step as above but at a height h O = 130 cm. Note that the spatial data are slightly greater than the previous case but still under-sampled [15].
First, we show the reconstruction of a shallowly buried target (the target depth and type are the same as above) by employing all the available frequency. These results are reported in the left column of Fig. 7. As can be appreciated, the target is clearly detected and the related hot spot indicator changes position accordingly to the investigation domain centre offset. On the same figure (right column), instead, we report the reconstructions obtained by processing background data, i.e. in absence of the target. Differently form the target case, now the reconstructions do not exhibit a clear hot spot. Moreover, the reconstruction changes as the investigation domain centre offset varies. This confirms previous discussion that processing air/soil interface reflection does not return a focused hot spot. Also, the reconstruction corresponding to this contribution is in general different when the investigated spatial region changes. This is in particulate true here because roughness entails that the air/soil interface has different spatial details. Eventually, these results suggest a possible strategy to recognize actual targets against surface clutter. Indeed, comparing images obtained using different investigation domain offsets, the actual targets are those ones for which the reconstructions "move" coherently with the change of the investigation domain center.
In Fig. 8 the same case as in Fig. 7 by keeping fixed the investigation domain with o f f setx = 2.0m; o f f sety = 0.3m is considered. However, different numbers of frequencies are employed. It can be seen, that increasing the number of frequencies aliasing artefacts actually tend to disappear and the hot spots narrows. This is, of course, expected and perfetly consistent to the theoretical arguments discussed above.
The simple proposed strategy for reducing spatial data hence works very well. To further check this procedure we consider an even more challenging case by further reducing the spatial data employed to achieve the reconstructions. More in detail, we collect data over a 4 × 4 measurement grid, with twice the spatial step, that is 40cm. The height of the measurement aperture is still h r = 130 cm and the frequencies taken within the same band as above. Note that in this case the number of spatial data is even lower than the ones used in the first example reported at the beginning of this section.
It can be seen by Fig. 9, that even under this more challenging situation where the spatial data have been further reduced, the proposed method works very well in detecting and localizing the target and in counteracting aliasing artefacts, though the number of frequencies is actually low as well.
As a final example, we consider the case the target is buried more in depth. In particular, in this case, the plate target is approximately buried 10 cm below the air/sand interface. Fig. 10 shows the reconstructions obtained by considering the image plane at different depths. In all the examples, N k = 7 frequencies (that showed to be enough for mitigating aliasing artefacts) have been employed and the same spatial grid measurement as in Fig. 9 considered. These results show that the proposed method is still effective in detecting the target. Moreover, as discussed above, changing the depth at which reconstruction is achieved allows to get information about the target depth as the one for which the reconstruction is more focused. Of course, since the soil permittivity is not known (and hence not accounted for in the reconstruction algorithm), the estimated target depth does not coincide with the actual one. In particular, since the relative dielectric permittivity of a wet sand likely stands between 4 and 9, the "apparent" target depth can range from 20cm to 30cm, which is perfectly consistent to the result shown in Fig. 10.

Conclusions
In this paper we have experimentally validated a method for detecting and localizing shallowly buried targets from under-sampled multi-frequency data. The method relies on two main ingredients: a suitable scattering model and a simple procedure to process multi-frequency under-sampled data.  Figure 8. The same case as in Fig. 7 with o f f setx = 2.0m; o f f sety = 0.3m. Comparison of reconstructions obtained by employing different number of frequencies N k .
(a) N k = 2 (b) N k = 3 (c) N k = 7 Figure 9. Normalized reconstructions of a shallowly buried target corresponding to the case reported in Fig, 7, o f f setx = 2.0m, o f f sety = 0.3m but with only 4 × 4 measurement grid collected with a spatial step of 40cm. The actual target location is denoted as red rectangles. Comparison of reconstructions obtained for different values of N k . The scattering model is based on the equivalence theorem which allows to cast the detection as the reconstruction of equivalent sources supported over a reference plane. This way detection becomes a 2D problem. Also, the incident field is embodied into the unknown equivalent sources. Therefore, coherence between the TXs and Rxs is not necessary when single-frequency data are employed. Reconstructions obtained at different frequencies are then suitably combined. This allows to mitigate aliasing artefacts and hence to achieve the reconstructions with a very reduced number of spatial measurements. Experimental results confirm the feasibility of the method even under a rather complex scattering scenario, as the one addressed herein. Indeed, it is shown that the targets can be detected by using very few spatial measurement points and a number of properly selected frequencies. Also, the results confirm that robustness of the methods against the reflection occurring at the air/soil interface even when this exhibit some degree of roughness.
The obtained results can be meant as a proof of concept. However, they encourage towards the application of the proposed measurement configuration and of the related inversion algorithm to more challenging scenarios, where many targets can be present inside a not homogeneous terrain and the sensors are deployed at a large stand-off distance.