Abstract
Imaging of tunnel networks under irregular terrain using RF tomography is generalized to include the possibility of magnetic dipoles (i.e., electric loops) either as transmitting or receiving devices. Forward scattering models are presented, and a generalized method for computing numerical dyadic Green’s functions is detailed. Explicit formulas for fast numerical implementation are also presented. The paper is corroborated with numerical simulations aimed at validating formulas.
1. Introduction
This work extends the theory of RF tomography for underground imaging of tunnels [1] under irregular terrains [2] to include the possibility of using either electrically small dipoles or loops as transmitters or receivers. Prior derivations of RF tomography for underground imaging [1–3] assumed only the use of electrically small dipoles. Alas, in the frequency interval between 0.1 and 10 MHz, which is appropriate for the application of imaging of tunnels, electrically small dipoles exhibit very small radiation efficiency [4]. Therefore, the novelty of this paper is the extension of the forward scattering model of RF tomography to include electrically small loops, such as multiturn ferrite loaded coils, as either transmitters or receivers, because they are generally more efficient.
This paper is organized as follows. Section 2 derives a comprehensive linearized forward model used for RF tomography for the electromagnetic field generated from both electric and magnetic current sources; Section 3 describes a method to obtain numerical Green’s functions based upon a generalization of the method of moments; Section 4 shows some relevant simulations, and Section 5 draws conclusions. In Appendix, explicit formulas to compute the half-space Green’s functions, their evaluation through integral over cubic regions, and their implementation for fast processing are provided.
2. Forward Model
2.1. Notations and Conventions
In this paper, the time dependence of the fields is assumed. Bolded lower case letters represent vectors, while bolded capital case letters represent matrices. Underscored capital letters represent dyadics. Circumflexed bold letters represent unit vectors. To keep formulas succinct, expressions separated by a semicolon represent two separate equations, for instance, the expression corresponds to two separate equations: and .
2.2. Mathematical Description of the Terrain
Consider the three-dimensional geometry shown in Figure 1. The interface between air and ground is not planar but its shape is assumed known. This information can be used to build a numerical Green’s function that accounts for the terrain shape, leaving only the underground targets as unknown of the problem. The numerical Green’s function of the irregular terrain is computed by applying a generalized method of moment based on half-space Green’s functions, as described later. In other words, the final aim is the reconstruction (in the case at hand in terms of localization and shape estimation) and not the shape of the interface air/soil. To do this appropriately, it is convenient to solve an inverse problem where the buried target is an anomaly with respect to a known background scenario; in this context, the more a priori information is available for the background scenario (e.g., the shape of the interface), the better the target reconstruction will be. Hence, the accurate computation of the effects of the irregular surface (through the determination of the numerical Green’s function) is pivotal to provide more reliable images of the underground scene.
It is beneficial to express the electrical properties of the irregular terrain as the superposition of two contributions. The first contribution is associated with an ideal planar half-space with an interface at , whose equivalent relative dielectric permittivity is described by a scalar function , the half-space equivalent dielectric permittivity function, defined as where is the angular frequency, is the frequency-independent background relative dielectric permittivity of the ground, is the frequency-independent background conductivity of the ground, and is a position vector used to describe the terrain. The second contribution is a volumetric function with finite support, denoted as the background contrast function, which represents the deviation in dielectric permittivity and conductivity at point between the actual values and the values expected from the ideal half-space geometry . The background contrast function is mathematically described as where , are the actual relative dielectric permittivity and the actual conductivity of the terrain within the support of . Note that returns the actual electrical properties of the irregular terrain.
2.3. Investigation Domain
The targets (i.e., tunnels or voids) are assumed to reside within the investigation domain D, fully belonging to the lower medium. To ease the mathematical derivation, the region D shall not include any point in which the background contrast function differs from zero. A point in region D shall be represented with a separate position vector .
For a target-free scenario, the background electrical properties in D are fully described by the quantities . However, the presence of any subsurface structure creates a deviation in dielectric permittivity and conductivity with respect to the background. The actual relative dielectric permittivity distribution and conductivity distribution inside the investigation domain D are the unknowns of the inverse problem. The inverse problem is recast in terms of an object contrast function , which is more physically related to the presence/absence of a target, defined as
In fact, when , a dielectric or conducting anomaly is present at location .
2.4. Transmitters and Receivers
In the HF region, deployable transmitters and receivers can be either electrically small dipoles or loops: the radiation characteristic of these antennas can be easily incorporated in the forward model, since both electrically small dipoles and electrically small loops can be modeled as Hertzian electric/magnetic dipoles. At the particular n-th measurement of the field, the transmitter is located at position and directed along the unit-norm vector , while the receiver is located at position and direction . At the receiver side, the measured field is if the receiver is a magnetic dipole, or if an electric dipole acts as a receiver.
2.5. Derivation of the Forward Model
For a fixed frequency , the incident electric and magnetic fields observed at a generic point inside due to elementary electric/magnetic sources distributed arbitrarily above/below an irregular terrain at locations can be expressed in time-harmonic form as In (4), is the electric dyadic Green’s function due to an electric source, is the electric dyadic Green’s function due to a magnetic source, is the magnetic dyadic Green’s function due to an electric source, and is the magnetic dyadic Green’s function due to a magnetic source. Mathematically, these Green’s functions are the solutions of the following differential equations, with the appropriate radiation condition at infinity and the wavenumber can be evaluated using Note that (4) is valid for any arbitrary distribution of ; therefore, these Green’s functions represent the solutions of (5) for the most general case.
The presence of a nonzero contrast function inside generates an equivalent current density equal to The approximation in (7) is called “Born Approximation” [5], and it is generally valid when only a qualitative reconstruction (in terms of presence, location, and geometry) of the underground electrical properties is sought [6].
The approximated scattering integral equation for the problem becomes So far, the solutions of (5) are known in closed form only for very limited distributions of , such as the homogenous space, half-space, or planarly layered media [5, 7]. The general solution for any arbitrary distribution of , such as a half-space having irregular boundary, must be sought numerically.
Combining (7) with (8), the contrast function can be linearly related to the total (scalar) complex-valued (phasor) electric field collected at the receiver due to a single transmitter of length via the following relations.(1)Tx electric dipole, Rx electric dipole (2)Tx magnetic dipole, Rx magnetic dipole (3)Tx electric dipole, Rx magnetic dipole (4)Tx magnetic dipole, Rx electric dipole These formulas represent the forward scattering model for RF tomography, valid for any arbitrary distribution of the background electrical properties. In (9)–(12), T denotes the operation of transposition, and represents an (unknown) error vector accounting for high-order Born series terms, clutter, and noise. Because of the first-order Born approximation, the validity of the reconstruction is limited to low/moderate-contrast targets.
Discretization and inversion of the forward model can be accomplished in several ways. In this work, the method described in [2] will be followed and, due to limited space, it will not be repeated here.
By using the principle of superposition, the general Green’s functions in (9)–(12) can be separated into two contributions: where is one of . The functions represent the Green’s function specific for the half-space geometry (having planar interface at = 0), whereas represent the scattered contribution due to the presence of the background contrast function defined in (2). Analytical expressions of the half-space Green’s function do exist in several forms [5, 7–12]; however, in Appendix, we propose better formulas that are explicit, easy to be integrated, and computed with FFT, and having very few singularities (to improve numerical stability). To employ RF tomography for the imaging of any irregular surface, needs to be computed numerically, possibly in a fast and efficient manner.
3. Numerical Green’s Functions
The application of RF tomography to imaging below irregular terrain according to (9) requires the numerical computation of the Green’s function of the problem, which is discussed in this section. Since the topic is more general than its application to RF tomography, the problem of finding the value of the numerical Green’s function observed at point , due to a source located at , having wavenumber is described first.
The following sections extend the work proposed in [2], which was itself a generalization of [1, 7–16], to the general case of electric/magnetic sources and electric/magnetic probes.
3.1. Method of Moments
The scattering integral equations in vector form for the distributed dielectric anomaly can be easily derived from Maxwell’s equation and are stated below In this context, , are the scattered electric and magnetic fields, and is the total electric field. The next goal is to generalize (14) and (15) to the case of Green’s functions. This step is shown in [2], and accordingly we have for (14) and for (15) To generalize the result, the operation is disregarded. Substituting (16) in (14), the dyadic scattering integral equation becomes Similarly, In these formulas, is the support of . Since and the incident electric field is we obtain that Part of (18)-(19) can be computed analytically. Therefore, the left-hand side of (18)-(19) is divided into two dyadic functions: a known part and an unknown part Retrieving from (27)–(30) is a hard task, so that an approximation is needed. The computation is facilitated by discretizing the region R in P small cubes (voxels), representing an orthogonal basis expansion for . Each cube is centered at position , whose volume is defined as where is the half length of the edge of elementary cube. The dimension of the voxel is chosen so that and can be assumed constant within the boundaries of the voxel itself. Theoretically, could represent any kind of orthogonal basis, such tetrahedral voxels; however, the choice of “cubes” will be particularly useful for fast computation, as shown later.
Accordingly, (23) and (27) can be discretized in an unknown part and a known part where the index is used for the same purpose of the index , but is distinguished to emphasize that the summations are computed separately. In (32) and (33) the dyadic expression represents the integral of the half-space Green’s function, performed over a cubic region having volume , centered along the source point , and computed at the observation point . The exact evaluation of this integral is described in Appendix.
The integral equations (18)-(19) can be discretized as Only (34) and (35) can be solved directly. Equations (36) and (37) will be computed after the previous two equations are solved.
Consider (34) and (35), assuming , in these equations, would be determined when the value of is known for any value of . To find for all , (34) or (35) are tested exactly at the locations , leading to the following set of equations: Due to the dyadic nature of (38), equations with unknowns can be defined, leading to the creation of the following matrix equation: where, By solving for , the knowledge of for all follows. Direct substitution of this information into (38) yields the full reconstruction of the scattering Green’s function at any given point , that is,
The newly acquired knowledge of is also useful to determine the remaining types of scattering Green’s functions The final numerical Green’s dyadic can be computed by substituting (41), (42), (43), or (44) in (13). This procedure of computing numerical Green’s functions can be interpreted as a generalization of the works in [17–19] to the case of 3D half-space geometry.
4. Simulations
A large-scale numerical simulation has been performed to prove the validity of this novel method. The conditions of the experiments are similar to the ones reported in [2]. Accordingly, a half-space geometry is considered and, to simulate irregular terrain, a parallelepiped of size m with the same properties of the ground is placed on top of the flat surface at the center of the scene (see [2] Figure 2). Five transmitters and receivers are distributed uniformly along a circle of radius 15 m, at 0.5 m above the box, to avoid numerical instabilities. The target is an inverse-L shaped tunnel of equal sides and radius 1 m, located at depth of 10 m, as shown in Figure 2. The frequency range is 3–6 MHz, at constant intervals of 0.5 MHz. The electrical properties of the soil are chosen to be similar to dry rocks; accordingly, the background dielectric permittivity is and the background dielectric conductivity is S/m. To avoid the so-called “inverse-problem crime,” the exact electric field at the receiver side has been computed numerically using the finite difference time domain (FDTD) code GPRMAX [20]. For the FDTD solver, the discretization step is m, and the box size containing the whole scene is .
For each Tx-Rx-frequency combination, two simulations were performed: one containing the target (i.e., the total field) and one without (i.e., the background). The scattered field in complex form can be obtained by subtracting these two time-domain data sets, followed by an FFT (see [1] for details).
To fill , the algorithm proposed in this paper was applied. The irregular surface (i.e., the flat surface with the box on top) was discretized in cubes having m, corresponding to , which is moderately coarse but manageable in terms of memory allocation. Better performances is expected using smaller , but more sophisticated and memory-saving algorithms required to fill and invert . For m, the approximation (A.43) holds true within an error less than 5% occurring only to the cubes adjacent to the singularity; this value can be considered satisfactory. The region of interest D is a horizontal slice located at depth −10 m, discretized using a cubic mesh of half size of m. The direct path was eliminated a priori since its value was accurately estimated using (9)–(12). The data is noiseless; however, discrepancies between the numerical FDTD solver and the proposed MoM-based method can be assumed as unknown perturbation of the measured field, thus implicitly testing the robustness of our algorithm with respect to disturbances.
The reconstruction has been performed using a modified FISTA algorithm; details can be found in [2–16, 21].
Assuming electric dipoles as sources and receivers, the reconstructed image of the belowground surface using the proposed algorithm is shown in Figure 3. The same irregular geometry, parameters, and numerical algorithm are used for the results in Figures 4–6; however, in Figure 4, electric dipoles are used as Tx, and magnetic dipoles are used as Rx, in Figure 5 magnetic dipoles are used as Tx, and electric dipoles are used as Rx, and in Figure 6, both Tx and Rx are magnetic dipoles. To reconstruct these images, all formulas presented in this work had to be used. By induction, the validity of these formulas has been verified.
5. Conclusions
In summary, a complete forward model for RF tomography with applications to imaging under irregular terrains has been presented. The novelties of this approach can be summarized as follows.(i)The method includes the possibility of using both dipoles and loops as transmitters or receivers, at arbitrary locations (either above or below surface) and directions.(ii)The model accurately accounts for the radiation pattern of the antenna.(iii)The model is derived from the vector wave equation, thus preserving the vector nature of the and fields, as well as accounting for changes of polarization.(iv)The method fully accounts for the irregular shape of the terrain.(v)The method accurately describes scattering processes occurring both in the near field and in the far field.(vi)The method is intrinsically parallel, and most of the computation can be done using FFT routines.For completeness, limitations and challenges are reported as well.(i)The dimension of matrix and vector grows fast when the number of voxels in or increases, thus dramatically limiting the performance for highly irregular terrains or large investigation domains.(ii)The terrain is assumed flat after a certain distance from the investigation domain; the model does not account for leveled grounds.(iii)The model assumes homogeneous electrical properties of the ground. Realistic soils are hardly homogeneous, and improved Green’s functions are desired.(iv)The discretization process is performed using equal-size cubes, to ease the integration of the Green’s functions. Other basis functions that adequately approximate the irregular shape of the terrain might be used.These challenges are currently under investigation.
Appendix
A. Half-Space Green’s Functions
Chew and Cui presented an approach to compute half-space dyadic Green’s function based on vector wave functions expressed in Cartesian form [7–12, 22–34] and several manipulations to obtain less singular solutions. In this work, their spirit is used to explicitly calculate the half-space dyadic Green’s functions needed in RF tomography.
The half-space dyadic Green’s functions in terms of electric and magnetic sources shall be defined as The actual expressions of the dyadic Green’s functions depend upon the location of the observation and source point; assuming that the source point is at the layer , and the observation point is at the layer , the Green’s functions are computed as Physically, the subscript represents the “secondary” field contribution, accounting for either reflected or transmitted fields through layers, while the subscript represents the “primary” field contribution, occurring only within the same layer. The secondary contribution to the electric layered Green’s function of electric type can be computed using the following expression [12]: In the above, , and implies that the curl operator is applied to the primed variables. Equation (A.3) is valid only for . In (A.3), the and functions represent the generalized Fresnel coefficients for the TE and TM waves, respectively, expressed in spectral form. Explicitly, The functions describe the variation of (A.4) and are solutions of the following differential equation [12]: Although the general solution of (A.5) is not trivial, four simple expressions are obtained when the number of layers is exactly two. For , the magnetic layered Green’s function of electric type can be easily computed from the electric layered Green’s function of electric type as follows: The secondary contribution of the magnetic layered Green’s function of magnetic type can be calculated by invoking duality. Explicitly and the half-space electric Green’s function of magnetic type follows as When the source and observation points are in the same layer , the primary field contributions need to be added: By following the steps (A.3)–(A.9), the resulting Green’s functions are still relatively singular, hence difficult to be computed numerically. However, for the particular case of half-space geometry, there exists a strategy to reduce the singularity of the Green’s functions (and its integral over a cubic region), so that their numerical computation become fast and stable.
For half-space geometry, when the source and observation points belong to different layers, a straightforward mathematical manipulation is sufficient to isolate and simplify the singularities occurring in the expressions. Conversely, when source and observation points belong to the same layer, the singularity reduction can be achieved by rewriting the corresponding Green’s function as follows: In the above, is the closed-form expression of the homogeneous Green’s function due to the image source: its value can be calculated by substituting in , where is the location of the image source. The function is the spectral representation of the homogeneous Green’s function due to the image source: the homogeneous spectral Green’s function can be calculated from (A.3), (A.6), (A.7), and (A.8) as limit case that , and substituting . In (A.10), the quantity represents the Green’s function contribution due to the superposition of the fields generated at the interface and the fields due to the image source. Mathematically, is less singular than . Therefore, the expression is less singular than (A.2). Using multiple substitutions and some tedious math, singular-free and simplified solutions can be derived, as shown in the following subsections.
Here, the steps described in this appendix have been performed, and the resulting equations have been simplified, uniformed, and tested using a specifically designed version of the freeware FDTD code GPRMAX [20] that accepts magnetic dipoles as sources. Intermediate derivations are omitted, but the final formulas are presented in explicit form, ready to be implemented in routines.
For each Green’s function, the superscript indicates that the source point is on the layer , and the observation point is on the layer ; being a half-space problem, the two possible layers are either “a” (air) or “e” (earth).
The following notations will be used:
A.1. Electric Green’s Function of Electric Type
The electric Green’s function of electric type relates electric fields with electric current sources in a half-space geometry, according to (A.1). The homogeneous electric Green’s function of electric type, both in air or earth, is The homogeneous electric Green’s function of electric type due to the image source is The electric half-space Green’s functions of electric type for source and observation points both in the air (or in the earth) can be written as follows: where The electric half-space Green’s function of electric type for a source point in the earth and observation point in the air is The electric half-space Green’s function of electric type for a source point in the air and observation point in the earth is
A.2. Magnetic Green’s Function of Electric Type
The magnetic Green’s function of electric type relates magnetic field to electric currents in a half-space geometry, according to (A.1).
The homogeneous magnetic Green’s function of electric type for the air and earth cases is given The homogeneous magnetic Green’s function of electric type for the air and earth cases due to the image source is given by The air-air and earth-earth magnetic half-space Green’s function of electric type can be written as follows: where The magnetic half-space Green’s function of electric type for source point in the earth and observation point in the air is The magnetic half-space Green’s function of electric type for source point in the air and observation point in the earth is
A.3. Magnetic Green’s Function of Magnetic Type
The homogeneous magnetic Green’s function of magnetic type for the air case and the earth case are both for the real and image sources The magnetic half-space Green’s functions of magnetic type for source and observation points both in air (or earth) are where The magnetic half-space Green’s function of magnetic type for source point in the earth and observation point in the air is The magnetic half-space Green’s function of magnetic type for source point in the air and observation point in the earth is
A.4. Electric Green’s Function of Magnetic Type
The electric Green’s function of magnetic type relates electric fields with magnetic currents in a half-space geometry, according to (A.1). The homogeneous electric Green’s function of magnetic type is The homogeneous electric Green’s function of magnetic type due to the image source is: The electric half-space Green’s functions of magnetic type for source and observation points both in air (or in earth) are where The electric half-space Green’s function of magnetic type for source point in earth and observation point in air is The electric half-space Green’s function of magnetic for source point in air and observation point in earth is
A.5. Integral of Green’s Functions
Formulas for the integral of half-space Green’s functions along a cubic region are provided. These formulas are necessary to compute the functions , whenever they appear in (34)–(37). The integral of the half-space Green’s function is intentionally expressed in a 2D spectral-like form, so that 2D FFT algorithms are applicable. The intermediate steps to obtain these formulas are omitted, although relevant citations have been properly distributed throughout this subsection. The integral to be computed is of the form This integral has a form that is independent from the type, but depends upon the location of the source and observation points.
The simplest cases are the air-earth and earth-air cases, having, respectively, the following form: In these expressions, the only undetermined point can be explicitly computed using De L’Hôpital rule When sources are both in air (or earth) the integration becomes too difficult to be computed in an exact form, and some approximations are needed. The general form is The spectral contribution is easily computed The integral of the homogeneous part has two main problems.(i)Formulas contain the term which is difficult to be integrated over a cubic region, but easier to be integrated over a spherical region.(ii)When , the integration contains a singular point.To address these issues, the work done by Gao and Torres-Verdin [35] is helpful. Their main idea is that the integral over a cubic region of the homogeneous Green’s function can be approximated to the integral over a spherical region of the homogeneous Green’s function, having volume equal to the one of the cube. Using this approximation, one can show that the primary field contribution is where .
Care is needed in applying this approximation, since its validity fails for . The range of validity is a nonlinear function of both and . It is advised to test the error of this approximation by comparing (A.43) with a thorough numerical integration of the half-space Green’s function over a cube including the singularity.
If cube, then the integration is performed in a volume that contains the source (i.e., the singular point). In this case, the result obtained by Yaghjian in [36], together with the duality theorem, shall be used, yielding Within the validity of the forward model, the image source can be considered always outside of the region of integration; hence, the integral of the homogeneous Green’s function due to the image source is simply
A.6. Application to RF Tomography
In the procedure described in [1, 2], it is generally unlikely that Green’s functions need to be computed at sporadic points; in general, Green’s functions need to be calculated along horizontal layers, or the entire 3D space. To this aim, rather than computing the double integration in variables using classical numerical solvers, a fast FFT-based method can be used. Using 2D FFT, Green’s functions and their integral can be computed along predetermined horizontal layers in “one shot,” after properly filling the corresponding 2D Fourier space.
More exactly, by keeping the variable and fixed, the spectral part of the Green’s functions or their integral can be recast in the following form: The generic dyadic function corresponds to either or functions, and corresponds to either or deprived from the term .
Physically, this quantity represents the value of the dyadic function computed at the horizontal plane given that the source is fixed. Numerically, each component of this integral can be computed extremely fast by using 2D FFT. Note that each component is also independent, so 9 FFT can be executed in parallel. This particular function shall be called “Dyadic plane.” Similarly, by keeping the observation point and the variable fixed, we obtain the following dyadic plane: Physically, this operation represents the computation of the dyadic expression spanning all source points located at the horizontal plane at height . Numerically, this can be easily computed using 2D IFFT.
Due to the properties of Fourier integral, an important symmetric property holds Hence, where These properties will be extremely useful for the fast computation of numerical Green’s functions in RF tomography. In fact, to numerically compute the proposed forward model in [2], only the following dyadic planes need to be computed
By using the shifting property in (A.49), the computation of the FFT or IFFT needs to be performed basically only once per plane. All these planes can be computed independently and in parallel, thus dramatically accelerating the computational time.
Acknowledgments
The authors are thankful to Timothy Poth, Gary Scalzi, Braham Himed, and Johnathan Scanlan, Air Force Research Laboratory, for sponsoring and funding this research. This work has been sponsored by the US Air Force Research Laboratory, under Contract no. FA-8650-10-2-7028.