Abstract
A two-way parabolic equation (2W-PE) method based on complex-valued neural networks (CVNNs) and a physics-informed neural network (PINN) is proposed to predict the spatial field in the environment of dielectric obstacles with high efficiency and accuracy. In the framework of the 2W-PE method, CVNNs are used to calculate the lumped transmission field and the lumped reflection field considering the influence of the obstacle, thus avoiding the long calculation time caused by the internal multilayered reflection processes. The incident directions and field strength of the waves on the regional boundaries vary greatly with the propagation environment, so coefficients of the boundary conditions are obtained by using the PINN. Next, the training results are applied to the examples using the continuation method and compared with the numerical results of the method of moments (MoM). The proposed 2W-PE method has high computational accuracy and efficiency, which reflects the applicability of machine learning in solving the computational efficiency problem of radio wave propagation. Therefore, this study provides a very effective and reliable method for solving the spatial field in the obstacle environment.
1. Introduction
In recent years, with the rapid development of information technology, the demand for the detection of communication signals and the reliability of their transmission in both civilian and military areas is increasing. Therefore, it is becoming more and more important to accurately predict the propagation performance of radio waves. Scholars have provided many methods for establishing effective and reliable radio wave propagation models, mainly including the Fresnel integral method [1], GO + UTD method [2], finite difference (FD) method [3], and method of moments (MoM) [4], but these methods are unsuitable for calculating the spatial field in real-world, large-scale, complex environments because of the large amount of calculations. The parabolic equation (PE) method, which is obtained by approximate simplification of the wave equation, has become the most effective model for estimating the propagation characteristics of tropospheric waves due to its high computational efficiency and accuracy, and it has become the subject of extensive research [5, 6].
Ozgun first proposed the traditional two-way parabolic equation (2W-PE) [7–9], which combines the stepped terrain model and impedance boundary conditions to calculate the radio wave propagation problem on complex irregular terrain. Since then, the traditional 2W-PE [10, 11] and its solution method [12–14] have been continuously developed, but traditional 2W-PE ignores the calculation of the spatial field inside the obstacle, which will cause a large error in the case of a low-lossy obstacle. In response to the abovementioned problems, in our previous research [15–17], we proposed a 2W-PE that considers multiple reflections inside the obstacle based on the principle of domain decomposition. Then, we used the split-step Fourier transform (SSFT) and FD methods to solve the 2W-PE with boundary conditions to determine the field values of the areas above and inside the obstacle, respectively. This approach obtained more accurate simulation results than the traditional 2W-PE.
However, in the process of solving the field inside the obstacle, we found that due to the multiple reflections inside the obstacle, it takes a long time to solve the 2W-PE step by step using the FD method. Additionally, the calculation burden inside the obstacle is greater for more complex terrain. Because of the difference of incident waves at the boundaries of obstacles at different heights and distances, the determination of the coefficients of the impedance boundary conditions is a difficult problem, and the boundary conditions will also greatly affect the spatial field accuracy around the obstacle. Therefore, new methods are needed to solve the abovementioned problems.
In recent years, machine learning has been widely used in face recognition, image processing, and other fields, including the field of radio wave propagation computing. For example, artificial neural networks based on the back-propagation algorithm have been successfully applied to calculate radio wave propagation loss in rural [18], suburban [19], urban [20], and campus environments [21] owing to their advantages in accuracy, complexity, and prediction time. Currently, though, these neural networks have limited generalization ability to complex environments with different topographies. Since the 2W-PE has provided a calculation framework with stable performance and strong generalization ability, we consider using neural networks to calculate the influence of the internal field of obstacles on the surrounding spatial field and determine boundary condition coefficients, which can improve the calculation speed and accuracy of the entire spatial field.
Therefore, in this study, we propose an efficient and accurate 2W-PE method based on complex-valued neural networks (CVNNs) and a physics-informed neural network (PINN), which is used to quickly and accurately solve the spatial field value. The rest of this article is organized as follows: in Section 2, the theoretical framework of the 2W-PE method is described and the computational ideas are analyzed. In Section 3, the structure and continuation method of the neural network are introduced. In Section 4, the training process and experimental results are discussed, and in Section 5, the conclusions of this study are presented.
2. Framework of the Problem
The use of the 2W-PE method to predict the spatial field in a scene with an obstacle is shown in Figure 1. The time specification exp (jωt) is adopted in this study, where ω is the angular frequency. Field satisfies the two-dimensional (2D) wave equation.where is the wave number in free space, meaning it is also the refractive index in the medium. In the case of different polarization, is defined as
We suppose that the electromagnetic waves propagate axially and the positive axis direction is . After factoring and demodulating the wave equation, we obtain 2W-PE corresponding to forward propagation and backward propagation as follows:
In the calculation, we choose the 2D current density distribution [16] to represent the initial field in vertical polarization and the 2D magnetic current density distribution to represent the initial field in horizontal polarization. To solve the 2W-PE in the area with the rectangular obstacle to obtain the spatial field, we introduce the domain decomposition method in [15] to discuss the field and the propagation process inside the obstacle.
2.1. The Algorithm Framework Based on the Domain Decomposition Method
In the case of a homogeneous medium ground with a rectangular obstacle as shown in Figure 2, the domain decomposition method specifically refers to decomposing the obstacle area for calculating the propagation of radio waves into the areas above and inside the obstacle, which are, respectively, denoted as . The flat terrain area in front of and behind the obstacle object can be recorded and is denoted as . After decomposing the computational area, we solve the 2W-PE separately in these subdomains in combination with the boundary conditions.
It can be seen from [15] that, for the semi-infinite areas and , after domain decomposition, the upper boundary is replaced by the absorption boundary in the calculation, and the lower boundary is the medium ground, which can be replaced by the impedance boundary. Therefore, when solving the 2W-PE in and step by step, Q in equation (3) is expanded according to the Feit–Fleck model, and with the given impedance boundary conditions, discrete mixed Fourier transform (DMFT) [15] is used to solve the field value in the region. The lower boundary of is a homogeneous medium ground, and its upper boundary is connected to the lower boundary of , so its upper and lower boundaries are calculated using impedance boundaries. As the area is a finite area, the FD method is used to calculate the field value inside the obstacle [15]. However, the calculation of the field inside the obstacle increases the computational burden and prolongs the computational time, so we consider using neural networks to optimize the algorithm.
2.2. The Algorithm Framework Based on Machine Learning
According to the calculation idea of using the FD method to solve the 2W-PE inside the obstacle [15], there should be multiple reflection and transmission processes between the front and rear edges of the rectangular obstacle, as shown in Figure 3, and the calculation process is performed until the error between propagation processes converges, which usually takes a long time. The reason we analyze the propagation process inside the obstacle is mainly to consider the influence of the internal field value of the obstacle on the field value of the space surrounding the obstacle. To speed up the calculation, we do not consider the internal field value of the obstacle in practical application but completely separate the calculation in the outer area of the obstacle and that in the inner area and replace the external influence of the field values of the internal area by setting the appropriate boundary conditions.
As shown in Figure 3, according to the derivation of equation (3.30) in [15], we can obtain the relationship between the field and the field at after one-step forward propagation. Propagation at each step is a nonlinear calculation process. Thus, assuming that the initial field at the starting position of the obstacle is , the total backward propagation field penetrating from the end position of the obstacle after multiple propagations is , and the total forward propagation field penetrating from the starting position of the obstacle is , then we obtainwhere and are weight matrices with boundary condition coefficients as variables, and are the corresponding bias matrices, and and are the nonlinear transformations existing in multiple transmission reflections. In this way, we can replace the complex nonlinear process of multiple propagations and superpositions with neural networks. The input of each equation in (4) is the initial field , and the outputs are the fields , which are transmitted along the front and rear edges, respectively, after multiple forward and reverse propagation processes. After training, the multilayered calculation process of the spatial field value in the case of a single obstacle can be simplified to the lumped reflection of the front edge of the obstacle and the lumped transmission of the rear edge, as shown in Figure 3.
3. Calculation Process with Neural Networks
The objective of this study is to develop reasonable neural networks to greatly speed up the calculation process and improve the accuracy of 2W-PE. In addition, the neural networks are designed to have good continuation ability, meaning the results can be extended to rectangular obstacles of any width and even undulating terrain.
3.1. The Choice of Training Basis
First, since the size of the obstacles is uncertain, to make the training results applicable to an arbitrary situation of terrain with a rectangular obstacle, we choose the training basis shown in Figure 4(a), which is taken from the wide rectangle in Figure 4(b), and the training results of which can be extended to any width obstacle by the continuation method proposed next. The width of each basis is two steps of the +x direction in the calculation of 2W-PE. For the training basis in Figure 4(a), we set the initial field of each training basis as the input of the corresponding neural network, and the outputs are the lumped reflection field at the start of an obstacle, the intermediate field at the midpoint, the lumped transmission field at the end, and the coefficients of the boundary conditions , as shown in Figure 4(a).
(a)
(b)
3.2. The Continuation Method
The previous section discussed the training situation of a single training basis, which is summarized in Figure 4(a). The analysis process of extending this result to obstacles of arbitrary width is summarized in Figure 4(b). Three field components at different positions of the training basis are obtained according to the network training. As shown in Figure 4(a), considering the multiple reflections in the obstacle, we introduce several variables to help us solve for the lumped field values: the step transfer function , the reflection coefficient of the front edge , and the proportion B of the backward wave in . As shown in Figure 3, includes the reflection field and several transmission fields at the front edge during multiple propagations. Thus, can be expressed aswhere represents the transmission coefficient when waves propagate from the obstacle to the free space.
Similarly, includes several transmission fields at the back edge and can be expressed as
According to the multilayered reflections, the forward and backward waves in can be expressed aswhere represents the transmission coefficient when waves propagate from the free space to the obstacle.
Dividing the two equations in (7), we can obtain
Substituting (8) into (5) and (6) and simplifying, we obtain
Then, we convert (9) into equation form as
Thus, we can solve (10), choose a suitable root from the solution to obtain the value of , and then solve the other parameters R and S.
After obtaining the key parameters in this way, we can carry out the continuation into rectangular obstacles of any width. As shown in Figure 3, assuming that the width of the rectangular obstacle is N grids, similar to the abovementioned process, we extract the intermediate field value of the first training base , distinguish the forward and backward waves by the ratio , use the transfer function to recursively extend to the back edge, and then reversely extend to the front edge. Finally, the lumped reflection field and the lumped transmission field can be expressed as
Therefore, as long as we accurately obtain the corresponding for the incident field given by the training basis, we can obtain the important parameters in the extension process according to (10) and then obtain the total incident field and the total reflection field after the extension. Meanwhile, we can retain the boundary conditions corresponding to each training base in Figure 4(b) and calculate the spatial field.
Now, what we need to do is to set up reasonable neural networks to obtain the parameters in Figure 4(a) for any given incident field .
3.3. The Design of the Network Structure and the Choice of the Cost Function
The above mentioned analysis shows that the input of each neural network is the incident field, and the outputs of the networks are the fields at the three distance points and the boundary condition coefficients. Therefore, we consider using parallel neural networks with similar structures for calculation. Since the purpose of the first three networks is to fit the output field and the input and output arrays are of the same size, we apply the fully connected CVNN with one hidden layer. Since the output of the fourth network is the boundary condition coefficients and there are no corresponding label values to verify the accuracy of the boundary condition coefficients, we consider adding the physical equation as a constraint to the neural network, so our network becomes a PINN. Based on the neural networks, physical constraints are added to the loss function by adding operators such as partial derivatives representing physical information. The overall structure of our designed networks is shown in Figure 5.
We compare the output fields of the neural networks with the label data obtained by MoM, and we use the mean square error as a cost value. For the first three neural networks, the outputs are , the cost function is the mean square error between them, and their corresponding label data are .
By feeding this cost value back to the network for training, we find that the field results obtained by training at the three distance points are close to the numerical results obtained by MoM through iterations.
For the fourth neural network PINN, whose output is the boundary condition coefficients , we use non-end-to-end training. The outputs are imported into the wave equation, and the difference between the label values and iteration results of the wave equation represented by FD is taken as the cost value, which is returned to the network for training. The specific PINN structure is shown in Figure 6.
We derive the setting of the cost function as follows: in Figure 7, red marks the fields on the boundaries, green marks the fields at the midpoints of the obstacles, and blue marks the fields at the edges of the obstacle. The meshing is shown in the enlarged circle in Figure 7. We denote the vertical grid points as , the horizontal grid points as , and the field at the position as .
Therefore, by substituting the difference expressions of the second-order partial differential, that is,into (1), we can obtain the wave equation expression at 1∼N−1 grid points inside the obstacle except at the upper and lower boundaries as follows:
The upper and lower impedance boundary conditions arewhere are the impedance boundary condition coefficients, which are the parameters of the output of the neural network. Thus, these boundary conditions are treated by FD and substituted into the wave equation, and then we obtain
|
In this way, a complete wave equation in difference format is developed, and the field value at the intermediate position is obtained by the iterative solution as shown in Algorithm 1. Then, the cost value calculated by using the physical equation is returned to the PINN for the training of weights and biases, so the result is close to the original physical equation law under the constraint of physical information.
3.4. The Complex-Valued Back-Propagation Algorithm
The artificial neural networks consist of a large number of connected artificial neurons, where each connection is assigned a weight, and then the weighted sum is passed on to an activation function. In the solution process, the input and output are both field values, and the spatial field values containing phase information are complex numbers, so the weights, biases, and neuron parameters of our network should be complex numbers, and the entire process is calculated through a complex-valued back-propagation algorithm.
3.4.1. Forward-Propagation Process
In the fully connected network layer, the output of each layer can be written as follows:where k represents the order number of neurons in the layer, is the output of the jth neuron in the (l+1)th layer, and is the weight of the kth neuron in the lth layer to the jth neuron in the (l+1)th layer. Furthermore, the superscripts represent the order numbers of the network layer, the subscripts represent neuron positions, and represents the nonlinear activation function, where the tanh function is adopted and is defined asthe derivative of which is
3.4.2. Back-Propagation Process
We want to find the minimum value of the loss function through training and iteratively adjusting the weights and biases. First, the input change of the jth neuron in the lth layer is defined as the change of the cost function and denoted as the error function .
The result of the derivation of the weight by the cost function C is a complex number. The real part of the complex number is the result of the derivation of the real part of the weight by the function C, and the imaginary part of the complex number is the result of the derivation of the imaginary part of the weight. The process of postadjusting weights and biases in the (t+1)th iteration of training iswherewhere denotes the learning rate and denotes the complex conjugate of . The abovementioned formula shows that the increment of the updated weight is the inverse of the product of the learning rate , the error function , and the conjugate of the output value of the neuron in the previous layer .
The iterative calculation formulas of from the back to the front are given, the initial value in the reverse calculation process is the output, and is the error of the output layer. The superscript L refers to the Lth layer, and the subscript k refers to the kth neuron. Thus, can be expressed as
Similarly, we can derive
Then, the complex error of the (L−1)th layer is
Thus, we can obtain
Therefore,
From (23) and (27), we can obtain the recursive relationship of the error of the front and rear layers, and then we can derive the error expression of any layer. By substituting these errors to update the weights and biases, we can complete the derivation process of the back-propagation algorithm in a CVNN.
4. Results and Discussion
According to our designed neural network, we give details on the training process and verify the training results.
4.1. The Training Process
As shown in Figure 5, the CVNNs created in MATLAB adopt a two-layer feed-forward neural network architecture, consisting of an input layer, a hidden layer with a tanh activation function, and an output layer. To balance the convergence rate and the computational complexity, we choose the number of hidden neurons to be 10. We train the network following the back-propagation flow, and we evaluate the prediction performance by the cost function given in Section 3.3, with a maximum number of iterations of 1000.
We use MoM to generate a data set of containing 3000 data points with different positions of the obstacle (600–900 m) and different heights of the transmitting antenna (20–120 m). Assuming that the upper half of the space is free space, the other parameters remain unchanged, namely, the ground is a homogeneous medium, the conductivity of which is 0.00001 S/m, the magnetic permeability is H/m, and the relative permittivity is 6. The obstacle as training basis is 40 m high and 2 m wide, with a magnetic permeability of H/m, a relative permittivity of 4, and an electrical conductivity of 0.0001 S/m.
Considering the training time and complexity, we split the 3000 data samples generated by MoM into two subdatasets, with 75% of the samples being used for training and 25% for validation. The training is conducted according to the complex-valued back-propagation algorithm in Section 3.4, and the network structure is trained and optimized using the stochastic gradient descent algorithm. The training process stops when the cost value is observed to no longer drop within 5 epochs or the total number of training iterations reaches 1000. The cost value changes on the training set and validation set of networks for during the training process are presented in Figure 8.
The total training duration of the network corresponding to is 74.3 s, 77.2 s, and 73.6 s, and the training duration of the PINN corresponding to the boundary condition coefficient is 243 s. With an existing data set, the training of obstacles with a single height can be conducted very quickly and with high training efficiency.
The relative mean square errors of the training set and the validation set at the end of training are shown in Table 1.
4.2. The Simulation Results
We apply the trained neural networks to different examples to test the accuracy of our proposed model and its generalization ability. Two examples are discussed as follows.
Example 1. The environment of a single rectangular obstacle and the medium parameters of the obstacle are shown in Figure 9. We assume that the upper half space is the free space, the ground is a homogeneous medium, the electrical conductivity is 0.00001 S/m, the magnetic permeability is H/m, and the relative permittivity is 6. Furthermore, the magnetic permeability of the medium obstacle is H/m, the relative permittivity is 4, the conductivity is 0.0001 S/m, and the width is 22 m. Additionally, the operating frequency of the emission source is 30 MHz, and the equivalent radiation power of the radiation source is W when is the free space wave impedance. Therefore, we change the values of the characteristic parameters and compare them with the results of MoM and the 2W-PE in [15]. The distance from the antenna to the obstacle is 700 m, the antenna height is 25 m, and the transmission frequency is 30 MHz. In this single rectangular obstacle environment, according to (27), the trained network performs 21 continuation calculations in the forward direction from the second grid of the obstacle and 22 times in the reverse direction after one reflection. Finally, the lumped reflection and transmission fields of obstacles and boundary condition coefficients are obtained. The simulation results of the three methods at the observation height of 25 m are shown in Figure 10.
From Figure 10, it can be seen that an electromagnetic wave interferes in front of the obstacle, and the amplitude of the field value oscillation is large. As the field value obtained through the obstacle transmission and diffraction is small behind the obstacle and there is no obstacle reflection, the amplitude of the field value changes very little. Our new 2W-PE method and the original 2W-PE are in good agreement before and after the obstacle, and there is only a small error in the area close to the obstacle. The enlarged image shows that the MoM and the new 2W-PE are in good agreement in both phase and amplitude. In terms of computing time, MoM takes 670s, the original 2W-PE takes 10.9s, and the new 2W-PE takes only 3.2s, which greatly improves the computing efficiency.
Example 2. We test the generalization ability of the training results in the double-rectangular obstacle environment shown in Figure 11. The two obstacles are both 22 m wide and 40m high, located at 700 m and 900 m, respectively, and other parameters are the same as in Example 1. The computing process in this environment is performed under the original 2W-PE framework [15]. When an electromagnetic wave propagates in space, the multilayered reflection process of the electromagnetic wave between the front and back edges of obstacles and between obstacles is considered. The neural networks are only used in the calculation within the single obstacle, and the field generated by multiple reflections between two obstacles is obtained by solving the 2W-PE.
The simulation results obtained from the observation height of 39 m are shown in Figure 12. The results in Figure 12 show that the new 2W-PE and the original 2W-PE are only in good agreement before the obstacles, and there are relatively large errors between the obstacles and after the obstacles. According to [17], this is caused by the inaccurate boundary conditions during the calculation of the internal field value of the obstacle. However, from Figure 12, we can see that our new 2W-PE and MoM always match better whether the observation point is close to the source or between two obstacles.
Figure 13 presents the results of the spatial field distribution obtained by the three methods. The original 2W-PE has a large error in predicting the field above the obstacles due to inaccurate boundary conditions, but the results of new 2W-PE and MoM are very close in the whole space, indicating our method can generalize to environments with more rectangular obstacles. The time required to calculate the entire spatial field is 4657.5 s for MoM, 478.2 s for original 2W-PE takes, and only 6.8 s for new 2W-PE.
(a)
(b)
(c)
5. Conclusions
In this study, we propose a new 2W-PE method, which applies CVNNs and a PINN to obtain the transmission field and boundary conditions of the obstacle after multiple propagations inside the obstacle so that the spatial field distribution can be efficiently and accurately obtained without calculating the propagation paths inside the obstacle. Through simulations, we prove that the proposed method has several times higher computational efficiency and accuracy to predict the spatial field value than the original 2W-PE method. Moreover, we extend the training results in the case of a single rectangular obstacle to the case of two rectangular obstacles, and the method has higher accuracy and faster calculation speed in the more complex double rectangular obstacle terrain. Finally, we optimize the structure of our neural network and demonstrate the potential to extend the training results to arbitrarily complex undulating terrain. Thus, we have provided a new and reliable method for spatial field prediction.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.