An Effective Method for Borehole Imaging of Buried Tunnels

Detection and imaging of buried tunnels is a challenging problem which is relevant to both geophysical surveys and security monitoring. To comply with the need of exploring large portions of the underground, electromagnetic measurements carried out under a borehole configuration are usually exploited. Since this requires to drill holes in the soil wherein the transmitting and receiving antennas have to be positioned, low complexity of the involved apparatus is important. On the other hand, to effectively image the surveyed area, there is the need for adopting efficient and reliable imaging methods. To address these issues, in this paper we investigate the feasibility of the linear sampling method (LSM), as this inverse scattering method is capable to provide almost real-time results even when 3D images of very large domains are built, while not requiring approximations of the underlying physics. In particular, the results of the reported numerical analysis show that the LSM is capable of performing the required imaging task while using a quite simple measurement configuration consisting of two boreholes and a few number of multiview-multistatic acquisitions.


Introduction
Detection and geometric characterization of underground tunnels and cavities are classical problems in geophysical monitoring that, nowadays, also represent challenging issues for both military and homeland security.As a matter of fact, manmade or natural cavities may engender surface subsidence and mine the stability of populated area, whereas clandestine tunnels provide hidden connections to move across countries borders unauthorized people or illicit goods, such as drugs, weapons, and explosives [1].
Due to the relevance of this topic, many active and passive geophysical techniques, using gravity, seismic and electromagnetic waves have been developed and are currently used [2][3][4][5][6][7][8].However, the need for achieving accurate images by means of low cost and not time-consuming sensing devices as well as the intrinsic limits of the existing imaging approaches still make tunnel detection a challenging issue.
In this framework, ground penetrating radar (GPR) is worth to be considered due to its flexibility and its successful use in many other subsurface prospecting applications.However, with respect to the specific problem of underground tunnel detection, surface-operated GPR suffers of several issues, such as low-resolution limits, interference due to external electromagnetic sources, and complexity of the measurement apparatus [9,10].Hence, borehole GPR configurations are commonly preferred, and several crossborehole imaging systems exploiting radar approaches have been proposed [11,12].
Moving in this direction and focusing the attention on standard borehole measurement configurations, this paper aims at investigating the performances offered by the LSM [13] to tackle this imaging problem.The LSM is an inverse scattering method for shape reconstruction of unknown targets that belongs to the class of full-wave approaches, thus being suitable to tackle "strong" scatterers having electric properties very different from those of soil, such as indeed hidden cavities and tunnels.On the other hand, the LSM does not suffer from the drawbacks which affect other nonlinear inverse scattering approaches (especially when dealing International Journal of Antennas and Propagation with large investigated domains).As a matter of fact, by restricting the aim to the detection and the morphological characterization of the target, the LSM can tackle the imaging task through the solution of an exactly linear inverse problem, where the knowledge or the computation of the total field into the investigated domain is not needed [14].Therefore, the LSM is simple, computationally effective and suitable to rapidly achieve full three-dimensional reconstructions of the investigated regions even of large extent.In addition, the capability of the LSM of efficiently delivering full three-dimensional images of the investigated scenario is particularly convenient in this framework, since the peculiar geometrical features of tunnels makes it possible to discriminate them from other spurious targets present in the surveyed area.
However, the aforementioned reduced complexity of the LSM is traded with some requirements on the measurement device.As a matter of fact, while the LSM can properly work with monochromatic and single polarization data [15,16], it requires that data are gathered under a multiview and multistatic measurement configuration, in which each antenna acts both as transmitter and receiver, and for each transmitting position the backscattered signal is measured at all receiving locations.Owing to this requirement, a known limitation of the LSM is the need of acquiring data with enough "diversity" [17,18], that is considering a sufficiently large number of transmitting/receving pairs.Accordingly, given the importance of reducing the system's complexity, one of the goal of this study is to assess the method's performance when only two boreholes and a reduced number of multistatic pairs are considered.
It is worth to remark that the LSM also favorably compares with usual approaches for GPR imaging based on linear approximations, such as migration and diffraction tomography [19][20][21].As a matter of fact, while opposite to these latter, the LSM cannot be applied for multi-monostatic multifrequency data (that would allow a further simplification of the measurement configuration), its not approximated nature makes it more effective to handle complex scenarios in which multiple scattering occurs.Moreover, the fact that the LSM adopts monochromatic data entails that it does not require to take into account the frequency dispersion of the involved media.In addition, the computational burden of the LSM is so low that real-time imaging is viable even for 3D problems.
To the best of our knowledge, the LSM has not been applied to this kind of problem, so that it is certainly interesting to verify its imaging performances taking into account the above considerations.
The paper is organized as follows.Section 2 describes the borehole configurations we consider in this study.In Section 3 the general formulation of the LSM in 3D fullvectorial case is recalled, while the peculiar aspects relevant to tunnel detection are outlined in Section 4. In Section 5, the performance of the imaging strategy is tested with respect to showing the method's capability of actually performing the desired imaging task even using very simplified measurement configurations.Conclusions follows.

Geometry of the Problem and Borehole Imaging Configurations
The considered geometry is sketched in Figure 1.The tunnel, whose electromagnetic features are those of air, is assumed to be embedded into a lossy background medium, that is, the soil, characterized by a permittivity s and a conductivity σ s .Assuming that the distance between the antennas and the air-soil interface is large, we model the background as a homogeneous medium, that is, we neglect the effect of the air-soil interface.Ω denotes the investigated volume.Two measurement setup are considered, which are typical of borehole imaging applications.In particular, Figure 1(a) shows the common cross-hole measurement configuration in which the imaged region is "enclosed" between the boreholes, while Figure 1(b) depicts the reflection setup, in which the two boreholes lie on the same side of the investigated region.
In both configurations, transmission and reflection data are collected by means of probes working under a multiviewmultistatic configuration moved along each borehole.
Accordingly, each probe acts as transmitter and receiver, that is, for each position of the transmitter, the field is collected by all other probes.

Formulation of the LSM for 3D Problems
The LSM is a qualitative imaging approach, whose aim is to reconstruct the morphological features of unknown targets by simply processing the scattered field data and without requiring any a priori knowledge on the electromagnetic nature of the objects.This goal is pursued by illuminating the region of interest by a set of known incident fields and by collecting in a set of measurement locations the scattered fields.Without loss of generality for the method applicability and by addressing our formulation for the imaging problem at hand, we assume that the number of transmitters and receivers is the same (V = M).Then, the imaging task is achieved by partitioning the investigated region Ω into an arbitrary grid of sampling points and by solving for each of them a linear equation.
By referring the reader to [14,15] for the LSM formulation in the 3D vectorial case in which multiview-multistaticmultipolarization single frequency data are required, in the following we adopt the LSM formulation when the data are collected with a fixed polarization of the transmitting and receiving antennas.In this respect, by denoting with E s the M × V multistatic-multiview data matrix collected at a fixed working frequency, the LSM requires to solve for each sampling point r s = (x s , y s , z s ) ∈ Ω the matrix equation: wherein x is the V-dimensional unknown vector, and g p d (r s ) is the M-dimensional vector, whose elements are the values at the M receiver positions of the p-component of the field radiated by a test source, that is, an elementary electric dipole located in r s .d denotes the orientation of the electric dipole.
Due to the ill-conditioned nature of (1) [14], a regularization strategy has to be considered to get a reliable solution.By exploiting the singular-value decomposition (SVD) of E s and Tikhonov regularization, a stable solution of (1) in each r s can be achieved as wherein λ k are the singular values of the discretized matrix operator E s , u k and v k are the left and right singular vectors, respectively, K = max(M, V ), while • and •, • denote the L 2 -norm and the scalar product, respectively.The Tikhonov weighting coefficient α can be fixed according to [22] in order to avoid that x(r s ) 2 blows up due to the vanishing behavior of the singular values.Notably, this can be done even in the lack of a precise knowledge of the signal-to-noise ratio and once for all the sampling points [22].
According to the LSM theory [13], x(r s ) 2 assumes its lowest values (with respect to the overall behavior in Ω), when r s belongs to the scatterer support and larger values elsewhere.Hence, it can work as an indicator of the target's support, since its plot is a straightforward way to achieve morphological reconstructions.Interestingly, the above property can be explained from a physical point of view through the analogy existing between the LSM and the problem of focusing the electromagnetic field in presence of unknown targets [23].
Last but not least, a consideration of the method flexibility is worth to be done.In this respect, we point out that the main computational effort involved by the LSM is the evaluation of the SVD of the single-frequency data matrix E s , and that this task has to be carried out only once, since the kernel of ( 1) is the same for all the sampling points.As a result, the dimension of the involved matrix is dictated by the number of measurements, thus the processing time takes only few seconds for a full 3D reconstruction.

Applying the LSM to Borehole Imaging
In this section, we discuss some specific aspects which are relevant to the application of the LSM to the imaging configuration and problem we are considering.In particular, we will discuss how to match the requirement of reducing the complexity of the setup as well as how to improve the result of the imaging process.
As far as the setup is concerned, let us first note that in casting (1), the polarization of the various fields involved has not been specified.As a matter of fact, (1) actually corresponds to a family of equations in which the polarization of the incident field, that of the scattered one, as well as that of the "testing" field (the right-hand side of the equation) is, to some extent, a degree of freedom.
In this respect, in [16] it has been shown that satisfactory shape reconstructions can be achieved via LSM even when polarization diversity is not considered, that is, using transmitting and receiving antennas of the same kind and having the same polarization p, thus allowing a remarkable simplification of the imaging setup.
Taking this result into account and recalling that the most utilized antennas in borehole imaging applications are wire antennas, as for instance, short-or half-wavelength dipoles [24,25], we consider two possible cases for the probes to be adopted: (i) the probes are vertical dipoles, that corresponds to assumption p = z; (ii) the probes are horizontal dipoles, that corresponds to assumption p = y.
Then, exploiting the guidelines suggested in [17,26], the probes are assumed to be evenly spaced of λ b /2 along the borehole length, λ b being the background medium's wavelength.By relying on the same considerations, the spacing between the two boreholes in the "reflection" configuration is also assumed to be λ b /2.
As far as the imaging is concerned, we can assume that the test elementary dipole is oriented in the same way of the measured field component (p = d).Note this choice is expected to maximize the sought fitting between the two International Journal of Antennas and Propagation sides of the LSM equation when using the same polarization for transmitters and receivers [16].Thanks to such a simplification, we are able to build a full 3D reconstruction of the imaged scenario while simply solving a single linear and scalar problem, which obviously entails a remarkable measuring and processing efficiency.Then, we take advantage of a normalization of the indicator function [16].As a matter of fact, since the probes are not located all around the investigated domain, the indicator x(r s ) 2 can lead to poor reconstructions of those parts of the imaged domain located far from the boreholes.This drawback can be overcome by considering the compensated indicator: wherein the weight at the denominator is used to take into account the different location of the sampling point with respect to the receivers.Moreover, to improve the readability of the results, as usual in the LSM imaging, we adopt a logarithmic plot of such indicator.Finally, since the survey's aim is to obtain an image capable of informing the operator of the tunnel's presence, the problem of how visualizing the result has to be faced.Actually, this task is accomplished by plotting a binarized version of the compensated indicator function (3) in the whole volume adopting some cutoff value, so to achieve an isosurface representation of the target.In this respect, the important issue is how to choose the cutoff value.To this end, we adopt a simple strategy that consists in inspecting only the most meaningful single cross-section of the indicator map.In particular, by relying on the results in [18], we consider the cross-section in which the borehole lie, in the cross-hole case, and the cross-section equidistant to the two boreholes, in the reflection setup.Even if the cutoff is determined heuristically, we rely on the observation that, again, the presence of the tunnel is expected to produce strong effects in terms of scattered fields, so that, especially in the selected crosssections, its presence should be easily detectable.Moreover, a slight error in the cutoff value is not expected to affect the information content of the overall result.

Performance Assessment
In the following, we assume the tunnel embedded into a homogeneous soil having b = 5 and σ b = 0.5 mS/m.The tunnel is simulated as a rectilinear vacuum with a square section of side 3 m, oriented parallel to the y-axis and lying at depth of 20 m below the air-soil interface.This test-bed scenario is depicted in Figure 2. The investigated domain is a cubic region having side of 40 m, and we suppose that the holes are drilled in z-direction.The total length of the array in the boreholes is 30 m starting from 5 m below the airsoil interface.According to the rules recalled in the previous section, a total number of M = V = 18 probes is positioned in the two boreholes.
We process synthetic data collected at 10 MHz in both the cross-hole and reflection setup configuration.In particular, the scattered field is simulated by means of a forward solver based on the method of moment (MoM), and it is disturbed with a Gaussian noise of SNR = 25 dB.Moreover, it is worth to notice, that the data elaboration we are considering is not affected by "inverse crime" (even with noise-free data), since two completely different models underlie the forward problem (i.e., the MoM) and the inversion procedure (i.e., the LSM).

Cross-Hole Measurement Setup.
First of all, we have assessed the reconstruction capabilities of the LSM against the usually adopted measurement configuration in which the region between the two boreholes is investigated.The boreholes are located at (x = −30 m; y = 0 m) and (x = 30 m; y = 0 m).First, we have considered the case in which the probes are vertical dipoles.
Once the LSM equation has been solved for all sampling points in the investigated volume, we exploit the simple strategy discussed in the previous section.Hence, we first observe the cross-section of the LSM indicator in the plane y = 0, that is, the plane wherein the boreholes are positioned.Figure 3(a) shows this result, while Figure 3(b) depicts the corresponding binarized image (the contour plot denotes the actual geometry of the tunnel).As we can see, a good estimation of both the location and shape of the tunnel is retrieved.Then, on the basis of the so-achieved threshold value, we build the 3D isosurface representation of the target, Figure 3(c).As expected, the low resolution achievable in the longitudinal direction, that is, along the orientation of the tunnel, prevents to achieve a satisfactory reconstruction in the lateral regions of the investigated domain located far from the position of the boreholes.However, by considering such a plot in the region closer to the boreholes' position, the presence of a target and its geometry can be clearly associated with that of a tunnel owing to its length.
The results observed when using horizontal dipoles are shown in Figures 3(d)-3(f).Also in this case, the results are quite satisfying both in the transverse section and in the isosurface 3D representation.

Reflection Measurement Setup.
In this configuration, the two boreholes are placed on the same side of the investigated domain, and, according to the criteria outlined in the previous section, their spacing is 7 m (x = 50 m; y = ±3.5 m).Figures 4(a), 4(b), 4(d) and 4(e) show that also in this case it is possible to successfully perform the imaging of the tunnel in the cross-section equidistant from the two boreholes.In particular, the position and the transverse dimension of the tunnel can be retrieved by looking at the projection of the LSM indicator over the z-x plane.On the other hand, the three-dimensional isosurface representation provides an underestimation of the tunnel length, see Figures 4(c) and 4(f).This result can be explained recalling the limited imaging capabilities of the LSM with depth when using aspect-limited configurations [18], which entails that a smaller portion of the probed area can be reliably imaged as compared to the previous configuration.

Cluttered Scenario.
As a further proof of the performance achievable by the LSM in the imaging of hidden tunnels, we consider as a last example the case wherein the tunnel is not simply embedded in a homogeneous background, but some other targets (representing clutter for the survey's aim) are present in the investigated region.In particular, in the considered scenario, several spurious targets lie all around the tunnel, see Figure 5(a).The permittivity values of these targets assume random values between those of the air and the soil.
For the sake of brevity, we consider only the crosshole measurement configuration with vertical dipoles, for which, on the base of the above-discussed results, better performances are also expected.
Moreover, in order to give a further proof of the achievable performance, we report also a comparison with the LSM imaging based on a common surface GPR measurement configuration.In such an arrangement, an array of 9 equispaced antennas is moved along the y-direction with an acquisition step of 5 m.Note that the acquisition step has been chosen according to the above outlined criteria concerning the spatial sampling of the scattered field (6 m λ b /2), see Figure 5(b).
When the cross-hole borehole configuration is used, the LSM imaging is still capable of imaging the transverse section of the tunnel, even if the sharpness of the imaging is lower as compared to the unperturbed scenario, that is, when no clutter is present, see Figure 5(c).This notwithstanding, by exploiting the threshold value resulting from this image, it is possible to obtain a satisfactory three-dimensional image, wherein the detection of the tunnel is still possible (Figures 5(d) and 5(e)).On the other hand, in the surface measurement configuration, the achieved results significantly deteriorate.As a matter of fact, the LSM processing gives back significantly worsen results as neither the transverse section of the tunnel is well localized and estimated, nor the achieved 3D isosurface rendering resembles the geometry of the tunnel (Figures 5(f)-5(h)).Finally, a general comment concerning the imaging of hidden tunnels in more realistic instances it is worth to be done.In particular, when the tunnel will not be parallel to the boreholes, we can expect that LSM will be still able to provide a fair estimate of the tunnel location and shape by considering the 2D plot of the indicator in the transverse plane.On the other hand, 3D reconstructions are expected to slightly deteriorate depending on both the inclination angle and the distance between the tunnel branch and the borehole location.

Conclusion
In this paper, we have investigated the suitability of the LSM in the framework of GPR borehole imaging for tunnels detection.The possibility to obtain a qualitative reconstruction of hidden tunnels has been corroborated through synthetic data.In particular, by taking into account the lowcomplexity measurement configurations required by the GPR surveys for such a kind of applications, we tackled the microwave imaging task with a low number of transmitting/receiving probes and by adopting only two boreholes.Satisfactory reconstructions from 3D data processing can be achieved even when the tunnel is far extended (in terms of wavelengths) in the dimension transverse to that of the boreholes axes.The high flexibility of the LSM in performing the imaging task by processing only one component of the scattered field (single-polarization data) as well as the possibility to handle the three-dimensional full-vectorial problem as a scalar one reveals the suitability of the method to tackle tunnel detection with low-complexity apparatus, very low computational burden and not time-consuming processing step.
This preliminary assessment encourages further analysis involving more realistic scenarios, wherein to take into account the antenna model and different soil/background medium.Moreover, the capability of the proposed method could be also tested with possible different configuration of the probes, as for instance, by positioning the antennas into nonrectilinear boreholes or by considering data collected by means of directional borehole arrays [27].Finally, the imaging capabilities on experimental data will be assessed as well.

Figure 1 :
Figure 1: Reference scenario and the adopted measurement configurations.The volume Ω is the imaged region which is probed by antennas along the two boreholes: (a) cross-hole measurement setup and (b) reflection measurement setup.

Figure 3 :
Figure 3: LSM imaging of the tunnel in the cross-hole measurement setup.Normalized logarithmic plot of the LSM indicator, binarized LSM indicator with threshold 0.75, and binarized three-dimensional LSM indicator in the whole region: (a)-(c) vertical dipoles and (d)-(f) horizontal dipoles.

Figure 4 :
Figure 4: LSM imaging of the tunnel in the reflection measurement setup.Normalized logarithmic plot of the LSM indicator, binarized LSM indicator with threshold 0.75, and binarized three-dimensional LSM indicator in the whole region: (a)-(c) vertical dipoles and (d)-(f) horizontal dipoles.

Figure 5 :
Figure 5: LSM imaging of the tunnel in cluttered scenario.(a) Cross-hole measurement configuration; (b) surface moving array measurement configuration.Results achieved with borehole measurement configuration: (c) normalized logarithmic plot of the LSM indicator; (d) binarized LSM indicator with threshold 0.75; (e) binarized LSM 3D indicator in the whole region.Results achieved with surface measurement configuration: (f) normalized logarithmic plot of the LSM indicator; (g) binarized LSM indicator with threshold 0.75; (h) binarized LSM 3D indicator in the whole region.