# Parabolic Equation Modeling of Propagation over Terrain Using Digital Elevation Model

Authors: Xiao-Wei Guan; Li-Xin Guo; Ya-Jiao Wang; Qing-Liang Li
Journal: International Journal of Antennas and Propagation
 (2018)

---

## Abstract

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.


---

## Body

## 1. 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–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 discretely-sampled 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.


## 2. The Parabolic Equation Method

### 2.1. Wide-Angle Parabolic Equation

In the ensuing theory, a time conventione−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
(1)ψx,z=Eyx,zfor horizontal polarization and
(2)ψx,z=Hyx,zfor vertical polarization. It satisfies the 2D scalar wave equation
(3)∂2ψ∂x2+∂2ψ∂z2+k2n2ψ=0,where k is the wave number in free space and n is the relative index of refraction of the medium.Figure 1
Paraxial propagation over irregular terrain.We introduce the reduced function associated with the paraxial directionx
(4)ux,z=e−ikxψx,z,and then the scalar wave equation in terms of u can be expressed by
(5)∂2u∂x2+2ik∂u∂x+∂2u∂z2+k2n2−1u=0.If the backward propagation is ignored, the PE corresponding to forward propagating waves will be(6)∂u∂x=−ik1−Qu,where the pseudo-differential operator Q is defined by
(7)Q=1k2∂2∂z2+n2.According to the operator splitting proposed by Feit and Fleck [9],
(8)1+a+b≈1+a+1+b−1.The error is of the order of the coupled termab. 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 writeQ as
(9)Q=1+A+B,where A and B are defined by
(10)A=1k2∂2∂z2,(11)B=n2−1,then Q can be approximated by
(12)Q≈1+A+1+B−1=1+1k2∂2∂z2+n−1.Thus, the Feit-Fleck WAPE is obtained.(13)∂u∂x=ik1+1k2∂2∂z2−1u+ikn−1u.


### 2.2. Split-Step Fourier Transform Solution

The above WAPE has the formal solution(14)ux+Δx,z=eiΔxk2+∂2/∂z2+kn−2ux,z,and if an appropriate reduced function v is introduced
(15)vx,z=eiΔxk2+∂2/∂z2ux,z,then we have
(16)ux+Δx,z=vx,zeikn−2Δx.The Fourier transform ofv is given by
(17)Vx,p=Fvx,z=∫−∞+∞vx,ze−ipzdz,and according to (15)
(18)Vx,p≅eiΔxk2−p2Ux,p,where U is the Fourier transform of u. Therefore,
(19)vx,z=F−1Vx,p=F−1eiΔxk2−p2Ux,p.Substitute (19) into (16), and we get
(20)ux+Δx,z=eikn−2ΔxF−1eiΔxk2−p2Ux,p.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.


### 2.3. 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 Figure2. The terrain elevation values of the study area are stored in a matrix whose size is determined by the resolution of the DEM.Figure 2
The storage format of the DEM data.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
(21)T′=tanβ,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 kcosβ and the WAPE in (13) now becomes
(22)∂u∂x=ik2cos2β+∂2∂z2⋅u+ikn2−sin2β⋅u,which can also be solved by using the SSFT algorithm. Besides, the phase discontinuity when propagating across a boundary x1,2 is accounted for by
(23)ux1,2,z=ux1,2,zeikzsinβ1−sinβ2.


## 2.1. Wide-Angle Parabolic Equation

In the ensuing theory, a time conventione−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
(1)ψx,z=Eyx,zfor horizontal polarization and
(2)ψx,z=Hyx,zfor vertical polarization. It satisfies the 2D scalar wave equation
(3)∂2ψ∂x2+∂2ψ∂z2+k2n2ψ=0,where k is the wave number in free space and n is the relative index of refraction of the medium.Figure 1
Paraxial propagation over irregular terrain.We introduce the reduced function associated with the paraxial directionx
(4)ux,z=e−ikxψx,z,and then the scalar wave equation in terms of u can be expressed by
(5)∂2u∂x2+2ik∂u∂x+∂2u∂z2+k2n2−1u=0.If the backward propagation is ignored, the PE corresponding to forward propagating waves will be(6)∂u∂x=−ik1−Qu,where the pseudo-differential operator Q is defined by
(7)Q=1k2∂2∂z2+n2.According to the operator splitting proposed by Feit and Fleck [9],
(8)1+a+b≈1+a+1+b−1.The error is of the order of the coupled termab. 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 writeQ as
(9)Q=1+A+B,where A and B are defined by
(10)A=1k2∂2∂z2,(11)B=n2−1,then Q can be approximated by
(12)Q≈1+A+1+B−1=1+1k2∂2∂z2+n−1.Thus, the Feit-Fleck WAPE is obtained.(13)∂u∂x=ik1+1k2∂2∂z2−1u+ikn−1u.


## 2.2. Split-Step Fourier Transform Solution

The above WAPE has the formal solution(14)ux+Δx,z=eiΔxk2+∂2/∂z2+kn−2ux,z,and if an appropriate reduced function v is introduced
(15)vx,z=eiΔxk2+∂2/∂z2ux,z,then we have
(16)ux+Δx,z=vx,zeikn−2Δx.The Fourier transform ofv is given by
(17)Vx,p=Fvx,z=∫−∞+∞vx,ze−ipzdz,and according to (15)
(18)Vx,p≅eiΔxk2−p2Ux,p,where U is the Fourier transform of u. Therefore,
(19)vx,z=F−1Vx,p=F−1eiΔxk2−p2Ux,p.Substitute (19) into (16), and we get
(20)ux+Δx,z=eikn−2ΔxF−1eiΔxk2−p2Ux,p.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.


## 2.3. 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 Figure2. The terrain elevation values of the study area are stored in a matrix whose size is determined by the resolution of the DEM.Figure 2
The storage format of the DEM data.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
(21)T′=tanβ,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 kcosβ and the WAPE in (13) now becomes
(22)∂u∂x=ik2cos2β+∂2∂z2⋅u+ikn2−sin2β⋅u,which can also be solved by using the SSFT algorithm. Besides, the phase discontinuity when propagating across a boundary x1,2 is accounted for by
(23)ux1,2,z=ux1,2,zeikzsinβ1−sinβ2.


## 3. 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 with3° beamwidth located 50 m above the origin. Assuming that f=1.0GHz, Δx=100m, Δy=Δz=λ, and Ny=Nz=1024. Besides, the terrain is perfectly conducting and a 3D Gaussian terrain profile is assumed, which has the analytical form
(24)zx,y=h1e−x−x12+y−y12/2σ12+h2e−x−x22+y−y22/2σ22,where 0m≤x,y≤50km, h1=30m, x1=15km, y1=25km, σ1=1.2×104, h2=40m, x2=40km, y2=25km, and σ2=6×103. Its 3D digital map is shown in Figure 3, and the DEM version with a resolution of 100 m is illustrated in Figure 4.Figure 3
The 3D digital map of the Gaussian terrain model.Figure 4
The DEM of the Gaussian terrain model with a resolution of 100 m.In order to simplify the calculation process and compare the DEM-based method with the shift map, Tx is set to (0 m, 25 km, 50 m) and Rx is set to (50 km, 25 km, 50 m); then the 2D terrain profile between Tx and Rx is given by(25)Tx=h1e−x−x12/2σ12+h2e−x−x22/2σ22,0m≤x≤50km,and thus, the shift map can be carried out.


### 3.1. Validation of the Proposed Method

As Figure5 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=100m. 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 ground-reflected 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.Figure 5
Propagation factor calculated by different methods: (a) the shit map; (b) the DEM-based method.
(a)
(b)The propagation factor at range 50 km is also compared, as shown in Figure6, 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.Figure 6
Propagation factor versus height at range 50 km.


### 3.2. 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 Table1, 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.Table 1
Average error and standard deviation.
Resolution (m)
Average error (dB)
Standard deviation (dB)
10
0.019
0.335
50
0.012
0.377
100
0.006
0.406
200
0.005
0.485
300
0.063
0.721
400
0.061
0.847
500
0.040
1.291
600
0.070
1.326
700
0.062
1.802
800
0.046
1.952Figure 7
Average error and standard deviation versu resolution.


## 3.1. Validation of the Proposed Method

As Figure5 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=100m. 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 ground-reflected 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.Figure 5
Propagation factor calculated by different methods: (a) the shit map; (b) the DEM-based method.
(a)
(b)The propagation factor at range 50 km is also compared, as shown in Figure6, 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.Figure 6
Propagation factor versus height at range 50 km.


## 3.2. 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 Table1, 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.Table 1
Average error and standard deviation.
Resolution (m)
Average error (dB)
Standard deviation (dB)
10
0.019
0.335
50
0.012
0.377
100
0.006
0.406
200
0.005
0.485
300
0.063
0.721
400
0.061
0.847
500
0.040
1.291
600
0.070
1.326
700
0.062
1.802
800
0.046
1.952Figure 7
Average error and standard deviation versu resolution.


## 4. 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 DEM-based 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.



---

*Source: 1878307-2018-03-21.xml*