# Air Combat Maneuver Trajectory Prediction Model of Target Based on Chaotic Theory and IGA-VNN

Authors: Zhi-fei Xi; An Xu; Ying-xin Kou; Zhan-wu Li; Ai-wu Yang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Target maneuver trajectory prediction is an important prerequisite for air combat situation awareness and threat assessment. Aiming at the problem of low prediction accuracy in traditional trajectory prediction methods, combined with the chaotic characteristics of the target maneuver trajectory time series, a target maneuver trajectory prediction model based on chaotic theory and improved genetic algorithm-Volterra neural network (IGA-VNN) model is proposed, mathematically deducing and analyzing the consistency between Volterra functional model and back propagation (BP) neural network in structure. Firstly, the C-C method is used to reconstruct the phase space of the target trajectory time series, and the maximum Lyapunov exponent of the time series of the target maneuver trajectory is calculated. It is proved that the time series of the target maneuver trajectory has chaotic characteristics, so the chaotic method can be used to predict the target trajectory time series. Then, the practicable Volterra functional model and BP neural network are combined together, learning the advantages of both and overcoming the difficulty in obtaining the high-order kernel function of the Volterra functional model. At the same time, an adaptive crossover mutation operator and a combination mutation operator based on the difference degree of gene segments are proposed to improve the traditional genetic algorithm; the improved genetic algorithm is used to optimize BP neural network, and the optimal initial weights and thresholds are obtained. Finally, the IGA-VNN model of chaotic time series is applied to the prediction of target maneuver trajectory time series, and the experimental results show that its estimated performance is obviously superior to other prediction algorithms.


---

## Body

## 1. Introduction

Maneuvering trajectory prediction is a process of learning and reasoning the inherent information contained in the target’s historical trajectory and then making a reasonable prediction of the target’s future trajectory. In modern air combat, it is of great significance to make a reasonable prediction of the maneuvering trajectory of enemy targets. According to the OODA (observation, orientation, decision, and action) cycle theory, the essence of winning an air is to form an OODA cycle prior to the enemy and accurately predict the target maneuver trajectory so as to achieve the purpose of preemption [1].In recent years, the research on target maneuver trajectory prediction methods can be divided into two categories. One is the traditional method based on Kalman filter algorithm,α/β filter algorithm, linear regression model, and particle motion model. In view of the complex and changeable characteristics of target motion patterns, a polynomial Kalman filter for the motion trajectory prediction algorithm is proposed for the change of the target motion mode [2]. Aiming at the problem of missing historical position information of the target, an improved Kalman filter algorithm with system noise estimation is proposed to predict the target maneuver trajectory [3]. Owing to general traditional fitting-based trajectory prediction algorithms cannot meet the requirements of high accuracy and real-time prediction, a dynamic Kalman filter for trajectory prediction approach is proposed [4]. In order to solve the problems of high nonlinear, difficult data processing, and low prediction accuracy in the high-order motion model of the target, an interactive multimodel trajectory prediction algorithm with control variables is proposed [5]. According to the above analysis, the traditional prediction method is only applicable to the relatively simple trajectory prediction of the target motion characteristics. However, in the process of air combat, the target motion is often a highly nonlinear and complex time sequence process, which is also affected by many uncertain factors. The traditional trajectory prediction algorithms cannot reflect the target’s maneuvering characteristics very well. In addition, the complexity of the trajectory prediction model is high, the adaptability of the prediction algorithm is poor, and the prediction accuracy cannot meet the requirements of air combat.The other algorithm is a machine learning prediction algorithm based on intelligent algorithms, which are mainly based on big data to establish a target maneuvering trajectory prediction model. In view of GRNN’s (generalized regression neural network) good nonlinear mapping ability, flexible network structure, and high fault tolerance and robustness, a trajectory prediction model based on GRNN neural network is proposed [6]; the target group trajectory is used to train the BP neural network, establishes a trajectory prediction model, and realizes the prediction of the flight trajectory [7]. However, the BP neural network has the disadvantages of being sensitive to the initial value and not having the global search ability. In order to solve these problems, the genetic algorithm and particle swarm optimization algorithm with global search ability are, respectively, used to optimize the neural network weights to improve the prediction accuracy [8, 9]. In addition, the essence of target maneuver trajectory prediction is the time series prediction problem, which is highly nonlinear and time-varying. The theory and practice show that the Volterra functional model can well represent the nonlinear system and be used to accurately predict the time series. However, the solution of the kernel function of the higher-order Volterra functional model is the bottleneck of its application. Due to the good adaptability, parallelism and fault tolerance of BP neural network, and the ability to approximate any nonlinear function, the Volterra functional model and BP neural network have the same structure [10–12]. Therefore, in this paper, a target maneuver trajectory prediction model based on BP neural network and Volterra series is proposed. In this paper, the theory of phase space reconstruction is used to reconstruct the time series of the target’s maneuvering trajectory, and the small data method is used to prove that the time series of the target’s maneuvering trajectory has chaotic characteristics. Then, the mathematical equivalence between the BP neural network model and the Volterra functional model is proved by theoretical analysis. Secondly, combining the chaotic characteristics of the target maneuver trajectory time series, an IGA-VNN target maneuver trajectory time series prediction model based on chaos theory is established. The model combines the advantages of the Volterra functional model and BP neural network and overcomes the difficulty of the Volterra series model in determining higher-order kernel function and the blindness of BP neural network modeling, and an improved genetic algorithm is proposed to optimize the weights and thresholds of VNN networks simultaneously. Finally, the performance of the prediction model proposed in this paper is verified by simulation. The simulation results show that the model can accurately and rapidly predict the maneuvering trajectory of the target, which provides a new way to solve the problem of trajectory prediction.


## 2. Dynamics Theory of Chaotic System

### 2.1. Phase Space Reconstruction Theory

PSRT (phase space reconstruction theory) is an effective method to analyze nonlinear time series, and it is the basis for determining and predicting chaotic characteristics of time series. The basic idea of phase space reconstruction theory is to regard time series as the component generated by nonlinear dynamic series. According to Takens theorem [13, 14], for a chaotic time series, when the embedding dimension m and attractor dimension D of the phase space reconstruction meet the condition of m≥2D+1, the attractor of the original time series can be replaced in the reconstructed phase space. At the same time, the time series after phase space reconstruction is topologically equivalent to the original system on the basis of differential homeomorphism [15]. If there is a time series xi,i=1,2,…,N, the time series reconstructed by phase space can be expressed as follows:(1)X=X1X2⋯XMT=x1x2⋯xMx1+τx2+τ⋯xM+τ⋮⋮⋱⋮x1+m−1τx2+m−1τ⋯xM+m−1τ,where N is the length of the time series; m is the embedding dimension; τ is the delay time; and M=N−m−1τ.In the process of reconstructing phase space, the determination of time delay and embedding dimension plays a crucial role in the quality of phase space reconstruction. If the embedding dimension is too low, the self-interaction of attractors will appear; if the parameter is too high, the distance between points will be too far. If the time delay is too small, the correlation between the adjacent points of the reconstructed attractor is too strong, and the analysis of the attractor is easily interfered by noise; if the parameter is too large, the originally close vectors will also be far away, resulting in an uncertain system state [16]. Therefore, it is very important to determine reasonable time delay and embedding dimension for phase space reconstruction.


### 2.2. C-C Method

The traditional calculation method regards the embedding dimension and time delay as two independent parameters, which need to be calculated independently. The research in recent years shows that the embedding dimension and time delay are interrelated. The determination of time delay should not be independent of the embedding dimension but should be determined by combining the time windowτw=m−1τ.In this paper, the C-C method [15] is used to determine the time delay and embedding dimension. The C-C method uses the correlation integral of time series to form statistics, which reflects the correlation of nonlinear time series. Based on the relationship between time delay and statistics, the time window and time delay can be calculated simultaneously so that the embedding dimension can be determined. The correlation integral is defined as follows:(2)Cm,r,t=2MM−1∑1≤i≤j≤MHr−dij,dij=Xi−Xj,where M is the number of phase points; r is the size of neighborhood radius; dij is the Euclidean distance between two phase points in phase space; and Hz is the Heaviside unit function, which is defined as follows:(3)Hz=1,z>0,0,z≤0.Dividing time seriesxi,i=1,2,…,N into t independent subtime series, the length of the subtime series is INTN/t, where INT is an integer function. For t disjoint subtime series, the test statistic Sm,N,r,t and ΔSm,t can be expressed as follows:(4)Sm,N,r,t=1t∑s=1tCsm,Nt,r,t−Csm1,Nt,r,t,S¯=1M⋅J∑m=1M∑j=1JSm,rj,t,where J is the number of r.In order to express the relationship betweenSm,N,r,t and r, the difference ΔSm,t between the maximum and minimum values of the two radius is expressed as follows:(5)ΔSm,t=maxSm,rj,t−minSm,rj,t.According to Brock mathematical statistical conclusions [16], in general, we can obtain the range values of m and r as m=2,3,4,5,rj=jσ/2,andj=1,2,3,4, where σ is the standard deviation of the time series. The three test statistical variables are defined as follows:(6)S¯t=116∑m=25∑j=14Sm,rj,t,ΔSt=14∑m=25ΔSm,t,Scort=ΔS¯t+S¯t,where the time delay τ is the time corresponding to the first minimum of ΔS¯t or the time corresponding to the first zero of S¯t. When Scort obtains the minimum value, the corresponding time is the time window τw. The C-C method uses statistical theory to obtain the time delay and embedding dimension. The calculation is simple, the calculation amount is small, and it has strong anti-interference ability.Based on the theoretical analysis of PSR in Section2.1, the C-C method is used to reconstruct the phase space of the three-dimensional coordinate time series of the target maneuver trajectory. The test statistics ΔS¯t, S¯t, and Scort calculated by the C-C method are shown in Figures 1–3. It can be seen from Figure 1 that the first minimum point of ΔS¯t is τx=19, and the minimum point of Scort is τwx=126. Then, based on the embedded time window equation τw=m−1τ, the optimal embedded dimension is moptx=7. As can be seen from Figure 2, the first minimum point of ΔS¯t is τy=25, and the minimum point of Scort is τwy=181. Then, based on the embedded time window equation τw=m−1τ, the optimal embedded dimension is myopt=8. Similarly, as can be seen from Figure 3, the first minimum point of ΔS¯t is τz=21, and the minimum point of Scort is τwz=157. Then, based on the embedded time window equation τw=m−1τ, the optimal embedded dimension is moptz=8. The determination of the above parameters provides a basis for proving that the target maneuver trajectory prediction system has chaotic characteristics.Figure 1
Statistics curves of the C-C method to reconstruct theX-direction trajectory time series.Figure 2
Statistics curves of the C-C method to reconstruct theY-direction trajectory time series.Figure 3
Statistics curves of the C-C method to reconstruct theZ-direction trajectory time series.


### 2.3. Calculation of Maximum Lyapunov Exponent

Butterfly effect is the basic form of the chaos system; that is, chaos is sensitive to the initial value. The trajectories generated by two initial values with little difference will show exponential separation with time. Lyapunov exponent is the eigenvalue of the chaotic system after averaging on the trajectories of the whole attractor [14]. Among them, the largest Lyapunov index can quantitatively describe the separation of orbits generated by two very close initial values as time goes by. It is theoretically proved in [17] that as long as the maximum Lyapunov exponent λ satisfies λ>0, the system has chaotic characteristics.The small data is one of the most commonly used methods to calculate the maximum Lyapunov exponent of the chaotic system. The specific calculation steps are as follows:Step 1: perform fast Fourier transform on the time seriesxi,i=1,2,…,N to determine the average period P.Step 2: use the C-C method to determine time delayτ and embedding dimension m and reconstruct phase space of time series Yj|j=1,2,…,M.Step 3: find the closest phase pointYj^ to each phase point Yj in the phase space after time series reconstruction and limit the short-term separation, namely:(7)dj0=minYj−Yj^,j−j^>p.Step 4: for each phase pointYj in the phase space, calculate the distance dji after j discrete-time steps of the adjacent point pair:(8)dji=minYj+i−Yj^+i,i=1,2,…,minM−j,M−j^.Step 5: for each discrete time stepi, calculate the average yj of all lndji:(9)yj=1qh∑i=1qlndji,whereq is the number that dji is nonzero. The least square method is used to make the regression line, and the slope of the line is the largest Lyapunov exponent λ. This method makes full use of the spatial evolution information of time series, which is more reliable for small data groups, less computation, and easy to operate.Based on the theoretical analysis of the small data method in Section2.3, the small data method is used to analyze the chaotic characteristics of the three-dimensional coordinate time series of the target maneuver trajectory. The least square fitting curve of the target maneuver trajectory time series and the corresponding slope changes are shown in Figures 4–6. Based on the small data method, the maximum Lyapunov exponent of the target maneuver trajectory in X-direction time series, Y-direction time series, and Z-direction time series is calculated as LyX=0.0011, LyY=0.0042, and LyZ=0.004682, respectively. It can be seen from the above that the maximum Lyapunov exponents of the three-dimensional coordinate time series of the target’s maneuvering trajectory are all positive numbers, which shows that the time series of the target’s maneuvering trajectory has chaotic characteristics.Figure 4
Fitted curve by the least square method and the corresponding slope changes of theX-direction trajectory time series.Figure 5
Fitted curve by the least square method and the corresponding slope changes of theY-direction trajectory time series.Figure 6
Fitted curve by the least square method and the corresponding slope changes of theZ-direction trajectory time series.


## 2.1. Phase Space Reconstruction Theory

PSRT (phase space reconstruction theory) is an effective method to analyze nonlinear time series, and it is the basis for determining and predicting chaotic characteristics of time series. The basic idea of phase space reconstruction theory is to regard time series as the component generated by nonlinear dynamic series. According to Takens theorem [13, 14], for a chaotic time series, when the embedding dimension m and attractor dimension D of the phase space reconstruction meet the condition of m≥2D+1, the attractor of the original time series can be replaced in the reconstructed phase space. At the same time, the time series after phase space reconstruction is topologically equivalent to the original system on the basis of differential homeomorphism [15]. If there is a time series xi,i=1,2,…,N, the time series reconstructed by phase space can be expressed as follows:(1)X=X1X2⋯XMT=x1x2⋯xMx1+τx2+τ⋯xM+τ⋮⋮⋱⋮x1+m−1τx2+m−1τ⋯xM+m−1τ,where N is the length of the time series; m is the embedding dimension; τ is the delay time; and M=N−m−1τ.In the process of reconstructing phase space, the determination of time delay and embedding dimension plays a crucial role in the quality of phase space reconstruction. If the embedding dimension is too low, the self-interaction of attractors will appear; if the parameter is too high, the distance between points will be too far. If the time delay is too small, the correlation between the adjacent points of the reconstructed attractor is too strong, and the analysis of the attractor is easily interfered by noise; if the parameter is too large, the originally close vectors will also be far away, resulting in an uncertain system state [16]. Therefore, it is very important to determine reasonable time delay and embedding dimension for phase space reconstruction.


## 2.2. C-C Method

The traditional calculation method regards the embedding dimension and time delay as two independent parameters, which need to be calculated independently. The research in recent years shows that the embedding dimension and time delay are interrelated. The determination of time delay should not be independent of the embedding dimension but should be determined by combining the time windowτw=m−1τ.In this paper, the C-C method [15] is used to determine the time delay and embedding dimension. The C-C method uses the correlation integral of time series to form statistics, which reflects the correlation of nonlinear time series. Based on the relationship between time delay and statistics, the time window and time delay can be calculated simultaneously so that the embedding dimension can be determined. The correlation integral is defined as follows:(2)Cm,r,t=2MM−1∑1≤i≤j≤MHr−dij,dij=Xi−Xj,where M is the number of phase points; r is the size of neighborhood radius; dij is the Euclidean distance between two phase points in phase space; and Hz is the Heaviside unit function, which is defined as follows:(3)Hz=1,z>0,0,z≤0.Dividing time seriesxi,i=1,2,…,N into t independent subtime series, the length of the subtime series is INTN/t, where INT is an integer function. For t disjoint subtime series, the test statistic Sm,N,r,t and ΔSm,t can be expressed as follows:(4)Sm,N,r,t=1t∑s=1tCsm,Nt,r,t−Csm1,Nt,r,t,S¯=1M⋅J∑m=1M∑j=1JSm,rj,t,where J is the number of r.In order to express the relationship betweenSm,N,r,t and r, the difference ΔSm,t between the maximum and minimum values of the two radius is expressed as follows:(5)ΔSm,t=maxSm,rj,t−minSm,rj,t.According to Brock mathematical statistical conclusions [16], in general, we can obtain the range values of m and r as m=2,3,4,5,rj=jσ/2,andj=1,2,3,4, where σ is the standard deviation of the time series. The three test statistical variables are defined as follows:(6)S¯t=116∑m=25∑j=14Sm,rj,t,ΔSt=14∑m=25ΔSm,t,Scort=ΔS¯t+S¯t,where the time delay τ is the time corresponding to the first minimum of ΔS¯t or the time corresponding to the first zero of S¯t. When Scort obtains the minimum value, the corresponding time is the time window τw. The C-C method uses statistical theory to obtain the time delay and embedding dimension. The calculation is simple, the calculation amount is small, and it has strong anti-interference ability.Based on the theoretical analysis of PSR in Section2.1, the C-C method is used to reconstruct the phase space of the three-dimensional coordinate time series of the target maneuver trajectory. The test statistics ΔS¯t, S¯t, and Scort calculated by the C-C method are shown in Figures 1–3. It can be seen from Figure 1 that the first minimum point of ΔS¯t is τx=19, and the minimum point of Scort is τwx=126. Then, based on the embedded time window equation τw=m−1τ, the optimal embedded dimension is moptx=7. As can be seen from Figure 2, the first minimum point of ΔS¯t is τy=25, and the minimum point of Scort is τwy=181. Then, based on the embedded time window equation τw=m−1τ, the optimal embedded dimension is myopt=8. Similarly, as can be seen from Figure 3, the first minimum point of ΔS¯t is τz=21, and the minimum point of Scort is τwz=157. Then, based on the embedded time window equation τw=m−1τ, the optimal embedded dimension is moptz=8. The determination of the above parameters provides a basis for proving that the target maneuver trajectory prediction system has chaotic characteristics.Figure 1
Statistics curves of the C-C method to reconstruct theX-direction trajectory time series.Figure 2
Statistics curves of the C-C method to reconstruct theY-direction trajectory time series.Figure 3
Statistics curves of the C-C method to reconstruct theZ-direction trajectory time series.


## 2.3. Calculation of Maximum Lyapunov Exponent

Butterfly effect is the basic form of the chaos system; that is, chaos is sensitive to the initial value. The trajectories generated by two initial values with little difference will show exponential separation with time. Lyapunov exponent is the eigenvalue of the chaotic system after averaging on the trajectories of the whole attractor [14]. Among them, the largest Lyapunov index can quantitatively describe the separation of orbits generated by two very close initial values as time goes by. It is theoretically proved in [17] that as long as the maximum Lyapunov exponent λ satisfies λ>0, the system has chaotic characteristics.The small data is one of the most commonly used methods to calculate the maximum Lyapunov exponent of the chaotic system. The specific calculation steps are as follows:Step 1: perform fast Fourier transform on the time seriesxi,i=1,2,…,N to determine the average period P.Step 2: use the C-C method to determine time delayτ and embedding dimension m and reconstruct phase space of time series Yj|j=1,2,…,M.Step 3: find the closest phase pointYj^ to each phase point Yj in the phase space after time series reconstruction and limit the short-term separation, namely:(7)dj0=minYj−Yj^,j−j^>p.Step 4: for each phase pointYj in the phase space, calculate the distance dji after j discrete-time steps of the adjacent point pair:(8)dji=minYj+i−Yj^+i,i=1,2,…,minM−j,M−j^.Step 5: for each discrete time stepi, calculate the average yj of all lndji:(9)yj=1qh∑i=1qlndji,whereq is the number that dji is nonzero. The least square method is used to make the regression line, and the slope of the line is the largest Lyapunov exponent λ. This method makes full use of the spatial evolution information of time series, which is more reliable for small data groups, less computation, and easy to operate.Based on the theoretical analysis of the small data method in Section2.3, the small data method is used to analyze the chaotic characteristics of the three-dimensional coordinate time series of the target maneuver trajectory. The least square fitting curve of the target maneuver trajectory time series and the corresponding slope changes are shown in Figures 4–6. Based on the small data method, the maximum Lyapunov exponent of the target maneuver trajectory in X-direction time series, Y-direction time series, and Z-direction time series is calculated as LyX=0.0011, LyY=0.0042, and LyZ=0.004682, respectively. It can be seen from the above that the maximum Lyapunov exponents of the three-dimensional coordinate time series of the target’s maneuvering trajectory are all positive numbers, which shows that the time series of the target’s maneuvering trajectory has chaotic characteristics.Figure 4
Fitted curve by the least square method and the corresponding slope changes of theX-direction trajectory time series.Figure 5
Fitted curve by the least square method and the corresponding slope changes of theY-direction trajectory time series.Figure 6
Fitted curve by the least square method and the corresponding slope changes of theZ-direction trajectory time series.


## 3. Volterra Functional Model of Nonlinear System

The Volterra filter is usually used in signal processing and model recognition of the nonlinear system [18], and the nonlinear expression ability is good, so the Volterra functional model can realize accurate prediction of chaotic time series. For nonlinear systems, the discretized Volterra functional model can be expressed as follows:(10)yk=∑i=1k∑l1,l2,…,li∞hil1,l2,…,li⋅xk−l1⋅xk−l2⋯xk−li,where i,li,k∈R; xk−li is the input of the system; yk is the output of the system, which is the state of the system at the next moment; and hil1,l2,…,li is the Volterra kernel function of order i.In theory, the discrete Volterra functional model can accurately predict the nonlinear time series. When the Volterra filter is used to predict the time series, because it is difficult for the functional model to solve the higher-order kernel function, the second-order or third-order truncation form of the Volterra functional model will be used in general, but this will greatly reduce the prediction accuracy and performance of the model [19–21]. In addition, the Volterra functional model has limited memory ability. When time point k−k0 is far away from the time point k, the input xk−k0 will not affect the output. Based on the above analysis, the practical Volterra functional model adopted in this paper is as follows:(11)yk=h0+∑l1=0m−1h1l1xk−l1τ+∑l1=0m−1∑l2=0m−1h2l1,l2xk−l1τxk−l2τ+⋯+∑l1=0m−1∑l2=0m−1⋯+∑lm=0m−1hml1,l2,…,lmxk−l1τxk−l2τ⋯xk−lmτ,where m is the memory length of the Volterra functional system, that is, the minimum embedding space dimension of the phase space reconstruction of the target maneuver trajectory time series, and τ is the time delay.Based on equation (11), it is possible to accurately predict the target maneuver trajectory. Through the above analysis, we can know that different time series have different characteristics, and the time delay and embedding dimension determined by the C-C method are also different. Therefore, the Volterra functional series prediction model used in the prediction of the three-dimensional coordinates of the target’s maneuvering trajectory is different.


## 4. Prediction Model of Chaotic Time Series Based on VNN

### 4.1. Model Equivalence Analysis

Them dimensional input, single hidden layer, and single output BP neural network are established, and its structure is shown in Figure 7. The input vector of the BP neural network is XTk=xk,0,xk,1,…,xk,n,…,xk,m−1, and the input is obtained by xk through phase space reconstruction theory, where xk,m=xk−m. The output of the lth unit of the hidden layer is as follows:(12)Zl,k=SlUl,k,l=1,2,…,L.Figure 7
The structure of three-layer neural networks.The discrete convolutionUl,k of system input time series is expressed by(13)Ul,k=∑n=0m−1wl,nxk,n.The incentive functionSl⋅ of the hidden layer of BP neural network is Sigmoid function and is expressed by(14)SlUl,k=11−exp−λUl,k−θl,where wl,n is the weight between the input layer and the hidden layer; rl is the weight between the hidden layer and the output layer; and λ is a fixed value, and the parameter indicates the transfer slope from the input layer to the hidden layer.The Taylor series expansion of each output of the hidden layer at the thresholdθl can be expressed by(15)Zl,k=SlUl,k=∑i=0∞diθlUl,ki,where diθl is the coefficient after the output of the hidden layer is expanded into a series form. The coefficient is a function of the hidden layer threshold BB of the BP neural network and is related to the excitation function of the hidden layer. The output of BP neural network is expressed in the form of linear summation. Based on Taylor series, the output of BP neural network can be expressed as follows:(16)yk=∑l=1LrlZl,k=∑l=1Lrl∑i=0∞diθlUl,ki=∑l=1Lrl∑i=0∞diθl∑n1=0m−1⋯∑ni=0m−1wl,n1⋯wl,nixl,n1⋯xl,ni.On the other hand, if the Volterra kernel function is expanded in the form of basis functionbnz, the deformed Volterra functional model can be obtained from the practical functional model, as shown in Figure 8.Figure 8
Series Volterra model after deformation.In Figure8, Vlk is the weighted sum of the input time series of the Volterra series model, that is, discrete convolution, which can be expressed as follows:(17)Vlk=∑z=0m−1bl,zxk,z.At the same time, the output of the system can be expressed in the polynomial form as follows:(18)yk=fV1k,…,Vlk,…,VLk=c0+∑l0=0m−1c1l0V1k+∑l0=0m−1∑l2=0m−1c2l1,l2⋅Vl1kVl2k+⋯+∑l0=0m−1∑l2=0m−1⋯∑lm=0m−1cml1,l2,…,lm⋅Vl1kVl2k⋯Vlmk,where cml1,l2,…,lm is the expanded polynomial coefficient. Selecting the appropriate basis function bl,z can make the error between each coefficient of the polynomial fV1k,…,Vlk,…,VLk and the corresponding kernel function of the model in the Volterra series approach zero. Based on this, the Volterra series kernel function can be expressed as follows:(19)h1z=∑l=0m−1c1lbl,z,h2z1,z2=∑l1=0m−1∑l2=0m−1c2l1,l2bl1,z1bl2,z2.Through the analysis of equations (13) and (17), we can see that there are implicit variables Ul,k,Vlk in the expressions and they are all discrete convolutions of the input time series.If the basis functionbnz in the practical Volterra functional model can be found, equation (18) can be expressed as the linear weighting of the Sigmoid function:(20)fV1k,…,Vlk,…,VLk=∑l=1LrlSlVlk.Comparing equations (16) and (17), it can be seen that there is an equivalent relationship between Volterra functional series model and BP neural network model, and the two models are completely consistent in theory. The multivariable function is linearly weighted by the univariate function glVlk and can be expressed as follows:(21)fV1k,…,Vlk,…,VLk=∑l=1LrlglVlk,where the univariate function glVlk can be expressed in any form.


### 4.2. Prediction Model of Target Maneuver Trajectory Based on VNN

In Section4.1, the equivalence between BP neural network and practical Volterra functional model is proved theoretically. Based on the equivalence between the two models, combined with BP neural network and Volterra series functional model, a time series prediction model of target maneuvering trajectory based on VNN is proposed. The model structure is shown in Figure 9. In the prediction model, the input vector XT=xk,xk+τ,…,xk+m−1τ of the prediction model is obtained by PSR. The convolution Vlk of target maneuvering trajectory time series can be expressed as follows:(22)Vlk=∑j=0m−1wl,jxk+jτ.Figure 9
Volterra neural network model of chaotic time series.In the model, the polynomial functiongl⋅l=1,2,…,L is selected as the excitation function, and it can be expressed as follows:(23)gl⋅=a0,l+a1,lx+a2,lx2+⋯+ai,lxi+⋯,where ai,l is the coefficient of the polynomial expression. The output expression of the VNN model of target maneuvering trajectory time series can be obtained as follows:(24)yk=∑l=1LrlglVlk=∑l=1L∑i=0∞rlai,lVlki=∑l=1L∑i=0∞rlai,l∑j=0m−1wl,jxk+jτi,where yk is the output of the system, which is the maneuver trajectory of the target at the next moment, and x is the input of the system, which is the historical trajectory of the target.Based on the above theoretical analysis, theith order Volterra series kernel function can be expressed as follows:(25)hiz1,z2,…,zi=∑l=1Lrlai,lwl,z1wl,z2⋯wl,zi.The weightswl,n and rl and thresholds θl of the VNN model are obtained by training, and glVlk is expanded by Taylor series at the threshold. Finally, the kernel functions of each order of the model can be obtained according to equation (25).


## 4.1. Model Equivalence Analysis

Them dimensional input, single hidden layer, and single output BP neural network are established, and its structure is shown in Figure 7. The input vector of the BP neural network is XTk=xk,0,xk,1,…,xk,n,…,xk,m−1, and the input is obtained by xk through phase space reconstruction theory, where xk,m=xk−m. The output of the lth unit of the hidden layer is as follows:(12)Zl,k=SlUl,k,l=1,2,…,L.Figure 7
The structure of three-layer neural networks.The discrete convolutionUl,k of system input time series is expressed by(13)Ul,k=∑n=0m−1wl,nxk,n.The incentive functionSl⋅ of the hidden layer of BP neural network is Sigmoid function and is expressed by(14)SlUl,k=11−exp−λUl,k−θl,where wl,n is the weight between the input layer and the hidden layer; rl is the weight between the hidden layer and the output layer; and λ is a fixed value, and the parameter indicates the transfer slope from the input layer to the hidden layer.The Taylor series expansion of each output of the hidden layer at the thresholdθl can be expressed by(15)Zl,k=SlUl,k=∑i=0∞diθlUl,ki,where diθl is the coefficient after the output of the hidden layer is expanded into a series form. The coefficient is a function of the hidden layer threshold BB of the BP neural network and is related to the excitation function of the hidden layer. The output of BP neural network is expressed in the form of linear summation. Based on Taylor series, the output of BP neural network can be expressed as follows:(16)yk=∑l=1LrlZl,k=∑l=1Lrl∑i=0∞diθlUl,ki=∑l=1Lrl∑i=0∞diθl∑n1=0m−1⋯∑ni=0m−1wl,n1⋯wl,nixl,n1⋯xl,ni.On the other hand, if the Volterra kernel function is expanded in the form of basis functionbnz, the deformed Volterra functional model can be obtained from the practical functional model, as shown in Figure 8.Figure 8
Series Volterra model after deformation.In Figure8, Vlk is the weighted sum of the input time series of the Volterra series model, that is, discrete convolution, which can be expressed as follows:(17)Vlk=∑z=0m−1bl,zxk,z.At the same time, the output of the system can be expressed in the polynomial form as follows:(18)yk=fV1k,…,Vlk,…,VLk=c0+∑l0=0m−1c1l0V1k+∑l0=0m−1∑l2=0m−1c2l1,l2⋅Vl1kVl2k+⋯+∑l0=0m−1∑l2=0m−1⋯∑lm=0m−1cml1,l2,…,lm⋅Vl1kVl2k⋯Vlmk,where cml1,l2,…,lm is the expanded polynomial coefficient. Selecting the appropriate basis function bl,z can make the error between each coefficient of the polynomial fV1k,…,Vlk,…,VLk and the corresponding kernel function of the model in the Volterra series approach zero. Based on this, the Volterra series kernel function can be expressed as follows:(19)h1z=∑l=0m−1c1lbl,z,h2z1,z2=∑l1=0m−1∑l2=0m−1c2l1,l2bl1,z1bl2,z2.Through the analysis of equations (13) and (17), we can see that there are implicit variables Ul,k,Vlk in the expressions and they are all discrete convolutions of the input time series.If the basis functionbnz in the practical Volterra functional model can be found, equation (18) can be expressed as the linear weighting of the Sigmoid function:(20)fV1k,…,Vlk,…,VLk=∑l=1LrlSlVlk.Comparing equations (16) and (17), it can be seen that there is an equivalent relationship between Volterra functional series model and BP neural network model, and the two models are completely consistent in theory. The multivariable function is linearly weighted by the univariate function glVlk and can be expressed as follows:(21)fV1k,…,Vlk,…,VLk=∑l=1LrlglVlk,where the univariate function glVlk can be expressed in any form.


## 4.2. Prediction Model of Target Maneuver Trajectory Based on VNN

In Section4.1, the equivalence between BP neural network and practical Volterra functional model is proved theoretically. Based on the equivalence between the two models, combined with BP neural network and Volterra series functional model, a time series prediction model of target maneuvering trajectory based on VNN is proposed. The model structure is shown in Figure 9. In the prediction model, the input vector XT=xk,xk+τ,…,xk+m−1τ of the prediction model is obtained by PSR. The convolution Vlk of target maneuvering trajectory time series can be expressed as follows:(22)Vlk=∑j=0m−1wl,jxk+jτ.Figure 9
Volterra neural network model of chaotic time series.In the model, the polynomial functiongl⋅l=1,2,…,L is selected as the excitation function, and it can be expressed as follows:(23)gl⋅=a0,l+a1,lx+a2,lx2+⋯+ai,lxi+⋯,where ai,l is the coefficient of the polynomial expression. The output expression of the VNN model of target maneuvering trajectory time series can be obtained as follows:(24)yk=∑l=1LrlglVlk=∑l=1L∑i=0∞rlai,lVlki=∑l=1L∑i=0∞rlai,l∑j=0m−1wl,jxk+jτi,where yk is the output of the system, which is the maneuver trajectory of the target at the next moment, and x is the input of the system, which is the historical trajectory of the target.Based on the above theoretical analysis, theith order Volterra series kernel function can be expressed as follows:(25)hiz1,z2,…,zi=∑l=1Lrlai,lwl,z1wl,z2⋯wl,zi.The weightswl,n and rl and thresholds θl of the VNN model are obtained by training, and glVlk is expanded by Taylor series at the threshold. Finally, the kernel functions of each order of the model can be obtained according to equation (25).


## 5. Chaos Adaptive Genetic Algorithm

### 5.1. Fitness Function

Fitness function is used to guide the evolutionary search process of the adaptive algorithm, and the design of fitness function is closely related to the convergence speed and accuracy of the algorithm. In general, the fitness function is derived from the objective function. In order to reduce the prediction error of the VNN model, in this paper, an improved GA is proposed to optimize the model parameters, and the optimization problem belongs to the minimum optimization problem, so the fitness function of the algorithm can be defined as follows:(26)FCh=1n∑l=1ny^l−yl2,where n is the number of training samples of the algorithm; yl is the predicted value of the lth sample; y^l is the expected value of the lth sample; and FCh represents the average prediction error when the parameters decoded by chromosome Ch are used as model parameters.


### 5.2. Selection Operator

In the selection operation of the genetic algorithm, roulette and optimal individual retention strategies are used to select the mating group [22]. When the roulette method is used to perform the selection operation, the probability of each chromosome being selected is positively correlated with its function fitness. pj is the probability that the jth individual is selected, which is expressed by(27)fj=kFChj,pj=fj∑j=1Mfj,where fj is the function fitness value of the jth chromosome; M is the size of the population; k is the correlation coefficient; and FChj can be calculated based on the fitness function defined above.


### 5.3. Adaptive Crossover Operator Based on the Difference Degree of Chromosome Gene Fragments

Crossover operation is a crucial operator in traditional genetic algorithms, which directly affects the convergence speed and global convergence ability of genetic algorithms.The crossover operation of traditional genetic algorithms is based on random selection mechanism. When there is a high degree of similarity between the selected chromosomal gene fragments, according to biological genetics theory, inbreeding will greatly reduce the probability of new individuals appearing; that is, it is most likely to be a meaningless crossover operation, which causes the convergence speed of the algorithm to slow down and the phenomenon of premature emergence. In order to improve the convergence speed and global convergence ability of the algorithm, inspired by the principle of excellent gene fragment cloning, crossover operation is carried out based on the differences between gene fragments. Based on the degree of difference between chromosomal gene fragments, the crossover probability of chromosomal gene fragments is determined, and the gene fragments are selected for crossover operation so as to achieve the purpose of reduced inbreeding and invalid crossover operations.In this paper, a continuous gene segment with a lengthL in the chromosome is defined as a chromosome gene segment, and the starting position of the chromosome gene segment is defined as a cross point.During the crossover operation, the length of the chromosomal gene fragment needs to be determined. According to the length coefficient and the randomly determined length ratio of the gene fragment, the length of the gene fragment can be obtained:(28)Lfragment=u1LchromI,where I is the chromosomal gene segment length coefficient; Lchrom is the total length of the chromosomal gene segment; and u1 is the ratio of the chromosomal gene segment length, which is a random number in the range of 0,1. Chromosome gene segment coefficient is a parameter that controls the length of gene segment from the whole. At the same time, based on the overall control, the random length ratio u1 is introduced and effectively increases the diversity of crossover operations.In the crossover operation based on the difference degree of gene fragments, the gene fragment length coefficient directly affects the convergence of the algorithm. In the early stage of the algorithm, the larger value of the gene fragment length coefficient is helpful to expand the search space of the population and avoid the premature algorithm; in the latter part of the algorithm, the algorithm should focus on local search to speed up the convergence rate of the algorithm, so the value of the gene fragment length coefficient should be smaller. Based on this, an adaptive gene fragment length coefficient is constructed, and it can be expressed as follows:(29)I=Imin−ImaxtTmax2−2tTmax+1+Imax,where Imin is the minimum value of the gene fragment length coefficient; Imax is the maximum value of the gene fragment length coefficient; t is the current iteration number of the algorithm; and Tmax is the maximum iteration number.The coding length coefficientI of the gene fragment decreases nonlinearly with the number of iterations of the algorithm. The specific steps of the crossover operation based on the improvement of the difference degree of the gene fragment are as follows:Step 1: determine the maximumImax and minimum Imin length coefficients of gene fragments and generate a random number u1∈0,1. According to the number of iterations of the current algorithm and the maximum number of iterations initially set by the algorithm, the gene fragment length is calculated by equation (29).Step 2: using the roulette method to randomly select individual parentsx1=x11,x21,…,xn1andx2=x12,x22,…,xn2.Step 3: calculate the degree of difference between gene fragments of lengthLfragment corresponding to parents x1 and x2, calculate the probability of crossover by determining the degree of difference between gene fragments, and select the intersection of gene fragments according to the probability of crossover of gene fragments. The gene fragment difference degree and crossover probability can be expressed as follows:(30)Di=1L∑j=ii+1∗L−1xi1−xi22;i=0,1,…,INTnL,Pi=Di∑j=0INTn/LDj;i=0,1,…,INTnL,whereDi is the difference degree of the i-th gene fragment of the parent individual; Dj is the difference degree of the j-th gene fragment of the parent individual; xi1 and xi2 are the gene encoding at the i-th gene position of the parent individual; and INT is the rounding function.Step 4: generate a random numberu2∈0,1 and get the individual y1=y11,y21,…,yn1 and y2=y12,y22,…,yn2 of the offspring from the parent x1=x11,x21,…,xn1 and x2=x12,x22,…,xn2 through the cross, which can be expressed as follows:(31)yi1=u2xi2−xi1+xi1,i∈istart,istart+Lfragment−1,xi1,i∉istart,istart+Lfragment−1,yi2=u2xi1−xi2+xi2,i∈istart,istart+Lfragment−1,xi2,i∉istart,istart+Lfragment−1,wherei is the gene position encoded by the gene and istart is the gene position where the starting point on the gene fragment is located.


### 5.4. Mutation Operation

The mutation probabilitypm directly affects the mutation of population. The appropriate mutation of individuals can keep the population diversity and prevent the algorithm from falling into local optimum. However, if the mutation probability pm is too large, the algorithm is similar to random search and loses the genetic evolution characteristics. In this paper, the mutation probability is improved from two aspects of the genetic evolution algebra and the fitness function value of individual population, and the mutation probability pm is updated adaptively, which is computed by(32)pmi=pmax+pmax−pmint2Tmax+fi−favg2fmax−favg,fi≥favg,pmin,fi<favg,where pmax is the maximum mutation probability; pmin is the minimum mutation probability; fi is the fitness value of the ith individual; and favg is the average of individual fitness of the current population.The existing mutation operators are mainly divided into two categories: one is the mutation operator with strong local search ability [23, 24], such as Gauss mutation operator; the other is the mutation operator with strong global search ability [24], such as Cauchy mutation operator. For the optimization problem with less extreme points, mutation operators with strong local search ability should be used; for the optimization problem with more extreme points, mutation operators with strong global search ability should be used, and if mutation operators with strong local search ability are used, it is easy to fall into local optimal value. In summary, it is difficult for a single mutation operator to take into account both global search and local search capabilities. Therefore, a combination mutation operator is proposed in this paper.In this paper, three kinds of mutation operators are used: the local search ability of the first mutation operator increases with the increase in the number of iterations; the global search ability of the second mutation operator is strong; and the local search ability of the third mutation operator is strong. The combination of the three mutation operators can make the mutation operation take into account both the global and local search ability. The combined mutation operator method can be described as follows:(1)
Adaptive mutation operator:(33)Xi′=Xi+λ1∗bT−aTt,whereXi′ is the offspring individual after mutation; a and b are the upper and lower limits of the variable values; and λ1=λ11,λ12,…,λ1m, where λ11,λ12,…,λ1m is a random number in the range of [−1, 1].(2)
Cauchy mutation operator:(34)Xi′=Xi+Cauchy0,1,whereCauchy0,1 is the standard Cauchy distribution.(3)
Gaussian mutation operator:(35)Xi′=Xi+Nμ,δ2,whereμ is the mean value in the normal distribution, and generally the value is μ=Xi; δ2 is the variance of the normal distribution, and this parameter is calculated by(36)δ2=Xbest−Xi6,whereXbest is the best individual in the population.(4)
Switching method of the combination mutation operator:(37)modt=0,Adaptive mutation,1,Cauchy mutation,2,Gaussian mutation,wheremod is the function to find the remainder.The first mutation operator has a strong global search ability at the beginning of the iteration; with the increase in the number of iterations, the local search capability of the algorithm gradually increases. But, when the number of iterations increases to a certain value, it almost loses the function of mutation. When the number of iterations is not very large, the effect of this mutation operator is obvious. The second mutation operator is Cauchy mutation. Compared with the normal mutation operator, the Cauchy mutation operator will produce larger mutation step length, so it will make the algorithm have better global search ability. The third mutation operator is normal mutation. The normal mutation operator focuses on searching a local area near the original individual. The local search ability is better, but the ability to guide the individual to jump out of the local better solution is weak, which is not conducive to global convergence. Therefore, the combined mutation operator takes into account both the global exploration capability and the local search capability, which enables GA to quickly converge to the global optimal solution.


### 5.5. The IGA-VNN Model Learning Algorithm for Target Maneuver Trajectory Time Series

Firstly, based on the PSR theory, the C-C method is used to determine the optimal time delayt and time window τw, and then the optimal embedding dimension is determined by the relationship τw=m−1τ between time delay and time window. Secondly, the small data method is used to calculate the maximum Lyapunov exponents of the target maneuvering trajectory time series, the chaotic characteristics of target maneuver trajectory time series is judged, and the time series of the target maneuver trajectory is reconstructed according to the obtained embedding dimension τ and delay time m. Then, the VNN model structure is determined. According to the embedding dimension obtained by the C-C method, the number of neurons in the input layer is determined as m; according to the empirical formula [25], the number of nodes in the hidden layer is determined as L; according to the sample characteristics, the number of neurons in the output layer is determined as one, and a chaotic time series IGA-VNN prediction model with m−L−1 network structure is established. Finally, the IGA-VNN model is constructed to obtain the kernel functions of each order of the Volterra series, and the specific execution steps are as follows:Step 1: the C-C method is used to determine the time delayτ and embedding dimension m of the phase space reconstruction. Based on this, the target maneuvering trajectory time series is reconstructed and M=N−m−1τ phase space vectors are obtained.Step 2: according to Figure9, the input vector of the IGA-VNN prediction model with m−L−1 structure is XT, the output of prediction model is yk, and the weights of the hidden layer and output layer are W=wi,jL×mj=1,2,…,m and rl, respectively.Step 3: the improved genetic algorithm is used to optimize the network parameters in the IGA-VNN target maneuver trajectory time series prediction model. The specific optimization steps are as follows:Firstly, according to the characteristics of the genetic algorithm, the network parameters including the connection weights and thresholds between input layer and hidden layer and hidden layer and output layer are real coded to form chromosomes, and the population size of the genetic algorithm is determined simultaneously. Secondly, the fitness function is determined, and it can be expressed as follows:(38)F=kek+1−x˜k+1,whereek+1 is the expected output of the network and x˜k+1 is the predicted output of the network.Then, selection, crossover, and mutation operations are carried out to generate a new generation of population. Finally, the fitness value of the evolved population is calculated, and the above operations are repeated until the conditionFl+1−Fl<ε is met, and then the optimal initial weight and threshold value of the network are obtained; otherwise, the previous step is returned until the iteration end condition is met.Step 4: based on the weights and thresholds obtained in Steps 1, 2, and 3, combined with equation (10), the IGA-VNN prediction model is used to predict the target maneuver trajectory time series.Step 5: calculate the prediction errorE of the network, and it can be expressed as follows:(39)E=12∑k=1N′yk−y˜k2,whereyk is the real value and y˜k is the predicted value.In this paper, we set the maximum errorEmax of network prediction as 0.025. If the condition E<Emax is met, the calculation will be stopped, and the network weight matrix W=wl,jL×m and weight coefficient rl obtained after training will be stored. The kernel functions hil1,l2,…,li of each order in the Volterra model are calculated with polynomial coefficients ai,li=1,2,…,m; otherwise, the next step is taken.Step 6: based on the gradient descent method, the local gradientδlk and the correction of the network weight Δwl,jk of the IGA-VNN prediction model of chaotic time series are calculated, respectively. This can be calculated as follows:(40)δlk=−∂E∂ykgl′Vlk,Δwl,jk=αΔwl,jk−1+ηδlkyk,whereαΔwl,jk−1 is the introduced moving vector and 0<α<1 and η is the learning step.Step 7: correct the weights and thresholds of the network and train the prediction network again, calculate the error between the predicted output of the network and the true value of the sample, and repeat the training until the conditionE<Emax is satisfied. After the training is completed, the weight matrix of the network is stored, and the Taylor series coefficients are decomposed at the threshold θl simultaneously, and then the Taylor expansion coefficient diθl can be obtained. Because the output excitation function is in the form of polynomials, diθl=ai,l can be obtained through mathematical analysis. At this time, the kernel functions of each order of the Volterra series can be obtained by substituting the coefficients of the Taylor expansion into equation (11).


### 5.6. The Framework of the Target Maneuver Trajectory Forecast Mode

The flowchart of the proposed target maneuver trajectory prediction model in this paper is depicted in Figure10, which consists of four modules: (a) preprocessing module; (b) the IGA-VNN model learning algorithm module; (c) the whole prediction algorithm module; and (d) evaluation module. Furthermore, the detailed explanation is listed as follows:(a)
Preprocessing module: In order to determine if there is chaotic character of the target maneuver trajectory time series, the small data quantity method is used to calculate the maximum Lyapunov exponent. In addition, the C-C method is used to determine the input forms of the proposed target maneuver trajectory forecast model. Then, the normalization method is applied.(b)
The IGA-VNN model learning algorithm module: In view of the shortcomings of the genetic optimization algorithm, a strategy of the adaptive crossover mutation operator and combination mutation operator based on difference degree of gene fragments is proposed in this module. In addition, the improved genetic optimization algorithm is applied to optimize the parameter of the VNN model.(c)
The whole prediction algorithm module: combining phase space reconstruction and the Volterra adaptive filter, a novel target maneuver trajectory forecast model is proposed. The IGA-VNN model is developed in this paper to search for the optimal kernel functions of each order of the Volterra adaptive filter which is conducive to enhance the forecast accuracy.(D)
Evaluation module: four metrics which areRmse (root mean square error), Mae (mean absolute error), Mape (mean absolute percent error), and Cor (correlation coefficient) are used as evaluating indicators in this paper.Figure 10
The flowchart of the proposed target maneuver trajectory prediction model.
(a)(b)(c)(d)


## 5.1. Fitness Function

Fitness function is used to guide the evolutionary search process of the adaptive algorithm, and the design of fitness function is closely related to the convergence speed and accuracy of the algorithm. In general, the fitness function is derived from the objective function. In order to reduce the prediction error of the VNN model, in this paper, an improved GA is proposed to optimize the model parameters, and the optimization problem belongs to the minimum optimization problem, so the fitness function of the algorithm can be defined as follows:(26)FCh=1n∑l=1ny^l−yl2,where n is the number of training samples of the algorithm; yl is the predicted value of the lth sample; y^l is the expected value of the lth sample; and FCh represents the average prediction error when the parameters decoded by chromosome Ch are used as model parameters.


## 5.2. Selection Operator

In the selection operation of the genetic algorithm, roulette and optimal individual retention strategies are used to select the mating group [22]. When the roulette method is used to perform the selection operation, the probability of each chromosome being selected is positively correlated with its function fitness. pj is the probability that the jth individual is selected, which is expressed by(27)fj=kFChj,pj=fj∑j=1Mfj,where fj is the function fitness value of the jth chromosome; M is the size of the population; k is the correlation coefficient; and FChj can be calculated based on the fitness function defined above.


## 5.3. Adaptive Crossover Operator Based on the Difference Degree of Chromosome Gene Fragments

Crossover operation is a crucial operator in traditional genetic algorithms, which directly affects the convergence speed and global convergence ability of genetic algorithms.The crossover operation of traditional genetic algorithms is based on random selection mechanism. When there is a high degree of similarity between the selected chromosomal gene fragments, according to biological genetics theory, inbreeding will greatly reduce the probability of new individuals appearing; that is, it is most likely to be a meaningless crossover operation, which causes the convergence speed of the algorithm to slow down and the phenomenon of premature emergence. In order to improve the convergence speed and global convergence ability of the algorithm, inspired by the principle of excellent gene fragment cloning, crossover operation is carried out based on the differences between gene fragments. Based on the degree of difference between chromosomal gene fragments, the crossover probability of chromosomal gene fragments is determined, and the gene fragments are selected for crossover operation so as to achieve the purpose of reduced inbreeding and invalid crossover operations.In this paper, a continuous gene segment with a lengthL in the chromosome is defined as a chromosome gene segment, and the starting position of the chromosome gene segment is defined as a cross point.During the crossover operation, the length of the chromosomal gene fragment needs to be determined. According to the length coefficient and the randomly determined length ratio of the gene fragment, the length of the gene fragment can be obtained:(28)Lfragment=u1LchromI,where I is the chromosomal gene segment length coefficient; Lchrom is the total length of the chromosomal gene segment; and u1 is the ratio of the chromosomal gene segment length, which is a random number in the range of 0,1. Chromosome gene segment coefficient is a parameter that controls the length of gene segment from the whole. At the same time, based on the overall control, the random length ratio u1 is introduced and effectively increases the diversity of crossover operations.In the crossover operation based on the difference degree of gene fragments, the gene fragment length coefficient directly affects the convergence of the algorithm. In the early stage of the algorithm, the larger value of the gene fragment length coefficient is helpful to expand the search space of the population and avoid the premature algorithm; in the latter part of the algorithm, the algorithm should focus on local search to speed up the convergence rate of the algorithm, so the value of the gene fragment length coefficient should be smaller. Based on this, an adaptive gene fragment length coefficient is constructed, and it can be expressed as follows:(29)I=Imin−ImaxtTmax2−2tTmax+1+Imax,where Imin is the minimum value of the gene fragment length coefficient; Imax is the maximum value of the gene fragment length coefficient; t is the current iteration number of the algorithm; and Tmax is the maximum iteration number.The coding length coefficientI of the gene fragment decreases nonlinearly with the number of iterations of the algorithm. The specific steps of the crossover operation based on the improvement of the difference degree of the gene fragment are as follows:Step 1: determine the maximumImax and minimum Imin length coefficients of gene fragments and generate a random number u1∈0,1. According to the number of iterations of the current algorithm and the maximum number of iterations initially set by the algorithm, the gene fragment length is calculated by equation (29).Step 2: using the roulette method to randomly select individual parentsx1=x11,x21,…,xn1andx2=x12,x22,…,xn2.Step 3: calculate the degree of difference between gene fragments of lengthLfragment corresponding to parents x1 and x2, calculate the probability of crossover by determining the degree of difference between gene fragments, and select the intersection of gene fragments according to the probability of crossover of gene fragments. The gene fragment difference degree and crossover probability can be expressed as follows:(30)Di=1L∑j=ii+1∗L−1xi1−xi22;i=0,1,…,INTnL,Pi=Di∑j=0INTn/LDj;i=0,1,…,INTnL,whereDi is the difference degree of the i-th gene fragment of the parent individual; Dj is the difference degree of the j-th gene fragment of the parent individual; xi1 and xi2 are the gene encoding at the i-th gene position of the parent individual; and INT is the rounding function.Step 4: generate a random numberu2∈0,1 and get the individual y1=y11,y21,…,yn1 and y2=y12,y22,…,yn2 of the offspring from the parent x1=x11,x21,…,xn1 and x2=x12,x22,…,xn2 through the cross, which can be expressed as follows:(31)yi1=u2xi2−xi1+xi1,i∈istart,istart+Lfragment−1,xi1,i∉istart,istart+Lfragment−1,yi2=u2xi1−xi2+xi2,i∈istart,istart+Lfragment−1,xi2,i∉istart,istart+Lfragment−1,wherei is the gene position encoded by the gene and istart is the gene position where the starting point on the gene fragment is located.


## 5.4. Mutation Operation

The mutation probabilitypm directly affects the mutation of population. The appropriate mutation of individuals can keep the population diversity and prevent the algorithm from falling into local optimum. However, if the mutation probability pm is too large, the algorithm is similar to random search and loses the genetic evolution characteristics. In this paper, the mutation probability is improved from two aspects of the genetic evolution algebra and the fitness function value of individual population, and the mutation probability pm is updated adaptively, which is computed by(32)pmi=pmax+pmax−pmint2Tmax+fi−favg2fmax−favg,fi≥favg,pmin,fi<favg,where pmax is the maximum mutation probability; pmin is the minimum mutation probability; fi is the fitness value of the ith individual; and favg is the average of individual fitness of the current population.The existing mutation operators are mainly divided into two categories: one is the mutation operator with strong local search ability [23, 24], such as Gauss mutation operator; the other is the mutation operator with strong global search ability [24], such as Cauchy mutation operator. For the optimization problem with less extreme points, mutation operators with strong local search ability should be used; for the optimization problem with more extreme points, mutation operators with strong global search ability should be used, and if mutation operators with strong local search ability are used, it is easy to fall into local optimal value. In summary, it is difficult for a single mutation operator to take into account both global search and local search capabilities. Therefore, a combination mutation operator is proposed in this paper.In this paper, three kinds of mutation operators are used: the local search ability of the first mutation operator increases with the increase in the number of iterations; the global search ability of the second mutation operator is strong; and the local search ability of the third mutation operator is strong. The combination of the three mutation operators can make the mutation operation take into account both the global and local search ability. The combined mutation operator method can be described as follows:(1)
Adaptive mutation operator:(33)Xi′=Xi+λ1∗bT−aTt,whereXi′ is the offspring individual after mutation; a and b are the upper and lower limits of the variable values; and λ1=λ11,λ12,…,λ1m, where λ11,λ12,…,λ1m is a random number in the range of [−1, 1].(2)
Cauchy mutation operator:(34)Xi′=Xi+Cauchy0,1,whereCauchy0,1 is the standard Cauchy distribution.(3)
Gaussian mutation operator:(35)Xi′=Xi+Nμ,δ2,whereμ is the mean value in the normal distribution, and generally the value is μ=Xi; δ2 is the variance of the normal distribution, and this parameter is calculated by(36)δ2=Xbest−Xi6,whereXbest is the best individual in the population.(4)
Switching method of the combination mutation operator:(37)modt=0,Adaptive mutation,1,Cauchy mutation,2,Gaussian mutation,wheremod is the function to find the remainder.The first mutation operator has a strong global search ability at the beginning of the iteration; with the increase in the number of iterations, the local search capability of the algorithm gradually increases. But, when the number of iterations increases to a certain value, it almost loses the function of mutation. When the number of iterations is not very large, the effect of this mutation operator is obvious. The second mutation operator is Cauchy mutation. Compared with the normal mutation operator, the Cauchy mutation operator will produce larger mutation step length, so it will make the algorithm have better global search ability. The third mutation operator is normal mutation. The normal mutation operator focuses on searching a local area near the original individual. The local search ability is better, but the ability to guide the individual to jump out of the local better solution is weak, which is not conducive to global convergence. Therefore, the combined mutation operator takes into account both the global exploration capability and the local search capability, which enables GA to quickly converge to the global optimal solution.


## 5.5. The IGA-VNN Model Learning Algorithm for Target Maneuver Trajectory Time Series

Firstly, based on the PSR theory, the C-C method is used to determine the optimal time delayt and time window τw, and then the optimal embedding dimension is determined by the relationship τw=m−1τ between time delay and time window. Secondly, the small data method is used to calculate the maximum Lyapunov exponents of the target maneuvering trajectory time series, the chaotic characteristics of target maneuver trajectory time series is judged, and the time series of the target maneuver trajectory is reconstructed according to the obtained embedding dimension τ and delay time m. Then, the VNN model structure is determined. According to the embedding dimension obtained by the C-C method, the number of neurons in the input layer is determined as m; according to the empirical formula [25], the number of nodes in the hidden layer is determined as L; according to the sample characteristics, the number of neurons in the output layer is determined as one, and a chaotic time series IGA-VNN prediction model with m−L−1 network structure is established. Finally, the IGA-VNN model is constructed to obtain the kernel functions of each order of the Volterra series, and the specific execution steps are as follows:Step 1: the C-C method is used to determine the time delayτ and embedding dimension m of the phase space reconstruction. Based on this, the target maneuvering trajectory time series is reconstructed and M=N−m−1τ phase space vectors are obtained.Step 2: according to Figure9, the input vector of the IGA-VNN prediction model with m−L−1 structure is XT, the output of prediction model is yk, and the weights of the hidden layer and output layer are W=wi,jL×mj=1,2,…,m and rl, respectively.Step 3: the improved genetic algorithm is used to optimize the network parameters in the IGA-VNN target maneuver trajectory time series prediction model. The specific optimization steps are as follows:Firstly, according to the characteristics of the genetic algorithm, the network parameters including the connection weights and thresholds between input layer and hidden layer and hidden layer and output layer are real coded to form chromosomes, and the population size of the genetic algorithm is determined simultaneously. Secondly, the fitness function is determined, and it can be expressed as follows:(38)F=kek+1−x˜k+1,whereek+1 is the expected output of the network and x˜k+1 is the predicted output of the network.Then, selection, crossover, and mutation operations are carried out to generate a new generation of population. Finally, the fitness value of the evolved population is calculated, and the above operations are repeated until the conditionFl+1−Fl<ε is met, and then the optimal initial weight and threshold value of the network are obtained; otherwise, the previous step is returned until the iteration end condition is met.Step 4: based on the weights and thresholds obtained in Steps 1, 2, and 3, combined with equation (10), the IGA-VNN prediction model is used to predict the target maneuver trajectory time series.Step 5: calculate the prediction errorE of the network, and it can be expressed as follows:(39)E=12∑k=1N′yk−y˜k2,whereyk is the real value and y˜k is the predicted value.In this paper, we set the maximum errorEmax of network prediction as 0.025. If the condition E<Emax is met, the calculation will be stopped, and the network weight matrix W=wl,jL×m and weight coefficient rl obtained after training will be stored. The kernel functions hil1,l2,…,li of each order in the Volterra model are calculated with polynomial coefficients ai,li=1,2,…,m; otherwise, the next step is taken.Step 6: based on the gradient descent method, the local gradientδlk and the correction of the network weight Δwl,jk of the IGA-VNN prediction model of chaotic time series are calculated, respectively. This can be calculated as follows:(40)δlk=−∂E∂ykgl′Vlk,Δwl,jk=αΔwl,jk−1+ηδlkyk,whereαΔwl,jk−1 is the introduced moving vector and 0<α<1 and η is the learning step.Step 7: correct the weights and thresholds of the network and train the prediction network again, calculate the error between the predicted output of the network and the true value of the sample, and repeat the training until the conditionE<Emax is satisfied. After the training is completed, the weight matrix of the network is stored, and the Taylor series coefficients are decomposed at the threshold θl simultaneously, and then the Taylor expansion coefficient diθl can be obtained. Because the output excitation function is in the form of polynomials, diθl=ai,l can be obtained through mathematical analysis. At this time, the kernel functions of each order of the Volterra series can be obtained by substituting the coefficients of the Taylor expansion into equation (11).


## 5.6. The Framework of the Target Maneuver Trajectory Forecast Mode

The flowchart of the proposed target maneuver trajectory prediction model in this paper is depicted in Figure10, which consists of four modules: (a) preprocessing module; (b) the IGA-VNN model learning algorithm module; (c) the whole prediction algorithm module; and (d) evaluation module. Furthermore, the detailed explanation is listed as follows:(a)
Preprocessing module: In order to determine if there is chaotic character of the target maneuver trajectory time series, the small data quantity method is used to calculate the maximum Lyapunov exponent. In addition, the C-C method is used to determine the input forms of the proposed target maneuver trajectory forecast model. Then, the normalization method is applied.(b)
The IGA-VNN model learning algorithm module: In view of the shortcomings of the genetic optimization algorithm, a strategy of the adaptive crossover mutation operator and combination mutation operator based on difference degree of gene fragments is proposed in this module. In addition, the improved genetic optimization algorithm is applied to optimize the parameter of the VNN model.(c)
The whole prediction algorithm module: combining phase space reconstruction and the Volterra adaptive filter, a novel target maneuver trajectory forecast model is proposed. The IGA-VNN model is developed in this paper to search for the optimal kernel functions of each order of the Volterra adaptive filter which is conducive to enhance the forecast accuracy.(D)
Evaluation module: four metrics which areRmse (root mean square error), Mae (mean absolute error), Mape (mean absolute percent error), and Cor (correlation coefficient) are used as evaluating indicators in this paper.Figure 10
The flowchart of the proposed target maneuver trajectory prediction model.
(a)(b)(c)(d)


## 6. Simulation Experiment

### 6.1. Experimental Data Source and Simulation Scenarios

In order to verify the effectiveness and accuracy of the target maneuver trajectory prediction model, in this paper, 4000 consecutive sets of flight confrontation training data are intercepted from the air combat simulation system, and the data contain the target’s three-dimensional coordinates and pitch angle and yaw angle data. A complete air combat confrontation is shown in Figure11, the initial confrontation scenario is shown in Table 1, and the three-dimensional coordinate components of the target flight trajectory are shown in Figure 12.Figure 11
A complete air combat trajectory map.Table 1
The scenarios for maneuver trajectories.
Scenario [1]Initial conditionVelocity (m/s)Pitch (rad)Roll (rad)Yaw (rad)Distance (m)Fighter(0, 0, 8569)214.19−0.78372.59411.482564411.44Target(-64259, -4360, 9327)248.760.00170.47441.4971Figure 12
Sample 3D data graph.


### 6.2. Experimental Data Processing

When using the IGA-VNN prediction model to predict the chaotic target maneuver trajectory time series, in order to obtain better prediction performance, it is necessary to normalize and inverse normalize the historical data [26]. Assuming that X=x1,x2,…,xN, X∈RN is the original data of the target maneuver trajectory time series, and the normalization processing can be described as follows:(41)xt′=xt−xminxmax−xmin,t=1,2,…,N,where N is the length of the target maneuvering trajectory time series; xt is the original value of the historical trajectory; xt′ is the normalized value; xmin is the minimum values in the historical trajectory; and xmax is the maximum values in the historical trajectory.Assuming that the output of the IGA-VNN target maneuvering trajectory prediction model isY=y1,y2,…,yNandY∈0,1, the antinormalization processing can be expressed as follows:(42)yt′=ytymax−ymin+ymin,t=1,2,…,N,where yt is the model output; yt′ is the output after inverse normalization; ymin is the minimum values in the model output; and ymax is the maximum values in the model output.


### 6.3. Simulation Experiment Settings

In order to evaluate the performance of the target maneuvering trajectory prediction model, in this paper, the concept of multistep prediction is introduced [27]. Multistep prediction has a high use value in the aspects of perception situation and threat assessment. Multistep prediction of target maneuvering trajectory means that, based on the existing target historical trajectory data, the model can not only predict the target trajectory at the next moment but also predict the target trajectory at multiple moments in the future. The structure of the multistep prediction model is shown in Figure 13.Figure 13
Multistep prediction model.Based on the one-step prediction model, the single-step prediction resultxk+1 is used to update the input vector p of the multistep prediction model, and p=xk+1xk⋯xk−m+1T, and the prediction value at k+2 time can be obtained by calling the Volterra series prediction model. In this way, the input vector is updated repeatedly, and the Volterra one-step prediction model is called recursively, so the output value at any time in the future can be predicted theoretically.In this paper, four prediction models are used to predict the target maneuvering trajectory time series in one step, two steps, four steps, and six steps, respectively. The specific steps are as follows.Firstly,N−m−1τ phase space vectors are obtained through phase space reconstruction, and M vectors are selected as training samples from the phase space vectors. After the training, the network parameters obtained by the training are stored. Then, the chaotic time series IGA-VNN prediction model, GA-VNN prediction model, RBF prediction model, second-order Volterra functional prediction model, and BP neural network prediction model are used to make 1, 2, 4, and 6 step predictions on the target maneuver trajectory, respectively.Finally, in order to quantify and compare the prediction accuracy and effectiveness of the four prediction algorithms, this paper uses root mean square errorRmse, mean absolute error Mae, mean absolute percentage error Mape, and correlation coefficient Cor as evaluation indicators:(43)Rmse=1n∑t=1ny˜k−yk2,Mae=∑t=1ny˜k−ykn,Mape=1n∑t=1ny˜k−yky˜k,Cor=∑t=1ny˜k−y˜¯⋅yk−y¯n⋅σy˜⋅σy,where n is the number of the sample, y˜ is the actual value, that is, the true position of the target at the next moment; y is the predicted value, that is, the predicted position of the target at the next moment; σy is the standard deviation of y; y˜¯ is the average value of y˜; and σy˜ is the standard deviation of y˜.


### 6.4. Analysis of Simulation Experiment Results

#### 6.4.1. Predictive Performance Comparison of the Six Models

Based on PSR, phase space reconstruction is performed on the target maneuver trajectory samples to obtain the input and output vectors of the prediction model. In this paper, IGA-VNN prediction model, GA-VNN prediction model, RBF prediction model, second-order Volterra functional prediction model, and BP neural network prediction model are used to predict the three-dimensional coordinates of the target’s maneuvering trajectory in one step, two steps, four steps, and six steps, respectively. The prediction results are shown in Figures14–16, and the comparison results of target maneuvering trajectory prediction are shown in Figures 17–20, and the comparison results of target maneuvering trajectory prediction error are shown in Figures 21–23.Figure 14
The curves of (a) 1-step ,(b) 2-step , (c) 4-step, and (d) 6-step X-coordinate prediction results.
(a)(b)(c)(d)Figure 15
The curves of (a) 1-step, (b) 2-step, (c) 4-step, and (d) 6-stepY-coordinate prediction results.
(a)(b)(c)(d)Figure 16
The curves of (a) 1-step, (b) 2-step, (c) 4-step, and (d) 6-stepZ-coordinate prediction results.
(a)(b)(c)(d)Figure 17
Comparison of one-step prediction results of target maneuvering trajectory.Figure 18
Comparison of two-step prediction results of target maneuvering trajectory.Figure 19
Comparison of four-step prediction results of target maneuvering trajectory.Figure 20
Comparison of six-step prediction results of target maneuvering trajectory.Figure 21
The curves of (a) 1-step, (b) 2-step, (c) 4-step, and (d) 6-stepX-coordinate prediction error.
(a)(b)(c)(d)Figure 22
The curves of (a) 1-step, (b) 2-step, (c) 4-step, and (d) 6-stepY-coordinate prediction error.
(a)(b)(c)(d)Figure 23
The curves of (a) 1-step, (b) 2-step, (c) 4-step, and (d) 6-stepZ-coordinate prediction error.
(a)(b)(c)(d)In order to more intuitively see the prediction accuracy of each model, the forecasting evaluation indicators of each model are calculated as shown in Tables2–4.Table 2
TargetX-coordinate multistep prediction performance comparison of four different prediction methods.
Prediction stepAlgorithmRMSEMAEMAPECOR1BP31.20545.25150.00110.9987IGA-VNN25.05544.25947.3741∗e−040.9984GA-VNN46.61616.29590.00160.9983Volterra88.09499.24940.00350.9985RBF65.58837.40460.00230.99872BP58.20297.12550.00210.9972IGA-VNN43.11025.66940.00130.9969GA-VNN51.51076.65890.00180.9961Volterra103.47909.41920.00360.9920RBF103.93748.98000.00330.98884BP100.80759.01220.00320.9983IGA-VNN64.82787.25190.00210.9952GA-VNN130.017010.60180.00460.9935Volterra227.380712.13730.00590.9816RBF171.446611.39850.00530.96366BP151.814211.97050.00580.9978IGA-VNN115.11469.73680.00380.9866Volterra200.244013.05920.00700.9980RBF238.182613.93210.00790.9531Table 3
TargetY-coordinate multistep prediction performance comparison of four different prediction methods.
Prediction stepAlgorithmRMSEMAEMAPECOR1BP48.02173.98510.00690.9965IGA-VNN15.65343.47660.00220.9987GA-VNN30.76535.22960.00470.9970Volterra40.57595.92180.00610.9971RBF48.87126.34570.00730.99742BP99.36028.65270.01390.9913IGA-VNN23.05824.26620.00330.9980GA-VNN57.35266.66140.00780.9932Volterra59.88367.10300.00890.9948RBF70.91277.88410.01090.98994BP123.700110.44900.01880.9851IGA-VNN35.48395.58170.00560.9976GA-VNN119.90399.97380.01750.9752Volterra148.417310.90430.02210.9904RBF119.68359.23100.01590.97986BP161.138611.65050.02360.9685IGA-VNN44.42396.41550.00730.9952GA-VNN158.597311.49050.02420.9846Volterra189.801212.73700.02950.9880RBF159.791610.98160.02030.9618Table 4
TargetZ-coordinate multistep prediction performance comparison of four different prediction methods.
Prediction stepAlgorithmRMSEMAEMAPECOR1BP3.81681.80708.7264∗e−040.9959IGA-VNN1.49141.09433.2070∗e−040.9984GA-VNN3.52071.74378.1077∗e−040.9958Volterra2.06361.32944.7240∗e−040.9937RBF5.04112.06800.00110.97202BP7.68772.68370.00190.9881IGA-VNN4.98932.12060.00120.9968GA-VNN4.97182.20860.00130.9957Volterra6.26892.40250.00150.9479RBF7.99692.32840.00180.95724BP11.98213.25660.00280.8768IGA-VNN6.47402.38600.00150.9757GA-VNN8.04132.78890.00210.9978Volterra8.18162.54050.00170.9316RBF14.08473.3960.00320.90446BP14.12093.27620.00290.8273IGA-VNN9.84692.88520.00220.9677GA-VNN10.98082.97730.00240.9585Volterra10.13712.94630.00230.8627RBF20.96114.43270.00520.8633


#### 6.4.2. Results Analysis

By comparing and analyzing the prediction performance of different models in Tables2–4, the following conclusions can be drawn:(1)
When performing multistep prediction on the target maneuver trajectory time series, the prediction performance of the chaotic time series IGA-VNN model is the best with the same prediction steps. Because the Volterra functional model can represent most nonlinear systems, it can accurately predict chaotic time series. However, the solution of high-order Volterra kernel function is the bottleneck of its application. In this paper, we establish the IGA-VNN model to determine the kernel functions of Volterra series, so that Volterra series can predict the target maneuver trajectory well.(2)
RBF neural network prediction model, second-order Volterra functional prediction model, and BP neural network prediction model show different prediction performance in different prediction steps and different coordinate directions of the target, which shows that the adaptability and robustness of these three prediction algorithms are poor, and the performance of the IGA-VNN prediction model is superior than that of these three prediction algorithms in this respect. In this paper, the phase space reconstruction theory is adopted to reconstruct the time series of target maneuver trajectory, which makes full use of the information contained in the time series of target maneuver trajectory so that the prediction accuracy is improved.(3)
Combined with the prediction error comparison charts in Figures14–16, it can be seen that the same prediction algorithm’s single-step prediction results for the three coordinates of the target maneuver trajectory are better than the 2-step prediction results; the 2-step prediction results are better than the 4-step prediction results; and the 4-step prediction results are better than 6-step prediction results. Therefore, it can be seen that for the same prediction model, the prediction performance of the model will also decline with the increase in prediction steps. In the multistep prediction of the target maneuver trajectory, as the number of prediction steps increases, the cumulative error of the prediction will continue to increase, so the prediction accuracy will also decrease.By comparing and analyzing the differences between the four models in the multistep prediction principle of target maneuver trajectory time series, the root cause of these differences can be further analyzed. BP neural network is based on sample data to train the network so as to realize the prediction of the future state. The important factors that affect the precision and prediction accuracy of the BP neural network are the subjectivity and blindness of its variable selection, which leads to the addition of extra irrelevant variables or ignoring important variables. In addition, BP neural network is sensitive to initial value and easy to fall into local minimum. The Volterra functional model essentially uses the chaotic characteristics of target maneuvering trajectory time series to predict the nonlinear system by training and fitting the trajectory of the chaotic attractor. The model is implemented by the multiplication-coupled form of the second-order Volterra series, which avoids the difficulty of solving the kernel function caused by the higher-order Volterra series and improves the prediction accuracy of the model. In addition, the Volterra filter uses the linear adaptive algorithm to adjust the filter parameters, and this dynamic property can improve the adaptability and accuracy of prediction and reduce the training time.The IGA-VNN prediction model based on chaos theory combines the discrete Volterra functional model, BP neural network, and the chaotic characteristics of the target maneuvering trajectory time series and uses the chaotic characteristics, causality, and memory function [25] of the target maneuvering trajectory time series to determine the truncation order and truncation term number of the Volterra functional model so as to establish an accurate Volterra functional model. At the same time, the equivalence of BP neural network and Volterra functional model is used to solve the problem of high-order kernel function of the Volterra functional model, and the global searching ability of the improved genetic algorithm is used to avoid falling into local minimum when optimizing the VNN model. In conclusion, the IGA-VNN prediction model based on chaos theory reasonably uses the advantages of the Volterra functional model and BP neural network and effectively solves the problem of subjectivity in BP neural network modeling and the difficulty in the Volterra functional model solving high-order kernel function, making the performance of the algorithm in multistep prediction better than that of the single BP neural network and Volterra functional model.


## 6.1. Experimental Data Source and Simulation Scenarios

In order to verify the effectiveness and accuracy of the target maneuver trajectory prediction model, in this paper, 4000 consecutive sets of flight confrontation training data are intercepted from the air combat simulation system, and the data contain the target’s three-dimensional coordinates and pitch angle and yaw angle data. A complete air combat confrontation is shown in Figure11, the initial confrontation scenario is shown in Table 1, and the three-dimensional coordinate components of the target flight trajectory are shown in Figure 12.Figure 11
A complete air combat trajectory map.Table 1
The scenarios for maneuver trajectories.
Scenario [1]Initial conditionVelocity (m/s)Pitch (rad)Roll (rad)Yaw (rad)Distance (m)Fighter(0, 0, 8569)214.19−0.78372.59411.482564411.44Target(-64259, -4360, 9327)248.760.00170.47441.4971Figure 12
Sample 3D data graph.


## 6.2. Experimental Data Processing

When using the IGA-VNN prediction model to predict the chaotic target maneuver trajectory time series, in order to obtain better prediction performance, it is necessary to normalize and inverse normalize the historical data [26]. Assuming that X=x1,x2,…,xN, X∈RN is the original data of the target maneuver trajectory time series, and the normalization processing can be described as follows:(41)xt′=xt−xminxmax−xmin,t=1,2,…,N,where N is the length of the target maneuvering trajectory time series; xt is the original value of the historical trajectory; xt′ is the normalized value; xmin is the minimum values in the historical trajectory; and xmax is the maximum values in the historical trajectory.Assuming that the output of the IGA-VNN target maneuvering trajectory prediction model isY=y1,y2,…,yNandY∈0,1, the antinormalization processing can be expressed as follows:(42)yt′=ytymax−ymin+ymin,t=1,2,…,N,where yt is the model output; yt′ is the output after inverse normalization; ymin is the minimum values in the model output; and ymax is the maximum values in the model output.


## 6.3. Simulation Experiment Settings

In order to evaluate the performance of the target maneuvering trajectory prediction model, in this paper, the concept of multistep prediction is introduced [27]. Multistep prediction has a high use value in the aspects of perception situation and threat assessment. Multistep prediction of target maneuvering trajectory means that, based on the existing target historical trajectory data, the model can not only predict the target trajectory at the next moment but also predict the target trajectory at multiple moments in the future. The structure of the multistep prediction model is shown in Figure 13.Figure 13
Multistep prediction model.Based on the one-step prediction model, the single-step prediction resultxk+1 is used to update the input vector p of the multistep prediction model, and p=xk+1xk⋯xk−m+1T, and the prediction value at k+2 time can be obtained by calling the Volterra series prediction model. In this way, the input vector is updated repeatedly, and the Volterra one-step prediction model is called recursively, so the output value at any time in the future can be predicted theoretically.In this paper, four prediction models are used to predict the target maneuvering trajectory time series in one step, two steps, four steps, and six steps, respectively. The specific steps are as follows.Firstly,N−m−1τ phase space vectors are obtained through phase space reconstruction, and M vectors are selected as training samples from the phase space vectors. After the training, the network parameters obtained by the training are stored. Then, the chaotic time series IGA-VNN prediction model, GA-VNN prediction model, RBF prediction model, second-order Volterra functional prediction model, and BP neural network prediction model are used to make 1, 2, 4, and 6 step predictions on the target maneuver trajectory, respectively.Finally, in order to quantify and compare the prediction accuracy and effectiveness of the four prediction algorithms, this paper uses root mean square errorRmse, mean absolute error Mae, mean absolute percentage error Mape, and correlation coefficient Cor as evaluation indicators:(43)Rmse=1n∑t=1ny˜k−yk2,Mae=∑t=1ny˜k−ykn,Mape=1n∑t=1ny˜k−yky˜k,Cor=∑t=1ny˜k−y˜¯⋅yk−y¯n⋅σy˜⋅σy,where n is the number of the sample, y˜ is the actual value, that is, the true position of the target at the next moment; y is the predicted value, that is, the predicted position of the target at the next moment; σy is the standard deviation of y; y˜¯ is the average value of y˜; and σy˜ is the standard deviation of y˜.


## 6.4. Analysis of Simulation Experiment Results

### 6.4.1. Predictive Performance Comparison of the Six Models

Based on PSR, phase space reconstruction is performed on the target maneuver trajectory samples to obtain the input and output vectors of the prediction model. In this paper, IGA-VNN prediction model, GA-VNN prediction model, RBF prediction model, second-order Volterra functional prediction model, and BP neural network prediction model are used to predict the three-dimensional coordinates of the target’s maneuvering trajectory in one step, two steps, four steps, and six steps, respectively. The prediction results are shown in Figures14–16, and the comparison results of target maneuvering trajectory prediction are shown in Figures 17–20, and the comparison results of target maneuvering trajectory prediction error are shown in Figures 21–23.Figure 14
The curves of (a) 1-step ,(b) 2-step , (c) 4-step, and (d) 6-step X-coordinate prediction results.
(a)(b)(c)(d)Figure 15
The curves of (a) 1-step, (b) 2-step, (c) 4-step, and (d) 6-stepY-coordinate prediction results.
(a)(b)(c)(d)Figure 16
The curves of (a) 1-step, (b) 2-step, (c) 4-step, and (d) 6-stepZ-coordinate prediction results.
(a)(b)(c)(d)Figure 17
Comparison of one-step prediction results of target maneuvering trajectory.Figure 18
Comparison of two-step prediction results of target maneuvering trajectory.Figure 19
Comparison of four-step prediction results of target maneuvering trajectory.Figure 20
Comparison of six-step prediction results of target maneuvering trajectory.Figure 21
The curves of (a) 1-step, (b) 2-step, (c) 4-step, and (d) 6-stepX-coordinate prediction error.
(a)(b)(c)(d)Figure 22
The curves of (a) 1-step, (b) 2-step, (c) 4-step, and (d) 6-stepY-coordinate prediction error.
(a)(b)(c)(d)Figure 23
The curves of (a) 1-step, (b) 2-step, (c) 4-step, and (d) 6-stepZ-coordinate prediction error.
(a)(b)(c)(d)In order to more intuitively see the prediction accuracy of each model, the forecasting evaluation indicators of each model are calculated as shown in Tables2–4.Table 2
TargetX-coordinate multistep prediction performance comparison of four different prediction methods.
Prediction stepAlgorithmRMSEMAEMAPECOR1BP31.20545.25150.00110.9987IGA-VNN25.05544.25947.3741∗e−040.9984GA-VNN46.61616.29590.00160.9983Volterra88.09499.24940.00350.9985RBF65.58837.40460.00230.99872BP58.20297.12550.00210.9972IGA-VNN43.11025.66940.00130.9969GA-VNN51.51076.65890.00180.9961Volterra103.47909.41920.00360.9920RBF103.93748.98000.00330.98884BP100.80759.01220.00320.9983IGA-VNN64.82787.25190.00210.9952GA-VNN130.017010.60180.00460.9935Volterra227.380712.13730.00590.9816RBF171.446611.39850.00530.96366BP151.814211.97050.00580.9978IGA-VNN115.11469.73680.00380.9866Volterra200.244013.05920.00700.9980RBF238.182613.93210.00790.9531Table 3
TargetY-coordinate multistep prediction performance comparison of four different prediction methods.
Prediction stepAlgorithmRMSEMAEMAPECOR1BP48.02173.98510.00690.9965IGA-VNN15.65343.47660.00220.9987GA-VNN30.76535.22960.00470.9970Volterra40.57595.92180.00610.9971RBF48.87126.34570.00730.99742BP99.36028.65270.01390.9913IGA-VNN23.05824.26620.00330.9980GA-VNN57.35266.66140.00780.9932Volterra59.88367.10300.00890.9948RBF70.91277.88410.01090.98994BP123.700110.44900.01880.9851IGA-VNN35.48395.58170.00560.9976GA-VNN119.90399.97380.01750.9752Volterra148.417310.90430.02210.9904RBF119.68359.23100.01590.97986BP161.138611.65050.02360.9685IGA-VNN44.42396.41550.00730.9952GA-VNN158.597311.49050.02420.9846Volterra189.801212.73700.02950.9880RBF159.791610.98160.02030.9618Table 4
TargetZ-coordinate multistep prediction performance comparison of four different prediction methods.
Prediction stepAlgorithmRMSEMAEMAPECOR1BP3.81681.80708.7264∗e−040.9959IGA-VNN1.49141.09433.2070∗e−040.9984GA-VNN3.52071.74378.1077∗e−040.9958Volterra2.06361.32944.7240∗e−040.9937RBF5.04112.06800.00110.97202BP7.68772.68370.00190.9881IGA-VNN4.98932.12060.00120.9968GA-VNN4.97182.20860.00130.9957Volterra6.26892.40250.00150.9479RBF7.99692.32840.00180.95724BP11.98213.25660.00280.8768IGA-VNN6.47402.38600.00150.9757GA-VNN8.04132.78890.00210.9978Volterra8.18162.54050.00170.9316RBF14.08473.3960.00320.90446BP14.12093.27620.00290.8273IGA-VNN9.84692.88520.00220.9677GA-VNN10.98082.97730.00240.9585Volterra10.13712.94630.00230.8627RBF20.96114.43270.00520.8633


### 6.4.2. Results Analysis

By comparing and analyzing the prediction performance of different models in Tables2–4, the following conclusions can be drawn:(1)
When performing multistep prediction on the target maneuver trajectory time series, the prediction performance of the chaotic time series IGA-VNN model is the best with the same prediction steps. Because the Volterra functional model can represent most nonlinear systems, it can accurately predict chaotic time series. However, the solution of high-order Volterra kernel function is the bottleneck of its application. In this paper, we establish the IGA-VNN model to determine the kernel functions of Volterra series, so that Volterra series can predict the target maneuver trajectory well.(2)
RBF neural network prediction model, second-order Volterra functional prediction model, and BP neural network prediction model show different prediction performance in different prediction steps and different coordinate directions of the target, which shows that the adaptability and robustness of these three prediction algorithms are poor, and the performance of the IGA-VNN prediction model is superior than that of these three prediction algorithms in this respect. In this paper, the phase space reconstruction theory is adopted to reconstruct the time series of target maneuver trajectory, which makes full use of the information contained in the time series of target maneuver trajectory so that the prediction accuracy is improved.(3)
Combined with the prediction error comparison charts in Figures14–16, it can be seen that the same prediction algorithm’s single-step prediction results for the three coordinates of the target maneuver trajectory are better than the 2-step prediction results; the 2-step prediction results are better than the 4-step prediction results; and the 4-step prediction results are better than 6-step prediction results. Therefore, it can be seen that for the same prediction model, the prediction performance of the model will also decline with the increase in prediction steps. In the multistep prediction of the target maneuver trajectory, as the number of prediction steps increases, the cumulative error of the prediction will continue to increase, so the prediction accuracy will also decrease.By comparing and analyzing the differences between the four models in the multistep prediction principle of target maneuver trajectory time series, the root cause of these differences can be further analyzed. BP neural network is based on sample data to train the network so as to realize the prediction of the future state. The important factors that affect the precision and prediction accuracy of the BP neural network are the subjectivity and blindness of its variable selection, which leads to the addition of extra irrelevant variables or ignoring important variables. In addition, BP neural network is sensitive to initial value and easy to fall into local minimum. The Volterra functional model essentially uses the chaotic characteristics of target maneuvering trajectory time series to predict the nonlinear system by training and fitting the trajectory of the chaotic attractor. The model is implemented by the multiplication-coupled form of the second-order Volterra series, which avoids the difficulty of solving the kernel function caused by the higher-order Volterra series and improves the prediction accuracy of the model. In addition, the Volterra filter uses the linear adaptive algorithm to adjust the filter parameters, and this dynamic property can improve the adaptability and accuracy of prediction and reduce the training time.The IGA-VNN prediction model based on chaos theory combines the discrete Volterra functional model, BP neural network, and the chaotic characteristics of the target maneuvering trajectory time series and uses the chaotic characteristics, causality, and memory function [25] of the target maneuvering trajectory time series to determine the truncation order and truncation term number of the Volterra functional model so as to establish an accurate Volterra functional model. At the same time, the equivalence of BP neural network and Volterra functional model is used to solve the problem of high-order kernel function of the Volterra functional model, and the global searching ability of the improved genetic algorithm is used to avoid falling into local minimum when optimizing the VNN model. In conclusion, the IGA-VNN prediction model based on chaos theory reasonably uses the advantages of the Volterra functional model and BP neural network and effectively solves the problem of subjectivity in BP neural network modeling and the difficulty in the Volterra functional model solving high-order kernel function, making the performance of the algorithm in multistep prediction better than that of the single BP neural network and Volterra functional model.


## 6.4.1. Predictive Performance Comparison of the Six Models

Based on PSR, phase space reconstruction is performed on the target maneuver trajectory samples to obtain the input and output vectors of the prediction model. In this paper, IGA-VNN prediction model, GA-VNN prediction model, RBF prediction model, second-order Volterra functional prediction model, and BP neural network prediction model are used to predict the three-dimensional coordinates of the target’s maneuvering trajectory in one step, two steps, four steps, and six steps, respectively. The prediction results are shown in Figures14–16, and the comparison results of target maneuvering trajectory prediction are shown in Figures 17–20, and the comparison results of target maneuvering trajectory prediction error are shown in Figures 21–23.Figure 14
The curves of (a) 1-step ,(b) 2-step , (c) 4-step, and (d) 6-step X-coordinate prediction results.
(a)(b)(c)(d)Figure 15
The curves of (a) 1-step, (b) 2-step, (c) 4-step, and (d) 6-stepY-coordinate prediction results.
(a)(b)(c)(d)Figure 16
The curves of (a) 1-step, (b) 2-step, (c) 4-step, and (d) 6-stepZ-coordinate prediction results.
(a)(b)(c)(d)Figure 17
Comparison of one-step prediction results of target maneuvering trajectory.Figure 18
Comparison of two-step prediction results of target maneuvering trajectory.Figure 19
Comparison of four-step prediction results of target maneuvering trajectory.Figure 20
Comparison of six-step prediction results of target maneuvering trajectory.Figure 21
The curves of (a) 1-step, (b) 2-step, (c) 4-step, and (d) 6-stepX-coordinate prediction error.
(a)(b)(c)(d)Figure 22
The curves of (a) 1-step, (b) 2-step, (c) 4-step, and (d) 6-stepY-coordinate prediction error.
(a)(b)(c)(d)Figure 23
The curves of (a) 1-step, (b) 2-step, (c) 4-step, and (d) 6-stepZ-coordinate prediction error.
(a)(b)(c)(d)In order to more intuitively see the prediction accuracy of each model, the forecasting evaluation indicators of each model are calculated as shown in Tables2–4.Table 2
TargetX-coordinate multistep prediction performance comparison of four different prediction methods.
Prediction stepAlgorithmRMSEMAEMAPECOR1BP31.20545.25150.00110.9987IGA-VNN25.05544.25947.3741∗e−040.9984GA-VNN46.61616.29590.00160.9983Volterra88.09499.24940.00350.9985RBF65.58837.40460.00230.99872BP58.20297.12550.00210.9972IGA-VNN43.11025.66940.00130.9969GA-VNN51.51076.65890.00180.9961Volterra103.47909.41920.00360.9920RBF103.93748.98000.00330.98884BP100.80759.01220.00320.9983IGA-VNN64.82787.25190.00210.9952GA-VNN130.017010.60180.00460.9935Volterra227.380712.13730.00590.9816RBF171.446611.39850.00530.96366BP151.814211.97050.00580.9978IGA-VNN115.11469.73680.00380.9866Volterra200.244013.05920.00700.9980RBF238.182613.93210.00790.9531Table 3
TargetY-coordinate multistep prediction performance comparison of four different prediction methods.
Prediction stepAlgorithmRMSEMAEMAPECOR1BP48.02173.98510.00690.9965IGA-VNN15.65343.47660.00220.9987GA-VNN30.76535.22960.00470.9970Volterra40.57595.92180.00610.9971RBF48.87126.34570.00730.99742BP99.36028.65270.01390.9913IGA-VNN23.05824.26620.00330.9980GA-VNN57.35266.66140.00780.9932Volterra59.88367.10300.00890.9948RBF70.91277.88410.01090.98994BP123.700110.44900.01880.9851IGA-VNN35.48395.58170.00560.9976GA-VNN119.90399.97380.01750.9752Volterra148.417310.90430.02210.9904RBF119.68359.23100.01590.97986BP161.138611.65050.02360.9685IGA-VNN44.42396.41550.00730.9952GA-VNN158.597311.49050.02420.9846Volterra189.801212.73700.02950.9880RBF159.791610.98160.02030.9618Table 4
TargetZ-coordinate multistep prediction performance comparison of four different prediction methods.
Prediction stepAlgorithmRMSEMAEMAPECOR1BP3.81681.80708.7264∗e−040.9959IGA-VNN1.49141.09433.2070∗e−040.9984GA-VNN3.52071.74378.1077∗e−040.9958Volterra2.06361.32944.7240∗e−040.9937RBF5.04112.06800.00110.97202BP7.68772.68370.00190.9881IGA-VNN4.98932.12060.00120.9968GA-VNN4.97182.20860.00130.9957Volterra6.26892.40250.00150.9479RBF7.99692.32840.00180.95724BP11.98213.25660.00280.8768IGA-VNN6.47402.38600.00150.9757GA-VNN8.04132.78890.00210.9978Volterra8.18162.54050.00170.9316RBF14.08473.3960.00320.90446BP14.12093.27620.00290.8273IGA-VNN9.84692.88520.00220.9677GA-VNN10.98082.97730.00240.9585Volterra10.13712.94630.00230.8627RBF20.96114.43270.00520.8633


## 6.4.2. Results Analysis

By comparing and analyzing the prediction performance of different models in Tables2–4, the following conclusions can be drawn:(1)
When performing multistep prediction on the target maneuver trajectory time series, the prediction performance of the chaotic time series IGA-VNN model is the best with the same prediction steps. Because the Volterra functional model can represent most nonlinear systems, it can accurately predict chaotic time series. However, the solution of high-order Volterra kernel function is the bottleneck of its application. In this paper, we establish the IGA-VNN model to determine the kernel functions of Volterra series, so that Volterra series can predict the target maneuver trajectory well.(2)
RBF neural network prediction model, second-order Volterra functional prediction model, and BP neural network prediction model show different prediction performance in different prediction steps and different coordinate directions of the target, which shows that the adaptability and robustness of these three prediction algorithms are poor, and the performance of the IGA-VNN prediction model is superior than that of these three prediction algorithms in this respect. In this paper, the phase space reconstruction theory is adopted to reconstruct the time series of target maneuver trajectory, which makes full use of the information contained in the time series of target maneuver trajectory so that the prediction accuracy is improved.(3)
Combined with the prediction error comparison charts in Figures14–16, it can be seen that the same prediction algorithm’s single-step prediction results for the three coordinates of the target maneuver trajectory are better than the 2-step prediction results; the 2-step prediction results are better than the 4-step prediction results; and the 4-step prediction results are better than 6-step prediction results. Therefore, it can be seen that for the same prediction model, the prediction performance of the model will also decline with the increase in prediction steps. In the multistep prediction of the target maneuver trajectory, as the number of prediction steps increases, the cumulative error of the prediction will continue to increase, so the prediction accuracy will also decrease.By comparing and analyzing the differences between the four models in the multistep prediction principle of target maneuver trajectory time series, the root cause of these differences can be further analyzed. BP neural network is based on sample data to train the network so as to realize the prediction of the future state. The important factors that affect the precision and prediction accuracy of the BP neural network are the subjectivity and blindness of its variable selection, which leads to the addition of extra irrelevant variables or ignoring important variables. In addition, BP neural network is sensitive to initial value and easy to fall into local minimum. The Volterra functional model essentially uses the chaotic characteristics of target maneuvering trajectory time series to predict the nonlinear system by training and fitting the trajectory of the chaotic attractor. The model is implemented by the multiplication-coupled form of the second-order Volterra series, which avoids the difficulty of solving the kernel function caused by the higher-order Volterra series and improves the prediction accuracy of the model. In addition, the Volterra filter uses the linear adaptive algorithm to adjust the filter parameters, and this dynamic property can improve the adaptability and accuracy of prediction and reduce the training time.The IGA-VNN prediction model based on chaos theory combines the discrete Volterra functional model, BP neural network, and the chaotic characteristics of the target maneuvering trajectory time series and uses the chaotic characteristics, causality, and memory function [25] of the target maneuvering trajectory time series to determine the truncation order and truncation term number of the Volterra functional model so as to establish an accurate Volterra functional model. At the same time, the equivalence of BP neural network and Volterra functional model is used to solve the problem of high-order kernel function of the Volterra functional model, and the global searching ability of the improved genetic algorithm is used to avoid falling into local minimum when optimizing the VNN model. In conclusion, the IGA-VNN prediction model based on chaos theory reasonably uses the advantages of the Volterra functional model and BP neural network and effectively solves the problem of subjectivity in BP neural network modeling and the difficulty in the Volterra functional model solving high-order kernel function, making the performance of the algorithm in multistep prediction better than that of the single BP neural network and Volterra functional model.


## 7. Conclusions

In this paper, a VNN network model is designed by analyzing the chaotic characteristics of the target maneuver trajectory time series, and the IGA-VNN prediction model of target maneuver trajectory based on chaos theory is proposed, and then the target maneuver trajectory time series is predicted in multiple steps. The prediction method proposed in this paper combines the advantages of the Volterra functional model and BP neural network algorithm, overcomes the blindness in building BP neural network, and solves the difficulties in determining the higher-order kernel function of the Volterra functional model, which has certain theoretical application value.(1)
In this paper, the C-C method is used to reconstruct the phase space of the target maneuver trajectory, and the chaotic characteristic of the target maneuver trajectory time series is proved by the method of small data volume.(2)
In this paper, the IGA-VNN target trajectory prediction model based on chaos theory is constructed. The model solves the problem of solving the higher-order kernel function of the Volterra functional model and realizes the accurate prediction of target maneuvering trajectory. In multistep prediction, the prediction performance of the model is better than that of the single BP neural network and Volterra functional prediction model.(3)
The IGA-VNN prediction model based on chaos theory proposed in this paper combines the chaotic characteristics of the target maneuver trajectory time series to determine the dimension of the input vector of the prediction network. It makes the truncation term and truncation order of the Volterra functional model have clear physical meaning, avoids the blindness brought by subjective determination of parameters, and improves the practical value of the model.(4)
In addition, in order to realize the application value of the IGA-VNN prediction model of chaotic time series, the realization of its hardware is also the focus of future research.



---

*Source: 8325498-2020-11-29.xml*