**1. Introduction**

The fifth-generation (5G) wireless system is in the rapid development stage of China's mobile telecommunication standards. Moreover, 5G relies on unique advantages, such as high speed, low latency, large capacity, and wide connectivity to rapidly integrate applications in various fields of production and life, and its application scope will eventually range from mobile broadband services to next-generation car and device networking services [1]. Meanwhile, in the sixth-generation mobile communication technology planning, satellite internet communication will play an important role in realizing the transition from the interconnection of everything to the intelligent connection of everything. With the rapid development of coastal cities, the surrounding sea–land junction meteorological environment and the distribution of ground objects will change and become more complex. The modeling requirements for the forward-propagation of communication signals when mixed ducts and non-uniform turbulence occur are more accurate.

**Citation:** Yang, K.; Guo, X.; Wu, Z.; Wu, J.; Wu, T.; Zhao, K.; Qu, T.; Linghu, L. Using Multi-Source Real Landform Data to Predict and Analyze Intercity Remote Interference of 5G Communication with Ducting and Troposcatter Effects. *Remote Sens.* **2022**, *14*, 4515. https://doi.org/10.3390/rs14184515

Academic Editor: Gang Zheng

Received: 20 July 2022 Accepted: 1 September 2022 Published: 9 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The 5G network must coordinate high, medium, and low frequencies according to user service coverage and speed requirements to form an efficient and coordinated layered coverage. According to public data, by the end of 2020, operators have opened a total of 718,000 5G base stations, and the number of 5G mobile terminal connections has exceeded 200 million. The dense distribution of base stations has higher requirements for accurate reception and anti-interference of communication links. We selected three sub-6 GHz in the 5G plan for the simulation, and evaluated the weather-induced intercity CCI that may have occurred in different urban links when ducts and tropospheric turbulence existed.

An atmospheric duct is mainly caused by the increase of atmospheric temperature with the increase of altitude or the rapid decrease of water vapor density with the increase of altitude. Turton [2] proposed five synoptic processes that can form atmospheric ducts—radiative temperature inversion, subsidence temperature inversion on land surface, temperature advection inversion at the sea–land interface, topographic inversion in local valleys, and frontal temperature inversion are prone to occur. According to data statistics, communication disconnections and other phenomena often occur at night after rain in the plains of China. The monitoring found that it was related to the interference of ducts. This paper selected urban links in the East China Plain (where ducts are prone to occur) to analyze the CCI of over-the-horizon communications.

In 1946, Leontovich and Fock [3] proposed a forward full-wave analytical method—the parabolic equation (PE) method—to solve the problem of forward-propagation in vertically layered inhomogeneous tropospheric ducts. In 1973, Hardin and Tappert [4] developed a split-step Fourier transform (SSFT) algorithm suitable for underwater acoustic transmission to solve parabolic equations. Subsequently, the PE model has been theoretically deduced and optimized by a series of scholars. Based on its framework, the Feit–Fleck wide-angle parabolic equation (WAPE) [5] method that can calculate a larger range of grazing angles was proposed. The grazing angle can reach 30◦. In 1983, Thomson and Chapman [6] formally proposed a new WAPE form based on the Feit–Fleck method, and Kuttler [7] verified the good performance of WAPE at wide elevation angles.

The PE method takes into account the refraction and diffraction effects of radio wave propagation. Moreover, it can simply and accurately describe the dielectric parameters of the irregular surface and the refractive index distribution of non-uniform atmospheric structures to simulate the propagation characteristics of radio waves in a complex environment. The step-by-step iterative solution method can predict the propagation loss of point-topoint links, and can further realize regional propagation loss prediction. In 2000, Levy [8] made a relatively complete summary of the derivation, terrain, and boundaries of electromagnetic wave propagation calculated by the parabolic equation method, and preliminarily deduced the three-dimensional parabolic equation algorithm. In 2003, Janaswamy [9] used 3D vector parabolic equations to predict path loss on flat terrain with buildings. In 2005, Hu et al. [10] established a macrocell radio propagation model for mobile communication based on 3DPE. In 2019, Rasool et al. [11] deduced the solution of 3DPE, and used 3DPE for field strength prediction for flat and irregular forest environments. In this paper, the PE model was improved to simulate the possible CCI of the over-the-horizon transmission of the communication signal on the intercity link.

When electromagnetic waves encounter irregular terrain during transmission, such as wedge-shaped boundaries, more intense reflection will occur, so the modeling of complex terrain is also particularly important. In 1987, Dockery and Konstanzer [12] first used PE for tropospheric wave propagation problems. In 1996, Kuttler and Dockery [13,14] proposed the method of continuous hybrid Fourier transform under the impedance boundary to solve PE, which greatly improved the efficiency of the SSFT algorithm. In 2000, Donohue and Kuttler [15] applied the SSFT PE algorithm to irregular terrain, extended the terrain shift map algorithm proposed by Beilis and Tappert to WAPE, and deduced a new impedance boundary condition for electromagnetic wave incidents on the surface of a finite conductor, such that the WAPE solution could be adapted to the previous MFT algorithm. Numerical examples demonstrate its effectiveness on impedance surfaces with terrain slopes of 10–15◦

and discontinuous slope changes of 15–20◦. To accommodate terrains with larger slopes, it is recommended to combine the terrain shift map with terrain masking (knife-edge diffraction) approximations. In this paper, the WAPE-based terrain shift map algorithm and the new impedance formula are applied to the real inter-city link distribution to evaluate the atmospheric loss distribution.

Since the influence of troposcatter cannot be ignored, scholars have begun to turn their focus to the influence of atmospheric ducts on troposcatter communications. In 1993, Hitney [16] proposed the scattering parabolic equation method based on the radio physical optics (RPO) model. This method can simply and effectively calculate the tropospheric scattering caused by atmospheric turbulence on the basis of PE. In 2009, Ivanov et al. [17] used this method to further study the boundary layer leakage of tropospheric scattering and the propagation characteristics in evaporation ducts. In 2013, Wang et al. [18] carried out experiments with three unequal-distance cross-sea circuits to study the signal fading characteristics of tropospheric scattering in the sea area where evaporation ducts occur. The evaporation duct phenomenon was found to significantly alter the fading level of the signal.

In 2015, Song et al. [19] found that the atmospheric meteorological state has a grea<sup>t</sup> influence on the scattering loss through the study of the scattering parabolic equation. In 2016, Wanger et al. [20] evaluated the refractive index of the lower atmosphere and gave a propagation loss (PL) measurement; that is, using height-independent PL measurements over a range of 10–80 km to infer information about the existence and potential parameters of atmospheric ducts in the lowest 1 km of the atmosphere. Its main improvement includes turbulence-induced distance-dependent fluctuations in the forward-propagation model, where the inhomogeneous turbulence model uses parameters that are height-dependent on the mixing ducts. In this paper, Hitney's RPO basic scattering model was applied to the transmission over inter-city links with evaporation ducts, and the non-uniform turbulence model proposed by Wanger was used for transmission calculations over inter-city links with mixed ducts.

In this paper, combined with open source data from multiple sources, the probability of occurrence of ducts in the selected area was verified by WRF. The land types in inter-city links were extracted based on MCD12Q1 from MODIS. The empirical values of complex dielectric parameters suitable for these land types were found from various references or measured values. In the end, we substituted impedance correction, real terrain relief, and turbulence correction into intercity links to evaluate the probability of possible interference when ducts and turbulence occurred between different cities.

### **2. A Multi-Source Real Data Model for Intercity Link Communication Propagation Computing**

### *2.1. WRF Models Regional Duct Distribution*

The new generation of the mesoscale advanced research WRF model system is a flexible and advanced atmospheric simulation system. It is suitable for a wide range of electromagnetic applications from meters to thousands of kilometers. The WRF modeling system is shown in Figure 1 and mainly consists of the following main programs: 1. WRF pre-processing system (WPS); 2. WRF data assimilation (WRF-DA) assimilation system; 3. advanced research WRF (ARW) solver; 4. post-processing and visualization tools.

We used the WRF pre-processing system (WPS) to specify the static terrain data path of the processing area and perform interpolation on the terrain data (including terrain, land use, and soil type). The background field and boundary field data used in WPS were NCEP FNL (final) business global analysis data [21,22]. The global meteorological data in GRIB2 (Gridded Binary2) format were downloaded. The grid data resolution was 1◦ × 1◦ and the time interval was 6 h.

The latitude and longitude coordinate ranged of the selected area were (25.846◦N– 36.224◦N, 117.472◦E–123.528◦E). We set the number of coarse grids and fine grids, and their relative positions. There were several cities that we simulated in the fine grids,

of which, the latitude and longitude coordinates were Wuxi (120.3166◦E, 31.5072◦N), Hangzhou (120.1602◦E, 30.2856◦N), Nanjing (118.7988◦E, 32.0651◦N), Shanghai (121.4802◦E, 31.2389◦N), and Jiaxing (120.7635◦E, 30.7509◦N), respectively. Using the Lambert conformal projection, the simulation area was used as a coarse grid, and the key selected urban area was a fine grid. The schematic diagram of the simulation area is shown in Figure 2, and the parameter settings are shown in Table 1. After the initialization parameters were set, we executed *geogrid.exe* to define the model domain to be processed, and imported the static terrain data into the grid; secondly, we commanded *ungrib.exe* to extract the meteorological fields from the imported grid-formatted terrain fields; finally, we executed *metgrid.exe* to horizontally interpolate the meteorological fields extracted by ungrib.exe in the second step into the grid model defined by *geogrid.exe* in the first step. The vertical interpolation of the meteorological field to the WRF eta level was realized in the actual program.

**Options Settings Nested Levels Double Nesting Coarse grid settings** Longitude and latitude coordinate of grid center (120.5◦E, 31.035◦N) Number of grids 192 × 384 Horizontal resolution 3 km × 3 km **Fine grid settings** Position in the coarse grid (32, 20) Number of grids 123 × 246 Horizontal resolution 1 km × 1 km

**Figure 1.** Procedures of the WRF pre-processing system.

The background field data of 1 March 2018, 1 June 2018, 1 September 2018, and 1 December 2018UTC were selected, respectively, to simulate the duct distribution in the selected plain area by seasons. The duct distribution located in East China is shown in Figure 3.

As shown in Figure 3, it can be seen that the ducts were distributed in the sea area throughout the year, while the inland ducts were less distributed and scattered in March. The inland ducts were densely distributed in June and December, and the inland ducts were also widely distributed in September, but not continuous. Due to the large amount of data to be analyzed in the selected area, only the first day of each quarter was selected for numerical simulation. The actual occurrence of the duct was also related to sudden changes in weather, such as rainfall and snowfall, or the diurnal temperature variation. However, based on the numerical simulation results of WRF in the East China Plain, it is fully proven that the ducting phenomenon occurs in this area all year around. Therefore, the surface duct with the base layer is used to evaluated the CCI of the intercity in East China.

**Figure 2.** The nested grids of the simulated areas in WRF modes.

Figures 4 and 5—we calculated the duct height and duct intensity distributions on 1 June and 1 December 2018 using WRF mesoscale numerical modeling software.

As shown in Figure 4, the average distribution of duct strengths between cities in June is more than 60 m, the intercity link from Shanghai, Jiaxing, to Wuxi also has a hybrid duct phenomenon of up to 240 m, and its strength can reach up to 35 M-units or so, which is shown in Figure 5a. In December, the marine ducts in this area were densely distributed, and the heights and strengths of the ducts were significantly higher than those in coastal cities. In inland areas, the duct heights on the Nanjing–Wuxi link were basically distributed at about 40–80 m, while duct heights on the intercity link of Shanghai–Wuxi and Jiaxing–Wuxi could reach more than 200 m; the average duct strength was about 20 M-units. Considering the limitation of the installation height of the communication signal antenna, based on the duct distribution results simulated by the real data using WRF, this paper studied the over-the-horizon transmission of the 5G communication signal on the intercity link when the non-uniform surface-based duct with the trapping layer height was 40–80 m.

### *2.2. Extraction of Digital Elevation Data on Intercity Link*

This paper used DEM digital elevation data with a resolution of 90 m for the analysis and calculation. The data set was provided by geospatial data cloud site, Computer Network Information Center, Chinese Academy of Sciences (http://www.gscloud.cn). We accessed the data in 17 July of 2022. The data came from SRTM (Shuttle Radar Topography Mission) and were jointly measured by the National Aeronautics and Space Administration (NASA) and the National Imagery and Mapping Agency (NIMA).The data volume of radar images acquired by the SRTM system was about 9.8 trillion bytes. After various preprocessing steps, such as editing, gross error removal, determining water surface elevation, and defining coastlines, etc., the results were cut into 15,000 pieces according to geographic coordinates. Each map frame was stored in latitude and longitude, and a digital elevation model (DEM) was made. The data accuracy that could be used for scientific research was the STRM-3 version, which was sampled every 3 arc seconds, and the data of 3 arc seconds were collected and generated from 1 arc second of data. They were divided in 1◦ × 1◦ of the latitude and longitude grid in the geographic projection, and the horizontal resolution was about 90 m. SRTM terrain data obtained by CIAT (the International Center for Tropical Agriculture) using a new interpolation algorithm, this method better fills the data gap of SRTM 90. The interpolation algorithm is from Reuter et al. [23].

**Figure 3.** The distribution of ducts in selected plain areas by seasons, where the yellow area represents the surface duct and the red area represents tge elevated duct. (**a**) March. (**b**) June. (**c**) September. (**d**) December.

SRTM-3 data are divided into files every 5 degrees of longitude and latitude, with a total grid data of 24 rows (including the area that latitude ranges from 60◦N to 60◦S), 72 columns (longitude ranged from −180◦W to 180◦E). According to DEM data product, the elevation data of several intercity links from the area of interest (25.846◦N–36.224◦N, 117.472◦E–123.528◦E), could be determined with a resolution of 90 m. Interpolating to fill the null zone, the global elevation distribution of the selected area is obtained as shown in Figure 6.

We extracted the link elevations from four surrounding cities to city A. Considering that the minimum distance between cities was about 90 km, which was about 0.81◦ when it was represented by the longitude and latitude coordinates. To not reduce the accuracies of DEM data on intercity links, we extracted the required link DEM data using traversal and interpolation. The processing steps are in Figure 7.

**Figure 4.** Duct height (m) calculated by WRF on 1 June and 1 December 2018. (**a**) Duct height (m) on 1 June 2018. (**b**) Duct height (m) on 1 December 2018.

**Figure 5.** Duct Strength (M-units) calculated by WRF on 1 June and 1 December 2018. (**a**) Duct strength (M-units) on 1 June 2018. (**b**) Duct strength (M-units) on 1 December 2018.

### *2.3. Data Extraction of Land Cover Distribution Types of Intercity Links from MODIS* 2.3.1. Land Cover Type Distribution on Intercity Links of East China

The MODIS Land cover product MCD12Q1 provides a series of 13 sets of scientific data sets with 500 m of resolution global land cover projections, which were published in an annual span, including five different land cover classification schemes of IGBP/UMD/LAI/BGC/PFT, and a new three-layer classification standard based on the land cover classification system (LCCS) from the Food and Agriculture Organization, including a quality assurance (QA) layer, posterior probability of the LCCS three-layer standard, and binary land water cover. MCD12Q1 was secondary validation datum based on cross-validation of the training dataset generated from projections.

MCD12Q1 was derived from observations from the Terra and Aqua data inputs and was generated using MODIS reflectance data supervised classifications. The primary land cover program identified 17 land cover categories as defined by the International Geosphere– Biosphere Program (IGBP), including 11 natural vegetation categories, 3 developed and planted land categories, and 3 non-vegetated land categories, as listed in Table 2.

**Figure 6.** DEM of the selected area.

**Figure 7.** Flow chart of the extracting line-distributed DEM data from polygon-distributed DEM data.


**Table 2.** IGBP Land Cover Type.

In this paper, the vegetation distribution in the latitude and longitude range areas was converted into WGS84 projection by positive selection projection and expressed in Figure 8:

**Figure 8.** MCD12Q1-LC1 distribution of the selected region.

We converted the original sinusoidal projection of the MCD12Q1 to WGS84 projections. The original data of MCD12Q1 was tile distribution, which was similar to that of DEM. After converting the projection of MCD12Q1, the original square tile became a diamond distribution. Then the required intercity land type data were extracted from the adjacent converted diamond tile data. The value of the land type of the latitude and longitude of the intercity link was uniformly interpolated according to the equidistant latitude and longitude interval. Finally, we traversed the data coordinates in the shorter dimension direction (longitude direction or latitude direction), found the point closest to the longitude and latitude slope of the link, and performed a one-dimensional intercity link land type on

the standard of the IGBP arrangement. The land cover types classified in the IGBP of the links from the four surrounding cities to Wuxi, marked in Figure 8, are shown in Figure 9:

**Figure 9.** Land Cover types of IGBP on intercity links from four different cities to Wuxi. (**a**) LC distribution on Hangzhou–Wuxi. (**b**) LC distribution on Nanjing–Wuxi. (**c**) LC distribution on Shanghai–Wuxi. (**d**) LC distribution on Jiaxing–Wuxi.

2.3.2. Determination of Corresponding Dielectric Parameters in the Area Covered by the Intercity Link

The classification of the IGBP features was relatively fine, and there have been many studies on the empirical model and test results of several kinds of vegetation dielectric permittivities. After consulting a large number of literature studies, we chose to use the test results in the literature to solve the simulation in this paper. The land types involved in several intercity links selected in this paper mainly include: 9 Savannas, 10 grasslands, 11 permanent wetlands, 12 croplands, 13 urban and built-ups, 15 snow and ice, 16 barren or sparsely vegetated, 17 water bodies.

According to [24], we chose the following experimental results: the permittivity of 18.2 + 4.5i of sandy loam from a frequency of 3 GHz and water content of 0.3 g/cm<sup>3</sup> as the permittivity of barren or sparsely vegetated; the no. 16 classification in the IGBP land cover (LC) rule. Silty clay in Figure 10 was selected as the no. 11 classification in the IGBP LC rule when the frequency was 1.4 GHz and the water content *mv* = 0.45; that is, the permittivity of the permanent wetland WET was assumed as 26.5 + 1.2i.

By referring to Reference [25], the crops with large soil fertility contribution rates distributed in the selected area (north China region) were mainly rice. The dielectric permittivities of rice leaves of 50 rice growth days of were selected as the permittivities of crops (CRO) in the IGBP classification (i.e., 4.0 + i0.93). Moreover, the permittivity of

rice straw grown for 10 days was taken as the permittivity of the savanna (SAV) in IGBP classification, which was 2.85 + i0.33.

**Figure 10.** Dielectric permittivity of wet soil at different volumetric moisture or frequencies. (**a**) Dielectric constant vs. Volumetric Moisture. (**b**) Dielectric constant vs. frequency.

According to Reference [26], combined with the test results, the permittivity of asphalt in the oil-generating area was generally = 5∼15, = 0.070∼0.30. Then we chose the dielectric permittivity of the asphalt of the intermediate state as the permittivity of urban and built-up (URB) areas in the IGBP classification, which was 5 + i0.070.

In Reference [27], from the relationship that the complex dielectric permittivity of wet snow versus water content or frequency at natural conditions (as shown in Figure 11), we choose the snow and ice (SNO, 15th classification in IGBP LC) dielectric permittivity 15.8 − i1.8, when the frequency was 3 GHz, the water content was 0.3, and the temperature was −5 ◦C.

**Figure 11.** Dielectric permittivity of snow at different volumetric moisture/frequencies. (**a**) Dielectric constant vs. Volumetric moisture. (**b**) Dielectric constant vs. Frequency.

The reference dielectric permittivity values used for the land cover types ultimately involved in this paper are listed in Table 3.

This paper used the WRF database source, MCD12Q1 land cover classification from MODIS, terrain relief from the DEM to model the over-the-horizon propagation of communication signals at a sub-6 frequency band, comprehensively considering atmospheric ducts and tropospheric inhomogeneous turbulence. The schematic diagram of the over-thehorizon propagation of intercity links is shown in Figure 12.


**Table 3.** Dielectric constant used on simulated intercity links.

**Figure 12.** Schematic diagram of over-the-horizon forward-propagation of intercity links based on multi-source data fusion.

### **3. Scattering Parabolic Equation Algorithm for Irregular Terrain and Inhomogeneous Turbulent Atmosphere**

When 5G signals are transmitted between cities, irregular terrain relief and dielectric roughness exist on the lower boundary of the atmospheric duct. At the same time, there may be non-uniform turbulence in the duct layer and even its upper and lower boundaries. This paper considers the actual situation and studies the radio wave propagation based on the complex meteorological environment when ducting and turbulence occurs. Compared with many existing terrain models, after research and comparison, a wide-angle shiftmap (WASM) algorithm was determined in the terrain boundary processing of the WAPE model. The over-the-horizon propagation effects caused by duct considering the irregular terrain, different surface dielectric parameters at the lower boundary, and non-uniform tropospheric turbulence at the upper boundary of the duct layer were simulated on intercity links. The simulation results were compared with the results of numerical modeling software, AREPS (advanced refractive effects prediction system).

The tropospheric atmosphere uses the refraction index n to describe the refraction characteristics of electromagnetic waves. The propagation of electromagnetic waves on the earth's surface needs to consider the influence of the earth's curvature. The atmospheric refraction index *n* is expressed by the modified refraction index *m* as

$$m(h) = n(h) + h/a,\tag{1}$$

The modified refractive index *M* is expressed as

$$M = (m - 1) \times 10^6 = N + 10^6 \times h/a = N + 157h,\tag{2}$$

When simulating the tropospheric duct propagation, the refractive index profile used a hybrid four-parameter three-line model to represent a hybrid duct with a base layer:

$$M(z) = \begin{cases} M\_0 + k\_1 z & 0 \le z \le z\_1 \\\ M\_1 + k\_2 (z - z\_1) & z\_1 < z \le z\_1 + z\_2 \\\ M\_2 + 0.118 (z - z\_1 - z\_2) & z > z\_1 + z\_2 \end{cases} \tag{3}$$

where *z*1 is the height of the duct base layer and *z*2 is the duct thickness.

When tropospheric ducts appear over land due to changes in meteorological conditions, the electromagnetic wave wavelength and emission elevation angle satisfy the following conditions:

Maximum cut-off wavelength that can trap the wave in the duct

$$
\lambda\_{\text{max}} = \mathbb{C} \int\_{h\_{\text{L}}}^{h\_{\text{U}}} [M(h) - M\_{\text{U}}]^{\frac{1}{2}} dh \tag{4}
$$

where *M*U is the modified refractive index of the top of the duct trapping layer. When studying the propagation through the duct structure with the base layer, the above formula is simplified as

$$
\lambda\_{\text{max}} = \frac{2}{3} \text{Ct} \sqrt{\triangle M} \tag{5}
$$

where *t* is the duct thickness in *m* and is the modified refractive index difference corresponding to the duct layer height and the surface.

The critical trapping angle is simplified by derivation using Snell's law, as:

$$
\theta\_c \approx \sqrt{2 \times 10^{-6} \cdot (M\_s - M\_{\rm U})} \tag{6}
$$

When the wavelength of the electromagnetic wave was less than the cut-off wavelength, and the emission elevation angle was within the trapping angle, the electromagnetic wave was trapped in the duct and propagated over a long distance beyond the horizon.

In the two-dimensional Cartesian coordinate system, it was assumed that the electromagnetic wave propagated along the x-axis, the electromagnetic field component was independent of the y-axis, the refractive index *n* of the propagating medium was isotropic, and the field component satisfied the two-dimensional scalar wave equation.

$$\frac{\partial^2 \psi(\mathbf{x}, z)}{\partial \mathbf{x}^2} + \frac{\partial^2 \psi(\mathbf{x}, z)}{\partial z^2} + k\_0^2 + n^2 \psi(\mathbf{x}, z) = 0 \tag{7}$$

where *k*0 is the electromagnetic wavenumber in vacuum, *ψ* is the electronic or magnetic field. *ψ*(*<sup>x</sup>*, *z*) = *Ey*(*<sup>x</sup>*, *z*) for horizontally polarized, and *ψ*(*<sup>x</sup>*, *z*) = *Hy*(*<sup>x</sup>*, *z*) for vertically polarized.

We numerically approximate the wave equation to obtain a wide-angle parabolic equation (WAPE) [5] form suitable for computing irregular terrain

$$\frac{\partial u(\mathbf{x},z)}{\partial \mathbf{x}} = ik \left[ \sqrt{1 + \frac{1}{k^2} \frac{\partial^2}{\partial z^2} - 1} \right] u(\mathbf{x},z) + ik [u(\mathbf{x},z) - 1] u(\mathbf{x},z) \tag{8}$$

The piecewise linear terrain shift map (PLTSM) algorithm proposed by Donohue and Kuttler was used to model the terrain boundary of the PE propagation model. The D–K algorithm is based on the continuous terrain shift map (CTSM) algorithm of Beilis–Tappert, which improves the difficulty of calculating the terrain curvature data.

The core idea of the PLTSM algorithm is to establish a new coordinate system on the irregular terrain section, and translate the PE calculation of each step into a flat ground form for solution. The terrain translation coordinates are expressed as

$$\mathbf{x} = \mathbf{u}, \mathbf{z} = \mathbf{v} - T(\mathbf{u}) \tag{9}$$

Let the field components expressed in the new coordinate system be equal to the field components expressed in the irregular terrain coordinate system *φ*(*<sup>x</sup>*, *z*) = <sup>Φ</sup>(*<sup>u</sup>*, *<sup>v</sup>*). We introduce simplified functions *φ*(*<sup>x</sup>*, *z*) = *<sup>e</sup><sup>i</sup>θψ*(*<sup>x</sup>*, *<sup>z</sup>*). According to the coordinate map formula and wave equation, the forward-propagation PE considering irregular terrain is obtained

$$\left[ \left\{ \frac{\partial}{\partial \mathbf{x}} + i \frac{\partial \theta}{\partial \mathbf{x}} - T' \left( \frac{\partial}{\partial z} + i \frac{\partial \theta}{\partial z} \right) \right\} - i \sqrt{\left( \frac{\partial}{\partial z} + i \frac{\partial \theta}{\partial z} \right)^2 + k^2 n^2} \right] \psi = 0 \tag{10}$$

In the piecewise linear terrain, the phase function is *<sup>θ</sup>*(*<sup>x</sup>*, *z*) = *kzT*(*x*) + *f*(*x*). Let *f* (*x*) = *k*(1 + *T*)/2, we obtain the WAPE in the coordinate system expressed in terms of the terrain curvature:

$$\frac{\partial}{\partial \mathbf{x}} = i\sqrt{k^2 + \frac{\partial^2}{\partial z^2}} \cdot \boldsymbol{\psi} + ik(n - zT'')\boldsymbol{\psi} \tag{11}$$

In this paper, digital elevation model (DEM) data with 90 m resolution were used, which were downloaded from the National Geospatial data cloud. Taking into account the resolution, it was more accurate to use the terrain slope to calculate the long-distance propagation. To avoid dealing with the terrain curvature directly, the PLTSM method was used. The phase factor is expressed in the form of

$$\theta(\mathbf{x}, z) = \mathbf{K}\_0 z T^\prime(\mathbf{x}) + f(\mathbf{x}) \tag{12}$$

Let

$$f'(x) = \mathbb{K}\_0(T^{\prime 2} - 1) \tag{13}$$

where *K*0 is the undetermined constant. Substituting Equations (11)–(13) into the forwardpropagation PE (10) after the coordinate transformation, the simplified approximation of PE can be obtained:

$$\frac{\partial \psi}{\partial x} = i \sqrt{\frac{k^2}{1 + T'^2} + \frac{\partial^2}{\partial z^2}} \cdot \psi + ik \sqrt{n^2 - \frac{T'^2}{1 + T'^2}} \cdot \psi \tag{14}$$

The above equation is the WAPE corresponding to the PLTSM method. We abbreviate this method as TWPE (terrain wide-angle parabolic equation).

The atmospheric refraction profile is expressed using the modified refraction index *<sup>m</sup>*(*<sup>x</sup>*, *z*) = *<sup>n</sup>*(*<sup>x</sup>*, *z*) − 1 + *zT* in Equation (14); Equation (14) is solved by the split-step Fourier transform (SSFT) method. The relationship between the transforming field and the physical field could be expressed as *φ* = *<sup>e</sup><sup>i</sup>θψ*. For the segmented irregular terrain surface, the terrain inflection point changes due to different terrain slopes on adjacent segments, while the phase factor on adjacent segments need to be processed continuously to satisfy the realistic terrain simulation. The required phase correction item when the signal reflects from segmen<sup>t</sup> 1 to segmen<sup>t</sup> 2 is as follows:

$$\begin{split} \psi\_2(\mathbf{x}\_{1,2}, z) &= \psi\_1(\mathbf{x}\_{1,2}, z) \exp\left\{ ik \left[ \left( \frac{T\_1'}{1 + T\_1^2} \right) - \left( \frac{T\_2'}{1 + T\_2^2} \right) \right] \right\} \\ &= \psi\_1(\mathbf{x}\_{1,2}, z) \exp\left[ ikz (\sin \beta\_1 - \sin \beta\_2) \right] \end{split} \tag{15}$$

where *<sup>x</sup>*1,2 is the reflection point from segmen<sup>t</sup> 1 to segmen<sup>t</sup> 2; *T* = tan *β*, where *β* is the angle that the local terrain slop makes with the horizontal. In this modification of the phase factor, this phase item makes the controlling factor of height-dependent sequences of transforming filed expressed by *ψ*(*<sup>x</sup>*1,2, *jdz*). The controlling angle happens to be the physics angle corresponding to the slope difference across the reflection point *<sup>x</sup>*1,2. Then when controlling the phase of the transform field, in order to suppress the side lobes on the height sequence, the height step size must satisfy the following condition:

$$dz < \frac{2\pi}{k(1 + \sin \theta\_M)}\tag{16}$$

where *θM* is the maximum angle of the controlling factor. The limiting condition on the height step size above is much more strict than the Nyquist sampling theorem. Considering the impedance boundary condition suitable for irregular terrain, Donohue and Kuttler deduced the impedance boundary condition on the terrain surface, which PLTSM needs on the basis of the standard Leontovich impedance boundary condition. It is expressed as follows: 

$$\frac{\partial \psi}{\partial z} = -ik \cos \beta \left[ \sin \gamma \left( \frac{1 - \Gamma}{1 + \Gamma} \right) + \tan \beta (1 - \cos \gamma) \right] \psi \tag{17}$$

where Γ is the Fresnel reflection coefficient, which forms differently according to different polarizations on radio waves. Moreover, *γ* = *ξ* + *β*, where *γ* represents the grazing angle on the terrain slope, and *ξ* represents the angle that the propagation direction makes with the horizontal. Compared with the standard impedance boundary condition, the new impedance factor could be expressed as:

$$\alpha' = \cos \beta [\alpha + \tan \beta (1 - \cos \gamma)] \tag{18}$$

In addition, we used an inhomogeneous turbulence model on the upper boundary of the duct with the based layer to express the complex tropospheric turbulence. Wanger et al. [20] considered the turbulence-induced range-dependent fluctuation of the lowest 1 km from the surface in the forward-propagation model of the duct. The maximum likelihood (ML) estimation of the atmospheric refractive index had good accuracy, and with prior information about ducting, the maximum a prior (MAP) refractivity estimate could be found.

Based on the empirical model of the three-line four-parameter refractive index profile of the duct with a base layer, the influence of the non-uniform turbulence was superimposed. The refractive index fluctuation is expressed as follows, and is related to the height of the duct base layer *m*3 and the duct thickness *m*4:

$$m\_{ih}(z) = \bar{n}(z)\left\{1 + K \exp\left[\frac{-(z - (m\_3 + m\_4/2))^2}{(m\_4/2)^2}\right]\right\} \tag{19}$$

where *n*˜ is generated by formula *C*2*n* = *Cs*, for the relationship between the structure constant *C*2*n* and refractive index, see [20].

### **4. Propagation Loss Calculation Using the Improved TWPE Model**

In this section, we simulate intercity link ducting propagation with turbulence. The terrain relief on links were extracted from real DEM data. The LC types are classified based on the IGBP standards from MODIS data. Some dielectric permittivities of LC types relative to the intercity links are referenced from different literature studies, and measured or modeled on empirical formulae. The frequency was 3 GHz. The antenna height was 60 m, the antenna beam width was 3◦, and the antenna elevation angle was 0◦. The surface-based duct with turbulence in the upper boundary of duct is assumed here. We describe the surface-based duct with the tri-fold four-parameter model, of which the four-parameter combination is [40 40 0.12 −0.2].

In this intercity link from Hangzhou to Wuxi, the propagation distance was as far as 136 km. The relief from Hangzhou to Wuxi reached more than 200 m, which exceeded the radio-wave elevation angle limits of WPE (so was the intercity link from Nanjing to Wuxi). As the results show in Figures 13 and 14, the obstacles greatly increased the propagation loss. The radio-wave could not travel further with the presence of very high obstacles. The effects of ducting and turbulence were cut off. Therefore, no CCI would occur on the intercity link when there is high topography.

**Figure 13.** The (**a**) DEM and (**b**) propagation loss with tropospheric scattering, IGBP LC, and DEM considered on Hangzhou–Wuxi. (**a**) DEM form Hangzhou to Wuxi. (**b**) Propagation loss on Hangzhou–Wuxi.

**Figure 14.** The (**a**) DEM and (**b**) propagation loss with tropospheric scattering, IGBP LC, and DEM were considered on Nanjing–Wuxi. (**a**) DEM from Nanjing to Wuxi. (**b**) Propagation Loss on Nanjing–Wuxi.

Taking into account the above analysis, the Hangzhou–Wuxi and Nanjing–Wuxi intercity links had large obstacles regarding cutting-off the ducting propagation. Therefore, the following selected the relatively gentle Shanghai–Wuxi and Jiaxing–Wuxi links to the intercity links. When the duct occurred on the intercity link, it was simulated whether to consider the non-uniform tropospheric turbulence or the influence of the PL distribution— of the distribution of real land covers—on the possible intercity CCI.

### *4.1. Joint Effects of Tropospheric Turbulence and Duct*

The same parameters as the previous section are used in this section to simulate the situation of Shanghai–Wuxi and Jiaxing–Wuxi without considering tropospheric turbulence.

It can be seen from the comparison in Figures 15b and 16b of the range-dependent PL at the antenna level (60 m) that the non-uniform turbulence in the troposphere mainly affected the middle and rear sections of the intercity link. When considering the joint effect

of the tropospheric inhomogeneous turbulence and the duct, in the range of 0–40 km, the tropospheric turbulence made no obvious effects on ducting propagation. At the end of the intercity link of 40–70 km, the tropospheric turbulence effect reduced the PL of the duct. After 70 km of the intercity link, the joint effect of tropospheric inhomogeneous turbulence and ducting on PL was reduced relative to the middle part compared with the ducting-only propagation situation.

**Figure 15.** The PL on Jiaxing–Wuxi (**a**) in space (**b**) at a height of 60 m. (**a**) PL on Jiaxing–Wuxi without the troposcatter. (**b**) PL on J–W at a height of 60 m with or without the troposcatter.

**Figure 16.** The PL on Shanghai–Wuxi (**a**) in space (**b**) at a height of 60 m. (**a**) PL on Shanghai–Wuxi without the troposcatter. (**b**) PL on S–W at a height of 6 0m with or without the troposcatter.

Based on this, it can be concluded that the non-uniform turbulence model superimposed on the non-uniform tri-fold four-parameter hybrid duct has a boosting effect on the inter-city propagation at 40 km and further. This may lead to the fact that, on the urban link from surrounding cities to Wuxi, such as the Shanghai–Wuxi link, the downlink signal of the remote base station will likely interfere with the uplink signal at the first 20–60 km of the link, affecting the random uplink access user. After 60 km, or due to sudden changes in terrain, or being too far away, the PL of the two links at the antenna height of 60 m exceeds 160 dB, which means the sensitivity that can be distinguished by the receiver. Therefore, the possibility of remote interference is not considered, but simulated by the improved scattering PE model, the PL of this middle section of the path link is reduced, the signal of the transmitting section may propagate to the user area of this section, and remote interference may still occur. When considering the farthest distance that useful

signals can propagate between intercity links, the PL calculated by adding the influencing factors of tropospheric turbulence in this section also shows that it is possible to deploy base stations as precise as possible based on ensuring accurate communications and no remote interference. For example, a base station relay can be set up at 20 km where similar terrain and LCs are distributed, such as Shanghai–Wuxi, and a base station relay can be assumed at 60 km where similar terrain and LCs are distributed, such as Jiaxing–Wuxi link. In addition, because the two intercity links are close together, if Wuxi receives a useful signal from the Jiaxing link under the influence of a stronger duct, it may also receive a remote signal from Jiaxing at the same time, causing user communication interference. It could be concluded that the influence of tropospheric turbulence has an important reference value for deploying base stations and avoiding remote interference between base stations.
