International Journal of Antennas and Propagation
Volume 2021 (2021), Article ID 5519586, 13 pages
https://doi.org/10.1155/2021/5519586

Research Article

Research on Fast Algorithm of Radio Wave Propagation in Low-Lossy Obstacles Environment

An-qi Li, Cheng-you Yin, Qian-qian Zhang, and Yong-ji Gan

College of Electronic Engineering, National University of Defense Technology, Hefei, China

Correspondence should be addressed to Cheng-you Yin; cyouyin@sina.com

Received 20 January 2021; Revised 9 March 2021; Accepted 13 March 2021; Published 20 April 2021

Academic Editor: Giuseppe Castaldi

Copyright © 2021 An-qi Li et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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, as

where 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 obtain

where 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)

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 as

where 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 as

where . 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 as

where

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

where 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 as

where

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, are

where 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, are

where 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

Thus, SI can be simplified into the addition of multiple Hankel functions. In (45), and are related parameters of the complex exponential obtained by the generalized function bundle (GPOF) method [16] when calculating the integral kernel function in (44).

5. Results and Discussion

Here, we apply the new 2W-PE proposed in this paper to the obstacle environment we set and compare it with MoM and the traditional 2W-PE proposed in [5]. We assume that the single rectangular obstacle environment under the medium ground situation is shown in Figure 6. All examples here were run on a computer equipment environment of Core i7 processor (1.80 GHz), 64-bit operating system, and 8 GB memory.

Figure 6: Simulation diagram of single rectangular obstacle.

Example 1. In the case of horizontal polarization, we suppose the upper half of the space is free space, the ground is a homogeneous medium, of which the conductivity is  (S/m), the permeability is  (H/m), and the relative permittivity is 6; the rectangular obstacle is a medium, of which the permeability is  (H/m), and the relative dielectric constant is 4; the working frequency of the emission source is 30 MHz, when its height is 50 m and its width is 33 m. The equivalent radiation power of the radiation source is , where is the free space wave impedance. Changing the conductivity of the obstacle medium, the simulation results of three algorithms at an observation height of 50 m are shown in Figure 7.
It can be seen from the result of Figure 7(a), when the conductivity is large, the field values of the three have a higher degree of agreement. The fluctuation of the field value in front of the obstacle is caused by the superposition and interference of the incident wave and the reflected wave from the surface of the obstacle and the ground. It can be seen from the detailed diagram that the field value calculated by 2W-PE and MoM has only a small difference in amplitude when the phases are almost the same. The field value inside the obstacle is almost zero, which is due to the large conductivity, and the electromagnetic wave decays rapidly after entering the obstacle. Afterwards, due to the reduction of the direct wave and the reflected wave behind the obstacle without the obstacle, the fluctuation of the field value after the obstacle is small.
From Figure 7(b), it can be found that the new 2W-PE method we proposed is in good agreement with MoM, while the error of traditional 2W-PE is large in this situation. This is because when the conductivity of the obstacle is weak, the electromagnetic wave will not quickly decay to zero after entering the obstacle; on the contrary, there are multiple reflections at the front and back edges of the obstacle and the process of transmitting the obstacle, so the field value and phase around the obstacle have changed greatly compared with Figure 7(a). It can also be seen from Figure 7(b) that compared to the difference in field values in the front and rear areas of the obstacle, the difference of the new 2W-PE and MoM inside the obstacle is relatively larger. This may be due to the fact that the radiation source has a certain beamwidth which causes the wave to not be incident perpendicular to the surface of the obstacle. Therefore, it can be concluded from Figure 7 that the results of our new 2W-PE and MoM are consistent with each other regardless of whether the conductivity of the obstacle is large or small, while the original traditional 2W-PE method is more suitable for the situation where the obstacle is a high-lossy medium.
In the calculation process, MoM takes 3460 s on average, and the new 2W-PE algorithm takes 13 s on average, so the new 2W-PE algorithm has obvious advantages in calculation speed.

(a)

(b)

Figure 7: Comparison of total field values in the case of single rectangular medium obstacle of H Pol. (a)  (S/m). (b)  (S/m).

Example 2. In the case of vertical polarization, we suppose the upper half of the space is free space, the ground is a homogeneous medium of which the conductivity is  (S/m), the permeability is  (H/m), and the relative permittivity is 6; the conductivity of the obstacle is  (S/m), the permeability is  (H/m), and the relative permittivity is 4; other parameters are the same as in Example 1. When the obstacle width changes, the comparison of the simulation calculation results of the three algorithms is shown in Figure 8.
It can be seen from Figure 8(a) that when D = 33 m, since the electromagnetic wave propagation inside the obstacle is not considered, the amplitude result of the traditional 2W-PE has a large error. The phase deviation of traditional 2W-PE also occurs in the detailed picture, while the new 2W-PE and MoM are always in good agreement. In the same way, the field value fluctuates due to interference between the direct wave and the reflected wave in front of the obstacle; since there are more creeping waves on the surface of the obstacle of vertical polarization, small interference ripples occurs in the field value behind the obstacle in this case compared with the case of horizontal polarization. It can be seen from Figure 8(b) that since the width of the obstacle is an integer multiple of the wavelength in the obstacle medium at this time, a phenomenon similar to the “radome” occurs, and no transmitted wave inside the obstacle passes through the obstacle. Therefore, the internal field value of the obstacle increases, and the amplitude of the front and back field values decreases compared to Figure 8(a). In addition, comparing Figures 8(a) and 8(b), the traditional 2W-PE does not consider the reflection inside the obstacle, and the results obtained in the above different environment with varying widths are completely the same. Of course, there is no phenomenon of “radome” in the result of traditional 2W-PE, which further illustrates that the field value calculated by traditional 2W-PE under low-lossy conditions will have a large error.
In the calculation process, MoM takes an average of 2410 s, and the new 2W-PE algorithm takes 14.92 s. The speed advantage is also obvious.

(a)

(b)

Figure 8: Comparison of total field values in the case of single rectangular medium obstacle of . (a) D = 33 m. (b) D = 40 m.

Example 3. In the case of horizontal polarization, we assume that electromagnetic waves propagate in a double rectangular obstacles environment as shown in Figure 9.
We assume that the upper half of the space is a free space, the ground is a homogeneous medium, the conductivity is  (S/m), the permeability is  (H/m), and the relative permittivity is 8; the medium obstacle, the permeability is  (H/m), and relative permittivity is 4; the working frequency of the emission source is 30 MHz, and its height is 80 m, when the widths of obstacles are and . The total field distribution calculated by new 2W-PE and MoM is shown in Figure 10.
From Figure 10, we can clearly see the reflection and refraction process of electromagnetic waves in the double rectangular obstacles environment. Therefore, the consistency of the calculation results of the two in the entire space is good as shown in Figure 10, so we can consider the new 2W-PE to be a reliable and efficient calculation method.

Figure 9: Simulation diagram of double rectangular obstacles.

(a)

(b)

Figure 10: The field strength distribution in the entire space of H Pol. (a) new 2W-PE. (b) MoM.

6. Conclusions

This paper presents a new 2W-PE method for solving field value of electric wave propagation when there are obstacles on the medium ground. The surface impedance boundary conditions on the boundaries of the subdomains above and inside the obstacle are deduced in detail for different regional boundaries. After that, DMFT and FD are used, respectively, in solving 2W-PE in two subdomains. This method is compared with the previous traditional 2W-PE and the calculation results of MoM based on E-DCIM, and the numerical accuracy of 2W-PE in the case of medium ground is verified by examples in the rectangular obstacle environment, as well as the superiority in speed. Thus, this paper provides a novel and reliable fast algorithm to effectively solve the problem of rapid calculation of radio wave propagation especially when there are steep low-lossy obstacles in the medium ground situation and the electrical parameters of the obstacles are different from that of the ground. Combined with the staircase terrain model, it can be further extended to predict the field value in an environment with more complex undulating terrain.

Data Availability

The electromagnetic simulation experiment data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

References

P. Zhang, L. Bai, Z. Wu, and L. Guo, “Applying the parabolic equation to tropospheric groundwave propagation: a review of recent achievements and significant milestones,” IEEE Antennas and Propagation Magazine, vol. 58, no. 3, pp. 31–44, 2016.

D. A. Mcnamra, C. W. I. Pistorius, and J. A. G. Malherb, Introduction to the Uniform Geometrical Theory of Diffraction, Artech House, Norwood, MA, USA, 1990.

R. L. Roberts and J. J. Daniels, “Modeling near‐field GPR in three dimensions using the FDTD method,” Geophysics, vol. 62, no. 4, pp. 1114–1126, 1997.

O. Ozgun, “Recursive two-way parabolic equation approach for modeling terrain effects in tropospheric propagation,” IEEE Transactions on Antennas and Propagation, vol. 57, no. 9, pp. 2706–2714, 2009.

K. Wang and Y. Long, “Propagation modeling over irregular terrain by the improved two-way parabolic equation method,” IEEE Transactions on Antennas and Propagation, vol. 60, no. 9, pp. 4467–4471, 2012.

D. D. Wang, X. L. Xi, Y. R. Pu et al., “LF groundwave propagation modeling over irregular terrain by the hybrid two-way parabolic equation method,” IEEE Transactions on Antennas and Propagation, vol. 67, no. 10, pp. 90–98, 2019.

Q. F. Wei, C. Y. Yin, and W. Wu, “Research and verification for an improved two‐way parabolic equation method in obstacle environment,” IET Microwaves, Antennas & Propagation, vol. 12, no. 4, pp. 576–582, 2018.

M. Levy, Parabolic Equation Methods for Electromagnetic Wave Propagation, IEE Press, London, UK, 2000.

M. A. Leontovich and V. A. Fock, “Solution of the problem of propagation of electromagnetic waves along the Earth’s surface by method of parabolic wave equation,” Journal of Physics, vol. 10, pp. 13–23, 1946.

W. C. Gibson, The Method of Moments in Electromagnetics, Chapman & Hall/CRC, London, UK, 2008.

G. Apaydm and L. Sevgi, “Groundwave propagation at short ranges and accurate source modelling,” IEEE Transactions on Antennas and Propagation, vol. 55, no. 3, pp. 244–262, 2013.

D. Dockery and J. R. Kuttler, “An improved impedance-boundary algorithm for Fourier split-step solutions of the parabolic wave equation,” IEEE Transactions on Antennas and Propagation, vol. 44, no. 12, pp. 1592–1599, 1996.

Z. Y. Ge, G. Z. Lu, and R. D. Wang, “Using the finite difference parabolic equation with nonlocal boundary condition for complex environment radiowave propagation prediction,” in Proceedings of the 11th International Symposium on Antennas, Propagation and EM Theory (ISAPE), pp. 478–480, Guilin, China, October 2016.

C. Y. Yin, J. Zhu, and Q. F. Wei, “Study on calculation and verification of radiowave propagation using parabolic equation for the antenna near the ground,” in Proceedings of the 37th Progress in Electromagnetics, pp. 655–1659, Shanghai, China, August 2016.

W. W. Ding, K. Wang, and Y. L. Long, “Study of the electromagnetic waves propagation over the improved fractal sea surface based on parabolic equation method,” International Journal of Antennas and Propagation, vol. 2016, Article ID 8235752, 7 pages, 2016.

T. K. Sarkar and O. Pereira, “Using the matrix pencil method to estimate the parameters of a sum of complex exponentials,” IEEE Antennas and Propagation Magazine, vol. 37, no. 1, pp. 48–55, 1995.