Signal Propagation in Soil Medium: A Two Dimensional Finite Element Procedure

Abstract

A two-Dimensional Finite Element Method of electromagnetic (EM) wave propagation through the soil is presented in this chapter. The chapter employs a boundary value problem (BVP) to solve the Helmholtz time-harmonic electromagnetic model. An infinitely large dielectric object of an arbitrary cross-section is considered for scattering from a dielectric medium and illuminated by an incident wave. Since the domain extends to infinity, an artificial boundary, a perfectly matched layer (PML) is used to truncate the computational domain. The incident field, the scattered field, and the total field in terms of the z-component are expressed for the transverse magnetic (TM) and transverse electric (TE) modes. The radar cross-section (RCS), as a function of several other parameters, such as operating frequency, polarization, illumination angle, observation angle, geometry, and material properties of the medium, is computed to describe how a scatterer reflects an electromagnetic wave in a given direction. Simulation results obtained from MATLAB for the scattered field, the total field, and the radar cross-section are presented for three soil types – sand, loam, and clay.

1. Introduction

Signal propagation underground and in soil medium constitute the backbone of the Internet of underground things (IoUT), which power many applications such as precision agriculture, border monitoring for intrusion detection, pipeline monitoring, etc. [1, 2, 3]. In the recent past, the study of signal propagation in soil medium for underground wireless communication has focused mainly on empirical techniques [4, 5, 6, 7, 8]. The most commonly used modeling techniques for implementing IoUT include electromagnetic waves, magnetic induction, and acoustic waves [9, 10]. Electromagnetic field analysis is performed in this chapter using numerical modeling with the finite element method to examine the signal strength in the soil medium.

The process of numerical modeling of how electromagnetic fields propagate and interact with physical objects and the environment is usually referred to as computational electromagnetics (CEM), numerical electromagnetics. The primary motivation of this process is to develop efficient approximations to Maxwell’s equations through numerical schemes for cases where closed-form analytical solutions of Maxwell’s equation cannot be obtained due to the complexity of geometries, material parameters, and boundary conditions. Therefore, several real-life problems that are not analytically computable, such as electromagnetic scattering, antenna radiation, electromagnetic wave propagation, electromagnetic compatibility, etc., can effectively be solved by numerical techniques. The mathematical model of the electromagnetic problem is usually obtained in terms of partial differential equations, integral equations, or integro-differential equations derived from Maxwell’s equations and a set of a priori constraints of the problem such as boundary and initial conditions material parameters and geometry. The problem is ideally defined on an infinite-dimensional function space. Numerical methods apply a discretization to the continuum to reduce infinite degrees of freedom to a finite degree of freedom. In other words, the solution of an infinite dimension-dimensional problem is projected into a finite-dimensional space. Hence, the solution to the problem becomes amenable on a digital computer. The main philosophy in most of the numerical methods is to apply the divide-and-conquer strategy. The idea is to divide an intractable continuous problem into smaller pieces (divide), express the solution over each small piece (conquer), and then combine the piecewise solutions to obtain a global solution. In this chapter, FEM will be applied to the two-dimensional boundary value problem in EM wave propagation through the soil to evaluate the signal strength of the wave propagation in soil. This evaluation will be based on the incidence angle of the transmitted wave. The radar cross-section of the scatterer will be used to evaluate the direction of the wave.

2. Time-harmonic EM model

When deriving the wave equation for electromagnetic wave propagation through a vacuum or a dielectric medium such as soil or pure water, the free charge density ${\rho }_f=\mathit{0}$ and the free current density $J_f=\mathit{0}$. However, in the case of conductors like seawater or metals, we cannot control the flow of charges and, in general, $J_f$ is certainly not equal to zero. With this, Maxwell’s equations for linear media assume the form [11].

$$\left\{\begin{array}{c}\left(i\right)\nabla ·E=\frac{\mathit{1}}{\epsilon }{\rho }_f\left(\mathit{iii}\right)\nabla \mathrm{x}E=‐\frac{\partial B}{\partial t}\\ \left(\mathit{ii}\right)\nabla ·B=\mathit{0}\left(\mathit{iv}\right)\nabla \mathrm{x}B=\mathit{\mu \sigma E}+\mathit{\mu \epsilon }\frac{\partial E}{\partial t}\end{array}\right.$$E1

Apply the curl to (iii) and (iv), and we obtain modified wave equations for the electric field E and the magnetic field B:

$${\nabla }^{\mathit{2}}u=\mathit{\mu \epsilon }\frac{{\partial }^{\mathit{2}}u}{\partial t^{\mathit{2}}}+\mathit{\mu \sigma }\frac{\partial u}{\partial t}$$E2

where u represents the scalar component of the electric or magnetic field, $\epsilon$ is the permittivity of the medium $\mu ={\mu }_r{\mu }_{\mathit{0}}$, is the magnetic permeability (${\mu }_r$ is the relative permeability of the soil ${\mu }_r=\mathit{1}$, for non-magnetic soil, ${\mu }_{\mathit{0}}=\mathit{4}\pi \mathrm{x}{\mathit{10}}^{‐\mathit{7}}N/A^{\mathit{2}}$ is the magnetic constant in vacuum), and $\sigma$ is the electric conductivity of the soil medium. Eq. (2) admits the plane wave solution

$$u\left(x,t\right)={\mathit{ue}}^{i\left(\mathit{kx}‐\mathit{\omega t}\right)}$$E3

where $\omega =\mathit{2}\mathit{\pi f}$ is the angular frequency, t is the time, and k is the complex wavenumber or propagation constant, which can be derived from Eq. (3). The boundary value problem (BVP) used to solve the time-harmonic electromagnetic problem in 2-D, can be expressed in its generic form as

(4a)$$-\frac{\partial }{\partial x}\left(p_x\frac{\partial u}{\partial t}\right)-\frac{\partial }{\partial y}\left(p_y\frac{\partial u}{\partial y}\right)+\mathit{qu}=f\mathit{\text{for}}\left(x,y\right)\in \mathrm{\Omega }$$E4

$$u=u_{\mathit{0}}\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{D}}:\text{Dirichlet}\text{Boundary}\text{Condition}$$E5

$$\left(p_x\frac{\partial u}{\partial t}{\widehat{a}}_x+p_y\frac{\partial u}{\partial y}{\widehat{a}}_y\right).\widehat{n}=\beta \mathrm{on}{\mathrm{\Gamma }}_{\mathrm{N}}:\text{Neumann}\text{Boundary}\text{Condition}$$E6

$$\left(p_x\frac{\partial u}{\partial t}{\widehat{a}}_x+p_y\frac{\partial u}{\partial y}{\widehat{a}}_y\right).\widehat{n}+\mathit{\alpha u}=\beta \mathrm{on}{\mathrm{\Gamma }}_{\mathrm{M}}:\text{Mixed}\text{Boundary}\text{Condition}$$E7

where $u\left(x,y\right)$ is the unknown function to be determined, and $p_x\left(x,y\right)$$p_y\left(x,y\right)$$q\left(x,y\right)$$f\left(x,y\right)$ are given functions. $u_{\mathsf{0}}$$\alpha$, and $\beta$ are given functions in boundary conditions (BCs); ${\Gamma }_D$${\Gamma }_N$ and ${\Gamma }_M$ refers to boundaries where Dirichlet, Neumann, and mixed BCs are imposed, respectively; $\widehat{n}={\widehat{a}}_x{n}_x+{\widehat{a}}_y{n}_y$ is the unit vector normal to the boundary in the outward direction.

The weak form is used to provide the finite element solution. The weak form is written as [12, 13, 14].

$$\underset{\Omega }{\iint }\left[p_x\frac{\partial u}{\partial x}\frac{\partial \psi }{\partial x}+p_y\frac{\partial u}{\partial y}\frac{\partial \psi }{\partial y}\right]\mathit{ds}-{\oint }_{\Gamma }\psi \left[p_x\frac{\partial u}{\partial y}{n}_x+p_y\frac{\partial u}{\partial y}{n}_y\right]\mathit{dl}+\underset{\Omega }{\iint }\mathit{q\psi u}\mathit{ds}-\underset{\Omega }{\iint }\mathit{\psi f}\mathit{ds}=\mathit{0}$$E8

where $\psi \left(x,y\right)$ is the weight function.

The weak form is applied in each element domain, and element matrices are formed by expressing the unknown function as a weighted sum of nodal shape functions. The sum of line integrals of two neighboring elements cancels out while combining the element matrices. Therefore, the line integral can be omitted for interior elements and should be considered only for elements adjacent to the boundary. Due to its special form, the line integral makes easier the imposition of mixed types of BCs. The mesh generation will be discussed first, and the shape functions will be given, and then the finite element solution of the time-harmonic problem of electromagnetic scattering in a dielectric medium (soil) will be presented.

3. The Helmholtz time harmonic model

In 2-D, the geometry and boundary conditions do not vary along an axis (say the z-axis). Hence fields can be represented as a superposition of fields of two orthogonal polarizations using the linearity property. Any field can be decomposed into transverse magnetic (TM) and transverse electric (TE) parts for the z-variable. In the TM case (horizontal polarization), only the z-component of the electric field $\left(\mathbf{E}={\widehat{a}}_z{E}_z\left(x,y\right)\right)$ and the x-y-components magnetic field $\left(\mathbf{H}={\widehat{a}}_x{H}_x\left(x,y\right)+{\widehat{a}}_y{H}_y(x,y)\right)$ exist. However, in the TE case (vertical polarization), which is a dual of TM, only the z-component of the magnetic field $\left(\mathbf{H}={\widehat{a}}_z{H}_z\left(x,y\right)\right)$ and the x-y-components electric field $\left(\mathbf{E}={\widehat{a}}_x{E}_x\left(x,y\right)+{\widehat{a}}_y{E}_y(x,y)\right)$ exist.

The material tensors in this case are defined as

$${\mathit{\epsilon }}_{\mathit{rc}}=\left[\begin{array}{ccc}{\epsilon }_{\mathit{rc}}^{\mathit{xx}}& {\epsilon }_{\mathit{rc}}^{\mathit{xy}}& 0\\ {\epsilon }_{\mathit{rc}}^{\mathit{yx}}& {\epsilon }_{\mathit{rc}}^{\mathit{yy}}& 0\\ 0& 0& {\epsilon }_{\mathit{rc}}^{\mathit{zz}}\end{array}\right]\mathit{\text{with}}{\left({\mathit{\epsilon }}_{\mathit{rc}}\right)}_{\mathit{sub}}=\left(\begin{array}{cc}{\epsilon }_{\mathit{rc}}^{\mathit{xx}}& {\epsilon }_{\mathit{rc}}^{\mathit{xy}}\\ {\epsilon }_{\mathit{rc}}^{\mathit{yx}}& {\epsilon }_{\mathit{rc}}^{\mathit{yy}}\end{array}\right),$$E9

$${\mathit{\mu }}_r=\left[\begin{array}{ccc}{\mu }_r^{\mathit{xx}}& {\mu }_r^{\mathit{xy}}& 0\\ {\mu }_r^{\mathit{yx}}& {\mu }_c^{\mathit{yy}}& 0\\ 0& 0& {\mu }_r^{\mathit{zz}}\end{array}\right]\mathit{\text{with}}{\left({\mathit{\mu }}_r\right)}_{\mathit{sub}}=\left(\begin{array}{cc}{\mu }_r^{\mathit{xx}}& {\mu }_r^{\mathit{xy}}\\ {\mu }_r^{\mathit{yx}}& {\mu }_r^{\mathit{yy}}\end{array}\right),$$E10

The generalized homogeneous Helmholtz equation in TM and TE case can be written, respectively, in the following forms [15].

$$\nabla \cdot \left({\mathbf{\Lambda }}_{\mu }\cdot \nabla E_z\right)+k_{\mathit{0}}^{\mathit{2}}{\epsilon }_{\mathit{rc}}^{\mathit{zz}}{E}_z=\mathit{0},$$E11

$$\nabla \cdot \left({\mathbf{\Lambda }}_{\epsilon }\cdot \nabla H_z\right)+k_{\mathit{0}}^{\mathit{2}}{\mu }_r^{\mathit{zz}}{H}_z=\mathit{0},$$E12

Where

$${\mathbf{\Lambda }}_{\mu }=\frac{1}{\left|{\left({\mathit{\mu }}_r\right)}_{\mathit{sub}}^T\right|}{\left({\mathit{\mu }}_r\right)}_{\mathit{sub}}^T,$$E13

$${\mathbf{\Lambda }}_{\epsilon }=\frac{1}{\left|{\left({\mathit{\epsilon }}_{\mathit{rc}}\right)}_{\mathit{sub}}^T\right|}{\left({\mathit{\epsilon }}_{\mathit{rc}}\right)}_{\mathit{sub}}^T,$$E14

Where the superscript T indicates the transpose and $\left|.\right|$ refers to the determinant of the corresponding submatrix; $k_{\mathit{0}}$ is the free-space wavenumber; ${\epsilon }_{\mathit{rc}}={\epsilon }_r-\mathit{j\sigma }/\omega {\epsilon }_{\mathit{0}}$ is the complex relative permittivity ${\epsilon }_r$, and ${\mu }_r$ are the relative permittivity and permeability, respectively; and $\sigma$ is the conductivity. For isotropic mediums, Eq. (8) becomes

$$\nabla \cdot \left({\mu }_r^{‐1}\nabla E_z\right)+k_{\mathit{0}}^{\mathit{2}}{\epsilon }_{\mathit{rc}}{E}_z=\mathit{0},$$E15

Or

$$\frac{\partial }{\partial x}\left(\frac{1}{{\mu }_r}\frac{\partial E_z}{\partial x}\right)+\frac{\partial }{\partial x}\left(\frac{1}{{\mu }_r}\frac{\partial E_z}{\partial y}\right)+k_{\mathit{0}}^{\mathit{2}}{\epsilon }_{\mathit{rc}}{E}_z=\mathit{0},$$E16

The same way Eq. (11) becomes

$$\nabla \cdot \left({\epsilon }_{\mathit{rc}}^{‐1}\nabla H_z\right)+k_{\mathit{0}}^{\mathit{2}}{\mu }_r{H}_z=\mathit{0},$$E17

Or

$$\frac{\partial }{\partial x}\left(\frac{1}{{\epsilon }_{\mathit{rc}}}\frac{\partial H_z}{\partial x}\right)+\frac{\partial }{\partial x}\left(\frac{1}{{\epsilon }_{\mathit{rc}}}\frac{\partial H_z}{\partial y}\right)+k_{\mathit{0}}^{\mathit{2}}{\mu }_r{H}_z=\mathit{0},$$E18

4. Scattering from a dielectric medium (soil)

In [16], a new wave number model is proposed with the combination of the Peplinski principle and multiple scattering from particles in the soil medium. The new wave number is used in the computation of the path loss. In another recent work, sensitivity analysis of the Ku-band scattering coefficient to soil moisture was performed under single-polarized, dual-polarized, and dual-angular combinations [17]. Similarly, a model of parabolic equations for reflection and refraction in an environment with an obstacle where the area is decomposed into two different domains. The discrete mixed Fourier transformation is used to compute the field strength in the upper subdomain, and the finite difference method is used to calculate the field strength in the lower subdomain [18].

In our case, an infinitely large dielectric object of an arbitrary cross-section is considered and illuminated by an incident wave that is not a function of z. An illustration of the scattering problem for TM mode is shown in Figure 1 where a general FEM is depicted. The problem is defined in TE mode by replacing the electric field with a magnetic field.

Figure 1.

Electromagnetic scattering in soil: FEM modeling with PML.

4.1 The perfectly matched layer (PML)

Since the domain extends to infinity, an artificial boundary or layer is used to truncate the computational domain. That is an absorbing boundary condition (ABC) or perfectly matched layer (PML) [15].

Locally-Conformal PML (LC-PML) is a powerful PML method whose implementation is straightforward. It is implemented by replacing the real coordinates $\left(\mathbf{\rho }=\left(x,y\right)\right)$ inside the PML region with their corresponding complex coordinates $\left(\tilde{\mathbf{\rho }}=\left(\tilde{x},\tilde{y}\right)\right)$ [15]. The computational domain is the union of the PML region, the free-space region, and the scatterer region.

Mesh generation is the process of representing the domain of interest as a collection of elements. The two most commonly used elements in 2-D problems are triangular and quadrilateral. Figure 2 shows the mesh that is formed by triangular elements for a rectangular domain of the boundary value problem (BVP). The figure was generated in MATLAB.

Figure 2.

Mesh formed by triangular elements for the rectangular domain of the BVP.

For a computational domain with such curved boundaries, discretization errors will occur due to the inability to capture the exact geometry of the domain. Triangular elements are preferred due to their simplicity and the possibility of developing algorithms for automatic triangulation for a computational domain such as Delaunay triangulation [19, 20]. Although fewer elements are needed when quadrilateral elements are used, Triangular elements are well-suited for complex geometries and cause fewer numerical dispersion errors. Furthermore, the calculation of element matrices is easier in triangular elements. Discretization error might inevitably occur unless sufficient mesh density is used. This happens regardless of whether triangular or quadrilateral elements are used for meshing. One way to overcome this problem is to refine the mesh. Mesh refinement might be needed especially if the geometry has a curved boundary or some corners, sharp edges, small features, or discontinuities. This might be important especially if linear interpolation functions are used. Another approach that increases accuracy is to use high-order elements at the expense of increased computational load. To achieve this, we use extra nodes within an element and use high-order interpolation functions. Figure 3 shows a simple mesh with six linear triangular elements having eight nodes. Local node numbers must follow the anticlockwise orientation in all elements to guarantee that the area of each element is obtained as a positive quantity. During the mesh generation, certain data arrays must be created [19].

Figure 3.

Mesh of a 2-D domain using linear triangular elements.

An element connectivity matrix of size Mx3, where M is the number of elements vectors of node coordinates (x and y), each of size Nx1, where N is the number of nodes representation of the given functions ($p_x$, $p_y$, $q$, and $f$) in each element (each of which is of size Mx1 Arrays containing special nodes and/or element. The connectivity matrix belonging to Figure 3 is given in Table 1.

Element (e)	Node 1	Node 2	Node 3
1	1	2	8
2	2	3	8
3	3	7	8
4	3	4	7
5	4	6	7
6	4	5	6

Table 1.

Structure of element connectivity for triangular mesh.

The Delaunay function is used to create the connectivity matrix and automatically enumerates the nodes. The triangular mesh is generated for a rectangular domain. The element size (els) = $\mathit{\Delta h}={\lambda }_{\mathit{0}}/\mathit{5}$, $\mathit{\Delta h}={\lambda }_{\mathit{0}}/\mathit{10}$, $\mathit{\Delta h}={\lambda }_{\mathit{0}}/\mathit{20}$ and the like. The higher the denominator, the finer the element. Figure 4 shows the mesh generated in MATLAB for els = $\mathit{\Delta h}={\lambda }_{\mathit{0}}/\mathit{15}$.

Figure 4.

Triangular mesh generated with $\mathit{\Delta h}={\lambda }_{\mathit{0}}/\mathit{15}$.

4.2 The scattered field formulation

Fields in the presence of the scatterer can be decomposed into two parts:

• The first one is the incident field which is produced without the scatterer and

• The second one is the scattered field produced by an equal amount of current induced on the scatterer or the surface enclosing the scatterer.

We first assume the TM polarization case, where the incident field $\left({\mathbf{E}}^{\mathrm{inc}}={\widehat{a}}_z{E}_z^{\mathrm{inc}}\right)$, and the scattered field $\left({\mathbf{E}}^{\text{scat}}={\widehat{a}}_z{E}_z^{\text{scat}}\right)$, are expressed in terms of the z-component, the z-component of the total field can be expressed as $E_z=E_z^{\text{scat}}+E_z^{\mathrm{inc}}$. For an isotropic case, the total field satisfies the differential Eq. [9].

$$-\frac{\partial }{\partial x}\left(\frac{1}{{\mu }_r}\frac{\partial E_z^{\text{scat}}}{\partial x}\right)-\frac{\partial }{\partial x}\left(\frac{1}{{\mu }_r}\frac{\partial E_z^{\text{scat}}}{\partial y}\right)-k_{\mathit{0}}^{\mathit{2}}{\epsilon }_{\mathit{rc}}{E}_z^{\text{scat}}=f\left(x,y\right),$$E19

Where the source term is given by

$$f\left(x,y\right)=\frac{\partial }{\partial x}\left(\frac{1}{{\mu }_r}\frac{\partial E_z^{\mathrm{inc}}}{\partial x}\right)+\frac{\partial }{\partial x}\left(\frac{1}{{\mu }_r}\frac{\partial E_z^{\mathrm{inc}}}{\partial y}\right)+k_{\mathit{0}}^{\mathit{2}}{\epsilon }_{\mathit{rc}}{E}_z^{\mathrm{inc}},$$E20

This differential equation is a special form of Eq. (4a), where $u=E_z^{\text{scat}}\left(x,y\right)$, $p_x=p_y=p=\frac{1}{\mu }$ and $q=-k_{\mathit{0}}^{\mathit{2}}{\epsilon }_{\mathit{rc}}$.

The incident field can be arbitrary and usually chosen as a uniform plane wave since the incident field sources are sufficiently far away from the object. The incident field with the unit magnitude is given by

$${\mathbf{E}}^{\mathrm{inc}}={\widehat{a}}_z{E}_z^{\mathrm{inc}}={\widehat{a}}_{z}exp\left[\mathit{jk}\left(xcos{\phi }^{\mathrm{inc}}+ysin{\phi }^{\mathrm{inc}}\right)\right],$$E21

where ${\phi }^{\mathrm{inc}}$ is the angle of incidence for the x-axis in cartesian coordinates.

Similar calculations can be performed in TE polarization mode by replacing the electric field with the magnetic field and permittivity with permeability. Hence, the differential equation in terms of the scattered magnetic field is given by

$$-\frac{\partial }{\partial x}\left(\frac{1}{{\epsilon }_{\mathit{rc}}}\frac{\partial H_z^{\text{scat}}}{\partial x}\right)-\frac{\partial }{\partial x}\left(\frac{1}{{\epsilon }_{\mathit{rc}}}\frac{\partial H_z^{\text{scat}}}{\partial y}\right)-k_{\mathit{0}}^{\mathit{2}}{\mu }_r{H}_z^{\text{scat}}=f\left(x,y\right),$$E22

Where the source term is given by

$$f\left(x,y\right)=\frac{\partial }{\partial x}\left(\frac{1}{{\epsilon }_{\mathit{rc}}}\frac{\partial H_z^{\mathrm{inc}}}{\partial x}\right)+\frac{\partial }{\partial x}\left(\frac{1}{{\epsilon }_{\mathit{rc}}}\frac{\partial H_z^{\mathrm{inc}}}{\partial y}\right)+k_{\mathit{0}}^{\mathit{2}}{\mu }_r{H}_z^{\mathrm{inc}},$$E23

This differential equation is a special form of Eq. (4a), where $u=H_z^{\text{scat}}\left(x,y\right)$, $p_x=p_y=p=\frac{1}{{\epsilon }_{\mathit{rc}}}$ and $q=-k_{\mathit{0}}^{\mathit{2}}\mu$.

For dielectric objects, the right-hand side of the matrix equation can be obtained by using the source term f(x, y) in each element. An alternative simpler approach is that the source term in Eq. (20) is just the differential operator applied to the known incident field within the object. The left-hand side of Eq. (19) is the same differential operator being applied to the unknown scattered field, which yields the left-hand side of the matrix equation $\left(\mathbf{Au}\right)$, where $\mathbf{u}$ refers to the vector of nodal values of the scattered field. Therefore, $\mathbf{Au}\mathbf{=}-{\mathbf{Au}}^{\mathrm{inc}}$ must be satisfied within the scatterer region. After forming the global matrix, as usual, the right-hand side vector can be modified by just multiplying the incident field vector by the global matrix, i.e. $\left(\mathbf{b}\mathbf{=}-{\mathbf{Au}}^{\mathrm{inc}}\right)$, only for entries $\mathbf{b}$ corresponding to the nodes lying inside the object.

4.3 Radar cross section

The radar cross-section (RCS) of the scatterer is perhaps the most critical parameter that must be evaluated in the post-processing phase of FEM for the electromagnetic scattering problem [9]. RCS is the reflection of the scattering electromagnetic ware at an incident angle in a particular direction. In other words, it is the area capturing that amount of scattered power produced at the receiver in an isotropic medium. This is a density that is equal to that scattered by the actual target. It is a function of several parameters, such as operation frequency, polarization, illumination angle, observation angle, geometry, and material properties of the object. In 2-D, it is mathematically defined as

$${\mathrm{\sigma }}_{\mathit{2}D}=\underset{r\to \infty }{lim}\mathit{2}\mathit{\pi \rho }\frac{{\left|{\mathrm{u}}_{\mathrm{far}}^{\text{scat}}\right|}^{\mathit{2}}}{{\left|{\mathrm{u}}^{\mathrm{inc}}\right|}^{\mathit{2}}}$$E24

Where ${\mathrm{u}}_{\mathrm{far}}^{\text{scat}}=\frac{1}{\sqrt{\rho }}f\left(\phi \right)$ is the scattered electric or magnetic field at far-zone observed in a given direction, when $\rho$ satisfies the inequality $\rho \gg \mathit{2}{D}^{\mathit{2}}/\lambda$, where D is the largest dimension of the scatterer and $\lambda$ is the wavelength.

If the incident and observation directions are the same, the RCS is called monostatic or backscatter RCS; otherwise, it is referred to as bistatic RCS. In 2-D, RCS is usually normalized for $\lambda$ (wavelength) or m (meter). The unit for ${\sigma }_{\mathit{2}D}/\lambda$ is dBw, and for ${\sigma }_{\mathit{2}D}/\mathrm{m}$ is dBm.

5. Scattering inside the dielectric medium

For each such element, the far-zone field is computed and then superposed over all elements, as follows:

$$E_{\mathrm{far}}^{\text{scat}}=\eta \sqrt{\frac{k}{\mathit{8}\mathit{\pi \rho }}}{e}^{\raisebox{1ex}{$-j\mathit{3}\pi $}\!\left/ \!\raisebox{-1ex}{$\mathit{4}$}\right.}{e}^{-\mathit{jk\rho }}\sum _{i=\mathit{1}}^K{J}_z^{\left(i\right)}{e}^{\mathit{jk}\left(x_m^{\left(i\right)}\mathit{\text{cosφ}}+y_m^{\left(i\right)}\mathit{\text{sinφ}}\right)}{\mathit{\Delta l}}^{\left(i\right)},$$E25

where K is the number of elements adjacent to the boundary, $\left(x_m^{\left(i\right)},y_m^{\left(i\right)}\right)$ the i-th segment midpoint, and ${\mathit{\Delta l}}^{\left(i\right)}$ the i-th segment length.

Since FEM is formulated in terms of the scattered electric field, computation of the derivatives of the scattered field requires additional effort. This can simply be achieved by using the weighted sum of derivatives of shape functions in terms of nodal scattered fields.

5.1 Scattering inside the dielectric medium, TM case

For the dielectric object, the TM case, the interior part of the object should be included. Since the magnetic current density becomes nonzero, the integral containing the magnetic current density should be evaluated. The magnetic current density has x- and y-components, which is defined as $\mathrm{M}={\widehat{a}}_x{M}_x+{\widehat{a}}_y{M}_y$. Hence, it shows that $\widehat{\rho }\times \mathrm{M}=\left(cos\phi M_y-sin\phi M_x\right){\widehat{a}}_z$. The components of the magnetic current density can be determined in terms of the scattered electric field as follows:

$$\mathrm{M}=\mathrm{E}\times \widehat{n}=-n_y{E}_z{\widehat{a}}_x+n_x{E}_z{\widehat{a}}_y,$$E26

Finally, the scattered electric field can be obtained as follows:

$$\begin{array}{l}{E}_{\mathrm{far}}^{\text{scat}}=\sqrt{\frac{k}{\mathit{8}\mathit{\pi \rho }}}{e}^{\raisebox{1ex}{$-j\mathit{3}\pi $}\!\left/ \!\raisebox{-1ex}{$\mathit{4}$}\right.}{e}^{-\mathit{jk\rho }}\sum _{i=\mathit{1}}^K\left(n_{x}cos\phi +n_{y}sin\phi \right)E_z^{\left(i\right)}{e}^{\mathit{jk}\left(x_m^{\left(i\right)}\mathit{\text{cosφ}}+y_m^{\left(i\right)}\mathit{\text{sinφ}}\right)}{\mathit{\Delta l}}^{\left(i\right)}\\ +\eta \sqrt{\frac{k}{\mathit{8}\mathit{\pi \rho }}}{e}^{\raisebox{1ex}{$-j\mathit{3}\pi $}\!\left/ \!\raisebox{-1ex}{$\mathit{4}$}\right.}{e}^{-\mathit{jk\rho }}\sum _{i=\mathit{1}}^K{J}_z^{\left(i\right)}{e}^{\mathit{jk}\left(x_m^{\left(i\right)}\mathit{\text{cosφ}}+y_m^{\left(i\right)}\mathit{\text{sinφ}}\right)}{\mathit{\Delta l}}^{\left(i\right)}\end{array}$$E27

Here, $E_z^{\left(i\right)}$ can be determined as the average of nodal field values connected to the boundary.

5.2 Scattering inside the dielectric medium, TE case

For dielectric object, TE case: The magnetic current density is nonzero and has only z-component, $\mathrm{M}={\widehat{a}}_z{M}_z$ because $\widehat{\rho }\times \left(\widehat{\rho }\times \mathrm{M}\right)=-M_z{\widehat{a}}_z$. The scattered magnetic field is determined as:

$$\begin{array}{l}{H}_{\mathrm{far}}^{\text{scat}}=\sqrt{\frac{k}{\mathit{8}\mathit{\pi \rho }}}{e}^{\raisebox{1ex}{$-j\mathit{3}\pi $}\!\left/ \!\raisebox{-1ex}{$\mathit{4}$}\right.}{e}^{-\mathit{jk\rho }}\sum _{i=\mathit{1}}^K\left(J_x^{\left(i\right)}sin\phi -J_y^{\left(i\right)}cos\phi \right)e^{\mathit{jk}\left(x_m^{\left(i\right)}\mathit{\text{cosφ}}+y_m^{\left(i\right)}\mathit{\text{sinφ}}\right)}{\mathit{\Delta l}}^{\left(i\right)}\\ +\eta \sqrt{\frac{k}{\mathit{8}\mathit{\pi \rho }}}{e}^{\raisebox{1ex}{$-j\mathit{3}\pi $}\!\left/ \!\raisebox{-1ex}{$\mathit{4}$}\right.}{e}^{-\mathit{jk\rho }}\sum _{i=\mathit{1}}^K{M}_z^{\left(i\right)}{e}^{\mathit{jk}\left(x_m^{\left(i\right)}\mathit{\text{cosφ}}+y_m^{\left(i\right)}\mathit{\text{sinφ}}\right)}{\mathit{\Delta l}}^{\left(i\right)}\end{array}$$E28

6. Simulation results and discussion

The operating frequency used in the simulation is 300 MHz, and the free-space wavelength is ${\mathrm{\lambda }}_0=1\mathrm{m}$. The mesh size is approximatively $\mathit{\Delta h}={\lambda }_{\mathit{0}}/\mathit{40}$. The dielectric medium is a rectangular computational domain, the relative permittivity for sandy soil ${\epsilon }_{\mathit{rc}\mathit{1}}=\mathit{4}+i\mathit{0.1}$, the relative permittivity for loamy soil ${\epsilon }_{\mathit{rc}\mathit{2}}=\mathit{12}+i\mathit{1.8}$, and the relative permittivity for clay soil ${\epsilon }_{\mathit{rc}\mathit{3}}=\mathit{25.3}+i\mathit{5.7}$. The angle of incidence is ${\mathit{0}}^{\mathsf{o}}$, that is, the incident plane wave travels in the $-y$ direction. The wavelength inside the dielectric medium (soil) decreases in terms of the wavelength in free-space and relative permittivity, i.e. $\lambda ={\lambda }_{\mathit{0}}/\sqrt{{\epsilon }_{\mathit{rc}}}$, and hence, mesh size $\mathit{\Delta h}={\lambda }_{\mathit{0}}/\mathit{40}$ resolution corresponds to $\lambda \sqrt{{\epsilon }_{\mathit{rc}}}/\mathit{40}$ resolution inside the dielectric. Therefore, smaller elements are used in dielectric media to preserve the level of accuracy. All simulation results are obtained using MATLAB.

Figure 5 shows scattering from a dielectric medium for TE mode, with sandy soil characteristics. In (a) and (b), the scattered and total fields are shown respectively. We observe a minimal electric field scattering and an intense total electric field inside the dielectric medium of sandy soil. This is because this medium is porous and allows for better signal propagation.

Figure 5.

Scattering from a dielectric medium (Sandy soil): (a) scattered field (TE) (b) Total field (TE).

Figure 6 shows scattering from a dielectric medium for TM mode, with sandy soil characteristics. In this case, we observe a minimal magnetic field scattering in (a) and an intense total magnetic field in (b) inside the dielectric medium of sandy soil.

Figure 6.

Scattering from a dielectric medium (Sandy soil): (a) scattered field (TM) (b) Total field (TM).

Figure 7 shows scattering from a dielectric medium for TE mode, with loamy soil characteristics. In (a) and (b), the scattered and total fields are shown respectively. We observe an intense electric field scattering and a minimal total electric field inside the dielectric medium of loamy soil. This is because this medium is less porous and presents some challenges in signal propagation.

Figure 7.

Scattering from a dielectric medium (loamy soil): (a) scattered field (TE) (b) Total field (TE).

Similarly, Figure 8 shows scattering from a dielectric medium for TM mode, with loamy soil characteristics. In this case, we observe a high magnetic field scattering in (a) and a low total magnetic field in (b) inside the dielectric medium of loamy soil.

Figure 8.

Scattering from a dielectric medium (loamy soil): (a) scattered field (TM) (b) Total field (TM).

Figure 9 shows scattering from a dielectric medium for TE mode, with clay soil characteristics. In (a) and (b), the scattered and total fields are shown respectively. We observe very high electric field scattering and a low total electric field inside the dielectric medium of clay soil. This is because this medium is non-porous and presents a very poor signal propagation.

Figure 9.

Scattering from a dielectric medium (clay soil): (a) scattered field (TE) (b) Total field (TE).

Similarly, Figure 10 shows scattering from a dielectric medium for TM mode, with clay soil characteristics. In this case, we observe a very high magnetic field scattering in (a) and a low total magnetic field in (b) inside the dielectric medium of clay soil.

Figure 10.

Scattering from a dielectric medium (clay soil): (a) scattered field (TM) (b) Total field (TM).

Figures 11–13 show the bistatic RCS profiles to describe how scatterers reflect the incident electromagnetic wave in a given direction. This is the area intercepting that amount of power which, when scattered in the soil medium, produces at the receiver a density which is equal to that scattered by the actual target. The RCS is a function of several parameters, such as operation frequency, polarization, illumination angle, observation angle, geometry, and properties of the soil medium. It is shown in the TE modes (a) and TM modes (b) for sandy soil, loamy soil, and clay soil, respectively.

Figure 11.

Radar cross section in Sandy soil for (a) TE mode and (b) TM mode.

Figure 12.

Radar cross section in loam soil for (a) TE mode and (b) TM mode.

Figure 13.

Radar cross section in clay soil for (a) TE mode and (b) TM mode.

7. Conclusion

In a two-Dimensional Finite Element Analysis of EM wave Propagation through the soil, a boundary value problem (BVP) used to solve the time-harmonic electromagnetic problem in 2-D, has been expressed in its generic form. In TM and TE cases, the Helmholtz model has considered an infinitely large dielectric object of an arbitrary cross-section for scattering from a dielectric medium and illuminated by an incident wave. Since the domain extends to infinity, an artificial boundary, an absorbing boundary condition (ABC), or a perfectly matched layer (PML), has been used to truncate the computational domain. The incident field, the scattered field, and the total field in terms of the z-component are expressed for the TM and TE modes. The radar cross-section (RCS), a function of several parameters, such as operation frequency, polarization, illumination angle, observation angle, geometry, and material properties of the medium, has been computed to describe how a scatterer reflects an incident electromagnetic wave in a given direction. Simulation results for the scattered field, the total field, have been presented for soil types, and the radar cross-section for different element refinements have also been presented.