# A Comprehensive Forward Model for Imaging under Irregular Terrain Using RF Tomography

Authors: Lorenzo Lo Monte; Francesco Soldovieri; Danilo Erricolo; Antonios Giannopoulos; Michael C. Wicks
Journal: International Journal of Antennas and Propagation
 (2012)

---

## 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.


---

## Body

## 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. Section2 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, theexp(-iωt) 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 Ab;c=(+;-)x corresponds to two separate equations: Ab=+x and Ac=-x.


### 2.2. Mathematical Description of the Terrain

Consider the three-dimensional geometry shown in Figure1. 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.Figure 1
Sample geometry where a portion of irregular surface is shown above the planez=0. An anomaly, represented by two cylindrical structures, is located underground in the investigation region D.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 atz=0, whose equivalent relative dielectric permittivity is described by a scalar function H:ℝ3→ℂ, the half-space equivalent dielectric permittivity function,defined as
(1)H(r~)={εD+iσD(ωε0)z~<0,1z~>0,
where ω=2πf is the angular frequency, εD is the frequency-independent background relative dielectric permittivity of the ground, σD is the frequency-independent background conductivity of the ground, and r~=x~x^+y~y^+z~z^ is a position vector used to describe the terrain. The second contribution is a volumetric function Q:ℝ3→ℂ with finite support, denoted as the background contrast function, which represents the deviation in dielectric permittivity and conductivity at point r~ between the actual values and the values expected from the ideal half-space geometry H(r). The background contrast function is mathematically described as
(2)Q(r~)={εr(r~)-εD+i(σ(r~)-σD)(ωε0)z~<0,εr(r~)-1+iσ(r~)(ωε0)z~>0,
where εr(r~), σ(r~) are the actual relative dielectric permittivity and the actual conductivity of the terrain within the support of Q. Note that H(r~)+Q(r~)  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 domainD,fully belonging to the lower medium. To ease the mathematical derivation, the region Dshall not include any point in which the background contrast function Q(r~) differs from zero.A point in region Dshall be represented with a separate position vector r-=x-x^+y-y^+z-z^.For a target-free scenario, the background electrical properties inDare fully described by the quantities εD,σD. 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 εrO(r-) and conductivity distribution σO(r-) inside the investigation domain D are the unknowns of the inverse problem. The inverse problem is recast in terms of an object contrast functionV:ℝ3→ℂ, which is more physically related to the presence/absence of a target, defined as
(3)V(r-)=εrO(r-)-εD+i(σO(r-)-σD)(ωε0).In fact, whenV(r-)≠0, a dielectric or conducting anomaly is present at location r-.


### 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 particularn-th measurement of the field, the transmitter is located at position rna and directed along the unit-norm vector a^n, while the receiver is located at position rnb and direction b^n. At the receiver side, the measured field is H if the receiver is a magnetic dipole, or E if an electric dipole acts as a receiver.


### 2.5. Derivation of the Forward Model

For a fixed frequencyf, the incident electric and magnetic fields observed at a generic point r inside D due to elementary electric/magnetic sources distributed arbitrarily above/below an irregular terrain at locations r′ can be expressed in time-harmonic form as
(4)[EH](r)=∭D[iωμ0G_eeIrrG_emIrrG_meIrriωε0G_mmIrr]·[JM](r′)dr′.
In (4), G_eeIrr is the electric dyadic Green’s function due to an electric source, G_emIrr is the electric dyadic Green’s function due to a magnetic source, G_meIrr is the magnetic dyadic Green’s function due to an electric source, and G_mmIrr 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
(5)∇×∇×G_uIrr(r,r′)-k2(r)G_uIrr(r,r′)={I_δ(r-r′)u=ee;mm∇×I_δ(r-r′)u=em;me
and the wavenumber k(r) can be evaluated using
(6)k(r)=ωμ0(ε0εr(r,ω)+iσ(r,ω)(ωε0)).
Note that (4) is valid for any arbitrary distribution of k(r); therefore, these Green’s functions represent the solutions of (5) for the most general case.The presence of a nonzero contrast function insideD generates an equivalent current density equal to
(7)JInc(r-)=-iωε0V(r-)·[EInc(r-)+ESca(r-)]≅-iωε0V(r-)·EInc(r-).
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(8)E(rnb)=iωμ0∭DG_eeIrr(rnb,r-)·Jeq(r-)dr-,H(rnb)=∭DG_meIrr(rnb,r-)·Jeq(r-)dr-.
So far, the solutions of (5) are known in closed form only for very limited distributions of k(r), such as the homogenous space, half-space, or planarly layered media [5, 7]. The general solution for any arbitrary distribution of k(r), such as a half-space having irregular boundary, must be sought numerically.Combining (7) with (8), the contrast function V(r-) can be linearly related to the total (scalar) complex-valued (phasor) electric field collected at the receiver due to a single transmitter of length ln via the following relations.(1)
Tx electric dipole, Rx electric dipole(9)En(rna,a^n,rnb,b^n,ωn)=+iωnμ0k02Ineln×∭Db^nT·G_eeIrr(rnb,r-)·G_eeIrr(r-,rna)·a^nV(r-)dr-+iωnμ0Inelnb^nT·G_eeIrr(rnb,rna)·a^n+b^nT·Z.(2)
Tx magnetic dipole, Rx magnetic dipole(10)Hn(rna,a^n,rnb,b^n,ωn)H=-iωε0k0n2Inhln×∭Db^nT·G_meIrr(rnb,r-)·G_emIrr(r-,rna)·a^nV(r-)dr-+iωε0Inhlnb^nT·G_mmIrr(rnb,rna)·a^n+b^nT·Z.(3)
Tx electric dipole, Rx magnetic dipole(11)Hn(rna,a^n,rnb,b^n,ωn)=k0n2Ineln×∭Db^nT·G_meIrr(rnb,r-)·G_eeIrr(r-,rna)·a^nV(r-)dr-+Inelnb^mT·G_meIrr(rnb,rna)·a^n+b^nT·Z.(4)
Tx magnetic dipole, Rx electric dipole(12)En(rna,a^n,rmb,b^m,ωn)=k0n2Inhln×∭Db^nT·G_eeIrr(rnb,r-)·G_emIrr(r-,rna)·a^nV(r-)dr-+Inhlnb^nT·G_emIrr(rnb,rna)·a^n+b^nT·Z.These formulas represent the forward scattering modelfor RF tomography, valid for any arbitrary distribution of the background electrical properties. In (9)–(12), T denotes the operation of transposition, and Z 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:
(13)G_uIrr(r,r′)=G_u(r,r′)+G_uS(r,r′),
where u is one of ee;em;me;mm. The functions G_u(r,r') represent the Green’s function specific for the half-space geometry (having planar interface at z = 0), whereas G_uS(r,r') represent the scattered contribution due to the presence of the background contrast function Q defined in (2). Analytical expressions of the half-space Green’s function G_(r,r') 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, G_uS(r,r′) needs to be computed numerically, possibly in a fast and efficient manner.


## 2.1. Notations and Conventions

In this paper, theexp(-iωt) 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 Ab;c=(+;-)x corresponds to two separate equations: Ab=+x and Ac=-x.


## 2.2. Mathematical Description of the Terrain

Consider the three-dimensional geometry shown in Figure1. 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.Figure 1
Sample geometry where a portion of irregular surface is shown above the planez=0. An anomaly, represented by two cylindrical structures, is located underground in the investigation region D.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 atz=0, whose equivalent relative dielectric permittivity is described by a scalar function H:ℝ3→ℂ, the half-space equivalent dielectric permittivity function,defined as
(1)H(r~)={εD+iσD(ωε0)z~<0,1z~>0,
where ω=2πf is the angular frequency, εD is the frequency-independent background relative dielectric permittivity of the ground, σD is the frequency-independent background conductivity of the ground, and r~=x~x^+y~y^+z~z^ is a position vector used to describe the terrain. The second contribution is a volumetric function Q:ℝ3→ℂ with finite support, denoted as the background contrast function, which represents the deviation in dielectric permittivity and conductivity at point r~ between the actual values and the values expected from the ideal half-space geometry H(r). The background contrast function is mathematically described as
(2)Q(r~)={εr(r~)-εD+i(σ(r~)-σD)(ωε0)z~<0,εr(r~)-1+iσ(r~)(ωε0)z~>0,
where εr(r~), σ(r~) are the actual relative dielectric permittivity and the actual conductivity of the terrain within the support of Q. Note that H(r~)+Q(r~)  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 domainD,fully belonging to the lower medium. To ease the mathematical derivation, the region Dshall not include any point in which the background contrast function Q(r~) differs from zero.A point in region Dshall be represented with a separate position vector r-=x-x^+y-y^+z-z^.For a target-free scenario, the background electrical properties inDare fully described by the quantities εD,σD. 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 εrO(r-) and conductivity distribution σO(r-) inside the investigation domain D are the unknowns of the inverse problem. The inverse problem is recast in terms of an object contrast functionV:ℝ3→ℂ, which is more physically related to the presence/absence of a target, defined as
(3)V(r-)=εrO(r-)-εD+i(σO(r-)-σD)(ωε0).In fact, whenV(r-)≠0, a dielectric or conducting anomaly is present at location r-.


## 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 particularn-th measurement of the field, the transmitter is located at position rna and directed along the unit-norm vector a^n, while the receiver is located at position rnb and direction b^n. At the receiver side, the measured field is H if the receiver is a magnetic dipole, or E if an electric dipole acts as a receiver.


## 2.5. Derivation of the Forward Model

For a fixed frequencyf, the incident electric and magnetic fields observed at a generic point r inside D due to elementary electric/magnetic sources distributed arbitrarily above/below an irregular terrain at locations r′ can be expressed in time-harmonic form as
(4)[EH](r)=∭D[iωμ0G_eeIrrG_emIrrG_meIrriωε0G_mmIrr]·[JM](r′)dr′.
In (4), G_eeIrr is the electric dyadic Green’s function due to an electric source, G_emIrr is the electric dyadic Green’s function due to a magnetic source, G_meIrr is the magnetic dyadic Green’s function due to an electric source, and G_mmIrr 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
(5)∇×∇×G_uIrr(r,r′)-k2(r)G_uIrr(r,r′)={I_δ(r-r′)u=ee;mm∇×I_δ(r-r′)u=em;me
and the wavenumber k(r) can be evaluated using
(6)k(r)=ωμ0(ε0εr(r,ω)+iσ(r,ω)(ωε0)).
Note that (4) is valid for any arbitrary distribution of k(r); therefore, these Green’s functions represent the solutions of (5) for the most general case.The presence of a nonzero contrast function insideD generates an equivalent current density equal to
(7)JInc(r-)=-iωε0V(r-)·[EInc(r-)+ESca(r-)]≅-iωε0V(r-)·EInc(r-).
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(8)E(rnb)=iωμ0∭DG_eeIrr(rnb,r-)·Jeq(r-)dr-,H(rnb)=∭DG_meIrr(rnb,r-)·Jeq(r-)dr-.
So far, the solutions of (5) are known in closed form only for very limited distributions of k(r), such as the homogenous space, half-space, or planarly layered media [5, 7]. The general solution for any arbitrary distribution of k(r), such as a half-space having irregular boundary, must be sought numerically.Combining (7) with (8), the contrast function V(r-) can be linearly related to the total (scalar) complex-valued (phasor) electric field collected at the receiver due to a single transmitter of length ln via the following relations.(1)
Tx electric dipole, Rx electric dipole(9)En(rna,a^n,rnb,b^n,ωn)=+iωnμ0k02Ineln×∭Db^nT·G_eeIrr(rnb,r-)·G_eeIrr(r-,rna)·a^nV(r-)dr-+iωnμ0Inelnb^nT·G_eeIrr(rnb,rna)·a^n+b^nT·Z.(2)
Tx magnetic dipole, Rx magnetic dipole(10)Hn(rna,a^n,rnb,b^n,ωn)H=-iωε0k0n2Inhln×∭Db^nT·G_meIrr(rnb,r-)·G_emIrr(r-,rna)·a^nV(r-)dr-+iωε0Inhlnb^nT·G_mmIrr(rnb,rna)·a^n+b^nT·Z.(3)
Tx electric dipole, Rx magnetic dipole(11)Hn(rna,a^n,rnb,b^n,ωn)=k0n2Ineln×∭Db^nT·G_meIrr(rnb,r-)·G_eeIrr(r-,rna)·a^nV(r-)dr-+Inelnb^mT·G_meIrr(rnb,rna)·a^n+b^nT·Z.(4)
Tx magnetic dipole, Rx electric dipole(12)En(rna,a^n,rmb,b^m,ωn)=k0n2Inhln×∭Db^nT·G_eeIrr(rnb,r-)·G_emIrr(r-,rna)·a^nV(r-)dr-+Inhlnb^nT·G_emIrr(rnb,rna)·a^n+b^nT·Z.These formulas represent the forward scattering modelfor RF tomography, valid for any arbitrary distribution of the background electrical properties. In (9)–(12), T denotes the operation of transposition, and Z 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:
(13)G_uIrr(r,r′)=G_u(r,r′)+G_uS(r,r′),
where u is one of ee;em;me;mm. The functions G_u(r,r') represent the Green’s function specific for the half-space geometry (having planar interface at z = 0), whereas G_uS(r,r') represent the scattered contribution due to the presence of the background contrast function Q defined in (2). Analytical expressions of the half-space Green’s function G_(r,r') 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, G_uS(r,r′) 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 r, due to a source located at r′, having wavenumber k0 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 anomalyQ(r) can be easily derived from Maxwell’s equation and are stated below
(14)ES(r)=+k02∭DQ(r~)G_ee(r,r~)·ETot(r~)dr~,(15)HS(r)=-iωε0∭DQ(r~)G_me(r,r~)·ETot(r~)dr~.
In this context, ES, HS are the scattered electric and magnetic fields, and ETot 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)
(16)ES(r)=iωμ0G_eeS(r,r′)·a^,ETot(r~)=iωμ0G_eeTot(r~,r′)·a^,
and for (15)
(17)ETot(r~)=G_emTot(r~,r′)·a^,HS(r)=G_meS(r,r′)·a^.
To generalize the result, the ·a^ operation is disregarded. Substituting (16) in (14), the dyadic scattering integral equation becomes
(18)G_eeS(r,r′)=k02∭RQ(r~)G_ee(r,r~)·G_eeTot(r~,r′)dr~.
Similarly,
(19)G_emS(r,r′)=k02∭RQ(r~)G_ee(r,r~)·G_emTot(r~,r′)dr~,G_meS(r,r′)=k02∭RQ(r~)G_me(r,r~)·G_eeTot(r~,r′)dr~,G_mmS(r,r′)=-∭RQ(r~)G_me(r,r~)·G_emTot(r~,r′)dr~.
In these formulas, R is the support of Q. Since
(20)ETot(r~)=EInc(r~)+ES(r~),
and the incident electric field is
(21)EInc(r~)={iωμ0G_ee(r~,r′)·a^electric  currents,G_em(r~,r′)·a^magnetic  currents,
we obtain that
(22)G_ee;emTot(r~,r′)=G_ee;emS(r~,r′)+G_ee;em(r~,r′).
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
(23)G_eeK(r,r′)=k02∭RQ(r~)G_ee(r,r~)·G_ee(r~,r′)dr~,(24)G_meK(r,r′)=k02∭RQ(r~)G_me(r,r~)·G_ee(r~,r′)dr~,(25)G_emK(r,r′)=k02∭RQ(r~)G_ee(r,r~)·G_em(r~,r′)dr~,(26)G_mmK(r,r′)=-∭RQ(r~)G_me(r,r~)·G_em(r~,r′)dr~,
and an unknown part
(27)G_eeU(r,r′)=k02∭RQ(r~)G_ee(r,r~)·G_eeS(r~,r′)dr~,(28)G_meU(r,r′)=k02∭RQ(r~)G_me(r,r~)·G_eeS(r~,r′)dr~,(29)G_emU(r,r′)=k02∭RQ(r~)G_ee(r,r~)·G_emS(r~,r′)dr~,(30)G_mmU(r,r′)=-∭RQ(r~)G_me(r,r~)·G_emS(r~,r′)dr~.
Retrieving G_S(r~,r') from (27)–(30) is a hard task, so that an approximation is needed. The computation is facilitated by discretizing the region Rin Psmall cubes (voxels), representing an orthogonal basis expansion for Q. Each cube p is centered at position r~p, whose volume Rp is defined as
(31)Rp={x~p-Δ≤x~≤x~p+Δ,y~p-Δ≤y~≤y~p+Δ,z~p-Δ≤z~≤z~p+Δ,
where Δ is the half length of the edge of elementary cube. The dimension of the voxel is chosen so that G_u(r~p,r′) and G_uS(r~p,r′) can be assumed constant within the boundaries of the voxel itself. Theoretically, Rp 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
(32)G_eeU(r,r′)≅k02∑p=1PQ(r~p)S_ee(r,r~p)·G_eeS(r~p,r′),G_meU(r,r′)≅k02∑p=1PQ(r~p)S_me(r,r~p)·G_eeS(r~p,r′),G_emU(r,r′)≅k02∑p=1PQ(r~p)S_ee(r,r~p)·G_emS(r~p,r′),G_mmU(r,r′)≅-∑p=1PQ(r~p)S_me(r,r~p)·G_emS(r~p,r′),
and a known part
(33)G_eeK(r,r′)≅k02∑q=1PQ(r~q)S_ee(r,r~q)·G_ee(r~q,r′),G_meK(r,r′)≅k02∑q=1PQ(r~q)S_me(r,r~q)·G_ee(r~q,r′),G_emK(r,r′)≅k02∑q=1PQ(r~q)S_ee(r,r~q)·G_em(r~q,r′),G_mmK(r,r′)≅-∑q=1PQ(r~q)S_me(r,r~q)·G_em(r~q,r′),
where the index q is used for the same purpose of the index p,but is distinguished to emphasize that the summations are computed separately. In (32) and (33) the dyadic expression S_(r,r′) represents the integral of the half-space Green’s function, performed over a cubic region having volume Rp, centered along the source point r′, and computed at the observation point r. The exact evaluation of this integral is described in Appendix.The integral  equations (18)-(19) can be discretized as
(34)G_eeS(r,r′)-k02∑p=1PQ(r~p)S_ee(r,r~p)·G_eeS(r~p,r′)=G_eeK(r,r′),(35)G_emS(r,r′)-k02∑p=1PQ(r~p)S_ee(r,r~p)·G_emS(r~p,r′)=G_emK(r,r′),(36)G_meS(r,r′)-k02∑p=1PQ(r~p)S_me(r,r~p)·G_eeS(r~p,r′)=G_meK(r,r′),(37)G_mmS(r,r′)+∑p=1PQ(r~p)S_me(r,r~p)·G_emS(r~p,r′)=G_mmK(r,r′).
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 ς=e;m, in these equations, G_eςS(r,r′) would be determined when the value of G_eςS(r~p,r′) is known for any value of p. To find G_S(r~p,r′) for all p, (34) or (35) are tested exactly at the locations r~p, leading to the following set of P equations:
(38)G_eςS(r~p,r')-k02∑q=1PQ(r~q)S_ee(r~p,r~q)·G_eςS(r~q,r′)=G_eςK(r~p,r′),p=1,…,P,ς=e;m.
Due to the dyadic nature of (38), 9P equations with 9P unknowns can be defined, leading to the creation of the following matrix equation:
(39)Aeςxeς=geς,
where,
(40)Aeς=[I9×9-Qeς11-Qeς12⋯-Qeς1P-Qeς21I9×9-Qeς22⋮⋱-Qeς1PI9×9-QeςPP],Qeςpq=k02Q(r~q)[I3×3SeςXXpqI3×3SeςXYpqI3×3SeςXZpqI3×3SeςYXpqI3×3SeςYYpqI3×3SeςYZpqI3×3SeςZXpqI3×3SeςZYpqI3×3SeςZZpq],xeς=[GeςXXS(r~1,r')GeςXYS(r~1,r')⋮GeςZZS(r~q-1,r')GeςXXS(r~q,r')⋮GeςZZS(r~P,r')],geς=[GeςXXK(r~1,r')GeςXYK(r~1,r')⋮GeςZZK(r~p-1,r')GeςXXK(r~p,r')⋮GeςZZK(r~P,r')].
By solving for xeς, the knowledge of G_eςS(r~p,r′) for all p follows. Direct substitution of this information into (38)  yields the full reconstruction of the scattering Green’s function at any given point r, r′, that is,
(41)G_eeS(r,r′)=G_eeK(r,r′)+k02∑p=1PQ(r~p)S_ee(r,r~p)·G_eeS(r~p,r′),(42)G_emS(r,r′)=G_emK(r,r′)+k02∑p=1PQ(r~p)S_ee(r,r~p)·G_emS(r~p,r′).The newly acquired knowledge ofG_eςS(r~p,r′) is also useful to determine the remaining types of scattering Green’s functions
(43)G_meS(r,r′)=G_meK(r,r′)+k02∑p=1PQ(r~p)S_me(r,r~p)·G_eeS(r~p,r′),(44)G_mmS(r,r′)=G_mmK(r,r′)-∑p=1PQ(r~p)S_me(r,r~p)·G_emS(r~p,r′).
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.


## 3.1. Method of Moments

The scattering integral equations in vector form for the distributed dielectric anomalyQ(r) can be easily derived from Maxwell’s equation and are stated below
(14)ES(r)=+k02∭DQ(r~)G_ee(r,r~)·ETot(r~)dr~,(15)HS(r)=-iωε0∭DQ(r~)G_me(r,r~)·ETot(r~)dr~.
In this context, ES, HS are the scattered electric and magnetic fields, and ETot 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)
(16)ES(r)=iωμ0G_eeS(r,r′)·a^,ETot(r~)=iωμ0G_eeTot(r~,r′)·a^,
and for (15)
(17)ETot(r~)=G_emTot(r~,r′)·a^,HS(r)=G_meS(r,r′)·a^.
To generalize the result, the ·a^ operation is disregarded. Substituting (16) in (14), the dyadic scattering integral equation becomes
(18)G_eeS(r,r′)=k02∭RQ(r~)G_ee(r,r~)·G_eeTot(r~,r′)dr~.
Similarly,
(19)G_emS(r,r′)=k02∭RQ(r~)G_ee(r,r~)·G_emTot(r~,r′)dr~,G_meS(r,r′)=k02∭RQ(r~)G_me(r,r~)·G_eeTot(r~,r′)dr~,G_mmS(r,r′)=-∭RQ(r~)G_me(r,r~)·G_emTot(r~,r′)dr~.
In these formulas, R is the support of Q. Since
(20)ETot(r~)=EInc(r~)+ES(r~),
and the incident electric field is
(21)EInc(r~)={iωμ0G_ee(r~,r′)·a^electric  currents,G_em(r~,r′)·a^magnetic  currents,
we obtain that
(22)G_ee;emTot(r~,r′)=G_ee;emS(r~,r′)+G_ee;em(r~,r′).
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
(23)G_eeK(r,r′)=k02∭RQ(r~)G_ee(r,r~)·G_ee(r~,r′)dr~,(24)G_meK(r,r′)=k02∭RQ(r~)G_me(r,r~)·G_ee(r~,r′)dr~,(25)G_emK(r,r′)=k02∭RQ(r~)G_ee(r,r~)·G_em(r~,r′)dr~,(26)G_mmK(r,r′)=-∭RQ(r~)G_me(r,r~)·G_em(r~,r′)dr~,
and an unknown part
(27)G_eeU(r,r′)=k02∭RQ(r~)G_ee(r,r~)·G_eeS(r~,r′)dr~,(28)G_meU(r,r′)=k02∭RQ(r~)G_me(r,r~)·G_eeS(r~,r′)dr~,(29)G_emU(r,r′)=k02∭RQ(r~)G_ee(r,r~)·G_emS(r~,r′)dr~,(30)G_mmU(r,r′)=-∭RQ(r~)G_me(r,r~)·G_emS(r~,r′)dr~.
Retrieving G_S(r~,r') from (27)–(30) is a hard task, so that an approximation is needed. The computation is facilitated by discretizing the region Rin Psmall cubes (voxels), representing an orthogonal basis expansion for Q. Each cube p is centered at position r~p, whose volume Rp is defined as
(31)Rp={x~p-Δ≤x~≤x~p+Δ,y~p-Δ≤y~≤y~p+Δ,z~p-Δ≤z~≤z~p+Δ,
where Δ is the half length of the edge of elementary cube. The dimension of the voxel is chosen so that G_u(r~p,r′) and G_uS(r~p,r′) can be assumed constant within the boundaries of the voxel itself. Theoretically, Rp 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
(32)G_eeU(r,r′)≅k02∑p=1PQ(r~p)S_ee(r,r~p)·G_eeS(r~p,r′),G_meU(r,r′)≅k02∑p=1PQ(r~p)S_me(r,r~p)·G_eeS(r~p,r′),G_emU(r,r′)≅k02∑p=1PQ(r~p)S_ee(r,r~p)·G_emS(r~p,r′),G_mmU(r,r′)≅-∑p=1PQ(r~p)S_me(r,r~p)·G_emS(r~p,r′),
and a known part
(33)G_eeK(r,r′)≅k02∑q=1PQ(r~q)S_ee(r,r~q)·G_ee(r~q,r′),G_meK(r,r′)≅k02∑q=1PQ(r~q)S_me(r,r~q)·G_ee(r~q,r′),G_emK(r,r′)≅k02∑q=1PQ(r~q)S_ee(r,r~q)·G_em(r~q,r′),G_mmK(r,r′)≅-∑q=1PQ(r~q)S_me(r,r~q)·G_em(r~q,r′),
where the index q is used for the same purpose of the index p,but is distinguished to emphasize that the summations are computed separately. In (32) and (33) the dyadic expression S_(r,r′) represents the integral of the half-space Green’s function, performed over a cubic region having volume Rp, centered along the source point r′, and computed at the observation point r. The exact evaluation of this integral is described in Appendix.The integral  equations (18)-(19) can be discretized as
(34)G_eeS(r,r′)-k02∑p=1PQ(r~p)S_ee(r,r~p)·G_eeS(r~p,r′)=G_eeK(r,r′),(35)G_emS(r,r′)-k02∑p=1PQ(r~p)S_ee(r,r~p)·G_emS(r~p,r′)=G_emK(r,r′),(36)G_meS(r,r′)-k02∑p=1PQ(r~p)S_me(r,r~p)·G_eeS(r~p,r′)=G_meK(r,r′),(37)G_mmS(r,r′)+∑p=1PQ(r~p)S_me(r,r~p)·G_emS(r~p,r′)=G_mmK(r,r′).
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 ς=e;m, in these equations, G_eςS(r,r′) would be determined when the value of G_eςS(r~p,r′) is known for any value of p. To find G_S(r~p,r′) for all p, (34) or (35) are tested exactly at the locations r~p, leading to the following set of P equations:
(38)G_eςS(r~p,r')-k02∑q=1PQ(r~q)S_ee(r~p,r~q)·G_eςS(r~q,r′)=G_eςK(r~p,r′),p=1,…,P,ς=e;m.
Due to the dyadic nature of (38), 9P equations with 9P unknowns can be defined, leading to the creation of the following matrix equation:
(39)Aeςxeς=geς,
where,
(40)Aeς=[I9×9-Qeς11-Qeς12⋯-Qeς1P-Qeς21I9×9-Qeς22⋮⋱-Qeς1PI9×9-QeςPP],Qeςpq=k02Q(r~q)[I3×3SeςXXpqI3×3SeςXYpqI3×3SeςXZpqI3×3SeςYXpqI3×3SeςYYpqI3×3SeςYZpqI3×3SeςZXpqI3×3SeςZYpqI3×3SeςZZpq],xeς=[GeςXXS(r~1,r')GeςXYS(r~1,r')⋮GeςZZS(r~q-1,r')GeςXXS(r~q,r')⋮GeςZZS(r~P,r')],geς=[GeςXXK(r~1,r')GeςXYK(r~1,r')⋮GeςZZK(r~p-1,r')GeςXXK(r~p,r')⋮GeςZZK(r~P,r')].
By solving for xeς, the knowledge of G_eςS(r~p,r′) for all p follows. Direct substitution of this information into (38)  yields the full reconstruction of the scattering Green’s function at any given point r, r′, that is,
(41)G_eeS(r,r′)=G_eeK(r,r′)+k02∑p=1PQ(r~p)S_ee(r,r~p)·G_eeS(r~p,r′),(42)G_emS(r,r′)=G_emK(r,r′)+k02∑p=1PQ(r~p)S_ee(r,r~p)·G_emS(r~p,r′).The newly acquired knowledge ofG_eςS(r~p,r′) is also useful to determine the remaining types of scattering Green’s functions
(43)G_meS(r,r′)=G_meK(r,r′)+k02∑p=1PQ(r~p)S_me(r,r~p)·G_eeS(r~p,r′),(44)G_mmS(r,r′)=G_mmK(r,r′)-∑p=1PQ(r~p)S_me(r,r~p)·G_emS(r~p,r′).
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 32×32×8 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 εD=10 and the background dielectric conductivity is σD=5×10-4 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 d=0.25m, and the box size containing the whole scene is 150m×150m×120m.Figure 2
Tunnel shape at depth of −10 m.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 fillL, 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 Δ=2m, corresponding to P=128, 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 A. For Δ=2m, 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 δ=1m. 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 Figure3. 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.Figure 3
Reconstructed image: Tx Dipole and Rx Dipole.Figure 4
Reconstructed image: Tx Dipole and Rx Loop.Figure 5
Reconstructed image: Tx Loop and Rx Dipole.Figure 6
Reconstructed image: Tx Loop and Rx Loop.


## 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 theE and H 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 matrixA and vector G_K grows fast when the number of voxels in D or Q 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.



---

*Source: 735414-2012-10-02.xml*