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

Authors: An-Qi Li; Cheng-You Yin; Qian-Qian Zhang; Han Liu
Journal: International Journal of Antennas and Propagation
 (2023)

---

## 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.


---

## Body

## 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 Section2, 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 Figure1. The time specification exp (jωt) is adopted in this study, where ω is the angular frequency. Field Ψ satisfies the two-dimensional (2D) wave equation.(1)∂2Ψ∂x2+∂2Ψ∂z2+k2n2Ψ=0,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, Ψ is defined as(2)Ψx,z=Hyx,z,forVPol,Eyx,z,forHPol.Figure 1
Schematic diagram of radio wave propagation in the obstacle area.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:(3)∂uf∂x=jk1−Quf,∂ub∂x=jkQ−1ub,Q=1k2∂2∂z2+n2x,z.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 Figure2, 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 Ω1andΩ2. The flat terrain area in front of and behind the obstacle object can be recorded and is denoted as Ω3. After decomposing the computational area, we solve the 2W-PE separately in these subdomains in combination with the boundary conditions.Figure 2
Schematic diagram of domain decomposition.It can be seen from [15] that, for the semi-infinite areas Ω1 and Ω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 Ω1 and Ω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 Ω2 is a homogeneous medium ground, and its upper boundary is connected to the lower boundary of Ω1, so its upper and lower boundaries are calculated using impedance boundaries. As the area Ω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 3
The multilayered propagation process inside an obstacle.As shown in Figure3, according to the derivation of equation (3.30) in [15], we can obtain the relationship between the field Ψn and the field at Ψ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 Ψ0, the total backward propagation field penetrating from the end position of the obstacle after multiple propagations is Ψ3, and the total forward propagation field penetrating from the starting position of the obstacle is Ψ1, then we obtain(4)Ψ1=ϑw1αn·Ψ0+b1,Ψ3=Iw2αn·Ψ0+b2,wherew1αn and w2αn are weight matrices with boundary condition coefficients αn as variables, b1 and b2 are the corresponding bias matrices, and I⋅ 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 Ψ0, and the outputs are the fields Ψ1andΨ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.


## 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 Figure2, 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 Ω1andΩ2. The flat terrain area in front of and behind the obstacle object can be recorded and is denoted as Ω3. After decomposing the computational area, we solve the 2W-PE separately in these subdomains in combination with the boundary conditions.Figure 2
Schematic diagram of domain decomposition.It can be seen from [15] that, for the semi-infinite areas Ω1 and Ω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 Ω1 and Ω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 Ω2 is a homogeneous medium ground, and its upper boundary is connected to the lower boundary of Ω1, so its upper and lower boundaries are calculated using impedance boundaries. As the area Ω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 3
The multilayered propagation process inside an obstacle.As shown in Figure3, according to the derivation of equation (3.30) in [15], we can obtain the relationship between the field Ψn and the field at Ψ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 Ψ0, the total backward propagation field penetrating from the end position of the obstacle after multiple propagations is Ψ3, and the total forward propagation field penetrating from the starting position of the obstacle is Ψ1, then we obtain(4)Ψ1=ϑw1αn·Ψ0+b1,Ψ3=Iw2αn·Ψ0+b2,wherew1αn and w2αn are weight matrices with boundary condition coefficients αn as variables, b1 and b2 are the corresponding bias matrices, and I⋅ 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 Ψ0, and the outputs are the fields Ψ1andΨ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 Figure4(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 Ψ0 of each training basis as the input of the corresponding neural network, and the outputs are the lumped reflection field Ψ1 at the start of an obstacle, the intermediate field Ψ2 at the midpoint, the lumped transmission field Ψ3 at the end, and the coefficients of the boundary conditions α1andα2, as shown in Figure 4(a).Figure 4
Schematic diagram of training basis and the continuation method.
(a)(b)


### 3.2. The Continuation Method

The previous section discussed the training situation of a single training basis, which is summarized in Figure4(a). The analysis process of extending this result to obstacles of arbitrary width is summarized in Figure 4(b). Three field components Ψ1,Ψ2,andΨ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 Ψ2. As shown in Figure 3, Ψ1 includes the reflection field Ψ0⋅R and several transmission fields at the front edge during multiple propagations. Thus, Ψ1 can be expressed as(5)Ψ1=Ψ0⋅R+Ψ2⋅B⋅S⋅1−R,where 1−R represents the transmission coefficient when waves propagate from the obstacle to the free space.Similarly,Ψ3 includes several transmission fields at the back edge and can be expressed as(6)Ψ3=Ψ2⋅1−B⋅S⋅1−R.According to the multilayered reflections, the forward and backward waves inΨ2 can be expressed as(7)1−BΨ2=Ψ01+R⋅S1+S4R2+S8R4+⋯,BΨ2=Ψ01+R⋅S3−R1+S4R2+S8R4+⋯,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(8)1−BB=−1S2R.Substituting (8) into (5) and (6) and simplifying, we obtain(9)a11−B−a2Ψ32=−Ψ22B⋅1−a11−B+a22,a1=Ψ1+Ψ3Ψ0,a2=Ψ3Ψ0.Then, we convert (9) into equation form as(10)c1A3+c2A2+c3A+c4=0,A=1−B,c1=Ψ221−a12,c2=2Ψ221−a1a2−Ψ221−a12,c3=Ψ22a22−2Ψ221−a1a2−a1Ψ32,c4=−Ψ22a22+a2Ψ32.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 Figure3, 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 Ψ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 Ψ1 and the lumped transmission field Ψ3 can be expressed as(11)Ψ21−B⋅SN−1⋅1−R=Ψ3,Ψ21−B⋅S2N−1⋅−R⋅1−R+Ψ0⋅R=Ψ1.Therefore, as long as we accurately obtain the correspondingΨ1,Ψ2,andΨ3 for the incident field Ψ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 Figure4(a) for any given incident field Ψ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 Figure5.Figure 5
The structure of the neural networks.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Ψ1,Ψ2,andΨ3, the cost function is the mean square error between them, and their corresponding label data are Ψ1∗,Ψ2∗,andΨ2∗.(12)MSEψ=1N∑m=1NΨm−Ψm∗2.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α1andα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 6
Schematic diagram of the PINN structure.We derive the setting of the cost function as follows: in Figure7, 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 zj=j∆zj=0,1,…,N, the horizontal grid points as xm=m∆xm=0,1,…,M, and the field Ψxm,zj at the position xm,zj as xm=m∆xm=0,1,…,M.Figure 7
The meshing of the FD method.Therefore, by substituting the difference expressions of the second-order partial differential, that is,(13)∂2Ψjm∂z2≈1∆z2Ψj+1m+Ψj−1m−2Ψjm,∂2Ψjm∂x2≈1∆x2Ψjm+1+Ψjm−1−2Ψjm,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:(14)Ψjm∆x2+−2∆z2−2∆x2+k2n2Ψjm+Ψj+1m+Ψj−1m∆z2+Ψjm−1∆x2=0.The upper and lower impedance boundary conditions are(15)∂ΨNm∂z+α1⋅ΨNm=0,∂Ψ0m∂z+α2⋅Ψ0m=0,where α1andα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(16)ΨNm+1∆x2+−2+2α1∆z∆z2+k2n2ΨNm+2⋅ΨN−1m∆z2+ΨNm−1∆x2=0,Ψ0m+1∆x2+−2+2α2∆z∆z2+k2n2Ψ0m+2⋅Ψ1m∆z2+Ψ0m−1∆x2=0.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Ψ1: total field at the start of obstacle calculated by MoMΨ2/∆x+1: total field at the end of obstacle calculated by MoMthreshold: threshold of iterative convergenceN: maximum number of iterationsα1,α2: coefficients of upper/lower impedance boundary conditionOutput: Cost: the cost value for back propagationInitialize: Ψ2:2/∆x=0;      //initial fields to be solvedn=1;       //initial number of iterationserror=1000;      //initial error of iterationsΨ1/∆x+1∗=Ψ1/∆x+1      //intermediate variableWhileerror>thresholdorn<Ndoforj=2:lengthΨ1−1doform=2:2/∆xdoΨ0m=Ψ0m+1/∆x2+2Ψ1m/∆z2+Ψ0m−1/∆x2/2−2α1∆z/∆z2−k2n2Ψ0m=Ψ0m+1/∆x2+2Ψ1m/∆z2+Ψ0m−1/∆x2/2−2α2∆z/∆z2−k2n2Ψjm=Ψjm+1/∆x2+Ψj+1m+Ψj−1m/∆z2+Ψjm−1/∆x2/2/∆z2+2/∆x2−k2n2error=MSEΨ1/∆x+1∗−Ψ1/∆x+1Ψ1/∆x+1∗=Ψ1/∆x+1endendendCost=Ψ0m−Ψ0m+1/∆x2+2Ψ1m/∆z2+Ψ0m−1/∆x2/2−2α1∆z/∆z2−k2n2+Ψ0m−Ψ0m+1/∆x2+2Ψ1m/∆z2+Ψ0m−1/∆x2/2−2α2∆z/∆z2−k2n2+∑j=1N−1Ψjm=Ψjm+1/∆x2+Ψj+1m+Ψj−1m/∆z2+Ψjm−1/∆x2/2/∆z2+2/∆x2−k2n2returnCostIn 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 Algorithm1. 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:(17)ajl+1=σzjl+1,zkl+1=∑jwkjl+1ajl+bkl+1=∑jwkjl+1σzjl+bkl+1,where k represents the order number of neurons in the layer, ajl+1 is the output of the jth neuron in the (l+1)th layer, and wjkl+1 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 as(18)σz=expz−exp−zexpz+exp−z≜σzR+iσzI,the derivative of which is(19)∂σz∂z∗=∂σzR+iσzI∂zR+i∂σzR+iσzI∂zI=σ2zI−σ2zR.


#### 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 thejth neuron in the lth layer is defined as the change of the cost function and denoted as the error function δjl.(20)δjl=∇zjiC=∂C∂zjl∗=∂C∂zj,Rl+i⋅∂C∂zj,Il.The result of the derivation of the weight by the cost functionC 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(21)wjkl+1t+1=wjkl+1t+∆wjkl+1t=wjk.Rl+∆wjk.Rl+iwjk,Il+∆wjk,Il,bjl+1t+1=bjl+1t+∆bjl+1t=bjk.Rl+∆bjk.Rl+ibjk,Il+∆bjk,Il,where(22)∆wjkl=−η∇wC=∂C∂wjkl∗=∂C∂wjk,Rl+j∂C∂wjk,Il=∂C∂zjl∗∂zjl∗∂wjkl∗=−ηδjlakl−1∗,∆bjl=−η∇bjlC=−η∂C∂bjl∗=−ηδjl,where η denotes the learning rate and akl−1∗ denotes the complex conjugate of akl−1. The abovementioned formula shows that the increment of the updated weight is the inverse of the product of the learning rate η, the error function δjl, and the conjugate of the output value of the neuron in the previous layer akl−1∗.The iterative calculation formulas ofδl from the back to the front are given, the initial value in the reverse calculation process is the output, and δkL 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, δkL can be expressed as(23)δkL=∂C∂zkL∗=∂∑k=1KakL−Tk2∂zkL∗=∂akL−TkakL−Tk∗∂akL∗⋅∂akL∗∂zkL∗+∂akL−TkakL−Tk∗∂akL⋅∂akL∂zkL∗=ak,RL−Tk,R⋅1−σ2zk,RL+i⋅ak,IL−Tk,I⋅1−σ2zk,IL.Similarly, we can derive(24)δkL∗=∂C∂zkL=ReakL−Tk⋅1−σ2zk,RL−i⋅ImakL−Tk⋅1−σ2zk,IL.Then, the complex error of the (L−1)th layer is(25)δjL−1=∂C∂zjL−1∗=∂C∂ajL−1∗⋅∂ajL−1∗∂zjL−1∗+∂C∂ajL−1⋅∂ajL−1∂zjL−1∗.Thus, we can obtain(26)∂C∂ajL−1∗=∑k∂C∂zkL∗⋅∂zkL∗∂ajL−1∗+∂C∂zkL⋅∂zkL∂ajL−1∗=∑kδkL⋅wjkL∗,∂C∂ajL−1=∑k∂C∂zkL⋅∂zkL∂ajL−1+∂C∂zkL∗⋅∂zkL∗∂ajL−1=∑kδkL∗⋅wjkL.Therefore,(27)δjL−1=∂C∂zjL−1∗=∑kδkL⋅wjkL∗⋅∂ajL−1∗∂zjL−1∗+∑kδkL∗⋅wjkL⋅∂ajL−1∂zjL−1∗,=1−σ2zj,RL−1⋅Re∑kδkL⋅wjkL∗+i⋅1−σ2zj,IL−1⋅Im∑kδkL⋅wjkL∗.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.


## 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 Figure4(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 Ψ0 of each training basis as the input of the corresponding neural network, and the outputs are the lumped reflection field Ψ1 at the start of an obstacle, the intermediate field Ψ2 at the midpoint, the lumped transmission field Ψ3 at the end, and the coefficients of the boundary conditions α1andα2, as shown in Figure 4(a).Figure 4
Schematic diagram of training basis and the continuation method.
(a)(b)


## 3.2. The Continuation Method

The previous section discussed the training situation of a single training basis, which is summarized in Figure4(a). The analysis process of extending this result to obstacles of arbitrary width is summarized in Figure 4(b). Three field components Ψ1,Ψ2,andΨ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 Ψ2. As shown in Figure 3, Ψ1 includes the reflection field Ψ0⋅R and several transmission fields at the front edge during multiple propagations. Thus, Ψ1 can be expressed as(5)Ψ1=Ψ0⋅R+Ψ2⋅B⋅S⋅1−R,where 1−R represents the transmission coefficient when waves propagate from the obstacle to the free space.Similarly,Ψ3 includes several transmission fields at the back edge and can be expressed as(6)Ψ3=Ψ2⋅1−B⋅S⋅1−R.According to the multilayered reflections, the forward and backward waves inΨ2 can be expressed as(7)1−BΨ2=Ψ01+R⋅S1+S4R2+S8R4+⋯,BΨ2=Ψ01+R⋅S3−R1+S4R2+S8R4+⋯,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(8)1−BB=−1S2R.Substituting (8) into (5) and (6) and simplifying, we obtain(9)a11−B−a2Ψ32=−Ψ22B⋅1−a11−B+a22,a1=Ψ1+Ψ3Ψ0,a2=Ψ3Ψ0.Then, we convert (9) into equation form as(10)c1A3+c2A2+c3A+c4=0,A=1−B,c1=Ψ221−a12,c2=2Ψ221−a1a2−Ψ221−a12,c3=Ψ22a22−2Ψ221−a1a2−a1Ψ32,c4=−Ψ22a22+a2Ψ32.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 Figure3, 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 Ψ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 Ψ1 and the lumped transmission field Ψ3 can be expressed as(11)Ψ21−B⋅SN−1⋅1−R=Ψ3,Ψ21−B⋅S2N−1⋅−R⋅1−R+Ψ0⋅R=Ψ1.Therefore, as long as we accurately obtain the correspondingΨ1,Ψ2,andΨ3 for the incident field Ψ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 Figure4(a) for any given incident field Ψ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 Figure5.Figure 5
The structure of the neural networks.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Ψ1,Ψ2,andΨ3, the cost function is the mean square error between them, and their corresponding label data are Ψ1∗,Ψ2∗,andΨ2∗.(12)MSEψ=1N∑m=1NΨm−Ψm∗2.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α1andα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 6
Schematic diagram of the PINN structure.We derive the setting of the cost function as follows: in Figure7, 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 zj=j∆zj=0,1,…,N, the horizontal grid points as xm=m∆xm=0,1,…,M, and the field Ψxm,zj at the position xm,zj as xm=m∆xm=0,1,…,M.Figure 7
The meshing of the FD method.Therefore, by substituting the difference expressions of the second-order partial differential, that is,(13)∂2Ψjm∂z2≈1∆z2Ψj+1m+Ψj−1m−2Ψjm,∂2Ψjm∂x2≈1∆x2Ψjm+1+Ψjm−1−2Ψjm,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:(14)Ψjm∆x2+−2∆z2−2∆x2+k2n2Ψjm+Ψj+1m+Ψj−1m∆z2+Ψjm−1∆x2=0.The upper and lower impedance boundary conditions are(15)∂ΨNm∂z+α1⋅ΨNm=0,∂Ψ0m∂z+α2⋅Ψ0m=0,where α1andα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(16)ΨNm+1∆x2+−2+2α1∆z∆z2+k2n2ΨNm+2⋅ΨN−1m∆z2+ΨNm−1∆x2=0,Ψ0m+1∆x2+−2+2α2∆z∆z2+k2n2Ψ0m+2⋅Ψ1m∆z2+Ψ0m−1∆x2=0.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Ψ1: total field at the start of obstacle calculated by MoMΨ2/∆x+1: total field at the end of obstacle calculated by MoMthreshold: threshold of iterative convergenceN: maximum number of iterationsα1,α2: coefficients of upper/lower impedance boundary conditionOutput: Cost: the cost value for back propagationInitialize: Ψ2:2/∆x=0;      //initial fields to be solvedn=1;       //initial number of iterationserror=1000;      //initial error of iterationsΨ1/∆x+1∗=Ψ1/∆x+1      //intermediate variableWhileerror>thresholdorn<Ndoforj=2:lengthΨ1−1doform=2:2/∆xdoΨ0m=Ψ0m+1/∆x2+2Ψ1m/∆z2+Ψ0m−1/∆x2/2−2α1∆z/∆z2−k2n2Ψ0m=Ψ0m+1/∆x2+2Ψ1m/∆z2+Ψ0m−1/∆x2/2−2α2∆z/∆z2−k2n2Ψjm=Ψjm+1/∆x2+Ψj+1m+Ψj−1m/∆z2+Ψjm−1/∆x2/2/∆z2+2/∆x2−k2n2error=MSEΨ1/∆x+1∗−Ψ1/∆x+1Ψ1/∆x+1∗=Ψ1/∆x+1endendendCost=Ψ0m−Ψ0m+1/∆x2+2Ψ1m/∆z2+Ψ0m−1/∆x2/2−2α1∆z/∆z2−k2n2+Ψ0m−Ψ0m+1/∆x2+2Ψ1m/∆z2+Ψ0m−1/∆x2/2−2α2∆z/∆z2−k2n2+∑j=1N−1Ψjm=Ψjm+1/∆x2+Ψj+1m+Ψj−1m/∆z2+Ψjm−1/∆x2/2/∆z2+2/∆x2−k2n2returnCostIn 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 Algorithm1. 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:(17)ajl+1=σzjl+1,zkl+1=∑jwkjl+1ajl+bkl+1=∑jwkjl+1σzjl+bkl+1,where k represents the order number of neurons in the layer, ajl+1 is the output of the jth neuron in the (l+1)th layer, and wjkl+1 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 as(18)σz=expz−exp−zexpz+exp−z≜σzR+iσzI,the derivative of which is(19)∂σz∂z∗=∂σzR+iσzI∂zR+i∂σzR+iσzI∂zI=σ2zI−σ2zR.


### 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 thejth neuron in the lth layer is defined as the change of the cost function and denoted as the error function δjl.(20)δjl=∇zjiC=∂C∂zjl∗=∂C∂zj,Rl+i⋅∂C∂zj,Il.The result of the derivation of the weight by the cost functionC 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(21)wjkl+1t+1=wjkl+1t+∆wjkl+1t=wjk.Rl+∆wjk.Rl+iwjk,Il+∆wjk,Il,bjl+1t+1=bjl+1t+∆bjl+1t=bjk.Rl+∆bjk.Rl+ibjk,Il+∆bjk,Il,where(22)∆wjkl=−η∇wC=∂C∂wjkl∗=∂C∂wjk,Rl+j∂C∂wjk,Il=∂C∂zjl∗∂zjl∗∂wjkl∗=−ηδjlakl−1∗,∆bjl=−η∇bjlC=−η∂C∂bjl∗=−ηδjl,where η denotes the learning rate and akl−1∗ denotes the complex conjugate of akl−1. The abovementioned formula shows that the increment of the updated weight is the inverse of the product of the learning rate η, the error function δjl, and the conjugate of the output value of the neuron in the previous layer akl−1∗.The iterative calculation formulas ofδl from the back to the front are given, the initial value in the reverse calculation process is the output, and δkL 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, δkL can be expressed as(23)δkL=∂C∂zkL∗=∂∑k=1KakL−Tk2∂zkL∗=∂akL−TkakL−Tk∗∂akL∗⋅∂akL∗∂zkL∗+∂akL−TkakL−Tk∗∂akL⋅∂akL∂zkL∗=ak,RL−Tk,R⋅1−σ2zk,RL+i⋅ak,IL−Tk,I⋅1−σ2zk,IL.Similarly, we can derive(24)δkL∗=∂C∂zkL=ReakL−Tk⋅1−σ2zk,RL−i⋅ImakL−Tk⋅1−σ2zk,IL.Then, the complex error of the (L−1)th layer is(25)δjL−1=∂C∂zjL−1∗=∂C∂ajL−1∗⋅∂ajL−1∗∂zjL−1∗+∂C∂ajL−1⋅∂ajL−1∂zjL−1∗.Thus, we can obtain(26)∂C∂ajL−1∗=∑k∂C∂zkL∗⋅∂zkL∗∂ajL−1∗+∂C∂zkL⋅∂zkL∂ajL−1∗=∑kδkL⋅wjkL∗,∂C∂ajL−1=∑k∂C∂zkL⋅∂zkL∂ajL−1+∂C∂zkL∗⋅∂zkL∗∂ajL−1=∑kδkL∗⋅wjkL.Therefore,(27)δjL−1=∂C∂zjL−1∗=∑kδkL⋅wjkL∗⋅∂ajL−1∗∂zjL−1∗+∑kδkL∗⋅wjkL⋅∂ajL−1∂zjL−1∗,=1−σ2zj,RL−1⋅Re∑kδkL⋅wjkL∗+i⋅1−σ2zj,IL−1⋅Im∑kδkL⋅wjkL∗.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.


## 3.4.1. Forward-Propagation Process

In the fully connected network layer, the output of each layer can be written as follows:(17)ajl+1=σzjl+1,zkl+1=∑jwkjl+1ajl+bkl+1=∑jwkjl+1σzjl+bkl+1,where k represents the order number of neurons in the layer, ajl+1 is the output of the jth neuron in the (l+1)th layer, and wjkl+1 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 as(18)σz=expz−exp−zexpz+exp−z≜σzR+iσzI,the derivative of which is(19)∂σz∂z∗=∂σzR+iσzI∂zR+i∂σzR+iσzI∂zI=σ2zI−σ2zR.


## 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 thejth neuron in the lth layer is defined as the change of the cost function and denoted as the error function δjl.(20)δjl=∇zjiC=∂C∂zjl∗=∂C∂zj,Rl+i⋅∂C∂zj,Il.The result of the derivation of the weight by the cost functionC 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(21)wjkl+1t+1=wjkl+1t+∆wjkl+1t=wjk.Rl+∆wjk.Rl+iwjk,Il+∆wjk,Il,bjl+1t+1=bjl+1t+∆bjl+1t=bjk.Rl+∆bjk.Rl+ibjk,Il+∆bjk,Il,where(22)∆wjkl=−η∇wC=∂C∂wjkl∗=∂C∂wjk,Rl+j∂C∂wjk,Il=∂C∂zjl∗∂zjl∗∂wjkl∗=−ηδjlakl−1∗,∆bjl=−η∇bjlC=−η∂C∂bjl∗=−ηδjl,where η denotes the learning rate and akl−1∗ denotes the complex conjugate of akl−1. The abovementioned formula shows that the increment of the updated weight is the inverse of the product of the learning rate η, the error function δjl, and the conjugate of the output value of the neuron in the previous layer akl−1∗.The iterative calculation formulas ofδl from the back to the front are given, the initial value in the reverse calculation process is the output, and δkL 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, δkL can be expressed as(23)δkL=∂C∂zkL∗=∂∑k=1KakL−Tk2∂zkL∗=∂akL−TkakL−Tk∗∂akL∗⋅∂akL∗∂zkL∗+∂akL−TkakL−Tk∗∂akL⋅∂akL∂zkL∗=ak,RL−Tk,R⋅1−σ2zk,RL+i⋅ak,IL−Tk,I⋅1−σ2zk,IL.Similarly, we can derive(24)δkL∗=∂C∂zkL=ReakL−Tk⋅1−σ2zk,RL−i⋅ImakL−Tk⋅1−σ2zk,IL.Then, the complex error of the (L−1)th layer is(25)δjL−1=∂C∂zjL−1∗=∂C∂ajL−1∗⋅∂ajL−1∗∂zjL−1∗+∂C∂ajL−1⋅∂ajL−1∂zjL−1∗.Thus, we can obtain(26)∂C∂ajL−1∗=∑k∂C∂zkL∗⋅∂zkL∗∂ajL−1∗+∂C∂zkL⋅∂zkL∂ajL−1∗=∑kδkL⋅wjkL∗,∂C∂ajL−1=∑k∂C∂zkL⋅∂zkL∂ajL−1+∂C∂zkL∗⋅∂zkL∗∂ajL−1=∑kδkL∗⋅wjkL.Therefore,(27)δjL−1=∂C∂zjL−1∗=∑kδkL⋅wjkL∗⋅∂ajL−1∗∂zjL−1∗+∑kδkL∗⋅wjkL⋅∂ajL−1∂zjL−1∗,=1−σ2zj,RL−1⋅Re∑kδkL⋅wjkL∗+i⋅1−σ2zj,IL−1⋅Im∑kδkL⋅wjkL∗.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 Figure5, 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Ψ1,Ψ2,andΨ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π×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π×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 Section3.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 Ψ1,Ψ2,andΨ3 during the training process are presented in Figure 8.Figure 8
Change of cost in the training and validation process.The total training duration of the network corresponding toΨ1,Ψ2,andΨ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 Table1.Table 1
Training results.
OutputsMSE in the training setMSE in the validation setΨ11.4030e−64.4306e−6Ψ22.9402e−66.5087e−5Ψ32.0135e−63.1585e−5α1,α22.5353e−42.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 Figure9. 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π×10−7 H/m, and the relative permittivity is 6. Furthermore, the magnetic permeability of the medium obstacle is 4π×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 Z0−1 W when Z0 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 Figure10, 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 9
Schematic diagram of the propagation environment.Figure 10
The simulation results for Example1.Example 2.
We test the generalization ability of the training results in the double-rectangular obstacle environment shown in Figure11. 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 Figure12. 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.
Figure13 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 11
Schematic diagram of the double obstacle environment.Figure 12
The simulation results for Example2.Figure 13
The distribution of the spatial field corresponding to Example2. (a) The spatial field calculated by MoM. (b) The spatial field calculated by original 2W-PE. (c) The spatial field calculated by new 2W-PE.
(a)(b)(c)


## 4.1. The Training Process

As shown in Figure5, 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Ψ1,Ψ2,andΨ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π×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π×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 Section3.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 Ψ1,Ψ2,andΨ3 during the training process are presented in Figure 8.Figure 8
Change of cost in the training and validation process.The total training duration of the network corresponding toΨ1,Ψ2,andΨ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 Table1.Table 1
Training results.
OutputsMSE in the training setMSE in the validation setΨ11.4030e−64.4306e−6Ψ22.9402e−66.5087e−5Ψ32.0135e−63.1585e−5α1,α22.5353e−42.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 Figure9. 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π×10−7 H/m, and the relative permittivity is 6. Furthermore, the magnetic permeability of the medium obstacle is 4π×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 Z0−1 W when Z0 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 Figure10, 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 9
Schematic diagram of the propagation environment.Figure 10
The simulation results for Example1.Example 2.
We test the generalization ability of the training results in the double-rectangular obstacle environment shown in Figure11. 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 Figure12. 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.
Figure13 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 11
Schematic diagram of the double obstacle environment.Figure 12
The simulation results for Example2.Figure 13
The distribution of the spatial field corresponding to Example2. (a) The spatial field calculated by MoM. (b) The spatial field calculated by original 2W-PE. (c) The spatial field calculated by 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.



---

*Source: 9485036-2023-02-10.xml*