1. Introduction
With the rapid development of marine activities, there has been a growing demand for high-speed and cost-efficient maritime communications [
1,
2,
3,
4,
5]. To achieve space–air–ground–sea integrated communication networks for sixth-generation (6G) communications, it is essential to build a seamless maritime communication network to extend broadband coverage to the sea area [
6,
7]. The LTE-Maritime project proved the feasibility of high-speed wireless communication in coastal areas (up to 100 km from a shore) [
8]. Unlike terrestrial scenarios, the maritime wireless propagation environment has unique meteorological characteristics, such as sea wave fluctuations [
9] and the ducting effect over the sea surface [
10]. Current fifth-generation (5G) cellular networks designed for terrestrial scenarios will be heavily affected by these maritime environmental factors [
11]. Therefore, it is crucial to establish reliable maritime channel models for maritime wireless networks, especially for near-coast wireless networks [
12].
Maritime communications face many application scenarios, such as land-to-ship, ship-to-ship, and buoy-to-ship [
13]. Researchers have recently launched several measurement campaigns to investigate wireless maritime channel properties [
14]. In particular, recent measurements have shown that evaporation ducts and elevated ducts can greatly reduce the transmission loss [
15]. The atmospheric duct is caused by a rapid decrease in atmospheric refractivity over the sea surface [
16]. The trapped signal in the duct layer has an extra gain and can travel hundreds of kilometers [
17]. In addition, atmospheric ducts have the characteristics of wide areas and long durations. The ducting effect usually causes unexpected remote interference and seriously degrades the performance of coastal 5G networks and LTE networks [
18]. More importantly, atmospheric ducts can be utilized to improve wireless coverage and transmission efficiency, especially for land-to-ship communications [
19]. It can also improve unmanned aerial vehicle (UAV)-aided maritime communication when the UAV is deployed at the proper height in atmospheric ducts [
20].
Therefore, the atmospheric duct effect on maritime radio propagation and channel modeling has drawn growing attention from researchers [
21]. Gunashekar et al. measured the 2 GHz maritime channel and analyzed the influence of the antenna height on the ducting effect [
22]. The authors in [
23] measured a near-shore wireless channel at 5 GHz, where the evaporation duct layer is assumed to be horizontally homogeneous. The authors in [
24] conducted a 30 km maritime communication experiment and proposed a ray-tracing model to analyze large-scale fading. In [
25], an 87 km maritime microwave measurement demonstrated that atmospheric ducts can enhance the received signal by more than 40 dB and revealed that the evaporation duct communication is feasible. In [
26], a 107 km measurement at the C band was conducted in the Bohai Sea and results indicated that the distribution of fast fading is close to a Rayleigh distribution.
These empirical models are completed through a large number of actual measurements and useful for maritime communications. Note that the measurements of maritime radio channels require a great deal of financial and human resources. Such statistical methods are not sufficient for fast variations in dynamic maritime environments [
27]. Ray-tracing and parabolic equation methods have better applicability. The authors in [
28] proposed a novel ship–ship round earth loss model by incorporating the ray-tracing method and wind speed. A novel three-dimensional non-stationary UAV-to-ship channel model was proposed, which considered the multi-bounce components introduced by atmospheric ducts [
29].
Recently, parabolic equation methods have been widely used for radio propagation simulations in maritime environmental conditions [
30]. A parabolic equation simulation tool (PETOOL) for evaporation ducts was proposed [
31,
32], which can support a wide variety of evaporation duct models. Using the propagation loss obtained from PETOOL, a vertically separated MIMO ship-to-ship communication system was proposed [
33]. In [
34], the authors proposed a large-scale path loss model for the ducting channel. For rough sea surface conditions, an improved parabolic equation method was proposed based on a pade approximation form [
35]. The fading characteristics in evaporation ducts were studied with average duct heights in the South China Sea [
36]. A ship-to-ship channel model was proposed, which considered the influence of evaporation ducts on the scatterer distribution [
37]. A prediction map of radio propagation was proposed with help of atmospheric duct data in the Korean coastal area [
38]. The authors in [
39] proposed a log-distance maritime wireless channel at 8 GHz, and the meteorological data pointed out that the difference in the evaporation duct profiles at the Tx and Rx increases with increasing link distance. Machine learning has been widely used in the processing of wireless channel measurements [
40]. In [
41], a two-stage deep learning network was designed to characterize duct effects. In [
42], a duct channel model based on a neural network was proposed to predict propagation loss. It should be noted that the training data in [
41,
42] were generated in ideal evaporation ducts. To improve the maritime communication reliability, the authors in [
43] proposed a hybrid method to exploit the combination of cooperative diversity and statistical evaporation duct heights.
In contrast to the homogeneous assumptions in [
30,
31,
32,
33,
34,
35,
36,
37,
38,
39,
40,
41,
42,
43], the duct layers are range-dependent. The strong spatial–temporal variability of the lower atmosphere in near-shore areas can cause horizontal variability in atmospheric ducts. A measurement campaign in the Indian Ocean and the South China Sea has also shown the horizontal variability of atmospheric ducts [
44,
45]. It is obvious that disregarding the range-dependent duct environment will lead to inaccurate prediction results [
46].
To the best of the authors’ knowledge, most existing models cannot reflect the effect of horizontally inhomogeneous ducts. The range-dependent characteristics should be taken into consideration. The decision methods of maritime communication with duct conditions are also lacking.
With the development of Maritime Internet of Things (M-IoT), future maritime wireless communications need to utilize meteorological information to establish sensitive maritime channels that cope with changing environments [
47]. Moreover, the state-of-the-art numerical weather prediction can provide data of high spatial–temporal resolution. This also motivates us to explore a framework of meteorological information-aided maritime radio propagation characterization. Our preliminary study [
48] presents a basic environmental information-aided interference simulation method and verifies the feasibility of meteorological data.
In this paper, we aim to study the effect of horizontally inhomogeneous ducts on maritime radio propagation. Considering the limitations of the existing methods, we propose a sliced parabolic equation algorithm to predict radio propagation characteristics in atmospheric ducts. The major contributions are summarized as follows.
We investigate existing atmospheric duct models and their effects on maritime radio propagation. Subsequently, we propose a range-dependent duct model to support the actual duct environment. The spatial–temporal characteristics of atmospheric ducts are studied by using the meteorological reanalysis data.
We propose a sliced parabolic equation algorithm to improve the prediction accuracy of path loss. Moreover, we derive the corresponding numerical solution and the approximate error expression in inhomogeneous environments.
We validate the algorithm’s effectiveness by the 3.5 GHz field measurement with the actual meteorological data. Compared with the traditional method, the proposed algorithm has higher accuracy during the multiple duct periods.
We evaluate the proposed algorithm in various range-dependent ducts, including evaporation ducts and elevated ducts. Simulation results show that the duct’s horizontal characteristics have a significant impact on path loss.
The rest of the paper is organized as follows. In
Section 2, maritime radio propagation scenarios and three basic duct models are introduced. Based on a range-dependent duct model,
Section 3 presents a sliced parabolic equation algorithm. The numerical solution and feasibility derivation are also provided. In
Section 4, the effectiveness of the proposed algorithm is verified by the 3.5 GHz field measurements.
Section 5 presents the simulation results in evaporation ducts and elevated ducts with different horizontal characteristics. Finally, conclusions are given in
Section 4.
3. Sliced Parabolic Equation Algorithm in Horizontally Inhomogeneous Duct Environments
In this section, we present a range-dependent duct model to make better use of meteorological reanalysis data. Then, we propose a sliced parabolic equation algorithm to simulate electromagnetic wave propagation in the maritime environment. The numerical solution and feasibility derivation are also provided. In the remainder of this section, we will present the workflow and implementation details.
3.1. Range-Dependent Duct Model
In this subsection, we will introduce the range-dependent refractivity and duct environments.
In the ideal case, assume that the meteorological data at different heights and distances can be measured by microwave refractometers and lidar. The refractivity
can be defined as
However, these methods are high-cost. On the other hand, the fast developments of mesoscale numerical weather prediction can provide meteorological data with a high spatial and temporal resolution and make it possible to analyze the horizontal heterogeneity of ducts. It is possible to estimate the refractivity index and duct heights from standard meteorological data. Specifically, meteorological measurements can provide the temperature, humidity, wind speed, pressure, and sea surface temperature at certain heights. The environmental information is used as input parameters. Then, the profiles of temperature, humidity, and air pressure can be estimated with the Monin–Obukhov similarity theory [
55], and the refractivity profiles can be obtained by (7). Finally, duct heights can be estimated.
For example, the fifth-generation European Center for Medium-Range Weather Forecasts (ECMWF) reanalysis data (ERA5) provide hourly products of meteorological variables. Its horizontal resolution is . Therefore, we aim to use the up-to-date meteorological reanalysis product to study the spatial and temporal distributions of the duct environment.
Therefore, we adopt a range-dependent duct model, considering the actual horizontal resolution of meteorological data and computational accuracy. As shown in
Figure 3, we assume that there are
L slices of different atmospheric vertical structures in the horizontal direction. In each slice,
is almost constant in the x direction. The formulas of the range-dependent duct environments are as follows:
3.2. Sliced Parabolic Equation Algorithm
In the rectangular coordinate system, the wave equation satisfies the approximation of the Helmholtz equation,
where
denotes the wavenumber in a vacuum,
n is the refractive index of the tropospheric medium,
x is the propagation distance, and
z is the height. The refractive index
varies with distance and height.
Substituting the simplified function
, the field
is expressed as
The differential operator in (
10) can be factored in terms of two pseudo-differential operators [
35]. Then, the above equation can be rewritten as
where
Q is expressed as
The forward propagation electromagnetic waves are expressed as
The approximate forms of the operator
Q determine the scope of the elevation angle and the accuracy of calculation. Specifically, by defining
,
,
Q is given by
. The pseudo-differential operator of the Feit–Fleck type is given by
Then, (
13) can be rewritten as
The wide-angle parabolic equation of the Feit–Fleck type can be expressed as
The formal solution can be formulated as
The split-step Fourier transform technique is widely used for parabolic equation’s numerical solutions. It requires a large amount of computational power when varies with distance.
In the
L-th slice, we define
, and the temporary variable is expressed as
Next, we aim to solve the differentiation equation, which can be formulated as
By applying the integral mean-value theorem in the i-th slice, it has an approximate formula as
Therefore, the general solution of the differentiation equation is given by
Then, the pseudo-differential operator is formulated as
Define
. The field
can be expressed as follows:
By using Taylor expansion,
can be expressed as
The Fourier transform of
u is defined as
The inverse Fourier transform is given by
Then, we have .
We apply the Fourier transform to both sides of (
24) and have
The inverse Fourier transform of
v is obtained as
By multiplying the refractive index term, the algorithm solution in the current slice is expressed as
The initial field at the starting position (
) is determined by the parameters of the antenna pattern. The initial field can be obtained by fast Fourier transform based on the aperture field and beam pattern. We use the Fourier shift theorems to embody the antenna height and elevation angle. Therefore, the key parameters of the antenna pattern are the height
, the
beamwidth
, and the elevation angle
, respectively. The initial field in the
p domain can be defined as [
31]
The initial field in the spatial z-domain is computed by taking the inverse fast Fourier transform. The Gaussian antenna pattern can represent various antenna types.
In this paper, the horizontally polarized Gaussian antenna pattern can be defined as
where
. The elevation angle is introduced by shifting the antenna pattern, i.e.,
.
At each step, reflects the refraction effect, and reflects the diffraction effect. Note that is assumed in the derivation of the numerical solution. It is essential to analyze the errors caused by the change in refractive index over distance in dynamic maritime environments.
3.3. Feasibility Analysis
In this subsection, we derive the approximate error expressions by the following theorems. With the Taylor expansion, we derive the actual field
u, which can be expressed as
where
, and
.
According to the Monin–Obukhov similarity theory and meteorological accuracy, meteorological conditions in the horizontal direction are consistent in a certain distance. In each slice, we assume that the refractivity profile does not change with range,
. The field at
can be expressed as
By comparing
with
, the relative error caused by the operator
M can be written as
It can be seen that the prediction error of the wide-angle operator mainly depends on the step size and . For horizontally inhomogeneous environments, is almost constant in the x direction in each slice, i.e., and . Thus, . Therefore, the propagation simulation of the whole path can be completed through the calculation of each slice.
3.4. Workflow of the Proposed Algorithm
The workflow of the sliced parabolic equation algorithm is shown in
Figure 4.
In the first step, the input parameters consist of three parts. (1) We need to set the wireless system parameters, including the frequency, antenna height, antenna pattern, beamwidth, and link distance. (2) Algorithm parameters include the max height, max range, horizontal step size, vertical step size, and surface conditions. (3) For the range-dependent duct model, the number of slices depends on the horizontal resolution of the meteorological data. The refractivity profiles of each slice can be directly provided from the reanalysis data, such as the duct height and duct strength. The refractivity profile can also be obtained by a bulk method based on the Monin–Obukhov similarity theory and flux–profile relationship [
48]. The refractivity profile can be computed from the sea surface temperature, wind speed, air temperature, relative humidity, and atmospheric pressure at several sampling heights.
In the second step, we use the proposed algorithm to compute the field strength and path loss in each slice. In the first slice, the initial field strength is based on the antenna pattern. The initial field strength of the second slice is based on the end of the first slice. The full procedure is summarized in Algorithm 1.
Algorithm 1. Sliced Parabolic Equation Algorithm. |
Numerical solution parameters of the proposed algorithm, Wireless system parameters, Duct model parameters. Path loss distribution. atmospheric vertical structures of L slices. each slice Initialize and based on (30)–(31). each step size Compute by (27). Compute by (26). Compute in the L-th slice. Update by (29). Obtain the path loss distribution in the whole computational domain based on (6).
|
4. Field Measurement and Verification
In this section, the accuracy of the proposed algorithm will be verified by 3.5 GHz field measurement results. We first give a brief description, and then we discuss the distribution characteristics of the atmospheric ducts by analyzing the meteorological data. The comparison between the measured and simulated path loss is given in the remainder of this section.
4.1. Measurement Scenario
Two microwave communication links were set up in the near-shore region, and the radio signal strength was recorded from September 2013 to November 2016. The radio measurement campaign was conducted by the Radio Communications Agency Netherlands. The 3.5 GHz band is very important for the 5G wireless communication system. Therefore, the measurement frequency is set to the 3.5 GHz band.
Figure 5 illustrates the geographic information of the two wireless links. The receiver is located at
N,
E, in the Netherlands. The transmitter at Goes is located at
N,
E. The transmitter at Amsterdam is located at
N,
E. The two fixed wireless links are the 138 km Amsterdam–Burum path and the 253 km Goes–Burum path. These links are located in coastal areas and river deltas. The terrain is very flat in the experimental zone. The parameters of the antenna height and frequency are listed in
Table 1.
The path loss of the 253 km link is approximately 220 dB under normal propagation conditions. On the other hand, the measurement system was able to cope with the path loss (approximately 150 dB) due to atmospheric ducts. Therefore, the upper and lower limits of the signal strength have been considered in the measurement system realization. A more detailed setup can be found in [
16]. The RX/TX antenna gain is approximately 26.1 dBi. The vertical 3 dB beamwidth is
. This measurement work improved the 3.5 GHz propagation prediction models and can be applied for 5G deployment. Moreover, the measurement results indicated that short-term abnormal propagation caused by atmospheric ducts exists.
In this paper, we will focus on the effect of range-dependent atmospheric ducts on the radio wave propagation. Note that post-measurement corrections and calibrations were performed. To obtain accurate path loss data, the antenna gains and system losses were removed.
4.2. Distribution Characteristics of Atmospheric Ducts
First, we aimed to characterize the spatial variability of the atmospheric ducts in the experimental zone. We downloaded synchronous meteorological data from the ERA5 atmospheric reanalysis. It can provide the sea surface temperature, 2 m temperature, 10 m wind speed, 2 m relative humidity, and pressure. It provides the duct top height, the duct base height, and the minimum vertical gradient of refractivity inside the duct layer. The horizontal resolution of ERA5 is (approximately 30 km).
Figure 6 shows the distribution characteristics of the atmospheric ducts in the experimental zone for the period of 2013–2016.
Figure 6a illustrates the atmospheric duct probability of the experimental zone. The color represents the probability of duct occurrence. For the wetland-to-sea path (the black line), the coastal zone has a higher probability of duct occurrence. This is because the eastern part of the ocean and the western part of the continent in the trade wind belt are high-probability zones of atmospheric ducts.
The above analysis can be used to identify signal enhancements caused by atmospheric ducts that last for more than three hours. In fact, atmospheric ducts occur more than 50% of the time in many areas, including the Red Sea, the Persian Gulf, and the South China Sea.
Figure 6b shows the mean top heights of atmospheric ducts from September 2013 to November 2016. The color represents the mean top height. The horizontal inhomogeneity of atmospheric ducts is significant in coastal regions. The top heights of the atmospheric ducts in the coastline area are distributed between 200 and 400 m. In coastal areas, the mean top heights of the ducting layer vary significantly. This indicates that the strong spatial and temporal variability of the lower atmosphere in coastal and near-shore areas can cause horizontal variability in the ducting layer. Therefore, the effect of range-dependent atmospheric ducts should be considered carefully if the wireless link is along different latitudes.
4.3. Comparison with Measurements
To further verify the effectiveness of the sliced parabolic equation algorithm, the comparison between the measured and simulated path loss is given in this subsection.
In order to be consistent with the time resolutions of meteorological data, we select the hourly measurement data, such as the path loss at 02:00. The spatial resolutions of slices are based on the meteorological resolution. The refractivity conditions are based on the reanalysis meteorological data grid points along the propagation link. Therefore, the Goes–Burum path is divided into seven slices and the Amsterdam–Burum path is divided into four slices based on the meteorological resolution of the ERA5 data.
Figure 7a shows the Goes–Burum path loss data from August 26th to August 27th, 2016. The path loss of the 253 km link is approximately 220 dB under normal propagation conditions. Atmospheric ducts can occur for up to three hours and the signals were received with approximately 50 dB enhancement. Note that the signal enhancement lasted for approximately 3 hours. An example of the refractivity M-profiles of the Goes–Burum path is plotted in
Figure 7b. The Goes–Burum path is divided into seven slices. Each distance is approximately 31 km. This indicates the appearance of a range-dependent elevated duct. It is obvious that the atmospheric ducts increase with range.
In this paper, we evaluate the proposed algorithm with the following two baselines: (a) the measured path loss and (b) the widely accepted parabolic equation method presented in references [
32,
34,
35,
36,
48], which is based on an ideal homogeneous assumption (only using the refractivity profile near the TX/RX).
The dates of significant signal enhancements caused by atmospheric ducts are listed in
Table 2. The input parameters of the simulation are listed in
Table 3. These simulation parameters and hourly range-dependent refractivity profiles are provided to the algorithm to simulate path loss over the entire height and range. We mainly selected the time period of the occurrence of atmospheric ducts. The discriminant method indicates that the path loss is reduced by more than 40 dB compared to the average path loss. Atmospheric attenuation at 3.5 GHz is small and can be neglected.
Figure 8 shows an hourly path loss comparison of the 253 km Goes–Burum link. The transmitting antenna height is 75 m and the receiver height is 6 m. When an atmospheric duct occurs, the path loss is approximately 170 dB at 253 km. Compared with the multiple measured path loss, the average prediction error of the traditional parabolic equation method is 9.4427 dB. However, the average prediction error of our proposed algorithm is 4.4286 dB. The proposed algorithm shows a good match with the measured path loss. The reasons are as follows. The duct condition on the propagation path varies widely. The sliced parabolic equation algorithm can reflect the inhomogeneity of the propagation path. However, the homogeneous assumption exaggerates the duct effect. Therefore, the effect of range-dependent ducts on electromagnetic wave propagation is significant and should be considered in coastal areas.
Based on the above simulation results, we compute the standard deviation of the duct heights of seven segments along the Goes–Burum path.
Figure 9 shows the prediction error under different standard deviations. We can observe that the proposed algorithm achieves robust performance, while the performance of the traditional parabolic equation method varies widely. This is because the traditional method considering the homogeneous assumption is inaccurate and unstable in actual and dynamic atmospheric ducts.
Figure 10a compares the prediction performance on the 138 km Amsterdam–Burum path. The path loss under standard atmospheric conditions is approximately 202.31 dB. When the atmospheric duct occurs, the path loss is approximately 165 dB. We observe that the simulation results of the proposed algorithm follow the measured data closely in most cases. The proposed algorithm’s average error is 3.1161 dB, while the average prediction error of the traditional method is 6.1810 dB. Since this path is on the sea surface or land, the horizontal inhomogeneity of the duct layers on the Amsterdam–Burum path is obvious. Most of the predicted values of the traditional method are unstable since it neglects the range-dependant atmospheric conditions in near-shore areas.
Figure 10b shows the relationship between the prediction performance in
Figure 10a and the standard deviation of the duct heights. As one can see, the standard deviation of the duct height on the Amsterdam–Burum path is smaller than that on the Goes–Burum path. When the standard deviation is more than 15 m, the sliced parabolic equation algorithm can achieve better performance.