Parabolic Equation Modeling of Propagation over Terrain Using Digital Elevation Model

The parabolic equation method based on digital elevation model (DEM) is applied on propagation predictions over irregular terrains. Starting from a parabolic approximation to the Helmholtz equation, a wide-angle parabolic equation is deduced under the assumption of forward propagation and the split-step Fourier transform algorithm is used to solve it. The application of DEM is extended to the Cartesian coordinate system and expected to provide a precise representation of a three-dimensional surface with high efficiency. In order to validate the accuracy, a perfectly conducting Gaussian terrain profile is simulated and the results are compared with the shift map. As a consequence, a good agreement is observed. Besides, another example is given to provide a theoretical basis and reference for DEM selection. The simulation results demonstrate that the prediction errors will be obvious only when the resolution of the DEM used is much larger than the range step in the PE method.


Introduction
For many years, the parabolic equation (PE) method has been widely used to model the propagation of electromagnetic waves in complex atmospheric, hydrological, and geographical environments.As a forward full-wave analysis method, PE modeling has attracted a lot of interest since it was first proposed by Leontovich and Fock [1] in 1946 to analyze the diffraction of radio waves around the earth.After that, Hardin and Tappert [2] imported the split-step Fourier transform (SSFT) from underwater acoustic research.Kuttler and Dockery [3] introduced the mixed Fourier transform (MFT), and on this basis, they proposed the discrete mixed Fourier transform (DMFT) [4] in 1996.Since then, this method has been effectively applied.
So far, the three commonly used algorithms of PE are SSFT, the finite-difference method (FDM) [5], and the finite-element method (FEM) [6], among which SSFT is an analytical method, while FDM and FEM are numerical methods.In terms of the problem of propagation in the troposphere under large-scene, long-distance, and low-elevation angle conditions, SSFT has many advantages over FDM and FEM, such as good in stability, small in error, and high in efficiency.In addition, the elevation angle is another critical issue that should be considered, because the accuracy of the PE method could be guaranteed only when the elevation angle is within the scope [7].The approximation accuracy of the standard PE model is limited to be no more than 10 °, so it is also called the narrow-angle parabolic equation (NAPE) method.Therefore, many wide-angle parabolic equation (WAPE) models have been proposed, which include the Claerbout PE method [8], Feit-Fleck PE modeling [9], and the Padé PE method [10].Among these models, the Feit-Fleck PE modeling is the most popular one and its elevation angle can be expanded to approximately 30 °.
In the past decades, a number of scholars [11][12][13][14][15][16][17][18] have studied the problem of propagation over irregular terrain by using the PE method (see [7] for a comprehensive reference list).Ayasli [12] established a terrain-masking model that uniformly zeros the field below the altitude of the boundary at each step.In 1979, Beilis and Tappert [11] developed the terrain-shift transform method, also referred to as the shift map.This method performs a general coordinate transformation for NAPE that flattens the irregular surface, thereby solving the problem through SSFT.Barrios [13] adopted this method to present the terrain PE model, which was applied to actual navy propagation calculations.In 1998, Janaswamy [14] improved this technique and introduced grid interpolation to determine the field distribution in height grids.Donohue and Kuttler [15] subsequently extended the shift-map method to the WAPE, and due to the fact that curvature information may not be available when working with actual terrain problems, they proposed the piecewise linear shift map (PLSM), an alternative method verified as the most accurate approach to handle discretelysampled terrain data.
In recent years, the vectorial version of the PE method has been employed to analyze the field propagation processes in road and railway tunnels where the presence of dissipative and field depolarization processes cannot be neglected.Bernardi et al. [19] used a model based on a hybrid implicit/ explicit finite-difference numerical scheme to solve the vector parabolic equation governing the electromagnetic field propagation in straight and curved tunnels.Zhang et al. [20] presented a hybrid channel modeling technique, which combines the ray tracing and vector parabolic equation methods, to enable the modeling of realistic railway scenarios including stations and long guideways within a unified simulation framework.Compared to the two-dimensional (2D) one, the vectorial version of the PE method can take the presence of lateral propagation environment and field depolarization processes take into account.However, at present, there is no report about three-dimensional (3D) terrain processing method based on SSFT, while irregular terrain modeling is very mature in the 2D PE method.
In this paper, we address the terrain problems in PE modeling and a WAPE method based on the digital elevation model (DEM) is developed using the SSFT algorithm.In fact, as a more advanced modeling approach, the DEM has already been applied in some propagation software programs.However, in the previous models, all the DEMs used are in the geographic coordinate system while in our model, the application of DEM is extended to the 3D Cartesian coordinate system and the effect of the resolution of the DEM on propagation predictions is discussed in order to provide a theoretical basis and reference for DEM selection.

The Parabolic Equation Method
2.1.Wide-Angle Parabolic Equation.In the ensuing theory, a time convention e −iωt is assumed for the fields and is suppressed throughout.As Figure 1 shows, we are concerned with 2D electromagnetic problems where all the fields can be decomposed into horizontally and vertically polarized components propagating independently.The field component ψ is defined by for horizontal polarization and for vertical polarization.It satisfies the 2D scalar wave equation where k is the wave number in free space and n is the relative index of refraction of the medium.We introduce the reduced function associated with the paraxial direction x u x, z = e −ikx ψ x, z , 4 and then the scalar wave equation in terms of u can be expressed by If the backward propagation is ignored, the PE corresponding to forward propagating waves will be where the pseudo-differential operator Q is defined by According to the operator splitting proposed by Feit and Fleck [9], The error is of the order of the coupled term ab.In a vacuum, expression (8) is exact.It provides a good approximation for radio wave propagation application; however, the elevation angle is limited to be approximately 30 If we write Q as where A and B are defined by then Q can be approximated by Thus, the Feit-Fleck WAPE is obtained.
Step Fourier Transform Solution.The above WAPE has the formal solution and if an appropriate reduced function v is introduced then we have The Fourier transform of v is given by and according to (15) where U is the Fourier transform of u.Therefore, Substitute ( 19) into ( 16), and we get It means that the forward propagating field can be obtained at a given range from the field at a previous range, which is quite suitable for numerical computation.

Terrain Modeling Using Digital Elevation
Model.As a digital model was created from terrain elevation data, a DEM is able to provide a 3D representation of a terrain's surface with high efficiency, especially when the environment is complex enough.The storage format of the DEM data in the 3D Cartesian coordinate system is shown in Figure 2. The terrain elevation values of the study area are stored in a matrix whose size is determined by the resolution of the DEM.
Assuming that the transmitting antenna is located at Tx and the receiving antenna Rx, the terrain elevation data between Tx and Rx in a 2D vertical plane can be easily obtained from the DEM, and according to the PLSM [15], the slope over the piecewise linear terrain is given by where β is the angle that the local terrain slope makes with the horizontal.Due to the fact that the tilting of the coordinate system by an angle relative to the surface has locally contracted the horizontal range step, the effective wavenumber in this case may be written as k cos β and the WAPE in ( 13) now becomes which can also be solved by using the SSFT algorithm.Besides, the phase discontinuity when propagating across a boundary x 1,2 is accounted for by , z e ikz sin β 1 −sin β 2 23

Results and Discussions
Two examples are given in this section.In the first one, the DEM-based modeling approach is compared with the shift map to validate its accuracy while the second one deals with the effect of DEM's resolution on prediction errors.In the following simulations, the transmitting antenna is a horizontally polarized Gaussian antenna with 3 °beamwidth located 50 m above the origin.Assuming that f = 1 0 GHz, Δx = 100 m, Δy = Δz = λ, and N y = N z = 1024.Besides, the terrain is perfectly conducting and a 3D Gaussian terrain profile is assumed, which has the analytical form where 0 m ≤ x, y ≤ 50 km, h 1 = 30 m, x 1 = 15 km, y 1 = 25 km, σ 1 = 1 2 × 10 4 , h 2 = 40 m, x 2 = 40 km, y 2 = 25 km, and σ 2 = 6 × 10 3 .Its 3D digital map is shown in Figure 3, and the DEM version with a resolution of 100 m is illustrated in Figure 4.
In order to simplify the calculation process and compare the DEM-based method with the shift map, Tx is set to (0 m, and thus, the shift map can be carried out.
3.1.Validation of the Proposed Method.As Figure 5 shows, the propagation factor is calculated using the shift map and the DEM-based method, respectively.The resolution of the DEM used is 100 m, for the reason that Δx = 100 m.It is easy to find that only a tiny difference is observed.Besides, both of them are able to yield good results in the far zone and the interference patterns of the direct wave and the groundreflected wave can be clearly seen.However, the results are not that satisfactory when the propagation angle is too large, because according to the PE method, electromagnetic waves propagate in a conical region along the range axis.
The propagation factor at range 50 km is also compared, as shown in Figure 6, to make a quantificational evaluation.A good agreement is seen with a mean error of 0.021 dB and a standard deviation of 0.241 dB.

Effect of DEM's Resolution on Prediction Errors.
Another issue of interest is the dependence of the prediction errors on the resolution of the DEM used.As shown in Table 1, DEMs with different resolutions are simulated and the average error and standard deviation are computed.All other parameters are the same with the previous example.Figure 7 shows the results in an intuitive way, and we can see that with the increase of the resolution, the standard deviation increases gradually while the average error changes very slowly and maintains close to zero.Besides, a slow increase in the standard deviation, as expected, is observed when the resolution is less than 100 m, which may be explained by the fact that the prediction errors will be obvious only when the resolution is much larger than the range step in the PE method.

Conclusion
In this paper, under the assumption of forward propagation, a DEM-based PE method is applied on propagation predictions over irregular terrains using the SSFT algorithm.The terrain is assumed to be perfectly conducting, and a 3D Gaussian terrain profile is simulated by using the DEMbased method and the shift map, respectively.As a result, a good agreement is observed.Furthermore, the simulation results indicate that the prediction errors will be obvious only when the resolution of the DEM used is much larger than the range step in the PE method.

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

Figure 2 :
Figure 2: The storage format of the DEM data.

Figure 3 :
Figure 3: The 3D digital map of the Gaussian terrain model.

Figure 4 :Figure 5 :Figure 6 :
Figure 4: The DEM of the Gaussian terrain model with a resolution of 100 m.

Table 1 :
Average error and standard deviation.