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

Authors: An-qi Li; Cheng-you Yin; Qian-qian Zhang; Yong-ji Gan
Journal: International Journal of Antennas and Propagation
 (2021)

---

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


---

## Body

## 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: Section2 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 Figure1. In the case of horizontal polarization and vertical polarization, the field Ψ satisfies the 2-D wave equation:(1)∂2Ψ∂x2+∂2Ψ∂z2+k2n2Ψ=0,where, in the case of horizontal and vertical polarization, Ψ in the 2-D rectangular coordinate system, is defined, respectively, as(2)Ψx,z=Eyx,z=Eyey,Ψx,z=Hyx,z=Hyey,where k is the wave number in the uniform environment of the medium and n represents the modified atmospheric refractive index. It is assumed that the refractive index n does not change with the distance, and the beamwidth of electromagnetic wave propagation is limited, when +x is the propagation direction. Then, after factoring (1), we can get(3)∂∂x+jk1k2∂2∂z2+n2x,z∂∂x−jk1k2∂2∂z2+n2x,zΨ=0.Figure 1
Schematic diagram of wave propagation problems.Introducing two different phase demodulation factors, we can obtain(4)Ψ=Ψfx,z+Ψbx,z=ufx,ze−jkx+ubx,zejkx,ufx,z=ejkxΨfx,z,ubx,z=e−jkxΨbx,z,where ΨfandΨb 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:(5)∂uf∂x=jk1−Quf,∂ub∂x=jkQ−1ub,Q=1k2∂2∂z2+n2x,z.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]:(6)u0,z=12J0,z,VPol,−12M0,z,HPol,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 Figure2, when meeting an obstacle, the area for calculating radio wave propagation can be decomposed into the upper subdomain Ω1 above the obstacle and the lower subdomain Ω2 inside the obstacle. The flat terrain area without obstacles can be recorded as Ω3.Figure 2
Schematic diagram of domain decomposition.The upper boundary ofΩ1 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 Ω3 is similar to that of Ω1, the upper boundary adopts the absorption boundary, and the lower boundary is the medium ground adopting SIBC.The lower boundary ofΩ2 is the uniform medium ground. The upper boundary of Ω2 is connected with the lower boundary of Ω1, which is exactly the top of obstacles. The upper and lower boundaries are obtained in the form of SIBC. Since Ω2 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 Figure3, 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.Figure 3
Schematic diagram of reflection and refraction at the interface of two media ofH pol. (a) Propagation from medium 1 to medium 2. (b) Propagation from medium 2 to medium 1.
(a)(b)We define the relationship between the total tangential electric field and the total tangential magnetic field on the boundary as(7)Etanω=Zsωen×Htanω=ZsωJω,where Zsω is the surface impedance. Now, we derive the relationship between Zsω and Fresnel reflection coefficient.The parameters of medium 1 and 2 are shown in the Figure3, so the wave impedance of medium 1 and 2 can be expressed as(8)Zi=μiε˜i=Z0μriεri11−jσi/ωεi,i=1,2,where ε˜i=εi−jσi/ω,εi=ε0εri,μi=μ0μri,Z0=μ0/ε0. Similarly, by the definition of wave number, the wave numbers in medium 1 and 2 are, respectively,(9)ki=ωμiε˜i=ωμiεi1−jσiωεi,i=1,2,According to Snell’s Law of refraction at the interface of the medium, we can get(10)k1sinθ1=k2sinθ2.Thus, from (9) and (10), we can obtain(11)cosθ2=1−k1k2sinθ12=1−μ1ε˜1μ2ε˜2sin2θ1.According to the propagation of horizontally polarized waves shown in Figure3, the relationship between the reflected and incident fields of horizontal polarization can be obtained as(12)Eyω=Eyiω+Eyrω=Eyiω1+Γhω,Hxω=Hxiω+Hxrω=Hxiω1−Γhω,where(13)Γhω=ErωEiω=EyrωEyiω=Z2cosθ1−Z1cosθ2Z2cosθ1+Z1cosθ2,Zi=EiωHiω=EyiωHxiω/cosθ1.According to the definition of surface impedance in (7), the surface impedance of horizontal polarization can be obtained as(14)Zsh=EyHx=Z1cosθ12Z2cosθ12Z1cosθ2=Z2cosθ2.Thus, we can obtain the SIBC of horizontal polarization:(15)Ey=Z2cosθ2Hx.In the case of horizontally polarized waves, there is only one electric field componentΨ=Ey, and from the Maxwell curl equation∇×E=−jωμH, we can get(16)−jωμHx=∂Ez∂y−∂Ey∂z=−∂Ey∂z.Substituting (16) into (15), and according to (4), the boundary condition expression after demodulation can be obtained:(17)∂ψ∂z=jk2cosθ2ψ⟹∂u∂z=jk2cosθ2u.


#### 3.2.1. The Lower Boundary ofΩ1

Here, we consider an ideal situation, electromagnetic waves are approximately horizontally irradiatedθ1≈90°, which enters the obstacle (medium 2) from the free space (medium 1), as shown in Figure 3(a). Thus, from (11), we can obtain(18)cosθ2≈1−ε˜r1ε˜r2.Substituting (18) into (17), we get(19)∂u∂zx,N+−jk21−ε˜r1ε˜r2⋅ux,N+=0.(19) can be further written as(20)∂u∂zx,N+−jk1ε˜r2ε˜r1−1⋅ux,N+=0.This formula is consistent with the conventional Leontovich impedance boundary, so the lower boundary ofΩ1 adopts this boundary condition. In the area without obstacle Ω3, 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Ω2

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 Figure3(b). We assume that the electromagnetic waves are approximately horizontally irradiated, θ1≈90°, 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 Ω2 can still be written as(21)∂u∂zx,N−−jk2ε˜r1ε˜r2−1⋅ux,N−=0.(21) can be further written as(22)∂u∂zx,N−+k1ε˜r2ε˜r1−1⋅ux,N−=0.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Ω2

We assume that the parameters of the ground medium areε3,μ3,σ3, and the lower boundary of the Ω2 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:(23)∂u∂zx,0+−jk2ε˜r3ε˜r2−1⋅ux,0+=0.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Ω1 with DMFT

In this paper, when solving 2W-PE inΩ1 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:(24)ufx+Δx,z=e−jkn−1Δx⋅F−1e−jkΔx1−kz2/k2−1⋅Fufx,z,ubx−Δx,z=e−jkn−1Δx⋅F−1e−jkΔx1−kz2/k2−1⋅Fubx,z.For the flat areaΩ3 around obstacles, the calculation method is the same as that in Ω1.


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

Considering that the calculation intervalΩ2 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 Figure4. Denote(25)zp=pΔz,p=0,1,…,N,where zp is the vertical grid point; denote x0,x1,…,xm,…,xN as the horizontal grid point and the x direction is the stepping direction. Set the midpoint of xm−1,zp and xm,zp as ζm,zp; thus, the difference expression of each partial differential at the point ζm,zp can be obtained as(26)∂u∂xξm,zp≈uxm,zp−uxm−1,zpΔx,∂2u∂z2ξm,zp≈12Δz2uxm,zp+1+uxm−1,zp+1+uxm,zp−1+uxm−1,zp−1−2uxm,zp−2uxm−1,zp.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 ofQ in (5) is obtained by using Pade (1, 1); then, we rewrite the forward parabolic equation in (6) as [15](27)∂∂x∂2u∂z2+j2k∂2u∂z2+k2n2+3∂u∂x+j2k3n2−1u=0.Then, substituting the differential expression (26) into (27), we can get simplified forward propagation equation:(28)capm+b−2upm+cup+1m+up−1m=c˜apm+b˜−2upm−1+c˜up+1m−1+up−1m−1,where(29)c=1+jkΔx,b=4k2Δz2c,apm=k2Δz2n2−1,c˜=1−jkΔx,b˜=4k2Δz2c˜.Equation (28) only provides the pthp=1,…,N−1 equations. According to (23), we can get the lower boundary condition of Ω2:(30)∂u∂zx,0++κ1⋅ux,0+=0,κ1=−jk2ε˜r3ε˜r2−1⋅ε˜r2ε˜r3,VPol,−jk2ε˜r3ε˜r2−1,HPol.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:(31)a1+b1−2c1d1u0m+2d1u1m=a1−b1−2c1e1u0m−1+2e1u1m−1,where(32)a1=k2n2+3Δz2,b1=jk3n2−1Δz2Δx,c1=1−κ1Δz,d1=1+jkΔx,e1=1−jkΔx.In the same way, the upper boundary condition ofΩ2 is obtained from (21):(33)∂u∂zx,N−+κ2⋅ux,N−=0,κ2=−jk2ε˜r1ε˜r2−1⋅ε˜r2ε˜r1,VPol,−jk2ε˜r1ε˜r2−1,HPol.Thus, the corresponding discrete equation on the upper boundary can be written as(34)a2+b2−2c2d2uNm+2d2uN−1m=a2−b2−2c2e2uNm−1+2e2uN−1m−1,where(35)a2=k2n2+3Δz2,b2=jk3n2−1Δz2Δx,c2=1+κ2Δz,d2=1+jkΔx,e2=1−jkΔx.Thus, we form a 0-N forward difference parabolic equation combining (28), (31), and (34):(36)capm+b−2upm+cup+1m+up−1m=c˜apm+b˜−2upm−1+c˜up+1m−1+up−1m−1,a1+b1−2c1d1u0m+2d1u1m=a1−b1−2c1e1u0m−1+2e1u1m−1,a2+b2−2c2d2uNm+2d2uN−1m=a2−b2−2c2e2uNm−1+2e2uN−1m−1.The recursive structure of the difference decomposition of the parabolic equation is expressed as a matrix:(37)AmUm=BmUm−1,where Um is the field vector at xm(38)Um=u0m⋯uNmT.According to (33), we can obtain the matrixces(39)Am=ca1+b1−2c1d1/c2d1/c000⋯001α1m100⋯0001α2m10⋯00⋮⋮0⋯01αN−1m10⋯02d2/ca2+b2−2c2d2/c,αpm=−2+b+apm,Bm=c˜a1−b1−2c1e1/c˜2e1/c˜000⋯001β1m100⋯0001β2m10⋯00⋮⋮0⋯01βN−1m10⋯02e2/c˜a2−b2−2c2e2/c˜,βpm=−2+b˜+apm.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 Figure5, 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 Figure5, 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.


## 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 Figure2, when meeting an obstacle, the area for calculating radio wave propagation can be decomposed into the upper subdomain Ω1 above the obstacle and the lower subdomain Ω2 inside the obstacle. The flat terrain area without obstacles can be recorded as Ω3.Figure 2
Schematic diagram of domain decomposition.The upper boundary ofΩ1 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 Ω3 is similar to that of Ω1, the upper boundary adopts the absorption boundary, and the lower boundary is the medium ground adopting SIBC.The lower boundary ofΩ2 is the uniform medium ground. The upper boundary of Ω2 is connected with the lower boundary of Ω1, which is exactly the top of obstacles. The upper and lower boundaries are obtained in the form of SIBC. Since Ω2 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 Figure3, 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.Figure 3
Schematic diagram of reflection and refraction at the interface of two media ofH pol. (a) Propagation from medium 1 to medium 2. (b) Propagation from medium 2 to medium 1.
(a)(b)We define the relationship between the total tangential electric field and the total tangential magnetic field on the boundary as(7)Etanω=Zsωen×Htanω=ZsωJω,where Zsω is the surface impedance. Now, we derive the relationship between Zsω and Fresnel reflection coefficient.The parameters of medium 1 and 2 are shown in the Figure3, so the wave impedance of medium 1 and 2 can be expressed as(8)Zi=μiε˜i=Z0μriεri11−jσi/ωεi,i=1,2,where ε˜i=εi−jσi/ω,εi=ε0εri,μi=μ0μri,Z0=μ0/ε0. Similarly, by the definition of wave number, the wave numbers in medium 1 and 2 are, respectively,(9)ki=ωμiε˜i=ωμiεi1−jσiωεi,i=1,2,According to Snell’s Law of refraction at the interface of the medium, we can get(10)k1sinθ1=k2sinθ2.Thus, from (9) and (10), we can obtain(11)cosθ2=1−k1k2sinθ12=1−μ1ε˜1μ2ε˜2sin2θ1.According to the propagation of horizontally polarized waves shown in Figure3, the relationship between the reflected and incident fields of horizontal polarization can be obtained as(12)Eyω=Eyiω+Eyrω=Eyiω1+Γhω,Hxω=Hxiω+Hxrω=Hxiω1−Γhω,where(13)Γhω=ErωEiω=EyrωEyiω=Z2cosθ1−Z1cosθ2Z2cosθ1+Z1cosθ2,Zi=EiωHiω=EyiωHxiω/cosθ1.According to the definition of surface impedance in (7), the surface impedance of horizontal polarization can be obtained as(14)Zsh=EyHx=Z1cosθ12Z2cosθ12Z1cosθ2=Z2cosθ2.Thus, we can obtain the SIBC of horizontal polarization:(15)Ey=Z2cosθ2Hx.In the case of horizontally polarized waves, there is only one electric field componentΨ=Ey, and from the Maxwell curl equation∇×E=−jωμH, we can get(16)−jωμHx=∂Ez∂y−∂Ey∂z=−∂Ey∂z.Substituting (16) into (15), and according to (4), the boundary condition expression after demodulation can be obtained:(17)∂ψ∂z=jk2cosθ2ψ⟹∂u∂z=jk2cosθ2u.


### 3.2.1. The Lower Boundary ofΩ1

Here, we consider an ideal situation, electromagnetic waves are approximately horizontally irradiatedθ1≈90°, which enters the obstacle (medium 2) from the free space (medium 1), as shown in Figure 3(a). Thus, from (11), we can obtain(18)cosθ2≈1−ε˜r1ε˜r2.Substituting (18) into (17), we get(19)∂u∂zx,N+−jk21−ε˜r1ε˜r2⋅ux,N+=0.(19) can be further written as(20)∂u∂zx,N+−jk1ε˜r2ε˜r1−1⋅ux,N+=0.This formula is consistent with the conventional Leontovich impedance boundary, so the lower boundary ofΩ1 adopts this boundary condition. In the area without obstacle Ω3, 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Ω2

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 Figure3(b). We assume that the electromagnetic waves are approximately horizontally irradiated, θ1≈90°, 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 Ω2 can still be written as(21)∂u∂zx,N−−jk2ε˜r1ε˜r2−1⋅ux,N−=0.(21) can be further written as(22)∂u∂zx,N−+k1ε˜r2ε˜r1−1⋅ux,N−=0.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Ω2

We assume that the parameters of the ground medium areε3,μ3,σ3, and the lower boundary of the Ω2 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:(23)∂u∂zx,0+−jk2ε˜r3ε˜r2−1⋅ux,0+=0.The above is the case of horizontal polarization. Similarly, SIBCs on each boundary of vertical polarization can be obtained.


## 3.2.1. The Lower Boundary ofΩ1

Here, we consider an ideal situation, electromagnetic waves are approximately horizontally irradiatedθ1≈90°, which enters the obstacle (medium 2) from the free space (medium 1), as shown in Figure 3(a). Thus, from (11), we can obtain(18)cosθ2≈1−ε˜r1ε˜r2.Substituting (18) into (17), we get(19)∂u∂zx,N+−jk21−ε˜r1ε˜r2⋅ux,N+=0.(19) can be further written as(20)∂u∂zx,N+−jk1ε˜r2ε˜r1−1⋅ux,N+=0.This formula is consistent with the conventional Leontovich impedance boundary, so the lower boundary ofΩ1 adopts this boundary condition. In the area without obstacle Ω3, 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Ω2

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 Figure3(b). We assume that the electromagnetic waves are approximately horizontally irradiated, θ1≈90°, 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 Ω2 can still be written as(21)∂u∂zx,N−−jk2ε˜r1ε˜r2−1⋅ux,N−=0.(21) can be further written as(22)∂u∂zx,N−+k1ε˜r2ε˜r1−1⋅ux,N−=0.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Ω2

We assume that the parameters of the ground medium areε3,μ3,σ3, and the lower boundary of the Ω2 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:(23)∂u∂zx,0+−jk2ε˜r3ε˜r2−1⋅ux,0+=0.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Ω1 with DMFT

In this paper, when solving 2W-PE inΩ1 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:(24)ufx+Δx,z=e−jkn−1Δx⋅F−1e−jkΔx1−kz2/k2−1⋅Fufx,z,ubx−Δx,z=e−jkn−1Δx⋅F−1e−jkΔx1−kz2/k2−1⋅Fubx,z.For the flat areaΩ3 around obstacles, the calculation method is the same as that in Ω1.


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

Considering that the calculation intervalΩ2 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 Figure4. Denote(25)zp=pΔz,p=0,1,…,N,where zp is the vertical grid point; denote x0,x1,…,xm,…,xN as the horizontal grid point and the x direction is the stepping direction. Set the midpoint of xm−1,zp and xm,zp as ζm,zp; thus, the difference expression of each partial differential at the point ζm,zp can be obtained as(26)∂u∂xξm,zp≈uxm,zp−uxm−1,zpΔx,∂2u∂z2ξm,zp≈12Δz2uxm,zp+1+uxm−1,zp+1+uxm,zp−1+uxm−1,zp−1−2uxm,zp−2uxm−1,zp.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 ofQ in (5) is obtained by using Pade (1, 1); then, we rewrite the forward parabolic equation in (6) as [15](27)∂∂x∂2u∂z2+j2k∂2u∂z2+k2n2+3∂u∂x+j2k3n2−1u=0.Then, substituting the differential expression (26) into (27), we can get simplified forward propagation equation:(28)capm+b−2upm+cup+1m+up−1m=c˜apm+b˜−2upm−1+c˜up+1m−1+up−1m−1,where(29)c=1+jkΔx,b=4k2Δz2c,apm=k2Δz2n2−1,c˜=1−jkΔx,b˜=4k2Δz2c˜.Equation (28) only provides the pthp=1,…,N−1 equations. According to (23), we can get the lower boundary condition of Ω2:(30)∂u∂zx,0++κ1⋅ux,0+=0,κ1=−jk2ε˜r3ε˜r2−1⋅ε˜r2ε˜r3,VPol,−jk2ε˜r3ε˜r2−1,HPol.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:(31)a1+b1−2c1d1u0m+2d1u1m=a1−b1−2c1e1u0m−1+2e1u1m−1,where(32)a1=k2n2+3Δz2,b1=jk3n2−1Δz2Δx,c1=1−κ1Δz,d1=1+jkΔx,e1=1−jkΔx.In the same way, the upper boundary condition ofΩ2 is obtained from (21):(33)∂u∂zx,N−+κ2⋅ux,N−=0,κ2=−jk2ε˜r1ε˜r2−1⋅ε˜r2ε˜r1,VPol,−jk2ε˜r1ε˜r2−1,HPol.Thus, the corresponding discrete equation on the upper boundary can be written as(34)a2+b2−2c2d2uNm+2d2uN−1m=a2−b2−2c2e2uNm−1+2e2uN−1m−1,where(35)a2=k2n2+3Δz2,b2=jk3n2−1Δz2Δx,c2=1+κ2Δz,d2=1+jkΔx,e2=1−jkΔx.Thus, we form a 0-N forward difference parabolic equation combining (28), (31), and (34):(36)capm+b−2upm+cup+1m+up−1m=c˜apm+b˜−2upm−1+c˜up+1m−1+up−1m−1,a1+b1−2c1d1u0m+2d1u1m=a1−b1−2c1e1u0m−1+2e1u1m−1,a2+b2−2c2d2uNm+2d2uN−1m=a2−b2−2c2e2uNm−1+2e2uN−1m−1.The recursive structure of the difference decomposition of the parabolic equation is expressed as a matrix:(37)AmUm=BmUm−1,where Um is the field vector at xm(38)Um=u0m⋯uNmT.According to (33), we can obtain the matrixces(39)Am=ca1+b1−2c1d1/c2d1/c000⋯001α1m100⋯0001α2m10⋯00⋮⋮0⋯01αN−1m10⋯02d2/ca2+b2−2c2d2/c,αpm=−2+b+apm,Bm=c˜a1−b1−2c1e1/c˜2e1/c˜000⋯001β1m100⋯0001β2m10⋯00⋮⋮0⋯01βN−1m10⋯02e2/c˜a2−b2−2c2e2/c˜,βpm=−2+b˜+apm.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 Figure5, 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 Figure5, 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 Figure1. 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 k1, the wave number in the obstacle medium is k2, and the wave number in the ground medium is k3, the PMCHW equation is obtained as follows:(40)z2L2−J1−K2−M1−z1L1J1+K1M1=Ei|t,1z2L2−M1+K2−J1−1z1L1M1−K1J1=Hi|t,where M1 and J1, respectively, represent the equivalent current and equivalent magnetic current on the obstacle, Ei and Hi, 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 L1andK1 expressing the electric field generated by the equivalent magnetic current should be decomposed as the sum of operators L10andK10 expressing the direct field and operators L1∗andK1∗ expressing reflection field by the ground:(41)L1=L10+L1∗,K1=K10+K1∗,where L10 and K10, respectively, are(42)L10X=−k04∫τ1+1k02∇∇⋅X⋅H02k0ρdτ′,K10X=j4∫τX×∇H02k0ρdτ′,where X is the equivalent current or magnetic current, ρ=ρ−ρ′=x−x′2+z−z′2, x,z and x′,z′ are, respectively, the coordinates of the field point and the source point, so ρ is the distance between the source point and the field point.L1∗ and K1∗, respectively, are(43)L1∗X=−k04∫τ1+1k02∇∇⋅X⋅1π∫−∞∞e−jk1xx−x′k1zR12e−jk1zz+z′dkxdτ′,K1∗X=j4∫τX×∇1π∫−∞∞e−jk1xx−x′k1zR12e−jk1zz+z′dkxdτ′,where k1xandk2x are the transverse wave numbers in free space and the ground, respectively, k1zandk2z are the longitudinal wave numbers in free space and the ground, and R12 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 byL2 and K2. 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 L10andK10 is k1 in free space, and k2 in the obstacle medium is the wave number used in L2andK2. 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:(44)SI=1π∫−∞∞e−jk1xx−x′k1zR12e−jk1zz+z′dkx≈∑n=1N1a1nH02kρ1n+∑n=1N2a2nH02kρ2n,where(45)ρ1n=x−x′2+z−z′−jα1n2,ρ2n=x−x′2+z−z′−jα2n2.Thus, SI can be simplified into the addition of multiple Hankel functions. In (45), a1n,a2n and α1n,α2n are related parameters of the complex exponential obtained by the generalized function bundle (GPOF) method [16] when calculating the integral kernel function R12 in (44).


## 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 Figure1. 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 k1, the wave number in the obstacle medium is k2, and the wave number in the ground medium is k3, the PMCHW equation is obtained as follows:(40)z2L2−J1−K2−M1−z1L1J1+K1M1=Ei|t,1z2L2−M1+K2−J1−1z1L1M1−K1J1=Hi|t,where M1 and J1, respectively, represent the equivalent current and equivalent magnetic current on the obstacle, Ei and Hi, 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 L1andK1 expressing the electric field generated by the equivalent magnetic current should be decomposed as the sum of operators L10andK10 expressing the direct field and operators L1∗andK1∗ expressing reflection field by the ground:(41)L1=L10+L1∗,K1=K10+K1∗,where L10 and K10, respectively, are(42)L10X=−k04∫τ1+1k02∇∇⋅X⋅H02k0ρdτ′,K10X=j4∫τX×∇H02k0ρdτ′,where X is the equivalent current or magnetic current, ρ=ρ−ρ′=x−x′2+z−z′2, x,z and x′,z′ are, respectively, the coordinates of the field point and the source point, so ρ is the distance between the source point and the field point.L1∗ and K1∗, respectively, are(43)L1∗X=−k04∫τ1+1k02∇∇⋅X⋅1π∫−∞∞e−jk1xx−x′k1zR12e−jk1zz+z′dkxdτ′,K1∗X=j4∫τX×∇1π∫−∞∞e−jk1xx−x′k1zR12e−jk1zz+z′dkxdτ′,where k1xandk2x are the transverse wave numbers in free space and the ground, respectively, k1zandk2z are the longitudinal wave numbers in free space and the ground, and R12 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 byL2 and K2. 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 L10andK10 is k1 in free space, and k2 in the obstacle medium is the wave number used in L2andK2. 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:(44)SI=1π∫−∞∞e−jk1xx−x′k1zR12e−jk1zz+z′dkx≈∑n=1N1a1nH02kρ1n+∑n=1N2a2nH02kρ2n,where(45)ρ1n=x−x′2+z−z′−jα1n2,ρ2n=x−x′2+z−z′−jα2n2.Thus, SI can be simplified into the addition of multiple Hankel functions. In (45), a1n,a2n and α1n,α2n are related parameters of the complex exponential obtained by the generalized function bundle (GPOF) method [16] when calculating the integral kernel function R12 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 is10−5 (S/m), the permeability is 4π×10−7 (H/m), and the relative permittivity is 6; the rectangular obstacle is a medium, of which the permeability is 4π×10−7 (H/m), and the relative dielectric constant is 4; the working frequency of the emission source is 30 MHz, when its height H is 50 m and its width D is 33 m. The equivalent radiation power of the radiation source is Z0−1w, where Z0 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 Figure7(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 Figure7(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.Figure 7
Comparison of total field values in the case of single rectangular medium obstacle ofH Pol. (a) σobstacle=0.1 (S/m). (b) σobstacle=10−5 (S/m).
(a)(b)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 is10−5 (S/m), the permeability is 4π×10−7 (H/m), and the relative permittivity is 6; the conductivity of the obstacle is 10−5 (S/m), the permeability is 4π×10−7 (H/m), and the relative permittivity is 4; other parameters are the same as in Example 1. When the obstacle width D changes, the comparison of the simulation calculation results of the three algorithms is shown in Figure 8.
It can be seen from Figure8(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.Figure 8
Comparison of total field values in the case of single rectangular medium obstacle ofVPol. (a) D = 33 m. (b) D = 40 m.
(a)(b)Example 3.
In the case of horizontal polarization, we assume that electromagnetic waves propagate in a double rectangular obstacles environment as shown in Figure9.
We assume that the upper half of the space is a free space, the ground is a homogeneous medium, the conductivity is10−4 (S/m), the permeability is 4π×10−7 (H/m), and the relative permittivity is 8; the medium obstacle, the permeability is 4π×10−7 (H/m), and relative permittivity is 4; the working frequency of the emission source is 30 MHz, and its height H is 80 m, when the widths of obstacles are D1=33m and D2=27m. The total field distribution calculated by new 2W-PE and MoM is shown in Figure 10.
From Figure10, 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.Figure 10
The field strength distribution in the entire space ofH Pol. (a) new 2W-PE. (b) MoM.
(a)(b)


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



---

*Source: 5519586-2021-04-20.xml*