Radio Frequency Tomography for Tunnel Detection

Radio frequency (RF) tomography is proposed to detect underground voids, such as tunnels or caches, over relatively wide areas of regard. The RF tomography approach requires a set of low-cost transmitters and receivers arbitrarily deployed on the surface of the ground or slightly buried. Using the principles of inverse scattering and diffraction tomography, a simplified theory for below-ground imaging is developed. In this paper, the principles and motivations in support of RF tomography are introduced. Furthermore, several inversion schemes based on arbitrarily deployed sensors are devised. Then, limitations to performance and system considerations are discussed. Finally, the effectiveness of RF tomography is demonstrated by presenting images reconstructed via the processing of synthetic data.


I. INTRODUCTION
I N RECENT years, underground void detection and localization has become a critical task, motivated by the incumbent need for protecting national borders and for monitoring sensitive areas, such as prisons, banks, and power plants. Additionally, tunnel detection is imperative for civil applications, including mining safety, search and rescue in devastated areas, environmental engineering, geophysics, archaeology, and speleology.
GPR appears as the most versatile approach. However, it cannot yet be considered a reliable and practical method in tunnel detection for the following reasons. First, the resolution for classical GPRs is generally improved by using large bandwidth, often requiring high frequencies. However, due to soil properties, higher frequencies experience higher attenuations, and the increased frequency/bandwidth affects the signal-tonoise ratio (SNR) and intensifies dispersion effects. Second, systems using lower frequencies [13], [14] require electrically small wideband antennas, which results in complex wideband systems; yet, they may not provide adequate resolution to determine the geometry of targets. Third, the available frequency spectrum for some applications can be severely limited by interference with external sources (e.g., broadcasting stations); therefore, a reliable system must operate using a small, discrete, and selectable number of frequencies. Fourth, the interpretation of the raw data is affected by the operator's expertise, and a priori information is necessary to obtain reliable results [10].
To overcome attenuation losses (thus enabling detection at large investigation depths), GPR may be inserted into boreholes [12]- [24], which generally provides an image of a vertical section in the plane between the logs. However, for the important task of tracing tunnel pathways and localizing adits, the "horizontal" prospecting is more desirable, leaving the determination of depth as second priority. Furthermore, boreholes are expensive, subject to drilling misalignments, and, most important, impractical in inaccessible terrains.
The enhancement of resolution may be accomplished by implementing the principles of RF/microwave tomography in the framework of GPR [25]- [31], referred to as GPRT. Although an improvement in reconstructed images may be achieved, GPRT still fails to be reliable for tunnel detection for several reasons, including as follows: 1) GPRT is typically employed in a bistatic configuration, where receive and transmit antennas are separated by an electrically small distance, while operating in proximity to the air-earth interface; therefore, it suffers from limited view diversity (using finite observation domains) [32]. 2) To compensate for this limited view, an increase in the amount of information is achieved by using multitude of frequency tones, thus making the system wideband [28], [31], [33]. 3) Any GPR survey requires the operator to be above the area under investigation. To date, no system has been designed to work effectively in remote situations where one is disallowed to reach the area of regard.
This paper extends, further improves upon, and gives a unified framework to preliminary ideas in RF tomography proposed in [33], [37]- [42]. In particular, Section II describes the concept of RF tomography for belowground imaging; Section III illustrates the principles of RF tomography, presenting the related forward model. Section IV describes several inversion schemes. Finally, Section V shows several reconstructed images obtained via the inversion techniques described in Section IV.

II. RF TOMOGRAPHY
To address the unresolved problems in Section I, we introduce a new approach based upon a multiview/multistatic design, where the view (associated with the distributed transmitters) and the observation (associated with multiple receivers) diversities increase the information pertaining to the scene [34], [35] while reducing the necessity for a large spectral content (bandwidth) of the probing waveform: In principle, with just a single frequency, it is possible to obtain high-resolution images [34], [36].
The proposed approach considers two separate sets of N electromagnetic transmitters (Tx) and M electromagnetic receivers (Rx), commonly referred with the generic name of Transponder. These transponders are placed on (or in) the ground at an arbitrary position properly defined by the operator. The transmitter Tx radiates a known waveform using a suitable polarization. The probing wave impinges upon a buried dielectric or conductive anomaly, thus generating a scattered wave-field. The receivers collect samples of the scattered instantaneous electric field and estimate the complex-valued electric field phasor at their locations. Subsequently, this information is relayed to a data collector (Fig. 1). To ease the detection process, one transmitter and frequency of operation per time are activated, thus simplifying the receiver's capability to properly discern the origin of the incoming wave-field. For a given sampling time, the used spectrum is virtually restricted to a single frequency component, thus ensuring ultranarrowband system architecture, low noise, and affordable cost. To ease the set up and portability of the system, sensors are intended to be "dart" shaped, as shown in Fig. 2. Sensors may be equipped with built-in GPS for precision timing and positioning and an S-band communication link to transfer the collected data to the overhead base station. During the sensor-dropout process, some transponders may fall in obstructed or heavily cluttered regions (e.g., Tx and Rx lay in proximity or over a vegetation layer) or they may not survive the impact with the ground. In any case, This approach may also operate using a discrete set of monochromatic frequency components, selected according to environmental conditions. A suitable modulation is stepped FM; however, in this paper, the conclusions are derived independently from the type of modulation. The operating frequencies must allow the electromagnetic wave to penetrate deeply into the ground, while simultaneously provide acceptable detection and resolution capabilities. Higher frequencies lead to better resolution, but strong attenuation limits the range of operation [60]. Conversely, if the frequencies are low, the corresponding resolution may not be adequate to localize tunnels, and the field behavior becomes diffusive, thus reducing the backscattered field [29]. Based upon electrical parameters of rocks reported in [10], [61], and [62], a suitable range of frequencies for this application is the range of 1-15 MHz, but the final choice strictly depends on the expected target type and/or depth.
When antennas are located belowground, three different modes of propagation between transmitter and scatterer, and between scatterer and receiver, are excited: direct, reflected, and lateral waves [52]- [54]. Generally, the lateral mode of propagation is the most undesirable: although it could be estimated and removed, its contribution to the overall field may be so high that it may saturate the LNA at the receiver side, thus masking the weak signal coming from the scatterer.

III. FORWARD MODEL
The first step to obtain a general tomographic inversion procedure is to define a suitable model for the electric field that closely represents the actual scene while keeping enough simplicity to facilitate the subsequent inversion. Similar procedures are also discussed in [2], [44], and [46]. The 3-D geometry shown in Fig. 3 is considered. For simplicity, a single operating frequency f is assumed, but extension to the multifrequency operation is straightforward [40]. Under the monochromatic assumption, the ground is modeled as a homogeneous medium with relative dielectric permittivity ε D , conductivity σ D , and permeability μ 0 . The targets (i.e., tunnels or voids) are assumed to reside in the investigation domain D. The sources are N electrically small dipoles (of length Δl t ) or loops (of area A t ) fed with current I t and located at position r t n (view diversity) and with (electric or magnetic) dipole moment directed along the unit vectorâ t n . For each transmitting antenna, the scattered field E S is collected by M receivers (observation diversity), located at r r m points in space. We assume the relative dielectric permittivity profile ε r (r ) and the conductivity profile σ(r ) inside the investigation domain D as unknowns of the problem.
Accordingly, the inverse problem is recast in terms of the unknown dielectric permittivity contrast [35], [43] In this way, the wavenumber inside D can be expressed as [30] The function in (1) accounts for the difference between the unknown dielectric permittivity of the object and that of the host medium. We also accounted for the conductivity profile in (1) and (3), because, in general, the halo of a tunnel shows an increased apparent conductivity [12], primarily due to salt dissolved in wet surfaces, steel-reinforced concrete for ground control, power lines, rails, ventilation tubes, etc. For each point r in region D, the vector wave equation holds The scattered wave in a point r / ∈ D that is a solution of (4) can be written in terms of integral equation of the dyadic Green's function where E(r ) is the total field in the investigation domain D, given as the superposition of the incident field E I (r ) (i.e., the field in the investigated area when objects are absent) and the field E S (r) scattered by the targets. The inverse scattering problem in (5), i.e., finding ε δ (r ) based on the knowledge of E S (r), is nonlinear. Nevertheless, it can be recast as a linear problem by means of the Born approximation (BA). In fact, from the Lippman-Schwinger representation [57] of the scattering problem, the total field in region D can be expressed in operator form where I is the identity vector, Γ represents the operation of convolution, and Ψ represents the operation of multiplication by ε δ . If the operator norm [71] of ΓΨ is less than one, then we can expand (6) using Neumann series The first-order BA accounts only for the first term in (7), thus the total field inside the integrand of (5) can be approximated by the known incident field [64], yielding Therefore, the corresponding inverse problem is cast as the inversion of the linear integral equation connecting the permittivity contrast function to the scattered field data.
The use of BA is justified by considering the following reasons.

1) Tunnels (or other targets of interest) are isolated, limited
in number, and embedded in a lossy medium. Therefore, mutual interaction between tunnels, a phenomenon neglected by BA, can be assumed negligible. 2) In general, the inhomogeneities of the soil are electrically small, and their conductivity remains low. Therefore, their scattered fields are insignificant as compared to the RF signal reirradiated by tunnels.
3) The operator norm of ΓΨ is generally < 1, due to the losses in the ground (that introduce an exponential attenuation of the fields) [68] and to the spherical spreading factor of the incident and scattered waves. 4) The main goal is to detect, localize, and approximately determine the geometry of the targets. Toward this objective, BA-based inversion algorithms preserve the information on target localization, even when objects are strong scatterers (although the quantitative description of ε δ in D is generally prejudiced) [30], [43].
The incident field can be expressed in terms of Green's functions E I n r , r t n = QG r , r t n ·â t n (9) where Q = jωμ 0 Δl t I t for an electrically small dipole or Q = −jωμ 0 A t I t for an electrically small loop. Additionally, the field received by a dipole or loop with moment directionâ r m positioned at r r m due to an equivalent (in terms of E I n ) current distribution defined inside the investigation domain D can be expressed as [28] Substituting (9) in (10), we obtain the scalar forward model for the scattered field as Equation (11) shows a linear relation between scattered field and the dielectric profile, and the reconstruction procedures presented in Section IV are based on its inversion.
However, the actual electric field measured at the receiver is the superposition of the scattered field from targets and the direct coupling between Tx and Rx. From (11) and (9), we determine the total field experienced at the receiver in time domain · G r , r t n · a t n ε δ (r )dr + Qa r m · G r r m , r t n · a t n + H r r m , r t n + N (t).
The term a r m · G(r r m , r t n ) · a t n in (12) represents the direct coupling between Tx and Rx, and it can be considered a source of deterministic clutter. The term H(r r m , r t n ) represents the nonlinear contribution due to higher order Born series terms. In the reconstruction process, H(r r m , r t n ) is considered as unknown bias (clutter) to the desired signal. The random variable N (t) can be modeled as Gaussian process with zero mean and variance equal to the noise power where B is the bandwidth, K B is the Boltzmann constant, T 0 is the environmental temperature, and the external noise figure F a can be inferred by consulting [63]. Tomography is inherently suited for noise mitigation, since it is ultranarrowband (therefore, P N is intrinsically very low), and by simply averaging n samples of the same signal, we obtain a theoretical SNR increase of √ n. Furthermore, the static contributions in (12) as functions of r r m , r t n are generally low correlated with the value of the scattered field, meaning that the combination of view and observation diversities randomizes the clutter due to the direct path coupling.
To obtain the electric field in phasor form from (12), the received field can be mixed with two coherent oscillators to retrieve the in-phase and quadrature components, since E r t n , r r m , t = Re E r t n , r r m exp(−jωt) = Re E r t n , r r m cos(ωt) + Im E r t n , r r m sin(ωt).
Hence, the real and imaginary parts are

IV. INVERSION PROCEDURES
From a mathematical point of view, the problem of finding the dielectric profile is to compute the inverse of the linear operator L in (11), connecting the unknown dielectric profile and the scattered field data. Noise and static clutter contributions are assumed sufficiently small, so that they can be considered as perturbations on the measured data.
A way to compute L −1 is to perform a numerical inversion of L. Let us collect the sampled field data in an ordered NM vector E S = {E S (r t n , r r m )} and discretize the domain region D in K voxels, each one located at position r k : The contrast dielectric permittivity can be embodied in a column vector ε δ = {ε δ (r k )} of length K, and it represents the set of unknown parameters. After this discretization, (11) can be rewritten in matrix form where L is now a matrix with dimensions NM × K and E S , ε δ are column vectors. The problem is then to invert the relation (17). Since L is generally not a square matrix, we need to consider its pseudoinverse, which we will still indicate with L −1 . Due to the usual location of Tx, Rx, and targets, L is generally ill conditioned. A common way to quantify the behavior of L is by inspection of its condition number κ. For the operator L, it is quite common to obtain values of κ above 10 6 .
This leads to artifacts in the reconstruction process, particularly exacerbated when noise (thermal, external, quantization) or clutter is impinging upon receivers.
According to the accuracy required from the system, we present four inversion strategies.
1) Levenberg-Marquardt (LM) regularization procedure. This method is relatively accurate for any environmental condition, and it is robust in presence of noise. It requires a proper choice of a regularization parameter.
2) Truncated singular-value decomposition (TSVD). This method is also relatively accurate in any scenario and fairly resistant to noise interference. This method also offers deeper insight into the physics behind the reconstruction, and the output can be easily adjusted by properly selecting the number of meaningful singular values. In addition, the number of retained singular values in TSVD plays the same role of the LM regularization parameter. 3) Back-propagation approach. This method works properly only when the operator L is relatively well conditioned. This implies that it can be used only for particular configurations and when the SNR is relatively high. However, the computational time is drastically reduced. 4) Fourier-Bojarski approach. This is the fastest inversion scheme, and it is suited for far-field probing and quasilossless Green's function. It privileges speed instead of accuracy.

A. LM Regularization
An efficient method to compute the inverse of an illconditioned matrix is by using the LM regularization procedure [2]. In this way, the dielectric-permittivity contrast is estimated asε where L H denotes the adjoint of L and β is the regularization parameter in the Tikhonov sense, which needs to be advantageously selected. Since a proper choice of β may be a difficult task, our initial guess for β is the midpoint value of the singular-value dynamic range of L. Oftentimes, it is necessary to determine β through optimizations and iterations before a meaningful, sharp, and low-blurred image is reconstructed [2], [69]. This implies that a (computationally expensive) matrix inversion for each attempt may be necessary. To accelerate this process, the SVD may be used advantageously, which is described next.

B. Truncated SVD
A more efficient way to invert the ill-conditioned matrix L is proposed in [34], [35], and [43] and takes advantage of SVD. In fact, L can be decomposed as follows: where S is a diagonal matrix containing the ordered singular values s i of L. The pseudoinverse of L can be written as Singular values associated to S −1 that are considerably large as compared to 1/s 1 represents the sensitive directions of L: along these directions, small amounts of noise or clutter in the collected electric field will lead to a large (undesired) deviation of ε δ . A way to remove this sensitivity is to consider only the first k smaller singular values of S −1 and setting to zero the remaining large ones. This strategy is commonly referred to as "truncated SVD" [66]. Accordingly, the dielectric profile can be estimated asε The SVD can also be very useful to properly dimension the β parameter in the LM method. In fact, if we rewrite (18) in terms of (19), we obtain TSVD and the SVD representation of the LM method in (22) have a remarkable feature. In both cases, once the evaluation of the singular system of L is performed (19), the reconstructed image is simply computed by a (fast) matrix multiplication as in (20). This means that a new image is obtained by varying the number of singular values (for the TSVD method) or varying β (for the LM method), and it can be computed extremely fast (often in real time).

C. Back Propagation
The approximate number of computations to perform the SVD of L is 9K 3 + 12MNK 2 [47]. When the size of the matrix is large, the evaluation of (19) may become prohibitive, and an alternative efficient strategy has to be pursued.
In this case, the contrast permittivity function can be estimated using the following approximation: Therefore, the inversion is carried out by just using the adjoint of L. This particular inversion holds when where αI represents a scaled identity matrix. This result implies that, in principle, κ(L H L) = 1. However [47] κ( Therefore, the use of the adjoint for the inversion is acceptable only when the singular values of L exhibit a limited dynamic range. In this case, an explicit formula for the solution of (23) in terms of each r k can be formulated · G H r t n , r k · a t n E S r t n , r r m . (26) Equation (26) is commonly referred to as matched filtering, migration, or back propagation [67]. This technique is suited for parallel processing [33].

D. Fourier-Bojarski Approach
A simple and fast approach (albeit less accurate) is to take advantage of the Fourier relation arising between scattered field and object shape, as discussed in the literature under the topic of diffraction tomography.
In fact, if targets and sensors are distant enough so that the EM field can be approximated as a locally lossless plane wave (normally occurring when the fields are primarily propagating as 1/r), then the forward model can be expressed as explained as follows.
We define the unit norm direction of propagation vectors aŝ l t n = −x sin θ t n cos ϕ t n −ŷ sin θ t n sin ϕ t n −ẑ cos θ t n (27) l r m =x sin θ r m cos ϕ r m +ŷ sin θ r m sin ϕ r m +ẑ cos θ r m .
Using the paraxial approximation, the transmit Green's function at the generic position r inside region D can be simplified as G t n r t n , r ∼ = exp (+jk D r t n ) exp +jk Dl t n · r 4πr t n (29) while the receive Green's function can be expressed as Therefore, for a pair of transmitters and receivers, the scattered field can be rewritten as [36] The quantity k D (l t n −l r m ) can be represented by a 3-D vector Equation (31) is then rewritten as [33] It is useful to consider a normalized version of (33) This result can be interpreted in the following way: Each collected sample E S (k mn ) returns the value of the k mn spectral component of the contrast function ε δ (r ). Theoretically, if enough samples are available to fully populate the spectral representation of ε δ (r ), the discrete function E S (k mn ) can be approximated (in the limit) as a continuous function E S (k), and (34) can be interpreted as a 3-D inverse Fourier transform of the permittivity contrast function. Therefore, an image can be reconstructed via direct Fourier transformation of (33), i.e., where the domain of integration K is the support of E S (k).
By inspection of (32), we conclude that when the sensors completely encircle the target, K is a sphere of radius 2k D . In other words, the available information of the spectral content of ε δ (r ) is limited up to the spectral component contained inside a sphere of radius 2k D . Therefore, the reconstructed image of the dielectric profile will be a low-pass-filtered version of the true image. By studying the impulse response of (34), the minimum resolution achievable using Fourier-Bojarski approach is shown to be d ∼ = λ D /3 [39]. For half-space problems, the resolution is further reduced [49].
In real cases, where a finite number of sensors are deployed (i.e., the spectral domain is undersampled) and external noise affects the measurements, the resolution is lower, and artifacts in the reconstructed image are very common. Severe smearing and blurring effects originates mainly from the invalidity of the paraxial approximation. A way to overcome this limitation is to segment the region D in smaller analysis regions and consider an inverse problem for each subregion, and the resulting subimages are concatenated to form the final image [39]. The undersampling of the spectral domain can be corrected with several approaches, such as trilinear interpolation [50] of the available samples, or using projection on convex sets [51] and smoothing of the peaks in the spectral domain to estimate the missing samples [39].

V. SIMULATIONS
We present some simulation results to test the proposed belowground imaging system. A set of six transmitters and 26 receivers, operating at 5 MHz, is placed belowground at a depth d = 0.15 m (a t = a r =x), as shown in Fig. 4. Two empty cylindrical structures with radius ρ = 1 m (representing two tunnels) are assumed to be embedded in a host medium with relative dielectric permittivity ε D = 10 and conductivity σ D = 5 × 10 −4 S/m [62] at a depth h = 25 m. The corresponding attenuation can be computed using [65] a ∼ = 4.34σ The scattered field was synthesized using the FDTD simulator GPRMAX [70]: The instantaneous scattered electric field has been correlated using (15) and (16) to retrieve the electric field in phasor form. For the sake of speed and simplicity, in (11), we selected the Green's dyad for the homogeneous medium with the same properties of the ground, i.e., [44]  This assumption is reasonable when the sensors are deployed below the air/ground interface, the frequencies involved are relatively low, and targets are assumed to reside nearly perpendicular to sensors. The deeper the scatterer, the lower the lateral mode. Nonetheless, more accurate models can take into account the distortions due to half-space or layered geometry by simply selecting the proper Green's function under a spectral form [28], [29], [44]- [46] or by numerically computing the Green's function using method of moments or fast eikonal equation solvers (subject to future research).
Direct-path coupling was mitigated thanks to the exact knowledge of the dyadic Green's function in this problem, which enabled accurate determination and subsequent cancellation. Random Gaussian noise has been added on the data according to [63]. Nevertheless, the SNR can be reduced to a desired value by opportunely sampling and averaging the received field, as discussed in Section IV.
Upon completion of the 3-D tomographic inversion (using a mesh of 1 m 3 per voxel), the horizontal section at 25-m depth (constant depth slice) is plotted. For more complex scenarios, a full 3-D image might be necessary.
Dielectric perturbations of values less than 5% of the peak value in the domain D have been neglected in the final images.
The inversion procedure using Fourier-Bojarski method has a different (coarse) voxel size equal to λ D /4, because sampling at finer discretization does not provide more information (due to resolution limitations [39]), while it dramatically increases the sparsity of the Fourier domain, leading to severe artifacts. Fig. 5 shows the reconstructed image using the (fast) Fourier-Bojarski method: Although the resolution is coarse, basic features of the two tunnels can be discerned. However, the back-propagation method (see Fig. 6) is clearly showing higher resolution while keeping the computational cost to a minimum. For high-level image reconstruction, regularized methods are paramount. In Fig. 7, the image has been reconstructed using LM method (in its SVD variant), and the regularization parameter has been empirically selected in order to achieve the sharpest solution. In Fig. 9, an image reconstructed using TSVD is shown: We assume that 10% of the total singular   values represent the sensitive directions of L, and therefore, they will be not included in the reconstruction based on TSVD (see Fig. 8). This threshold has been chosen heuristically, and it may vary according to the geometry and the SNR.  Finally, Fig. 10 shows the reconstruction starting from the scattered field data computed using the forward model (11) and the Green's function in (37). Clearly, the TSVD provides very similar reconstruction results for the two cases of FDTD and Born field data. This result supports the accuracy of BA for belowground imaging of void targets.

VI. CONCLUSION
We proposed a practical method for tunnel detection that does not require boreholes, is easily deployable, and covers relatively wide areas.
We applied diffraction tomography and inverse-scattering principles to our geometry. We proposed four simple inversion methods to reconstruct images that are suited for the environment encountered in practical situations.
Finally, we showed several reconstructed images, and in all cases, the location of the two tunnels is discernible when noise is low, while tunnel detection amid high noise is only possible when LM and TSVD methods are used.
In particular, by comparing Figs. 6 and 7 with Figs. 9 and 10, it is clear that LM and TSVD methods yield higher quality images. Therefore, if the computational load is not a consideration, these two methods should be preferred and become imperative when environmental conditions are hostile. The proposed strategy offers the following advantages. 1) RF tomography is able to surveil from local/shallow to global/deep areas of regard, and rapidly focus on specified areas, by simply changing frequency of operation and the delimitation of the investigation domain D. 2) The system is suited for both cooperative and denied scenes, where the physical presence of the human operator is hazardous.
3) The ultranarrowband sensor design and the system architecture are fiscal and manpower affordable solutions. 4) The sensor deployment is arbitrary and modular (i.e., the addition or removal of sensors does not compromise the remainder of the system). 5) The resolution of the system can be subwavelength and range independent.
However, many aspects still need deeper investigations, such as more accurate inversion models, the use of different Green's functions, improved methods of direct-path cancellation, or more considerations on the actual soil and antennas behavior.
In particular, the issue of unwanted lateral waves must be addressed in two ways: by investigating clutter-suppression techniques or by defining a suitable Green's function that accounts for this effect. We are currently pursuing further research in both directions.