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

the


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. erefore, it is particularly important to study methods that can effectively predict the spatial field value under the influence of irregular terrain. e 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. e 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. ese 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. erefore, Wang [6] proposed a hybrid method that uses both the staircase terrain model and the conformal mapping model to solve the irregular terrain problem. e 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. erefore, 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. en, 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. us, 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. erefore, 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. e 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. en, 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.

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. erefore, 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 Ω 1 above the obstacle and the lower subdomain Ω 2 inside the obstacle. e flat terrain area without obstacles can be recorded as Ω 3 . e 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. erefore, the 2W-PE in this area can be solved step by step using DMFT combined with SIBC. e 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. e lower boundary of Ω 2 is the uniform medium ground. e upper boundary of Ω 2 is connected with the lower boundary of Ω 1 , which is exactly the top of obstacles. e 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.

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.
We define the relationship between the total tangential electric field and the total tangential magnetic field on the boundary as where Z s (ω) is the surface impedance. Now, we derive the relationship between Z s (ω) and Fresnel reflection coefficient. e 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 Above the obstacle Ω 1 (ε r1 , σ 1 , k 1 )  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 us, 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 us, we can obtain the SIBC of horizontal polarization: In the case of horizontally polarized waves, there is only one electric field component Ψ � E y , and from the Maxwell Substituting (16) into (15), and according to (4), the boundary condition expression after demodulation can be obtained: 3.2.1. e 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). us, from (11), we can obtain Substituting (18) into (17), we get (19) can be further written as is 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. erefore, the lower boundary of this area can be obtained in the same way.

e 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 Figure 3(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) 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. e electromagnetic wave decays rapidly as the distance of the incident into the free space increases. is contrasts with the actual situation where the electromagnetic wave is totally reflected on the boundary. erefore, (21) is used as the upper impedance boundary condition inside the obstacle.

e 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: e 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: For the flat area Ω 3 around obstacles, the calculation method is the same as that in Ω 1 .

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. e grid division is shown in Figure 4. Denote where z p is the vertical grid point; denote . , x N as the horizontal grid point and the x direction is the stepping direction. Set the midpoint of (x m−1 , z p ) and (x m , z p ) as (ζ m , z p ); thus, the difference International Journal of Antennas and Propagation expression of each partial differential at the point (ζ m , z p ) can be obtained as 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] z zx (27) en, substituting the differential expression (26) into (27), we can get simplified forward propagation equation: where Equation (28) only provides the pth(p � 1, . . . , N − 1) equations. According to (23), we can get the lower boundary condition of Ω 2 : 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 Ω 2 is obtained from (21): zu zz x, N − + κ 2 · u x, N − � 0, us, the corresponding discrete equation on the upper boundary can be written as where us, we form a 0-N forward difference parabolic equation combining (28), (31), and (34): e recursive structure of the difference decomposition of the parabolic equation is expressed as a matrix: where U m is the field vector at x m According to (33), we can obtain the matrixces e above is the case of forward propagation, and similarly, the method can be applied to the case of backward propagation.

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 International Journal of Antennas and Propagation 7 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. e 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. e algorithm needs to complete the calculation once as soon as a new wave source is generated in the calculation process of the same level. us, 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.
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. erefore, 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. at 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.
rough the above theoretical analysis, we can more accurately use 2W-PE to solve the 2D field in the obstacle environment.

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.

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. e tangential field formed on the surface satisfies the PMCHW equation. After the discretization by the basis Step 1 Step 2 Step 3 Step 4 Step 5 Step 6 Step 7 Step 8 Step 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 k 1 , the wave number in the obstacle medium is k 2 , and the wave number in the ground medium is k 3 , the PMCHW equation is obtained as follows: where M 1 and J 1 , respectively, represent the equivalent current and equivalent magnetic current on the obstacle, E i and H i , 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. e operators L 1 and K 1 expressing the electric field generated by the equivalent magnetic current should be decomposed as the sum of operators L 0 1 and K 0 1 expressing the direct field and operators L * 1 and K * 1 expressing reflection field by the ground: where L 0 1 and K 0 1 , respectively, are where X is the equivalent current or magnetic current, ) 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. L * 1 and K * 1 , respectively, are where k 1x and k 2x are the transverse wave numbers in free space and the ground, respectively, k 1z and k 2z are the longitudinal wave numbers in free space and the ground, and R 12 is the reflection coefficient of free space incident on the ground. e 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 L 2 and K 2 . Since the internal equivalent obstacle is a zero field outside, the ground does not need to be considered.
us, the only difference between them and (42) is that the wave number in the operator L 0 1 and K 0 1 is k 1 in free space, and k 2 in the obstacle medium is the wave number used in L 2 and K 2 . According to the principle of duality, the case in which magnetic current generates magnetic field and current generates magnetic field can be obtained.

Numerical Processing of DCIM.
e 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. erefore, we use E-DCIM to perform numerical approximation on the influence of the medium ground: where us, SI can be simplified into the addition of multiple Hankel functions. In (45), a 1n , a 2n 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 R 12 in (44).

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 International Journal of Antennas and Propagation 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.
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 10 − 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. e equivalent radiation power of the radiation source is Z −1 0 w, where Z 0 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. e 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. e 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. is 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. is 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. erefore, 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.
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 10 − 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 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. e 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. erefore, 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. e speed advantage is also obvious.
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 10 − 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 D 1 � 33 m and D 2 � 27 m. e 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. erefore, 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.

Conclusions
is paper presents a new 2W-PE method for solving field value of electric wave propagation when there are obstacles on the medium ground. e 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.
is method is compared with the previous traditional 2W- PE   environment, as well as the superiority in speed. us, 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 e electromagnetic simulation experiment data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.