An Efficient and Accurate 2W-PE Method for Solving Spatial Field Strength Based on CVNN and PINN

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.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 $\mathrm{\Psi }$ satisfies the two-dimensional (2D) wave equation.$$\begin{array}{}\text{(1)}& \frac{{\partial }^2\Psi }{\partial x^2}+\frac{{\partial }^2\Psi }{\partial z^2}+k^2{n}^2\Psi =0\text{,}\end{array}$$where $k$ is the wave number in free space, meaning it is also the refractive index in the medium. In the case of different polarization, $\mathrm{\Psi }$ is defined as$$\begin{array}{}\text{(2)}& \Psi x,z=\begin{array}{cc}{H}_{y}x,z,& \text{for}\hspace{1em}\mathrm{V}\text{Pol},\\ E_{y}x,z,& \text{for}\hspace{1em}\mathrm{H}\text{Pol}.\end{array}\end{array}$$

[figure(s) omitted; refer to PDF]

We suppose that the electromagnetic waves propagate axially and the positive axis direction is $+x$. After factoring and demodulating the wave equation, we obtain 2W-PE corresponding to forward propagation and backward propagation as follows:$$\begin{array}{}\text{(3)}& \begin{array}{c}\frac{\partial u_f}{\partial x}=jk1-Q{u}_f\text{,}\\ \frac{\partial u_b}{\partial x}=jkQ-1u_b\text{,}\\ Q=\sqrt{\frac{1}{k^2}\frac{{\partial }^2}{\partial z^2}+n^{2}x,z}\text{.}\end{array}\end{array}$$

In the calculation, we choose the 2D current density distribution [16] $J0,z$ to represent the initial field in vertical polarization and the 2D magnetic current density distribution $M0,z$ 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 ${\mathrm{\Omega }}_1\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{\Omega }}_2$. The flat terrain area in front of and behind the obstacle object can be recorded and is denoted as ${\mathrm{\Omega }}_3$. After decomposing the computational area, we solve the 2W-PE separately in these subdomains in combination with the boundary conditions.

[figure(s) omitted; refer to PDF]

It can be seen from [15] that, for the semi-infinite areas ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_3$, 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 ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_3$ 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 ${\mathrm{\Omega }}_2$ is a homogeneous medium ground, and its upper boundary is connected to the lower boundary of ${\mathrm{\Omega }}_1$, so its upper and lower boundaries are calculated using impedance boundaries. As the area ${\mathrm{\Omega }}_2$ 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.

[figure(s) omitted; refer to PDF]

As shown in Figure 3, according to the derivation of equation (3.30) in [15], we can obtain the relationship between the field ${\mathrm{\Psi }}_n$ and the field at ${\mathrm{\Psi }}_{n+1}$ 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 ${\mathrm{\Psi }}_0$, the total backward propagation field penetrating from the end position of the obstacle after multiple propagations is ${\mathrm{\Psi }}_3$, and the total forward propagation field penetrating from the starting position of the obstacle is ${\mathrm{\Psi }}_1$, then we obtain$$\begin{array}{}\text{(4)}& \begin{array}{c}{\Psi }_1=\vartheta {\mathbf{w}}_1{\alpha }_n·{\Psi }_0+{\mathbf{b}}_1\text{,}\\ {\Psi }_3=\mathfrak{I}{\mathbf{w}}_2{\alpha }_n·{\Psi }_0+{\mathbf{b}}_2\text{,}\end{array}\end{array}$$where${\mathbf{w}}_1{\alpha }_n$ and ${\mathbf{w}}_2{\alpha }_n$ are weight matrices with boundary condition coefficients ${\alpha }_n$ as variables, ${\mathbf{b}}_1$ and ${\mathbf{b}}_2$ are the corresponding bias matrices, and $\mathfrak{I}\cdot$ and $\vartheta \cdot$ 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 ${\mathrm{\Psi }}_0$, and the outputs are the fields ${\mathrm{\Psi }}_1\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{\Psi }}_3$, 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 ${\mathrm{\Psi }}_0$ of each training basis as the input of the corresponding neural network, and the outputs are the lumped reflection field ${\mathrm{\Psi }}_1$ at the start of an obstacle, the intermediate field ${\mathrm{\Psi }}_2$ at the midpoint, the lumped transmission field ${\mathrm{\Psi }}_3$ at the end, and the coefficients of the boundary conditions ${\alpha }_1\mathrm{a}\mathrm{n}\mathrm{d}{\alpha }_2$, as shown in Figure 4(a).

[figure(s) omitted; refer to PDF]

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 ${\mathrm{\Psi }}_1,{\mathrm{\Psi }}_2,\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{\Psi }}_3$ 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 $S$, the reflection coefficient of the front edge $R$, and the proportion B of the backward wave in ${\mathrm{\Psi }}_2$. As shown in Figure 3, ${\mathrm{\Psi }}_1$ includes the reflection field ${\mathrm{\Psi }}_0\cdot R$ and several transmission fields at the front edge during multiple propagations. Thus, ${\mathrm{\Psi }}_1$ can be expressed as$$\begin{array}{}\text{(5)}& {\Psi }_1={\Psi }_0\cdot R+{\Psi }_2\cdot B\cdot S\cdot 1-R\text{,}\end{array}$$where $1-R$ represents the transmission coefficient when waves propagate from the obstacle to the free space.

Similarly, ${\mathrm{\Psi }}_3$ includes several transmission fields at the back edge and can be expressed as$$\begin{array}{}\text{(6)}& {\Psi }_3={\Psi }_2\cdot 1-B\cdot S\cdot 1-R\text{.}\end{array}$$

According to the multilayered reflections, the forward and backward waves in ${\mathrm{\Psi }}_2$ can be expressed as$$\begin{array}{}\text{(7)}& \begin{array}{c}1-B{\Psi }_2={\Psi }_{0}1+R\cdot S1+S^4{R}^2+S^8{R}^4+\cdots \text{,}\\ B{\Psi }_2={\Psi }_{0}1+R\cdot S^3-R1+S^4{R}^2+S^8{R}^4+\cdots \text{,}\end{array}\end{array}$$where $1+R$ represents the transmission coefficient when waves propagate from the free space to the obstacle.

Dividing the two equations in (7), we can obtain$$\begin{array}{}\text{(8)}& \frac{1-B}{B}=-\frac{1}{S^{2}R}\text{.}\end{array}$$

Substituting (8) into (5) and (6) and simplifying, we obtain$$\begin{array}{}\text{(9)}& \begin{array}{c}{a}_{1}1-B-a_2{\Psi }_3^2=-{\Psi }_2^{2}B\cdot {1-a_{1}1-B+a_2}^2\text{,}\\ a_1=\frac{{\Psi }_1+{\Psi }_3}{{\Psi }_0},a_2=\frac{{\Psi }_3}{{\Psi }_0}\text{.}\end{array}\end{array}$$

Then, we convert (9) into equation form as$$\begin{array}{}\text{(10)}& \begin{array}{c}{c}_1{A}^3+c_2{A}^2+c_{3}A+c_4=0\text{,}\\ A=1-B\text{,}\\ c_1={\Psi }_2^2{1-a_1}^2\text{,}\\ c_2=2{\Psi }_2^{2}1-a_1{a}_2-{\Psi }_2^2{1-a_1}^2\text{,}\\ c_3={\Psi }_2^2{a}_2^2-2{\Psi }_2^{2}1-a_1{a}_2-a_1{\Psi }_3^2\text{,}\\ c_4=-{\Psi }_2^2{a}_2^2+a_2{\Psi }_3^2\text{.}\end{array}\end{array}$$

Thus, we can solve (10), choose a suitable root from the solution to obtain the value of $B$, 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 ${\mathrm{\Psi }}_2$, distinguish the forward and backward waves by the ratio $B$, use the transfer function $S$ to recursively extend to the back edge, and then reversely extend to the front edge. Finally, the lumped reflection field ${\mathrm{\Psi }}_1$ and the lumped transmission field ${\mathrm{\Psi }}_3$ can be expressed as$$\begin{array}{}\text{(11)}& \begin{array}{c}{\Psi }_{2}1-B\cdot S^{N-1}\cdot 1-R={\Psi }_3\text{,}\\ {\Psi }_{2}1-B\cdot S^{2N-1}\cdot -R\cdot 1-R+{\Psi }_0\cdot R={\Psi }_1\text{.}\end{array}\end{array}$$

Therefore, as long as we accurately obtain the corresponding ${\mathrm{\Psi }}_1,{\mathrm{\Psi }}_2,\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{\Psi }}_3$ for the incident field ${\mathrm{\Psi }}_0$ 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 ${\mathrm{\Psi }}_0$.

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.

[figure(s) omitted; refer to PDF]

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 ${\mathrm{\Psi }}_1,{\mathrm{\Psi }}_2,\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{\Psi }}_3$, the cost function is the mean square error between them, and their corresponding label data are ${\mathrm{\Psi }}_1^{\ast },{\mathrm{\Psi }}_2^{\ast },{\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{\Psi }}_2^{\ast }$.$$\begin{array}{}\text{(12)}& MS{E}_{\psi }=\frac{1}{N}{\displaystyle \sum }_{m=1}^N{{\Psi }_{\mathbf{m}}-{\Psi }_{\mathbf{m}}^{\ast }}^2\text{.}\end{array}$$

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 ${\alpha }_1\mathrm{a}\mathrm{n}\mathrm{d}{\alpha }_2$, 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.

[figure(s) omitted; refer to PDF]

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 $z_j=j∆zj=\mathrm{0,1},\dots ,N$, the horizontal grid points as $x_m=m∆xm=\mathrm{0,1},\dots ,M$, and the field $\mathrm{\Psi }{x}_m,z_j$ at the position $x_m,z_j$ as $x_m=m∆xm=\mathrm{0,1},\dots ,M$.

[figure(s) omitted; refer to PDF]

Therefore, by substituting the difference expressions of the second-order partial differential, that is,$$\begin{array}{}\text{(13)}& \begin{array}{c}\frac{{\partial }^2{\Psi }_j^m}{\partial z^2}\approx \frac{1}{∆z^2}{\Psi }_{j+1}^m+{\Psi }_{j-1}^m-2{\Psi }_j^m\text{,}\\ \frac{{\partial }^2{\Psi }_j^m}{\partial x^2}\approx \frac{1}{∆x^2}{\Psi }_j^{m+1}+{\Psi }_j^{m-1}-2{\Psi }_j^m\text{,}\end{array}\end{array}$$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:$$\begin{array}{}\text{(14)}& \frac{{\Psi }_j^m}{∆x^2}+-\frac{2}{∆z^2}-\frac{2}{∆x^2}+k^2{n}^2{\Psi }_j^m+\frac{{\Psi }_{j+1}^m+{\Psi }_{j-1}^m}{∆z^2}+\frac{{\Psi }_j^{m-1}}{∆x^2}=0\text{.}\end{array}$$

The upper and lower impedance boundary conditions are$$\begin{array}{}\text{(15)}& \begin{array}{c}\frac{\partial {\Psi }_N^m}{\partial z}+{\alpha }_1\cdot {\Psi }_N^m=0\text{,}\\ \frac{\partial {\Psi }_0^m}{\partial z}+{\alpha }_2\cdot {\Psi }_0^m=0\text{,}\end{array}\end{array}$$where ${\alpha }_1\mathrm{a}\mathrm{n}\mathrm{d}{\alpha }_2$ 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$$\begin{array}{}\text{(16)}& \begin{array}{c}\frac{{\Psi }_N^{m+1}}{∆x^2}+-\frac{2+2{\alpha }_1∆z}{∆z^2}+k^2{n}^2{\Psi }_N^m+2\cdot \frac{{\Psi }_{N-1}^m}{∆z^2}+\frac{{\Psi }_N^{m-1}}{∆x^2}=0\text{,}\\ \frac{{\Psi }_0^{m+1}}{∆x^2}+-\frac{2+2{\alpha }_2∆z}{∆z^2}+k^2{n}^2{\Psi }_0^m+2\cdot \frac{{\Psi }_1^m}{∆z^2}+\frac{{\Psi }_0^{m-1}}{∆x^2}=0\text{.}\end{array}\end{array}$$

Algorithm 1: The calculation method of the cost function value of the PINN.

  Input: $∆x,∆$z: accuracy of meshing in the x-direction/z-direction

      ${\mathrm{\Psi }}_1$: total field at the start of obstacle calculated by MoM

     ${\mathrm{\Psi }}_{2/∆x+1}$: total field at the end of obstacle calculated by MoM

     $\mathrm{t}\mathrm{h}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{l}\mathrm{d}$: threshold of iterative convergence

      $N$: maximum number of iterations

      ${\alpha }_1,{\alpha }_2$: coefficients of upper/lower impedance boundary condition

  Output: $\mathrm{C}\mathrm{o}\mathrm{s}\mathrm{t}$: the cost value for back propagation

  Initialize: ${\mathrm{\Psi }}_{2:2/∆x}=\mathbf{0}$;      //initial fields to be solved

      $n=1$;       //initial number of iterations

      $\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}=1000$;      //initial error of iterations

      ${\mathrm{\Psi }}_{1/∆x+1}^{\ast }={\mathrm{\Psi }}_{1/∆x+1}$      //intermediate variable

  While $\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}>\mathrm{t}\mathrm{h}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{l}\mathrm{d}$ or $n<N$ do

   for$j=2:\mathrm{l}\mathrm{e}\mathrm{n}\mathrm{g}\mathrm{t}\mathrm{h}{\mathrm{\Psi }}_1-1$do

    for$m=2:2/∆x$do

     ${\mathrm{\Psi }}_0^m={\mathrm{\Psi }}_0^{m+1}/{∆x}^2+{2\mathrm{\Psi }}_1^m/{∆z}^2+{\mathrm{\Psi }}_0^{m-1}/{∆x}^2/2-2{\alpha }_1∆z/{∆z}^2-k^2{n}^2$

     ${\mathrm{\Psi }}_0^m={\mathrm{\Psi }}_0^{m+1}/{∆x}^2+{2\mathrm{\Psi }}_1^m/{∆z}^2+{\mathrm{\Psi }}_0^{m-1}/{∆x}^2/2-2{\alpha }_2∆z/{∆z}^2-k^2{n}^2$

     ${\mathrm{\Psi }}_j^m={\mathrm{\Psi }}_j^{m+1}/{∆x}^2+{\mathrm{\Psi }}_{j+1}^m+{\mathrm{\Psi }}_{j-1}^m/{∆z}^2+{\mathrm{\Psi }}_j^{m-1}/{∆x}^2/2/{∆z}^2+2/{∆x}^2-k^2{n}^2$

     $\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}=\mathrm{M}\mathrm{S}\mathrm{E}{\mathrm{\Psi }}_{1/∆x+1}^{\ast }-{\mathrm{\Psi }}_{1/∆x+1}$

     ${\mathrm{\Psi }}_{1/∆x+1}^{\ast }={\mathrm{\Psi }}_{1/∆x+1}$

    end

   end

  end

  $\mathrm{C}\mathrm{o}\mathrm{s}\mathrm{t}={\mathrm{\Psi }}_0^m-{\mathrm{\Psi }}_0^{m+1}/{∆x}^2+{2\mathrm{\Psi }}_1^m/{∆z}^2+{\mathrm{\Psi }}_0^{m-1}/{∆x}^2/2-2{\alpha }_1∆z/{∆z}^2-k^2{n}^2$

  $+{\mathrm{\Psi }}_0^m-{\mathrm{\Psi }}_0^{m+1}/{∆x}^2+{2\mathrm{\Psi }}_1^m/{∆z}^2+{\mathrm{\Psi }}_0^{m-1}/{∆x}^2/2-2{\alpha }_2∆z/{∆z}^2-k^2{n}^2$

  $+{{\displaystyle \sum }}_{j=1}^{N-1}{\mathrm{\Psi }}_j^m={\mathrm{\Psi }}_j^{m+1}/{∆x}^2+{\mathrm{\Psi }}_{j+1}^m+{\mathrm{\Psi }}_{j-1}^m/{∆z}^2+{\mathrm{\Psi }}_j^{m-1}/{∆x}^2/2/{∆z}^2+2/{∆x}^2-k^2{n}^2$

  return $\mathrm{C}\mathrm{o}\mathrm{s}\mathrm{t}$

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:$$\begin{array}{c}\text{(17)}& a_j^{l+1}=\sigma z_j^{l+1}\text{,}{z}_k^{l+1}={\displaystyle \sum }_j{w}_{kj}^{l+1}{a}_j^l+b_k^{l+1}={\displaystyle \sum }_j{w}_{kj}^{l+1}\sigma z_j^l+b_k^{l+1}\text{,}\end{array}$$where k represents the order number of neurons in the layer, $a_j^{l+1}$ is the output of the j^th neuron in the (l+1)^th layer, and $w_{jk}^{l+1}$ is the weight of the k^th neuron in the l^th layer to the j^th neuron in the (l+1)^th layer. Furthermore, the superscripts represent the order numbers of the network layer, the subscripts represent neuron positions, and $\sigma \cdot$ represents the nonlinear activation function, where the tanh function is adopted and is defined as$$\begin{array}{}\text{(18)}& \sigma z=\frac{\exp z-\exp -z}{\exp z+\exp -z}\triangleq \sigma z_R+i\sigma z_I\text{,}\end{array}$$the derivative of which is$$\begin{array}{}\text{(19)}& \frac{\partial \sigma z}{\partial z^{\ast }}=\frac{\partial \sigma z_R+i\sigma z_I}{\partial z_R}+i\frac{\partial \sigma z_R+i\sigma z_I}{\partial z_I}={\sigma }^2{z}_I-{\sigma }^2{z}_R\text{.}\end{array}$$

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 j^th neuron in the l^th layer is defined as the change of the cost function and denoted as the error function ${\delta }_j^l$.$$\begin{array}{}\text{(20)}& {\delta }_j^l={\nabla }_{z_j^i}C=\frac{\partial C}{\partial z_j^{l\ast }}=\frac{\partial C}{\partial z_{j,R}^l}+i\cdot \frac{\partial C}{\partial z_{j,I}^l}\text{.}\end{array}$$

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 is$$\begin{array}{}\text{(21)}& \begin{array}{c}{w}_{jk}^{l+1}t+1=w_{jk}^{l+1}t+∆w_{jk}^{l+1}t=w_{jk.R}^l+∆w_{jk.R}^l+i{w}_{jk,I}^l+∆w_{jk,I}^l\text{,}\\ b_j^{l+1}t+1=b_j^{l+1}t+∆b_j^{l+1}t=b_{jk.R}^l+∆b_{jk.R}^l+i{b}_{jk,I}^l+∆b_{jk,I}^l\text{,}\end{array}\end{array}$$where$$\begin{array}{}\text{(22)}& \begin{array}{c}∆w_{jk}^l=-\eta {\nabla }_{w}C=\frac{\partial C}{\partial {w_{jk}^l}^{\ast }}=\frac{\partial C}{\partial w_{jk,R}^l}+j\frac{\partial C}{\partial w_{jk,I}^l}=\frac{\partial C}{\partial z_j^{l\ast }}\frac{\partial z_j^{l\ast }}{\partial w_{jk}^{l\ast }}=-\eta {\delta }_j^l{a_k^{l-1}}^{\ast }\text{,}\\ ∆b_j^l=-\eta {\nabla }_{b_j^l}C=-\eta \frac{\partial C}{\partial b_j^{l\ast }}=-\eta {\delta }_j^l\text{,}\end{array}\end{array}$$where $\eta$ denotes the learning rate and ${a_k^{l-1}}^{\ast }$ denotes the complex conjugate of $a_k^{l-1}$. The abovementioned formula shows that the increment of the updated weight is the inverse of the product of the learning rate $\eta$, the error function ${\delta }_j^l$, and the conjugate of the output value of the neuron in the previous layer ${a_k^{l-1}}^{\ast }$.

The iterative calculation formulas of ${\delta }^l$ from the back to the front are given, the initial value in the reverse calculation process is the output, and ${\delta }_k^L$ is the error of the output layer. The superscript L refers to the L^th layer, and the subscript k refers to the k^th neuron. Thus, ${\delta }_k^L$ can be expressed as$$\begin{array}{c}\text{(23)}& {\delta }_k^L=\frac{\partial C}{\partial z_k^{L\ast }}=\frac{\partial {\displaystyle \sum }_{k=1}^K{a_k^L-T_k}^2}{\partial z_k^{L\ast }}=\frac{\partial a_k^L-T_k{a_k^L-T_k}^{\ast }}{\partial a_k^{L\ast }}\cdot \frac{\partial a_k^{L\ast }}{\partial z_k^{L\ast }}+\frac{\partial a_k^L-T_k{a_k^L-T_k}^{\ast }}{\partial a_k^L}\cdot \frac{\partial a_k^L}{\partial z_k^{L\ast }}=a_{k,R}^L-T_{k,R}\cdot 1-{\sigma }^2{z}_{k,R}^L+i\cdot a_{k,I}^L-T_{k,I}\cdot 1-{\sigma }^2{z}_{k,I}^L\text{.}\end{array}$$

Similarly, we can derive$$\begin{array}{}\text{(24)}& {\delta }_k^{L\ast }=\frac{\partial C}{\partial z_k^L}=\mathrm{R}\mathrm{e}{a}_k^L-T_k\cdot 1-{\sigma }^2{z}_{k,R}^L-i\cdot \mathrm{I}\mathrm{m}{a}_k^L-T_k\cdot 1-{\sigma }^2{z}_{k,I}^L\text{.}\end{array}$$

Then, the complex error of the (L−1)^th layer is$$\begin{array}{}\text{(25)}& {\delta }_j^{L-1}=\frac{\partial C}{\partial {z_j^{L-1}}^{\ast }}=\frac{\partial C}{\partial a_j^{L-1\ast }}\cdot \frac{\partial a_j^{L-1\ast }}{\partial z_j^{L-1\ast }}+\frac{\partial C}{\partial a_j^{L-1}}\cdot \frac{\partial a_j^{L-1}}{\partial z_j^{L-1\ast }}\text{.}\end{array}$$

Thus, we can obtain$$\begin{array}{}\text{(26)}& \begin{array}{c}\frac{\partial C}{\partial a_j^{L-1\ast }}={\displaystyle \sum }_k\frac{\partial C}{\partial z_k^{L\ast }}\cdot \frac{\partial z_k^{L\ast }}{\partial a_j^{L-1\ast }}+\frac{\partial C}{\partial z_k^L}\cdot \frac{\partial z_k^L}{\partial a_j^{L-1\ast }}={\displaystyle \sum }_k{\delta }_k^L\cdot w_{jk}^{L\ast }\text{,}\\ \frac{\partial C}{\partial a_j^{L-1}}={\displaystyle \sum }_k\frac{\partial C}{\partial z_k^L}\cdot \frac{\partial z_k^L}{\partial a_j^{L-1}}+\frac{\partial C}{\partial z_k^{L\ast }}\cdot \frac{\partial z_k^{L\ast }}{\partial a_j^{L-1}}={\displaystyle \sum }_k{\delta }_k^{L\ast }\cdot w_{jk}^L\text{.}\end{array}\end{array}$$

Therefore,$$\begin{array}{c}\text{(27)}& {\delta }_j^{L-1}=\frac{\partial C}{\partial {z_j^{L-1}}^{\ast }}={\displaystyle \sum }_k{\delta }_k^L\cdot w_{jk}^{L\ast }\cdot \frac{\partial a_j^{L-1\ast }}{\partial z_j^{L-1\ast }}+{\displaystyle \sum }_k{\delta }_k^{L\ast }\cdot w_{jk}^L\cdot \frac{\partial a_j^{L-1}}{\partial z_j^{L-1\ast }}\text{,}=1-{\sigma }^2{z}_{j,R}^{L-1}\cdot \mathrm{R}\mathrm{e}{\displaystyle \sum }_k{\delta }_k^L\cdot w_{jk}^{L\ast }+i\cdot 1-{\sigma }^2{z}_{j,I}^{L-1}\cdot \mathrm{I}\mathrm{m}{\displaystyle \sum }_k{\delta }_k^L\cdot w_{jk}^{L\ast }\text{.}\end{array}$$

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 ${\mathrm{\Psi }}_1,{\mathrm{\Psi }}_2,\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{\Psi }}_3$ 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 $4\pi \times {10}^{-7}$ 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 $4\pi \times {10}^{-7}$ 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 ${\mathrm{\Psi }}_1,{\mathrm{\Psi }}_2,{\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{\Psi }}_3$ during the training process are presented in Figure 8.

[figure(s) omitted; refer to PDF]

The total training duration of the network corresponding to ${\mathrm{\Psi }}_1,{\mathrm{\Psi }}_2,\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{\Psi }}_3$ 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.

Table 1

Training results.

Outputs	MSE in the training set	MSE in the validation set
${\mathrm{\Psi }}_1$	$1.4030e^{-6}$	$4.4306e^{-6}$
${\mathrm{\Psi }}_2$	$2.9402e^{-6}$	$6.5087e^{-5}$
${\mathrm{\Psi }}_3$	$2.0135e^{-6}$	$3.1585e^{-5}$
${\alpha }_1,{\alpha }_2$	$2.5353e^{-4}$	$2.6225e^{-4}$

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 $4\pi \times {10}^{-7}$ H/m, and the relative permittivity is 6. Furthermore, the magnetic permeability of the medium obstacle is $4\pi \times {10}^{-7}$ 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 $Z_0^{-1}$ W when $Z_0$ 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.

[figure(s) omitted; refer to PDF]

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.

[figure(s) omitted; refer to PDF]

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.

References

[1] R. L. Roberts, J. J. Daniels, "Modelling near-field GPR in three dimensions using the FDTD method," Geophysics, vol. 6 no. 24, pp. 1114-1126, 1997.

[2] M. Wang, Geometric Diffraction Theory, 1994.

[3] O. Ozgun, S. Levent, "Comparative study of analytic and numerical techniques in modelling electromagnetic scattering from single and double knife-edge in 2D ground wave propagation," Problemos, vol. 29 no. 5, pp. 376-388, 2012.

[4] Z. Su, J. D. Tan, J. Zhang, "Calculation of electromagnetic wave propagation based on high-order finite-difference time-domain algorithm," Journal of Radio Wave Science, vol. 29 no. 3, pp. 431-436, 2014.

[5] M. A. Leontovich, V. A. Fock, "Solution of the problem of propagation of electromagnetic waves along the Earth’s surface by method of parabolic wave equation," Journal of Physics-USSR, vol. 10, pp. 13-23, 1946.

[6] F. D. Tappert, "The parabolic approximation method," Wave Propagation and Underwater Acoustics, pp. 22-287, 1977.

[7] G. D. Dockery, G. C. Konstanzer, "Recent advances in prediction of tropospheric propagation using the parabolic equation," Johns Hopkins APL Technical Digest, vol. 8 no. 4, pp. 404-412, 1987.

[8] M. Levy, Parabolic Equation Methods for Electromagnetic Wave Propagation, 2000.

[9] O. Ozgun, G. Apaydin, M. Kuzuoglu, "Two-way Fourier split step algorithm over variable terrain with narrow and wide-angle propagators," Proceedings of the 2010 IEEE International Symposium on Antennas and Propagation, .

[10] O. Ozgun, "Two-way propagation modelling in waveguides with three-dimensional finite-element and split-step fourier-based PE approaches," IEEE Trans. Antennas and Wireless Propagation Letters, vol. 10, pp. 975-978, 2011.

[11] O. Ozgun, V. Sahin, M. E. Erguden, G. Apaydin, A. E. Yilmaz, M. Kuzuoglu, L. Sevgi, "PETOOL v2.0: parabolic Equation Toolbox with evaporation duct models and real environment data," Computer Physics Communications, vol. 256,DOI: 10.1016/j.cpc.2020.107454, 2020.

[12] D. D. Wang, X. L. Xi, Y. R. Pu, J. F. Liu, L. L. Zhou, "LF groundwave propagation modelling over irregular terrain by the hybrid two-way parabolic equation method," IEEE Transactions on Antennas and Propagation, vol. 67 no. 10, pp. 90-98, 2019.

[13] Q. Guo, Y. L. Long, "Two-way parabolic equation method for radio propagation over rough sea surface," IEEE Transactions on Antennas and Propagation, vol. 68 no. 6, pp. 4839-4847, DOI: 10.1109/tap.2020.2975288, 2020.

[14] P. D. Holm, "Wide-angle shift-map PE for a piecewise linear terrain: a finite-difference approach," IEEE Transactions on Antennas and Propagation, vol. 55 no. 10, pp. 2773-2789, DOI: 10.1109/tap.2007.905865, 2007.

[15] A. Q. Li, C. Y. Yin, Q. q. Zhang, Y. Gan, "Research on fast algorithm of radio wave propagation in low-lossy obstacles environment," International Journal of Antennas and Propagation, vol. 2021,DOI: 10.1155/2021/5519586, 2021.

[16] A. Q. Li, C. Y. Yin, Y. Gan, W. Wu, "Accurate prediction of radio wave propagation in an environment of dielectric ground and obstacles based on the principle of domain decomposition," IET Microwaves, Antennas & Propagation, vol. 15 no. 11, pp. 1473-1489, DOI: 10.1049/mia2.12130, 2021.

[17] A. Q. Li, C. Y. Yin, Q. Zhang, "Using 2W-PE method based on machine learning to accurately predict field strength distribution in flat-top obstacle environment," AEU - International Journal of Electronics and Communications, vol. 144,DOI: 10.1016/j.aeue.2021.154037, 2022.

[18] E. Ostlin, H. J. Zepernick, H. Suzuki, "Macrocell path-loss prediction using artificial neural networks," IEEE Transactions on Vehicular Technology, vol. 59 no. 6, pp. 2735-2747, DOI: 10.1109/tvt.2010.2050502, 2010.

[19] I. Popescu, D. Nikitopoulos, P. Constantinou, I. Nafornita, "ANN prediction models for outdoor environment," Proceedings of the 2006 IEEE 17th International Symposium on Personal, Indoor and Mobile Radio Communications, .

[20] Z. Stankovic, B. Milovanovic, M. Veljkovic, A. Dordevic, "The hybrid-neural empirical model for the electromagnetic field level prediction in urban environments," Proceedings of the 7th Seminar on Neural Network Applications in Electrical Engineering, pp. 189-192, .

[21] H. Singh, S. Gupta, C. Dhawan, "Path loss prediction in smart campus environment: machine learning-based approaches," Proceedings of the 2020 IEEE 91st Vehicular Technology Conference (VTC2020-Spring), .