Combination of the Improved Diffraction Nonlocal Boundary Condition and Three-Dimensional Wide-Angle Parabolic Equation Decomposition Model for Predicting Radio Wave Propagation

Diffraction nonlocal boundary condition (BC) is one kind of the transparent boundary condition which is used in the finitedifference (FD) parabolic equation (PE).The greatest advantage of the diffraction nonlocal boundary condition is that it can absorb the wave completely by using one layer of grid. However, the speed of computation is low because of the time-consuming spatial convolution integrals. To solve this problem, we introduce the recursive convolution (RC) with vector fitting (VF) method to accelerate the computational speed. Through combining the diffraction nonlocal boundary with RC, we achieve the improved diffraction nonlocal BC. Then we propose a wide-angle three-dimensional parabolic equation (WA-3DPE) decomposition algorithm in which the improved diffraction nonlocal BC is applied and we utilize it to predict the wave propagation problems in the complex environment. Numeric computation and measurement results demonstrate the computational accuracy and speed of the WA-3DPE decomposition model with the improved diffraction nonlocal BC.


Introduction
Electromagnetic lateral propagation, which includes lateral scattering, lateral diffraction, and the depolarization effect, affects the actual propagation in irregular terrains, especially in situations with steep transverse gradients.In this situation, the two-dimensional (2D) parabolic equation (PE) is no longer accurate.In order to get more accurate results, the three-dimensional parabolic equation (3DPE) is adopted.The 3DPE has two forms of scalar and vector, which comes from the scalar and vector wave equation, respectively.The scalar 3DPE has only one equation that can just reflect the propagation characteristics of an electromagnetic field in a particular component in 3D space.By contrast, vector 3DPE is composed of three equations, which can include all field components, with their boundary conditions.Levy [1] first used scalar 3DPE to calculate the diffraction problem behind urban buildings.Saini and Casiragh [2] studied the 3D split-step Fourier transformation method, in which a 2D plane Fourier transformation was used, and the plane was perpendicular to the propagation direction.Compared with the 3D finite-difference method, 3D split-step Fourier transformation method decreased computational complexity, but the computational amount was still very large for largescale computing problems.Zelley and Constantinou [3] established a 3DPE model based on implicit finite-difference method.The implicit finite-difference method can deal with the boundary conditions flexibly; therefore, it is relatively simple to solve the radio wave propagation problem on different terrains.However, due to the large-scale matrix operation, the computation amounts are still very large and the demand for resource is high.In recent years, using 3DPE to predict the radio wave propagation in a small area such as in the urban district environment has become a hot topic and has achieved many results [4][5][6].
The gradual improvement of the 3D geographic information system in recent years has caused the large-scale 3D electromagnetic propagation prediction to become a trend.Similar to the case of the 2DPE development, the computational amount of the 3DPE is the main obstacle, especially for radio wave propagation on irregular terrains.The convergence rate depends on the height and lateral width of the computational domain.For large-scale and real-time computation, there still exists a considerable challenge.To accelerate the computational speed of the scalar wide-angle 3DPE (WA-3DPE), we propose a decomposition algorithm.According to this algorithm, the total field at the receiving point is approximately the sum of the straight wave and diffraction wave with the shortest propagation path.In our proposed algorithm, the WA-3DPE can be decomposed into two wideangle 2DPEs which is called a WA-3DPE decomposition model.The WA-3DPE decomposition model can deal with the horizontal diffraction and vertical diffraction caused by irregular terrain obstructions simultaneously.Furthermore this approximation can largely decrease the computational complexity of the WA-3DPE and can accelerate the speed of computation.
Boundary condition (BC) is critical to the use of PE for accurate prediction of radio wave propagation.The 2DPE prediction of radio wave propagation involves two BCs, which are the upper BC and the lower BC, respectively.For the lower BC, before the 1990s, the smooth and perfect conducting BC was adopted.However, the ideal conducting BC caused large errors in practice.Until 1991, the introduction of hybrid Fourier transform method [7] made it possible to apply impedance BC in the solution PE.For non-PEC surfaces and buildings, the fields on them generally meet the impedance BC.For the upper boundary condition, the most used BCs include window function BC [8], perfect matching layers [9], and transparent BC [10,11].For the 3DPE, it contains four BCs.It needs to consider the left and right BCs as well as the upper and lower BCs like the 2DPE.And the most used left and right BCs of 3DPE are the window function BCs [12].
The nonlocal BC is an exact boundary that does not require the additional computational space in the form of absorbing layers.In the previous studies [13], we used the traditional nonlocal boundary conditions in wide-angle two-dimensional parabolic equation model to calculate the propagation loss.However, it has the disadvantage that the 2DPE ignored the lateral scattering and diffraction, and, due to the presence of a spatial convolution integral, the computation is time-consuming.This computation is also memory demanding as it requires the storage and use of all previous values of the field along the boundary.To solve this problem, [14] proposed the recursive convolution (RC) formulation to reduce the computational burden.The RC was achieved using the vector fitting (VF) method.The nonlocal BC combined with RC formulation was applied in a wide-angle 2DPE based on the FD solving method.The computational speed and the accuracy of the nonlocal BC combined with RC formulation were demonstrated.For the actual radio wave propagation problem, the initial source is often fixed at the start position, and there is no need to use the nonlocal BC.Using the simpler diffractive nonlocal BC [15] is enough.Therefore, we adopt the RC formulation to improve the computational speed of the diffraction nonlocal BC.As the diffraction nonlocal BC is a simplified form of the nonlocal BC, the RC formulation can be applied directly from nonlocal BC to diffraction nonlocal BC.The RC formulation is applicable to the nonlocal BC as well as to the corresponding simpler diffractive nonlocal BC.We call the diffraction nonlocal BC combined with RC as improved diffraction nonlocal BC.We will apply this improved diffraction nonlocal BC in the WA-3DPE decomposition model (described in Section 2).
For numerical implementation, we need a suitable approximation of the  operator.With the help of Padé approximation [15], a WA-3DPE is obtained.
(3) Equation ( 3) is called the Claerbout WA-3DPE, which keeps the propagation angle up to 40 degrees.Although the WA-3DPE can model radio wave propagation in 3D space, the computing efficiency is low because of large-scale matrix operations.
In order to accelerate the WA-3DPE speed of computation, similar to the ray tracing method, the main consideration is given to the straight wave and diffraction wave with the shortest propagation path, which can largely decrease the computational complexity of the WA-3DPE.It is assumed that the 3D function (, , ) can be expressed as the sum of two 2D functions: where  ℎ (, ) is propagation waves in (, ) plane and  V (, ) is propagation waves in (, ) plane.By applying (4), (3) can be reduced to two 2DPEs: International Journal of Antennas and Propagation 3 Equations ( 5) and ( 6) are called WA-3DPE decomposition model.The WA-3DPE decomposition model actually contains two wide-angle 2DPEs.The total field at the receiving point is the sum of  ℎ (, ) and  V (, ).

Finite-Difference Method of the WA-3DPE Decomposition Model
Since ( 5) and ( 6) are symmetrical in form, here we only take (5) as an example to introduce finite-difference implementation, and the finite-difference implementation of ( 6) can be obtained in the same way.We start by defining the integration grid, which is fixed in the  direction, but not in the  direction, so that it can adapt to the terrain shape.We let  ℎ = ℎ⋅Δ, ℎ = 0, . . ., ℎ be the grid point at  direction and   =  ⋅ Δ,  = 0, . . ., , . . .be the successive integration ranges.Using the central difference approximation, the derivatives of ( 5) are discretized as follows: In order to advance the solution from range  −1 to range   , ( 7) and ( 8) give the midpoint Equation ( 9) gives the central finite-difference approximation of the derivative in range Equation ( 10) gives the central finite-difference approximation of the second-order derivative in : Equation (11) gives the central finite-difference approximation of the third-order mixed derivative.Substituting ( 8)-( 11) in ( 5) yields where we have used the following notations: Equation (12) has expressed values at range   as a function of values at range  −1 in the form of a linear system.When  = 1, (12) can be written as For the Gaussian source, where  is a normalization constant;  is half-power beamwidth;  is elevation angle; and   is transmitting antenna height.The finite-difference implementation of ( 6) is similar to (5).We let  V = V ⋅ Δ, ( V = 0, . . ., V) be the grid point in  direction and   =  ⋅ Δ, ( = 0, . . ., , . ..) be the successive integration ranges.The difference equation of ( 6) can be written as follows: where we have used the following notations: (17) Equation ( 16) has expressed values at range   as a function of values at range  −1 in the form of a linear system.When  = 1, the values of  0 V,V can be obtained by the same way as  0 ℎ,ℎ .Equations ( 12) and ( 16) are the difference equations of the WA-3DPE decomposition model.In order to complete the equation system, we need to include equations at boundaries.For the WA-3DPE decomposition model, it contains 4 boundaries in ± and ± directions, respectively.Here we employ the 2nd-order accurate discretization of the derivatives  ℎ / and  V /.Thus, at the boundaries, Equations ( 18) give the 2nd-order accurate finite-difference approximation of the boundary conditions at ±: Equations ( 19) give the 2nd-order accurate finite-difference approximation of the boundary conditions at ±.
According to [14,15], the diffraction nonlocal BC in (, ) plane can be written as where the convolution kernel  is given by where  = 1 for  =  max and  =  max ,  = −1 for  =  min and  = 0.  0 is the zero-order Bessel function.

Recursive Convolution.
According to [14], if the function in the convolution term in (20) can be expressed as   exp(  ⋅ ) then Therefore From ( 22) and (23), the following recursive formula can be derived: In general, () is approximately expressed as thus, where  0  = 0.
The integral term in (26) can be discretized as follows: From ( 26) and ( 27), it follows that where we have used the following notations: In (28), the integral term in (20) is expressed by recursive convolution method.

The VF Approximation.
As indicated in (25), the function () must be curve fitted with a function ω() which is a summation of exponential terms.In   exp(  ⋅ ),   and   are complex numbers.To achieve this we use first the VF method [14] to express the function () =  0 () as a sum of exponential terms.In the VF method, the Laplace transform of (), that is, (), is approximated with a sum of fractions   /( −   ): thus Therefore, (21) can be written as From ( 25) and (32), it follows that   =   ,   =  0 (  − ).
According to [14], () is curve fitted over the argument region 0 ≤  ≤ 65000.It was found that  = 20 was sufficient for the computations.The values of   and   can be found in [14].
International Journal of Antennas and Propagation

Free Space Propagation. The biggest advantage of diffraction nonlocal
BCs is that it can be used in conjunction with any initial source.In the following computation, we choose Gaussian source [15] as the initial field.We consider radio wave propagation along the positive direction of -axis in free space, which means that the improved diffraction nonlocal BC can be applied in ± and ± four boundaries.The sizes of computational region are that  max = 5 m,  min = 0 m,  max = 20m,  min = 0m,  max = 100 m, and  min = 0 m.The transmitting antenna is fixed at the point ( = 0,  = 10,  = 2.5).The transmitting antenna frequency is 3 GHz; the 3 dB beamwidth is 5 deg and elevation angle is 0 deg.The refractive index  = 1.In order to verify the computational accuracy and efficiency of the WA-3DPE decomposition model with the improved diffraction nonlocal BCs, we choose the computation results of the standard 3DPE decomposition model [16] as a comparison.As it is accurate when the elevation angle is 0 deg, the standard 3DPE decomposition model computation results can be chosen as a comparison.The standard 3DPE decomposition model was solved by the Split-Step Fourier Transform (SSFT) method with the window function BCs.Figures 1(a Table 1 gives some specific parameters involved in the computation of the WA-3DPE decomposition model and the standard 3DPE decomposition model.The computer CPU is the Intel(R) Core(TM) i5-2400 @3.1 GHz.For the standard 3DPE decomposition model solved by SSFT, four window functions are necessary in ± and ± directions.The window function in the computation is Hanning window function [8].Theoretically, the greater the thickness of the absorbing layer is, the better the absorbing effect is; but for large-scale radio propagation problems in complex environments, the thickness of the absorbing layer cannot be large enough due to the restriction of computation amount.In our computation, the thickness of the absorbing layer is 20 times larger than the computation area in  and  directions.The size of the actual computation area of the standard 3DPE decomposition   1 that the computational speed of the standard 3DPE decomposition model is faster than the WA-3DPE decomposition model.But it is easy to find in Figure 1 that the computational accuracy of the WA-3DPE decomposition model is higher than the standard 3DPE decomposition model.Besides, it can be found that although the discretization size of FD is more restrictive than SSFT, the improved nonlocal BCs can absorb wave completely without increasing extra computation space.For the window function BCs, it is difficult to absorb the wave completely even if the computation space is further increased.Paper [14] had already investigated the improvement of computation speed brought by RC formulation with VF method, so it does not need to be analyzed again.For the above-mentioned 3D problem, the computing time of Dalrymple and Martin method [14] is 18.9054 s.

Knife Edge Propagation.
To verify the effectiveness of the WA-3DPE decomposition model with the improved diffraction nonlocal BCs, knife edge propagation problems are computed and measured.As depicted in Figure 2, two knife edges are center symmetrically situated on a metal plane.The parameters of their size and location are given in Table 2.
The measurements are made by transmitting a continuous wave signal which are generated by the Agilent Signal Generator.The transmitter output power is 0 dBm.The transmitting antenna for the measurements is a quarter wavelength monopole antenna, and the gain of it is ignored.The Agilent Spectrum Analyzer is used as a receiver with very narrow bandwidth.The receiving antenna is a waveguide antenna, having a gain of 0 dBi and working at the frequencies from 8 GHz to 10 GHz.According to the aforementioned properties of experiment equipment, the working frequency is chosen at 9.375 GHz.Both the transmitting antenna and receiving antenna are set on an aluminous plane which is surrounded by wave-absorbing materials to simulate free space propagation environment.
Both of the single knife edge propagation and double knife edges propagation are simulated and measured.In the simulation and experiment of single knife edge propagation, we only keep the knife edge O1 and take the knife edge O2 away.The curves of propagation loss varying with range are given in Figure 3. Figure 3(a) gives the simulation and experiment results of the single knife edge propagation.Then we put on the knife edge O2 and measured again.Figure 3(b) gives the simulation and experiment results of the double knife edges propagation.In order to increase the comparability, we also give the MoM (Method of Moment) and WA-2DPE (wide-angle two-dimensional parabolic equation) computation results in Figure 3.In the simulations, the surface BC is set as PEC, which means the improved diffraction nonlocal BC is applied in up, left, and right boundaries.The results of MoM are obtained by the commercial solver FEKO 7.0.
It can be seen from Figure 3 that the results of the WA-3DPE decomposition model with the improved diffraction nonlocal BCs, MoM, measurement, and WA-2PDE are basically the same no matter for the single knife edge propagation or the double knife edge propagation.However, there exists a big difference in which MoM computation presents fluctuant results while the WA-3DPE decomposition model with the improved diffraction nonlocal BCs and measurement does not.This is because PE only deals with the forward propagation wave, and the signal fluctuation caused by forward propagation and backward propagation interference is not modeled in front of the edges.Theoretically, the measurement should also present fluctuant results in front of the edges, but they are not be shown in Figure 3 in that the receiving antenna in our experiment is a waveguide antenna, which can only receive forward propagation waves.The computation results of WA-3DPE are more closer to measurement results than those of WA-2PDE, as the WA-3DPE considers the horizontal diffraction propagation.

Conclusion
In order to accelerate the computational speed, we proposed to apply the RC formulation with VF method in the diffrac-

Figure 1 :
Figure 1: Pseudocolor maps of field amplitude.Field amplitude of  ℎ in (, ) plane and at  = 2.5 m, computed by (a) the standard 3DPE decomposition model with the window function BCs; (b) the WA-3DPE decomposition model with the improved diffraction nonlocal BCs.Field amplitude of  V in (, ) plane and at  = 10 m, computed by (c) the standard 3DPE decomposition model with the window function BCs; (d) the WA-3DPE decomposition model with the improved diffraction nonlocal BCs.
) and 1(b) are the pseudocolor map of field amplitude of  ℎ in (, ) plane and at  = 2.5 m. Figure 1(a) gives the standard 3DPE decomposition model with window function BC computation results and Figure 1(b) shows the WA-3DPE decomposition model with the improved diffraction nonlocal BC computation results.Figures 1(c) and 1(d) are the pseudocolor map of field amplitude of  V in (, ) plane and at  = 10 m. Figure 1(c) gives the standard 3DPE decomposition model with window function BC computation results and Figure 1(d) shows the WA-3DPE decomposition model with the improved diffraction nonlocal BC computation results.

Table 1 :
Computation parameters for the WA-3DPE decomposition model and the standard 3DPE decomposition model.

Figure 3 :
Figure 3: Comparison of the results of the WA-3DPE decomposition model with the improved diffraction nonlocal BCs with the measurement, MoM, and WA-2PDE results: (a) propagation loss versus range for single knife edge and (b) propagation loss versus range for double knife edges.

Table 2 :
Parameters of knife edge propagation. direction reaches  max = 400 m.It means that the total thickness of the absorbing layer is 380 m and 190 m in the + and − direction, respectively.The size of the actual computation area of the standard 3DPE decomposition model in  direction reaches  max = 100 m.It means that the total thickness of the absorbing layer is 95 m and 47.5 m in the + and − direction, respectively.It can be seen from Table ℎ,  ℎ,  ℎ, 0 0 0 ⋅ ⋅ ⋅  ℎ,  ℎ,  ℎ,