Abstract

To predict the propagation of radio waves in the environment of dielectric ground and dielectric obstacles, a new two-way parabolic equation (2W-PE) method based on the domain decomposition principle and surface impedance boundary conditions (SIBC) is proposed. First, we decompose the obstacle area into different subdomains and derive the SIBC in each subdomain in detail; then, the discrete hybrid Fourier transform (DMFT) in the upper subdomain and finite difference (FD) algorithm in the lower subdomain is used to solve 2W-PE combined with SIBC, respectively. After that, we explain the algorithm steps in the process of calculating the total field, compared with the traditional 2W-PE, and then finally introduce the method of moments (MoM) combined with the enhanced discrete complex image (E-DCIM) method for accuracy verification of the new 2W-PE algorithm. The simulation results show that no matter how the obstacle medium parameters change, the results of 2W-PE method proposed in this paper and MoM are always in good agreement, which proves the accuracy of 2W-PE and its superiority in speed. Therefore, this paper provides a reliable and efficient method for solving the problem of radio wave propagation in the presence of obstacles, especially in the case of low-lossy obstacles.

1. Introduction

In many fields of contemporary society, especially the fields of communication and national defense construction, it is becoming more and more important to predict the field value quickly and accurately. In actual scenes, the terrain is usually undulating which reflects and transmits electromagnetic waves in a complex manner, and the undulations of the ground usually have a great influence on the propagation of radio waves in the troposphere. Therefore, it is particularly important to study methods that can effectively predict the spatial field value under the influence of irregular terrain. The parabolic equation (PE) method [1] is a numerical calculation approach widely used in radio wave propagation modelling and tropospheric propagation prediction in recent years. Compared with several other commonly used methods, such as finite difference method [2] and optical approximation method [3], it has the advantages of fast calculation speed and good numerical stability.

The traditional 2W-PE was first proposed by Ozgun [4]. It uses the staircase terrain model to perform equivalent treatment of irregular terrain. Combined with impedance boundary conditions, it is able to calculate the radio wave propagation problems on more complex irregular terrain [5]. However, the traditional 2W-PE does not consider the propagation paths inside the obstacle. These paths do not need to be considered when the electrical loss of obstacle is large because the electromagnetic wave inside the obstacle attenuates very quickly, but it cannot be ignored when the obstacle is low-lossy. Furthermore, the traditional 2W-PE based on the staircase terrain model is not able to deal with the cases of inclined surface. Therefore, Wang [6] proposed a hybrid method that uses both the staircase terrain model and the conformal mapping model to solve the irregular terrain problem. The 2W-PE method models low-frequency (LF) ground wave propagation, which improves the prediction accuracy of forward and backward propagation. However, it is still difficult to deal with the impact of steep obstacles.

Therefore, in response to above problems, Wei [7] divides the field value calculation of the obstacle area into two calculation areas above and inside the obstacle based on the principle of domain decomposition. Then, the upper area is downward adding a section of calculation range and the lower area is extended up to infinity, in which 2W-PE can be solved, respectively. Wei et al. [7] takes the propagation inside the obstacle into consideration and solves the problem of terrain inclination, but currently only discusses the situation of PEC ground.

When using 2W-PE to solve the radio wave propagation problem, we only care about the field above the ground, so the impedance boundary condition is used to deal with the influence of the field below the radio wave propagation in total space as the field value relationship is obtained on the boundary. Thus, we only need to calculate the field strength of the upper half space, without calculating the field distribution of the lower space, which greatly saves the amount of calculation and storage space. In [8], the generalized impedance boundary conditions are derived and summarized from the flat rectangular coordinate system and the flat Earth model, respectively. However, in this method, the angle of incidence plays a leading role in the calculation of impedance boundary conditions. To solve this problem, Fock and Leontovich [9] proposed the Leontovich impedance boundary condition, also known as the surface impedance boundary condition (SIBC), which is widely used in electromagnetic field problems in the frequency domain. However, the above boundary conditions are usually applied in 2W-PE when the electromagnetic wave is incident from the low-lossy medium above to the high-lossy medium below where electromagnetic waves are bounded, such as propagating from free space to obstacles or the ground [5]. But when the propagation process is from the inside of the obstacle to the free space, the conventional Leontovich impedance boundary condition is not satisfied.

Finally, in the verification of the accuracy of the 2W-PE algorithm, high-frequency approximation methods are usually applied, such as geometric optics method (GO) and geometric diffraction theory (UTD) methods, but this method ignores the reflected field generated by the edges of obstacles, which can only be used to verify the calculation accuracy of the 1W-PE. MoM is a more rigorous calculation method that uses the electric field integral equation (EFIE) and the magnetic field integral equation (MFIE) to solve the field distribution in space [10]. On this basis, Apaydm and Sevgi [11] proposed to use MoM to divide the ground, calculate the segmented current generated by PEC ground, and obtain the scattered field by superimposing the contribution of the segmented current through the Green function. Finally, by adding the incident field, the field value of total space can be obtained. However, this method only discusses the situation of PEC and flat terrain. Therefore, Wei et al. [7] derives MoM model in the case of undulating terrain and accurately solves the field distribution in the obstacle environment space, but it is also only applicable to PEC ground.

Aiming at the problem above, this paper proposes a 2W-PE method based on the principle of domain decomposition and SIBC, which can solve the total field distribution in complex environment with medium obstacles and overcome the shortcomings of 2W-PE that are limited by the terrain inclination and neglect field inside obstacle terrain. The paper is structured as follows: Section 2 describes the problem we need to solve and briefly explains the theoretical basis of 2W-PE. Section 3 specifically introduces the process of domain decomposition, using reflection and transmission coefficients to derive the SIBCs on each boundary. Then, discrete mixed Fourier transform (DMFT) [12] is used to solve the field value of the area above the obstacle and flat medium ground combined with the SIBC; the finite difference algorithm (FD) [13] is used to calculate the internal field value of the obstacle in combination with the upper and lower SIBCs, and the calculation process of the total field value is introduced. Section 4 introduces MoM under the condition of the medium ground as a verification method of the parabolic equation method. Section 5 shows the simulation results, and finally, Section 6 gives a comprehensive summary of the article.

2. Description of the Problem

We study the application of 2W-PE in field value prediction in the radio wave propagation scene as shown in Figure 1. In the case of horizontal polarization and vertical polarization, the field satisfies the 2-D wave equation:where, in the case of horizontal and vertical polarization, in the 2-D rectangular coordinate system, is defined, respectively, aswhere is the wave number in the uniform environment of the medium and represents the modified atmospheric refractive index. It is assumed that the refractive index does not change with the distance, and the beamwidth of electromagnetic wave propagation is limited, when is the propagation direction. Then, after factoring (1), we can get

Figure 1 

Schematic diagram of wave propagation problems.

Introducing two different phase demodulation factors, we can obtainwhere represent the forward field and the backward field, respectively. Substituting equation (4) into the equations in (3) after factoring, we obtain the forward and backward propagation parabolic equations as follows:

In the calculation, the initial field of 2W-PE is directly solved by the actual current distribution on the antenna. The 2-D current density distribution is equivalent to the initial field of vertical polarization, and the 2-D magnetic current density distribution is equivalent to the initial field of horizontal polarization [14]:which unifies the incident source environment of subsequent analysis of MoM and 2W-PE.

3. Calculation Process of 2W-PE Based on Domain Decomposition Principle

3.1. Conception of Domain Decomposition Algorithm

In the traditional 2W-PE calculation of undulating terrain, obstacles and the ground are considered as a whole, but this process may cause the terrain inclination to exceed the range of elevation angle in which the results of 2W-PE are calculated accurately. Besides, it is impossible to discuss the case when parameters of obstacle are different from those of the ground. Therefore, here, we introduce the domain decomposition method, focusing on the discussion of the field strength and the propagation paths inside the obstacle. After that, the calculation angle can be better controlled within the applicable elevation angle range of the parabolic equation, and the total field values can be calculated more accurately considering the internal field of the obstacle.

In the case of a homogeneous medium ground as shown in Figure 2, when meeting an obstacle, the area for calculating radio wave propagation can be decomposed into the upper subdomain above the obstacle and the lower subdomain inside the obstacle. The flat terrain area without obstacles can be recorded as .

Figure 2 

Schematic diagram of domain decomposition.

The upper boundary of adopts the absorbing boundary, which can be considered as infinity, and the lower boundary is the top of the obstacle, where SIBC is used for calculation. Therefore, the 2W-PE in this area can be solved step by step using DMFT combined with SIBC. The situation of is similar to that of , the upper boundary adopts the absorption boundary, and the lower boundary is the medium ground adopting SIBC.

The lower boundary of is the uniform medium ground. The upper boundary of is connected with the lower boundary of , which is exactly the top of obstacles. The upper and lower boundaries are obtained in the form of SIBC. Since is a limited area, FD is chosen to solve the 2W-PE inside the area.

3.2. Solving Impedance Boundary Conditions

When we solve the 2W-PE in each subdomain, in addition to the choice of calculation method, the most important thing is the choice of boundary conditions between subdomains. As shown in Figure 3, when there are medium 1 and 2 in the space, the lower boundary of the medium 1 area and the upper boundary of the medium 2 area are the same boundary in the physical sense, but in the actual domain solution process, the boundary conditions of the calculation in these two areas may be different. Next, we will discuss the surface impedance and corresponding SIBCs generated by electromagnetic waves incident on the boundary of different polarizations.

(a)

(b)

(a)
(b)

Figure 3 

Schematic diagram of reflection and refraction at the interface of two media of H pol. (a) Propagation from medium 1 to medium 2. (b) Propagation from medium 2 to medium 1.

We define the relationship between the total tangential electric field and the total tangential magnetic field on the boundary aswhere is the surface impedance. Now, we derive the relationship between and Fresnel reflection coefficient.

The parameters of medium 1 and 2 are shown in the Figure 3, so the wave impedance of medium 1 and 2 can be expressed aswhere . Similarly, by the definition of wave number, the wave numbers in medium 1 and 2 are, respectively,

According to Snell’s Law of refraction at the interface of the medium, we can get

Thus, from (9) and (10), we can obtain

According to the propagation of horizontally polarized waves shown in Figure 3, the relationship between the reflected and incident fields of horizontal polarization can be obtained aswhere

According to the definition of surface impedance in (7), the surface impedance of horizontal polarization can be obtained as

Thus, we can obtain the SIBC of horizontal polarization:

In the case of horizontally polarized waves, there is only one electric field component , and from the Maxwell curl equation, we can get

Substituting (16) into (15), and according to (4), the boundary condition expression after demodulation can be obtained:

3.2.1. The Lower Boundary of

Here, we consider an ideal situation, electromagnetic waves are approximately horizontally irradiated , which enters the obstacle (medium 2) from the free space (medium 1), as shown in Figure 3(a). Thus, from (11), we can obtain

Substituting (18) into (17), we get

(19) can be further written as

This formula is consistent with the conventional Leontovich impedance boundary, so the lower boundary of adopts this boundary condition. In the area without obstacle , electromagnetic waves enter the ground medium from free space; medium 1 is free space, and medium 2 is the ground. Therefore, the lower boundary of this area can be obtained in the same way.

3.2.2. The Upper Boundary of

If we consider the upper boundary condition inside the obstacle, the electromagnetic wave enters the free space (medium 1) from the obstacle (medium 2), as shown in Figure 3(b). We assume that the electromagnetic waves are approximately horizontally irradiated, , as the electromagnetic wave is incident from the optically dense medium (obstacles) to the optically thin medium (free space). Due to the large incident angle, the electromagnetic wave will be totally reflected at the boundary. Similarly, in this case, according to the reflection and transmission coefficients, the upper boundary conditions of can still be written as

(21) can be further written as

It can be seen from (22) that the electromagnetic wave enters the free space from the obstacles and takes the form of surface wave. The electromagnetic wave decays rapidly as the distance of the incident into the free space increases. This contrasts with the actual situation where the electromagnetic wave is totally reflected on the boundary. Therefore, (21) is used as the upper impedance boundary condition inside the obstacle.

3.2.3. The Lower Boundary of

We assume that the parameters of the ground medium are , and the lower boundary of the lies between the obstacle medium and the ground medium. Electromagnetic waves enter the ground of medium 1 from the obstacle of medium 1 and enter the ground of medium 2. If the dielectric constant of the ground is greater than the dielectric constant of the obstacle, the boundary condition is consistent with the conventional Leontovich boundary condition in (20). If the dielectric constant of the ground is greater than that of the obstacle, total reflection will occur, and the boundary conditions are similar to (21). Both boundaries of the above situations can be written uniformly:

The above is the case of horizontal polarization. Similarly, SIBCs on each boundary of vertical polarization can be obtained.

3.3. Solving 2W-PE in with DMFT

In this paper, when solving 2W-PE in of the domain decomposition step by step, the lower boundary is a homogeneous medium obstacle, so the iterative solution of the upper half space is equivalent to solving 2W-PE under the impedance boundary. First, we expand equation (5) according to the Feit–Fleck model and then combine the impedance boundary conditions given in (20) and the DMFT [12] method to calculate the field value:

For the flat area around obstacles, the calculation method is the same as that in .

3.4. Solving 2W-PE in with Finite Difference Method

Considering that the calculation interval is a finite area, here taking forward propagation as an example, we apply FD method to calculate the internal field value of the obstacle.

FD method needs to rewrite both the first-order partial differential and the second-order partial differential in the 2W-PE in a differential form, so that the step-by-step method is used to solve the problem to obtain the field value in the entire calculation area. The grid division is shown in Figure 4. Denotewhere is the vertical grid point; denote as the horizontal grid point and the x direction is the stepping direction. Set the midpoint of and as ; thus, the difference expression of each partial differential at the point can be obtained as

Figure 4 

Schematic diagram of finite difference mesh generation.

In order to realize 2W-PE method of the radio wave propagation calculation in the wide-angle case, the approximate expression of Q in (5) is obtained by using Pade (1, 1); then, we rewrite the forward parabolic equation in (6) as [15]

Then, substituting the differential expression (26) into (27), we can get simplified forward propagation equation:where

Equation (28) only provides the equations. According to (23), we can get the lower boundary condition of :

Substituting (30) into the difference expression of partial differential and then discretizing the forward parabolic equation (27), the equation on the lower boundary can be obtained:where

In the same way, the upper boundary condition of is obtained from (21):

Thus, the corresponding discrete equation on the upper boundary can be written aswhere

Thus, we form a 0-N forward difference parabolic equation combining (28), (31), and (34):

The recursive structure of the difference decomposition of the parabolic equation is expressed as a matrix:where is the field vector at

According to (33), we can obtain the matrixces

The above is the case of forward propagation, and similarly, the method can be applied to the case of backward propagation.

3.5. Calculation of Total Field Strength

When the traditional 2W-PE encounters undulating terrain, the obstacle and the ground are usually regarded as a whole, and the field value below the undulating terrain is set to 0, where the field value and propagation process inside the obstacle are not considered. Our method considers the propagation paths inside the obstacle and divides the propagation paths into forward and backward paths according to whether the final point is the end point or back to the starting point. Each type of paths is divided into different levels according to the number of reflections, in which, the first level means that the wave source propagates directly from the starting point to the end point, and the second level means that the new wave source generated by reflection propagates directly back to the starting point. The higher the level is, the greater the calculation accuracy is. In the calculation process of each level, the propagation paths of all wave sources under the specified reflection conditions can be completely calculated, and finally, the total field in space can be obtained through phase correction and boundary correction. The algorithm needs to complete the calculation once as soon as a new wave source is generated in the calculation process of the same level. Thus, it can record the historical paths of all wave sources during the propagation process and further accurately solve the field distribution in the space in the case of the undulating terrain. We assume the obstacle environment and incident source as shown in Figure 5, where ①, ②, ③, and ④ are reflection edges, respectively. Now, we describe the first 9 steps of propagation in detail as follows.

Figure 5 

Schematic diagram of the propagation paths inside the obstacle.

As can be seen from Figure 5, there are 3 forward paths F1∼F3 and 4 backward paths B1∼B4 in the first 9 steps. Therefore, we summarize the algorithm calculation ideas as follows:(1)Before the calculation starts, the reflection points are sorted and stored according to the obstacles, and two modules are set to record the path propagation trajectory and hierarchical matrix of the wave source.(2)Each time a reflecting surface is encountered during step calculation, a new source is set and its position and direction are recorded.(3)When adding a new source at each reflecting surface, it is necessary to determine whether to add it based on the relative size of the total field values before and after the previous source increase. That is, we need to calculate differences norm between each column before and after the previous source increase. If the sum of the differences norm is greater than a certain threshold, we choose to add the new source; otherwise, discard it. Moreover, we also give up adding a new source if the number of reflections exceeds a certain number.(4)Record the paths history of each new source for the final calculation of the phase correction and amplitude correction of the total field.(5)When calculating, we need to judge the end for each new source. If the path reaches the end or starting point, the calculation is over.

Through the above theoretical analysis, we can more accurately use 2W-PE to solve the 2D field in the obstacle environment.

4. The Analysis of Method of Moments

We use MoM to verify the calculation accuracy of the field distribution obtained by the proposed 2W-PEover the medium ground. Due to the large number of unknowns and large amount of calculation, MoM is inefficient as a prediction method of actual radio wave propagation. However, due to its analytical nature, it can be effectively compared and verified with 2W-PE in the case of small scales.

4.1. Establishment of PMCHW Equation

Here, we choose the PMCHW equation obtained by the combination of the principle of internal and external space plane equivalence and boundary conditions to solve the spatial field intensity distribution in the environment in Figure 1. Meanwhile, on the surface of the obstacle, there are the field generated by the direct wave of the original source, the field generated by the reflected wave from the direct wave reflected by the ground, the scattered field directly generated by the obstacle, and the reflected field generated by the direct scattered field through the ground. The tangential field formed on the surface satisfies the PMCHW equation. After the discretization by the basis function and the test function in MoM, the total field distribution in space can be calculated by the equivalent current, magnetic current, electric field, and magnetic field obtained by solving the equation. Assuming that the wave number in free space is , the wave number in the obstacle medium is , and the wave number in the ground medium is , the PMCHW equation is obtained as follows:where and , respectively, represent the equivalent current and equivalent magnetic current on the obstacle, and , respectively, represent the sum of the incident field and the field reflected by the ground from the incident field generated by the radiation source in (6). When the external equivalent is established, the equivalent internal field is 0, and the external field remains unchanged, so it means that the equivalent current generates an electric field during the external equivalent process of MoM. The operators expressing the electric field generated by the equivalent magnetic current should be decomposed as the sum of operators expressing the direct field and operators expressing reflection field by the ground:where and , respectively, arewhere is the equivalent current or magnetic current, , and are, respectively, the coordinates of the field point and the source point, so is the distance between the source point and the field point.

and , respectively, arewhere are the transverse wave numbers in free space and the ground, respectively, are the longitudinal wave numbers in free space and the ground, and is the reflection coefficient of free space incident on the ground.

The equivalent process of an electric field generated by the equivalent current and the operators of the equivalent magnetic current to generate the electric field are represented by and . Since the internal equivalent obstacle is a zero field outside, the ground does not need to be considered. Thus, the only difference between them and (42) is that the wave number in the operator is in free space, and in the obstacle medium is the wave number used in . According to the principle of duality, the case in which magnetic current generates magnetic field and current generates magnetic field can be obtained.

4.2. Numerical Processing of DCIM

The integral term in (43) conforms to the form of 1D Sommerfeld integral (SI). Although the simple numerical integration method is accurate to solve the SI, it is very time-consuming, so it is rarely used. Therefore, we use E-DCIM to perform numerical approximation on the influence of the medium ground:where