# Tracking Air-to-Air Missile Using Proportional Navigation Model with Genetic Algorithm Particle Filter

Authors: Hongqiang Liu; Lei Yu; Chenwei Ruan; Zhongliang Zhou
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

The purpose of this paper is to track the air-to-air missile. Here we put forward the PN-GAPF (Proportional Navigation motion model and Genetic Algorithm Particle Filter) method to solve the problem. The main jobs we have done can be listed as follows: firstly, we establish the missile state space model named as the Proportional Navigation (PN) motion model to simulate the real motion of the air-to-air missile; secondly, the PN-EKF and PN-PF methods are proposed to track the missile, through combining PN motion model with EKF and PF; thirdly, in order to solve the particle degeneracy and diversity loss, we introduce the intercross and variation in GA to the particles resampling step and then the PN-GAPF method is put forward. The simulation results show that the PN motion model is better than the CV and CA motion models for tracking the air-to-air missile and that the PN-GAPF method is more efficient than the PN-EKF and PN-PF.


---

## Body

## 1. Introduction

The air-to-air missile is the main weapon in the air combat. The aircraft been chased should maneuver to avoid the attacking missile, after the air-to-air missile was fired by an opposition fighter. The aircraft guidance method for maneuver evasion is based on knowing the missile location and the real-time track [1]. However, we cannot acquire the exact location information of the missile because the measurement for the missile location has great error. We need an on-line filter method to eliminate the error and track the air-to-air missile. How to effectively track the missile is the research content in this paper.This problem is the domain of the single target tracking. The common target motion models such as the CV, CA, and Current Statistical (CS) models cannot well and truly describe the air-to-air missile maneuver because the missile has a good maneuverability and a supersonic speed [2]. However, we can research from the navigation law point to establish a smarter motion model because the missile maneuver obeys some navigation law [3]. The Proportional Navigation (PN) law is the most common in the air-to-air missiles [4]. Therefore we establish a new PN motion model in 3d Cartesian coordinate system by analysis of the PN mechanism.The state space model for tracking the air-to-air missile is obtained further, through combining with the nonlinear measurement equation in the radar spherical coordinate system. The standard Kalman Filter (KF) cannot be used to track the missile because of the nonlinear measurement equation. The Extended Kalman Filter (EKF) and Unscented Kalman Filter (UKF) are used to solve the nonlinear problem, but the methods are not efficient enough in practice [5]. Recently, the Particle Filter (PF) or Monte Carlo (MC) method is widely applied in the nonlinear and non-Gaussian filter problem because PF is to be able to represent the required unknown probability density function by a scatter of particles sampled from a known probability density function [6]. However, the resampling of standard PF has a disadvantage of the loss of particle diversity [7]. We are inspired by the evolution idea of Genetic Algorithm (GA) to improve the performance of the PF at resampling step [8]. We take the particles and weights in the PF as the chromosomes and adaptability in the GA, respectively. Then the intercross and variation steps in the GA can be adopted to deal with the particles in resampling step for avoiding the particle degeneracy and loss of diversity. In this paper, Genetic Algorithm Particle Filter (GAPF) combined with the PN motion model (PN-GAPF) is used to track the air-to-air missile. PN-GAPF has better the estimate accuracy and the tracking stability, in comparison with other algorithms and models in computer simulation.This paper is divided into five sections. Section1 is introduction. Section 2 presents the state space model for tracking the air-to-air missile and the method of establishing the PN motion model. Section 3 describes the nonlinear filter algorithms such as the EKF, PF, and GAPF proposed. Section 4 discusses the performance of the PN motion model and GAPF algorithm on the base of the simulation results. Section 5 summarizes the main research content.


## 2. Problem Formulation

### 2.1. The State Space Model

The state space model for tracking the air-to-air missile consists of the missile state equation and the measurement equation. In order to simplify the problem, the missile state equation is usually modeled linearly in Cartesian coordinate system and the measurement equation is expressed nonlinearly in spherical coordinate system because of the air-to-air missile measurements given in the radar spherical coordinate system [9].


#### 2.1.1. The Linear Missile State Equation

In general, we estimate the state of the moving air-to-air missile with the discrete-time linear state space in 3d Cartesian coordinate system [10]:(1)Xk+1=AXk+vk,where Xk+1 and Xk are the state vector at time t=k+1T and t=kT respectively and T is a sampling period. A is the state transition matrix. vk is the process noise which is modeled as a zero-mean white Gaussian process with covariance matrix Q(k).The CV model and the CA model are the most common motion models for tracking air-to-air missile [11]. The state vector XCVk and the state transition matrix ACV of CV model in 3d Cartesian coordinate system are(2)XCVk=xk,vxk,yk,vyk,zk,vzk,ACV=1T0000010000001T0000010000001T000001.And the state vector XCAk and the state transition matrix ACA of CA model are(3)
X
C
A
k
=
x
k
,
v
x
k
,
a
x
k
,
y
k
,
v
y
k
,
a
y
k
,
z
k
,
v
z
k
,
a
z
k
,
A
C
A
=
E
3
×
3
dig
T
3
×
3
dig
1
2
·
T
2
3
×
3
0
3
×
3
E
3
×
3
dig
T
3
×
3
0
3
×
3
0
3
×
3
E
3
×
3
,
dig
T
3
×
3
=
T
0
0
0
T
0
0
0
T
,
dig
T
2
2
3
×
3
=
T
2
2
0
0
0
T
2
2
0
0
0
T
2
2
,where xk, yk, and zk are the missile positions along each axis of coordinates. Correspondingly, vxk, vyk, and vzk are the velocities. axk, ayk, and azk are the accelerators. If sampling the measurement in a very short time, we can consider the missiles maneuvering as an approximate CV model or CA model [11].


#### 2.1.2. The Nonlinear Measurement Equation

In the real air combat, the measurement of air-to-air missile is as follows in the radar spherical coordinate system [12]:(4)Zk=rk,θk,φk,where Zk is the measurement vector. rk is the distance between the aircraft attacked and the air-to-air missile which can be detected by the airborne laser range finder. θk and φk are the pitching and azimuth angles, respectively, which can be detected by the aircraft radar warning device [13].The nonlinear measurement equation is(5)Zk=hXk+wk,where h(·) is the nonlinear measurement function which will be given at the next section. w(k) is the measurement error which is zero-mean white Gaussian noise with covariance matrix R(k) [14].


### 2.2. The Proportional Navigation Motion Model

If the air-to-air missile movement is modeled just simply as the CV or CA motion model, it will produce great errors of the missile state estimation [4]. It is known that the air-to-air missile is designed to maneuver by the Proportional Navigation law. Therefore, the motion model based on the Proportional Navigation law is a smarter model to reflect the states in the attack process. We describe the relative movement between the aircraft attacked and the air-to-air missile in Figure 1.Figure 1
The relative movement and variables.


#### 2.2.1. The Geographic Coordinate System

We take “north-sky-west” geographic coordinate system [15] as the 3d Cartesian coordinate system to describe the relative movement in Figure 1. The origin of coordinates “O” is set at the firing point of the air-to-air missile. The “OX” axis is tangent with longitude and points to the north. The “OY” axis points to the sky and the “OZ” axis points to the west.


#### 2.2.2. The Relative Movement

We assume that the attacked aircraft position and the missile position arexf,yf,zf and xm,ym,zm in “OXYZ” coordinate system, respectively. r,θr,φr are the relative distance and the relative pitching and relative azimuth angles between them in the spherical coordinate system. The relationship of those variables is [16](6)r=xf-xm2+yf-ym2+zf-zm2,φr=arctan⁡zf-zmxf-xm,θr=arctan⁡yf-ymxf-xm2+zf-zm2.The velocity vectors of the aircraft and the missile arevfx,vfy,vfz and vmx,vmy,vmz in the “OXYZ” coordinate system; then vm,θm,φm and vf,θf,φf are in the spherical coordinate system. The relationships of those variables are(7)dxmdt=vmx=vmsin⁡θmcos⁡φm,dymdt=vmy=vmcos⁡θm,dzmdt=vmz=vmsin⁡θmsin⁡φm,(8)dxfdt=vfx=vfsin⁡θfcos⁡φf,dyfdt=vfy=vfcos⁡θf,dzfdt=vfz=vfsin⁡θfsin⁡φf.Popularly, the Proportional Navigation law can be decoupled in the “OXY” plane and “OXZ” plane [17]:(9)θ˙m=Kθθ˙r,φ˙m=Kφφ˙r,where ⁡Kθ and Kφ are the Proportional Navigation coefficients, 2≤K≤6. θ˙mφ˙mθ˙r and φ˙r are the derivations of θmφmθr and φr, respectively.


#### 2.2.3. Modeling the Proportional Navigation State Space Model

We assume that the velocityvm is a fixed value [18] because the roll booster of air-to-air missile works in a short time. Then we take derivative of (7):(10)v˙mx=θ˙mvmcos⁡θmcos⁡φm-φ˙mvmsin⁡θmsin⁡φm,v˙my=θ˙mvmsin⁡θm,v˙mz=-θ˙mvmcos⁡θmsin⁡φm+φ˙mvmsin⁡θmcos⁡φm.Combining (7) with (8), (10) will be changed to(11)v˙mx=vmyKθθ˙rcos⁡φm-vmzKφφ˙r,v˙my=vmyKθθ˙rtan⁡θm,v˙mz=vmxKφφ˙r+vmyKθθ˙rsin⁡φm.Now we take the missile state vector asXPN=xm,vmx,ym,vmx,zm,vmz. According to (11), the continuous-time linear state equation [19] can be established as(12)X˙=FPNX+v,(13)FPN=010000000Kθθ˙rcos⁡φm0Kφφ˙r000100000Kθθ˙rtan⁡θm000000010Kφφ˙r0Kθθ˙rsin⁡φm00.Then we convert (12) to the discrete linear state equation as(14)Xk+1=APNXk+vk,APN=eFPNT.The precise analytic expression ofAPN cannot be acquired, but we use (15) to get an approximation by taking k=3 [20]:(15)
e
F
P
N
T
=
∑
k
=
0
∞
F
P
N
T
k
k
!
≈
I
+
F
P
N
T
+
F
P
N
2
T
2
2
!
+
F
P
N
3
T
3
3
!
,
(16)
A
P
N
=
1
t
-
1
6
t
3
b
2
0
1
2
t
2
a
+
1
6
t
3
a
c
+
b
d
0
1
2
t
2
b
0
-
1
2
t
2
b
2
+
1
0
t
a
+
1
2
t
2
a
c
+
b
d
+
1
6
t
3
-
b
2
a
+
a
c
+
b
d
c
0
t
b
-
1
6
t
3
b
3
0
0
1
t
+
1
2
t
2
c
+
1
6
t
3
c
2
0
0
0
0
0
t
c
+
1
2
t
2
c
2
+
1
6
t
3
c
3
+
1
0
0
0
-
1
2
t
2
b
0
1
2
t
2
d
+
1
6
t
3
-
a
b
+
d
c
1
t
-
1
6
t
3
b
2
0
-
t
b
+
1
6
t
3
b
3
0
t
d
+
1
2
t
2
-
a
b
+
d
c
+
1
6
t
3
-
a
b
+
d
c
c
-
b
2
d
0
-
1
2
t
2
b
2
+
1
,
(17)
a
=
-
K
θ
θ
˙
r
cos
⁡
φ
m
,
b
=
-
K
φ
φ
˙
r
,
c
=
-
K
θ
θ
˙
r
tan
⁡
θ
m
,
d
=
K
θ
θ
˙
r
sin
⁡
φ
r
,
t
=
T
.The measurement equation can directly be modeled in discrete form like (5) because the measurements of the missile are detected discretely. According to (6), hXk is(18)hXk=xfk-xmk2+yfk-ymk2+zfk-zmk2arctan⁡yfk-ymkxfk-xmk2+zfk-zmk2arctan⁡zfk-zmkxfk-xmk.


## 2.1. The State Space Model

The state space model for tracking the air-to-air missile consists of the missile state equation and the measurement equation. In order to simplify the problem, the missile state equation is usually modeled linearly in Cartesian coordinate system and the measurement equation is expressed nonlinearly in spherical coordinate system because of the air-to-air missile measurements given in the radar spherical coordinate system [9].


### 2.1.1. The Linear Missile State Equation

In general, we estimate the state of the moving air-to-air missile with the discrete-time linear state space in 3d Cartesian coordinate system [10]:(1)Xk+1=AXk+vk,where Xk+1 and Xk are the state vector at time t=k+1T and t=kT respectively and T is a sampling period. A is the state transition matrix. vk is the process noise which is modeled as a zero-mean white Gaussian process with covariance matrix Q(k).The CV model and the CA model are the most common motion models for tracking air-to-air missile [11]. The state vector XCVk and the state transition matrix ACV of CV model in 3d Cartesian coordinate system are(2)XCVk=xk,vxk,yk,vyk,zk,vzk,ACV=1T0000010000001T0000010000001T000001.And the state vector XCAk and the state transition matrix ACA of CA model are(3)
X
C
A
k
=
x
k
,
v
x
k
,
a
x
k
,
y
k
,
v
y
k
,
a
y
k
,
z
k
,
v
z
k
,
a
z
k
,
A
C
A
=
E
3
×
3
dig
T
3
×
3
dig
1
2
·
T
2
3
×
3
0
3
×
3
E
3
×
3
dig
T
3
×
3
0
3
×
3
0
3
×
3
E
3
×
3
,
dig
T
3
×
3
=
T
0
0
0
T
0
0
0
T
,
dig
T
2
2
3
×
3
=
T
2
2
0
0
0
T
2
2
0
0
0
T
2
2
,where xk, yk, and zk are the missile positions along each axis of coordinates. Correspondingly, vxk, vyk, and vzk are the velocities. axk, ayk, and azk are the accelerators. If sampling the measurement in a very short time, we can consider the missiles maneuvering as an approximate CV model or CA model [11].


### 2.1.2. The Nonlinear Measurement Equation

In the real air combat, the measurement of air-to-air missile is as follows in the radar spherical coordinate system [12]:(4)Zk=rk,θk,φk,where Zk is the measurement vector. rk is the distance between the aircraft attacked and the air-to-air missile which can be detected by the airborne laser range finder. θk and φk are the pitching and azimuth angles, respectively, which can be detected by the aircraft radar warning device [13].The nonlinear measurement equation is(5)Zk=hXk+wk,where h(·) is the nonlinear measurement function which will be given at the next section. w(k) is the measurement error which is zero-mean white Gaussian noise with covariance matrix R(k) [14].


## 2.1.1. The Linear Missile State Equation

In general, we estimate the state of the moving air-to-air missile with the discrete-time linear state space in 3d Cartesian coordinate system [10]:(1)Xk+1=AXk+vk,where Xk+1 and Xk are the state vector at time t=k+1T and t=kT respectively and T is a sampling period. A is the state transition matrix. vk is the process noise which is modeled as a zero-mean white Gaussian process with covariance matrix Q(k).The CV model and the CA model are the most common motion models for tracking air-to-air missile [11]. The state vector XCVk and the state transition matrix ACV of CV model in 3d Cartesian coordinate system are(2)XCVk=xk,vxk,yk,vyk,zk,vzk,ACV=1T0000010000001T0000010000001T000001.And the state vector XCAk and the state transition matrix ACA of CA model are(3)
X
C
A
k
=
x
k
,
v
x
k
,
a
x
k
,
y
k
,
v
y
k
,
a
y
k
,
z
k
,
v
z
k
,
a
z
k
,
A
C
A
=
E
3
×
3
dig
T
3
×
3
dig
1
2
·
T
2
3
×
3
0
3
×
3
E
3
×
3
dig
T
3
×
3
0
3
×
3
0
3
×
3
E
3
×
3
,
dig
T
3
×
3
=
T
0
0
0
T
0
0
0
T
,
dig
T
2
2
3
×
3
=
T
2
2
0
0
0
T
2
2
0
0
0
T
2
2
,where xk, yk, and zk are the missile positions along each axis of coordinates. Correspondingly, vxk, vyk, and vzk are the velocities. axk, ayk, and azk are the accelerators. If sampling the measurement in a very short time, we can consider the missiles maneuvering as an approximate CV model or CA model [11].


## 2.1.2. The Nonlinear Measurement Equation

In the real air combat, the measurement of air-to-air missile is as follows in the radar spherical coordinate system [12]:(4)Zk=rk,θk,φk,where Zk is the measurement vector. rk is the distance between the aircraft attacked and the air-to-air missile which can be detected by the airborne laser range finder. θk and φk are the pitching and azimuth angles, respectively, which can be detected by the aircraft radar warning device [13].The nonlinear measurement equation is(5)Zk=hXk+wk,where h(·) is the nonlinear measurement function which will be given at the next section. w(k) is the measurement error which is zero-mean white Gaussian noise with covariance matrix R(k) [14].


## 2.2. The Proportional Navigation Motion Model

If the air-to-air missile movement is modeled just simply as the CV or CA motion model, it will produce great errors of the missile state estimation [4]. It is known that the air-to-air missile is designed to maneuver by the Proportional Navigation law. Therefore, the motion model based on the Proportional Navigation law is a smarter model to reflect the states in the attack process. We describe the relative movement between the aircraft attacked and the air-to-air missile in Figure 1.Figure 1
The relative movement and variables.


### 2.2.1. The Geographic Coordinate System

We take “north-sky-west” geographic coordinate system [15] as the 3d Cartesian coordinate system to describe the relative movement in Figure 1. The origin of coordinates “O” is set at the firing point of the air-to-air missile. The “OX” axis is tangent with longitude and points to the north. The “OY” axis points to the sky and the “OZ” axis points to the west.


### 2.2.2. The Relative Movement

We assume that the attacked aircraft position and the missile position arexf,yf,zf and xm,ym,zm in “OXYZ” coordinate system, respectively. r,θr,φr are the relative distance and the relative pitching and relative azimuth angles between them in the spherical coordinate system. The relationship of those variables is [16](6)r=xf-xm2+yf-ym2+zf-zm2,φr=arctan⁡zf-zmxf-xm,θr=arctan⁡yf-ymxf-xm2+zf-zm2.The velocity vectors of the aircraft and the missile arevfx,vfy,vfz and vmx,vmy,vmz in the “OXYZ” coordinate system; then vm,θm,φm and vf,θf,φf are in the spherical coordinate system. The relationships of those variables are(7)dxmdt=vmx=vmsin⁡θmcos⁡φm,dymdt=vmy=vmcos⁡θm,dzmdt=vmz=vmsin⁡θmsin⁡φm,(8)dxfdt=vfx=vfsin⁡θfcos⁡φf,dyfdt=vfy=vfcos⁡θf,dzfdt=vfz=vfsin⁡θfsin⁡φf.Popularly, the Proportional Navigation law can be decoupled in the “OXY” plane and “OXZ” plane [17]:(9)θ˙m=Kθθ˙r,φ˙m=Kφφ˙r,where ⁡Kθ and Kφ are the Proportional Navigation coefficients, 2≤K≤6. θ˙mφ˙mθ˙r and φ˙r are the derivations of θmφmθr and φr, respectively.


### 2.2.3. Modeling the Proportional Navigation State Space Model

We assume that the velocityvm is a fixed value [18] because the roll booster of air-to-air missile works in a short time. Then we take derivative of (7):(10)v˙mx=θ˙mvmcos⁡θmcos⁡φm-φ˙mvmsin⁡θmsin⁡φm,v˙my=θ˙mvmsin⁡θm,v˙mz=-θ˙mvmcos⁡θmsin⁡φm+φ˙mvmsin⁡θmcos⁡φm.Combining (7) with (8), (10) will be changed to(11)v˙mx=vmyKθθ˙rcos⁡φm-vmzKφφ˙r,v˙my=vmyKθθ˙rtan⁡θm,v˙mz=vmxKφφ˙r+vmyKθθ˙rsin⁡φm.Now we take the missile state vector asXPN=xm,vmx,ym,vmx,zm,vmz. According to (11), the continuous-time linear state equation [19] can be established as(12)X˙=FPNX+v,(13)FPN=010000000Kθθ˙rcos⁡φm0Kφφ˙r000100000Kθθ˙rtan⁡θm000000010Kφφ˙r0Kθθ˙rsin⁡φm00.Then we convert (12) to the discrete linear state equation as(14)Xk+1=APNXk+vk,APN=eFPNT.The precise analytic expression ofAPN cannot be acquired, but we use (15) to get an approximation by taking k=3 [20]:(15)
e
F
P
N
T
=
∑
k
=
0
∞
F
P
N
T
k
k
!
≈
I
+
F
P
N
T
+
F
P
N
2
T
2
2
!
+
F
P
N
3
T
3
3
!
,
(16)
A
P
N
=
1
t
-
1
6
t
3
b
2
0
1
2
t
2
a
+
1
6
t
3
a
c
+
b
d
0
1
2
t
2
b
0
-
1
2
t
2
b
2
+
1
0
t
a
+
1
2
t
2
a
c
+
b
d
+
1
6
t
3
-
b
2
a
+
a
c
+
b
d
c
0
t
b
-
1
6
t
3
b
3
0
0
1
t
+
1
2
t
2
c
+
1
6
t
3
c
2
0
0
0
0
0
t
c
+
1
2
t
2
c
2
+
1
6
t
3
c
3
+
1
0
0
0
-
1
2
t
2
b
0
1
2
t
2
d
+
1
6
t
3
-
a
b
+
d
c
1
t
-
1
6
t
3
b
2
0
-
t
b
+
1
6
t
3
b
3
0
t
d
+
1
2
t
2
-
a
b
+
d
c
+
1
6
t
3
-
a
b
+
d
c
c
-
b
2
d
0
-
1
2
t
2
b
2
+
1
,
(17)
a
=
-
K
θ
θ
˙
r
cos
⁡
φ
m
,
b
=
-
K
φ
φ
˙
r
,
c
=
-
K
θ
θ
˙
r
tan
⁡
θ
m
,
d
=
K
θ
θ
˙
r
sin
⁡
φ
r
,
t
=
T
.The measurement equation can directly be modeled in discrete form like (5) because the measurements of the missile are detected discretely. According to (6), hXk is(18)hXk=xfk-xmk2+yfk-ymk2+zfk-zmk2arctan⁡yfk-ymkxfk-xmk2+zfk-zmk2arctan⁡zfk-zmkxfk-xmk.


## 2.2.1. The Geographic Coordinate System

We take “north-sky-west” geographic coordinate system [15] as the 3d Cartesian coordinate system to describe the relative movement in Figure 1. The origin of coordinates “O” is set at the firing point of the air-to-air missile. The “OX” axis is tangent with longitude and points to the north. The “OY” axis points to the sky and the “OZ” axis points to the west.


## 2.2.2. The Relative Movement

We assume that the attacked aircraft position and the missile position arexf,yf,zf and xm,ym,zm in “OXYZ” coordinate system, respectively. r,θr,φr are the relative distance and the relative pitching and relative azimuth angles between them in the spherical coordinate system. The relationship of those variables is [16](6)r=xf-xm2+yf-ym2+zf-zm2,φr=arctan⁡zf-zmxf-xm,θr=arctan⁡yf-ymxf-xm2+zf-zm2.The velocity vectors of the aircraft and the missile arevfx,vfy,vfz and vmx,vmy,vmz in the “OXYZ” coordinate system; then vm,θm,φm and vf,θf,φf are in the spherical coordinate system. The relationships of those variables are(7)dxmdt=vmx=vmsin⁡θmcos⁡φm,dymdt=vmy=vmcos⁡θm,dzmdt=vmz=vmsin⁡θmsin⁡φm,(8)dxfdt=vfx=vfsin⁡θfcos⁡φf,dyfdt=vfy=vfcos⁡θf,dzfdt=vfz=vfsin⁡θfsin⁡φf.Popularly, the Proportional Navigation law can be decoupled in the “OXY” plane and “OXZ” plane [17]:(9)θ˙m=Kθθ˙r,φ˙m=Kφφ˙r,where ⁡Kθ and Kφ are the Proportional Navigation coefficients, 2≤K≤6. θ˙mφ˙mθ˙r and φ˙r are the derivations of θmφmθr and φr, respectively.


## 2.2.3. Modeling the Proportional Navigation State Space Model

We assume that the velocityvm is a fixed value [18] because the roll booster of air-to-air missile works in a short time. Then we take derivative of (7):(10)v˙mx=θ˙mvmcos⁡θmcos⁡φm-φ˙mvmsin⁡θmsin⁡φm,v˙my=θ˙mvmsin⁡θm,v˙mz=-θ˙mvmcos⁡θmsin⁡φm+φ˙mvmsin⁡θmcos⁡φm.Combining (7) with (8), (10) will be changed to(11)v˙mx=vmyKθθ˙rcos⁡φm-vmzKφφ˙r,v˙my=vmyKθθ˙rtan⁡θm,v˙mz=vmxKφφ˙r+vmyKθθ˙rsin⁡φm.Now we take the missile state vector asXPN=xm,vmx,ym,vmx,zm,vmz. According to (11), the continuous-time linear state equation [19] can be established as(12)X˙=FPNX+v,(13)FPN=010000000Kθθ˙rcos⁡φm0Kφφ˙r000100000Kθθ˙rtan⁡θm000000010Kφφ˙r0Kθθ˙rsin⁡φm00.Then we convert (12) to the discrete linear state equation as(14)Xk+1=APNXk+vk,APN=eFPNT.The precise analytic expression ofAPN cannot be acquired, but we use (15) to get an approximation by taking k=3 [20]:(15)
e
F
P
N
T
=
∑
k
=
0
∞
F
P
N
T
k
k
!
≈
I
+
F
P
N
T
+
F
P
N
2
T
2
2
!
+
F
P
N
3
T
3
3
!
,
(16)
A
P
N
=
1
t
-
1
6
t
3
b
2
0
1
2
t
2
a
+
1
6
t
3
a
c
+
b
d
0
1
2
t
2
b
0
-
1
2
t
2
b
2
+
1
0
t
a
+
1
2
t
2
a
c
+
b
d
+
1
6
t
3
-
b
2
a
+
a
c
+
b
d
c
0
t
b
-
1
6
t
3
b
3
0
0
1
t
+
1
2
t
2
c
+
1
6
t
3
c
2
0
0
0
0
0
t
c
+
1
2
t
2
c
2
+
1
6
t
3
c
3
+
1
0
0
0
-
1
2
t
2
b
0
1
2
t
2
d
+
1
6
t
3
-
a
b
+
d
c
1
t
-
1
6
t
3
b
2
0
-
t
b
+
1
6
t
3
b
3
0
t
d
+
1
2
t
2
-
a
b
+
d
c
+
1
6
t
3
-
a
b
+
d
c
c
-
b
2
d
0
-
1
2
t
2
b
2
+
1
,
(17)
a
=
-
K
θ
θ
˙
r
cos
⁡
φ
m
,
b
=
-
K
φ
φ
˙
r
,
c
=
-
K
θ
θ
˙
r
tan
⁡
θ
m
,
d
=
K
θ
θ
˙
r
sin
⁡
φ
r
,
t
=
T
.The measurement equation can directly be modeled in discrete form like (5) because the measurements of the missile are detected discretely. According to (6), hXk is(18)hXk=xfk-xmk2+yfk-ymk2+zfk-zmk2arctan⁡yfk-ymkxfk-xmk2+zfk-zmk2arctan⁡zfk-zmkxfk-xmk.


## 3. Nonlinear Filtering Algorithm

The Kalman filter assumptions do not hold because the measurement equation is a nonlinear function in the measurement equation [21]. Therefore, for the tracking missile, there are three approximate nonlinear filters to solve the nonlinear problem in the paper.


### 3.1. The Extended Kalman Filter with the PN Motion Model

The EKF approximates the nonlinear function through utilizing the first term in a Taylor expansion of the nonlinear function [22]:(19)Hk=dhXdXX=Xk∣k-1,(20)Hk=∂r∂xm∂r∂vmx∂r∂ym∂r∂vmy∂r∂zm∂r∂vmz∂θr∂xm∂θr∂vmx∂θr∂ym∂θr∂vmy∂θr∂zm∂θr∂vmz∂φr∂xm∂φr∂vmx∂φr∂ym∂φr∂vmy∂φr∂zm∂φr∂vmz,(21)∂r∂xm=-xf-xmxf-xm2+yf-ym2+zf-zm2,∂r∂ym=-yf-ymxf-xm2+yf-ym2+zf-zm2,∂r∂zm=-zf-zmxf-xm2+yf-ym2+zf-zm2,∂r∂vmx=∂r∂vmy=∂r∂vmz=0,∂θr∂xm=yf-ymxf-xmxf-xm2+zf-zm2·xf-xm2+yf-ym2+zf-zm2,∂θr∂ym=-xf-xm2+zf-zm2xf-xm2+zf-zm2·xf-xm2+yf-ym2+zf-zm2,∂θr∂zm=yf-ymzf-zmxf-xm2+zf-zm2·xf-xm2+yf-ym2+zf-zm2,∂θr∂vmx=∂θr∂vmy=∂θr∂vmz=0,∂φr∂xm=zf-zmxf-xm2+zf-zm2,∂φr∂zm=xf-xmxf-xm2+zf-zm2,∂qz∂vmx=∂qz∂vmy=∂qz∂vmz=∂qz∂ym=0.K
θ
K
φ
φ
m
θ
m
φ
˙
r and θ˙r should be estimated on the basis of the prior estimated information because they are unknown variables in state transport matrix:(22)θ^rk=arccos⁡y^mk-1-yfk-1x^mk-1-xfk-12+y^mk-1-yfk-12+z^mk-1-zfk-12,φ^rk=arcsin⁡z^mk-1-z^fk-1x^mk-1-xfk-12+z^mk-1-zfk-12,θ^mk=arccos⁡v^myk-1v^mx2k-1+v^my2k-1+v^mz2k-1,φ^mk=arcsin⁡v^mzk-1v^mx2k-1+v^mz2k-1,θ˙^rk=θ^rk-θ^rk-1T,φ˙^rk=φ^rk-φ^rk-1T,θ˙^mk=θ^mk-θ^mk-1T,φ˙^mk=φ^mk-φ^mk-1T,K^θk=θ˙^mkθ˙^rk,K^φk=φ˙^mkφ˙^rk.The EKF with the PN motion model to track air-to-air missile is shown in Algorithm1.Algorithm 1
NT is the total tracking time
Initialize the variables whenk=1 and k=2
FORk=3:N
Xk∣k-1=APNXk-1
Pk∣k-1=APNPk-1APNT+Qk-1
Kk=Pk∣k-1HkTHkPk∣k-1HkT+Rk-1-1
Z-k=Zk-hXk∣k-1
Xk=Xk∣k-1+KkZ-k
Pk=I-KkHkPk∣k-1
Calculate (22) to update APN
END FOR


### 3.2. The Particle Filter

#### 3.2.1. The Sequential Importance Sampling (SIS) Algorithm

The Particle Filter algorithm to deal with the nonlinear tracking problem has been used on a wide consensus [6]. Let Z1:k={Z1,Z2,…,Zk} denote a sequence of observations and pXk∣Z1:k denote the posterior density of Xk. We can acquire Bayesian estimation X^k and P^k by(23)X^k=∫XkpXk∣Z1:kdXk,P^k=∫Xk-X^kXk-X^kTpXk∣Z1:kdXk.The calculation of the posterior density requires integrals, but the normalizing integral cannot be evaluated analytically and numerical integration over possibly high-dimensional spaces is infeasible. The ideal of SIS algorithm is that the required posterior density function is represented by a set of random samples with associated weight{xik,ωik}i=1Ns and to compute state estimates in the basis of these samples and weights. According to SIS algorithm, the posterior density pXk∣Z1:k can be approximated as(24)pXk∣Z1:k≈∑i=1NsωikδXk-xik,where the weights ωik are defined in(25)ωik∝ωik-1pZk∣xikpxik∣xik-1qxik∣xik-1,Zk,where qxik∣xik-1,Zk is termed the importance sampling proposal distribution. It is easy to draw samples from the distribution defined by qxik∣xik-1,Zk, but it is hard to draw samples directly from the posterior density pXk∣Z1:k. In standard Particle Filter, we set qxik∣xik-1,Zk as pxik∣xik-1:(26)qxik∣xik-1,Zk=pxik∣xik-1.Then(27)ωik∝ωik-1pZk∣xik.Assuming thatvk and w(k) are mutually independent zero-mean white Gaussian processes [7] in (1) and (5), the prior density pxik∣xik-1 and the likelihood density pZk∣xik become(28)pxik∣xik-1=Nxik;APNxik-1,Qk,pZk∣xik=NZk;hxik,Rk.


#### 3.2.2. The Resampling Step

After a few iterations with SIS, the variance of the importance weights increases because most particles will have negligible weights. This is called particle degeneracy. Fixing on an importance sampling proposal distribution, one way to mitigate the particle degeneracy is to adopt the resampling step. The effective sample sizeNeff denotes an index variable of degeneracy introduced in [6]. Neff can be estimated approximately by(29)Neff≈1∑i=1Nsω~ik2.IfNeff falls below a threshold NT, the resampling algorithm will sample from (24) to generate a new set xi∗ki=1Ns and reset the weights ωi∗k to 1/Ns.The standard PF with the PN motion model to track the air-to-air missile is shown in Algorithm2.Algorithm 2
NT is the total tracking time
Initialize the variables whenk=1 and k=2
FORk=3:N
FORi=1:Ns
Drawxik~pxik∣xik-1
Assign the particle a weight according to (27)
END FOR
FORj=1:Ns
Normalize the weights by
ω~jk=ωjk/∑i=1Nsωik2
END FOR
EstimateX^k and P^k by
X^k=∑j=1Nsω~jkxjk
P^k=∑j=1Nsω~jkxjk-X^kxjk-X^kT
Calculate (22) to update APN
CalculateNeff using (29)
IFNeff<NT
Take resampling algorithm
End FOR


### 3.3. The Genetic Algorithm Particle Filter

The standard Particle Filter can solve the particle degeneracy, but the resampling algorithm relying on the high weight particles leads to the loss diversity [7]. The evolution idea of Genetic Algorithm (GA) [8] can be adopted in the resampling step to solve the problem. The chromosomes and adaptability in the GA are corresponding to the particle and weight in the PF, respectively. So we can use the selection, intercross, and variation steps to deal with the particle in resampling step to increase the diversity.GA for the resampling step is shown in Algorithm3.Algorithm 3
IFNeff<NT
FOR i=1:N1
Chose randomly two particlesxak and xbk from xik,ωiki=1Ns to intercross by
x~ak=εxak+1-εxbkx~bk=εxbk+1-εxakω~ak=εωak+1-εωbkω~bk=εωbk+1-εωak
whereε~U0,1 and a,b∈1,2,…Ns
END FOR
Create a new setx~ik,ω~iki=1Ns
FORi=1:N2
Chose randomly a particlex~ck from x~ik,ω~iki=1Ns to vary by
x⌢ck=x~ck+η
whereη~NX^k,P^k
IfpZk∣x⌢ck>pZk∣x~ck, accept x⌢ck
Else ifpZk∣x⌢ck/pZk∣x~ck>u  u~U0,1, accept x⌢ck
Otherwise acceptx~ck
END IF
END FOR
Create a new setx⌢ik,ω⌢iki=1Ns
Else leap over the resampling step
END IFThe GAPF uses the GA to replace the standard resampling step in the PN and increases the diversity, which is validated by simulations.


## 3.1. The Extended Kalman Filter with the PN Motion Model

The EKF approximates the nonlinear function through utilizing the first term in a Taylor expansion of the nonlinear function [22]:(19)Hk=dhXdXX=Xk∣k-1,(20)Hk=∂r∂xm∂r∂vmx∂r∂ym∂r∂vmy∂r∂zm∂r∂vmz∂θr∂xm∂θr∂vmx∂θr∂ym∂θr∂vmy∂θr∂zm∂θr∂vmz∂φr∂xm∂φr∂vmx∂φr∂ym∂φr∂vmy∂φr∂zm∂φr∂vmz,(21)∂r∂xm=-xf-xmxf-xm2+yf-ym2+zf-zm2,∂r∂ym=-yf-ymxf-xm2+yf-ym2+zf-zm2,∂r∂zm=-zf-zmxf-xm2+yf-ym2+zf-zm2,∂r∂vmx=∂r∂vmy=∂r∂vmz=0,∂θr∂xm=yf-ymxf-xmxf-xm2+zf-zm2·xf-xm2+yf-ym2+zf-zm2,∂θr∂ym=-xf-xm2+zf-zm2xf-xm2+zf-zm2·xf-xm2+yf-ym2+zf-zm2,∂θr∂zm=yf-ymzf-zmxf-xm2+zf-zm2·xf-xm2+yf-ym2+zf-zm2,∂θr∂vmx=∂θr∂vmy=∂θr∂vmz=0,∂φr∂xm=zf-zmxf-xm2+zf-zm2,∂φr∂zm=xf-xmxf-xm2+zf-zm2,∂qz∂vmx=∂qz∂vmy=∂qz∂vmz=∂qz∂ym=0.K
θ
K
φ
φ
m
θ
m
φ
˙
r and θ˙r should be estimated on the basis of the prior estimated information because they are unknown variables in state transport matrix:(22)θ^rk=arccos⁡y^mk-1-yfk-1x^mk-1-xfk-12+y^mk-1-yfk-12+z^mk-1-zfk-12,φ^rk=arcsin⁡z^mk-1-z^fk-1x^mk-1-xfk-12+z^mk-1-zfk-12,θ^mk=arccos⁡v^myk-1v^mx2k-1+v^my2k-1+v^mz2k-1,φ^mk=arcsin⁡v^mzk-1v^mx2k-1+v^mz2k-1,θ˙^rk=θ^rk-θ^rk-1T,φ˙^rk=φ^rk-φ^rk-1T,θ˙^mk=θ^mk-θ^mk-1T,φ˙^mk=φ^mk-φ^mk-1T,K^θk=θ˙^mkθ˙^rk,K^φk=φ˙^mkφ˙^rk.The EKF with the PN motion model to track air-to-air missile is shown in Algorithm1.Algorithm 1
NT is the total tracking time
Initialize the variables whenk=1 and k=2
FORk=3:N
Xk∣k-1=APNXk-1
Pk∣k-1=APNPk-1APNT+Qk-1
Kk=Pk∣k-1HkTHkPk∣k-1HkT+Rk-1-1
Z-k=Zk-hXk∣k-1
Xk=Xk∣k-1+KkZ-k
Pk=I-KkHkPk∣k-1
Calculate (22) to update APN
END FOR


## 3.2. The Particle Filter

### 3.2.1. The Sequential Importance Sampling (SIS) Algorithm

The Particle Filter algorithm to deal with the nonlinear tracking problem has been used on a wide consensus [6]. Let Z1:k={Z1,Z2,…,Zk} denote a sequence of observations and pXk∣Z1:k denote the posterior density of Xk. We can acquire Bayesian estimation X^k and P^k by(23)X^k=∫XkpXk∣Z1:kdXk,P^k=∫Xk-X^kXk-X^kTpXk∣Z1:kdXk.The calculation of the posterior density requires integrals, but the normalizing integral cannot be evaluated analytically and numerical integration over possibly high-dimensional spaces is infeasible. The ideal of SIS algorithm is that the required posterior density function is represented by a set of random samples with associated weight{xik,ωik}i=1Ns and to compute state estimates in the basis of these samples and weights. According to SIS algorithm, the posterior density pXk∣Z1:k can be approximated as(24)pXk∣Z1:k≈∑i=1NsωikδXk-xik,where the weights ωik are defined in(25)ωik∝ωik-1pZk∣xikpxik∣xik-1qxik∣xik-1,Zk,where qxik∣xik-1,Zk is termed the importance sampling proposal distribution. It is easy to draw samples from the distribution defined by qxik∣xik-1,Zk, but it is hard to draw samples directly from the posterior density pXk∣Z1:k. In standard Particle Filter, we set qxik∣xik-1,Zk as pxik∣xik-1:(26)qxik∣xik-1,Zk=pxik∣xik-1.Then(27)ωik∝ωik-1pZk∣xik.Assuming thatvk and w(k) are mutually independent zero-mean white Gaussian processes [7] in (1) and (5), the prior density pxik∣xik-1 and the likelihood density pZk∣xik become(28)pxik∣xik-1=Nxik;APNxik-1,Qk,pZk∣xik=NZk;hxik,Rk.


### 3.2.2. The Resampling Step

After a few iterations with SIS, the variance of the importance weights increases because most particles will have negligible weights. This is called particle degeneracy. Fixing on an importance sampling proposal distribution, one way to mitigate the particle degeneracy is to adopt the resampling step. The effective sample sizeNeff denotes an index variable of degeneracy introduced in [6]. Neff can be estimated approximately by(29)Neff≈1∑i=1Nsω~ik2.IfNeff falls below a threshold NT, the resampling algorithm will sample from (24) to generate a new set xi∗ki=1Ns and reset the weights ωi∗k to 1/Ns.The standard PF with the PN motion model to track the air-to-air missile is shown in Algorithm2.Algorithm 2
NT is the total tracking time
Initialize the variables whenk=1 and k=2
FORk=3:N
FORi=1:Ns
Drawxik~pxik∣xik-1
Assign the particle a weight according to (27)
END FOR
FORj=1:Ns
Normalize the weights by
ω~jk=ωjk/∑i=1Nsωik2
END FOR
EstimateX^k and P^k by
X^k=∑j=1Nsω~jkxjk
P^k=∑j=1Nsω~jkxjk-X^kxjk-X^kT
Calculate (22) to update APN
CalculateNeff using (29)
IFNeff<NT
Take resampling algorithm
End FOR


## 3.2.1. The Sequential Importance Sampling (SIS) Algorithm

The Particle Filter algorithm to deal with the nonlinear tracking problem has been used on a wide consensus [6]. Let Z1:k={Z1,Z2,…,Zk} denote a sequence of observations and pXk∣Z1:k denote the posterior density of Xk. We can acquire Bayesian estimation X^k and P^k by(23)X^k=∫XkpXk∣Z1:kdXk,P^k=∫Xk-X^kXk-X^kTpXk∣Z1:kdXk.The calculation of the posterior density requires integrals, but the normalizing integral cannot be evaluated analytically and numerical integration over possibly high-dimensional spaces is infeasible. The ideal of SIS algorithm is that the required posterior density function is represented by a set of random samples with associated weight{xik,ωik}i=1Ns and to compute state estimates in the basis of these samples and weights. According to SIS algorithm, the posterior density pXk∣Z1:k can be approximated as(24)pXk∣Z1:k≈∑i=1NsωikδXk-xik,where the weights ωik are defined in(25)ωik∝ωik-1pZk∣xikpxik∣xik-1qxik∣xik-1,Zk,where qxik∣xik-1,Zk is termed the importance sampling proposal distribution. It is easy to draw samples from the distribution defined by qxik∣xik-1,Zk, but it is hard to draw samples directly from the posterior density pXk∣Z1:k. In standard Particle Filter, we set qxik∣xik-1,Zk as pxik∣xik-1:(26)qxik∣xik-1,Zk=pxik∣xik-1.Then(27)ωik∝ωik-1pZk∣xik.Assuming thatvk and w(k) are mutually independent zero-mean white Gaussian processes [7] in (1) and (5), the prior density pxik∣xik-1 and the likelihood density pZk∣xik become(28)pxik∣xik-1=Nxik;APNxik-1,Qk,pZk∣xik=NZk;hxik,Rk.


## 3.2.2. The Resampling Step

After a few iterations with SIS, the variance of the importance weights increases because most particles will have negligible weights. This is called particle degeneracy. Fixing on an importance sampling proposal distribution, one way to mitigate the particle degeneracy is to adopt the resampling step. The effective sample sizeNeff denotes an index variable of degeneracy introduced in [6]. Neff can be estimated approximately by(29)Neff≈1∑i=1Nsω~ik2.IfNeff falls below a threshold NT, the resampling algorithm will sample from (24) to generate a new set xi∗ki=1Ns and reset the weights ωi∗k to 1/Ns.The standard PF with the PN motion model to track the air-to-air missile is shown in Algorithm2.Algorithm 2
NT is the total tracking time
Initialize the variables whenk=1 and k=2
FORk=3:N
FORi=1:Ns
Drawxik~pxik∣xik-1
Assign the particle a weight according to (27)
END FOR
FORj=1:Ns
Normalize the weights by
ω~jk=ωjk/∑i=1Nsωik2
END FOR
EstimateX^k and P^k by
X^k=∑j=1Nsω~jkxjk
P^k=∑j=1Nsω~jkxjk-X^kxjk-X^kT
Calculate (22) to update APN
CalculateNeff using (29)
IFNeff<NT
Take resampling algorithm
End FOR


## 3.3. The Genetic Algorithm Particle Filter

The standard Particle Filter can solve the particle degeneracy, but the resampling algorithm relying on the high weight particles leads to the loss diversity [7]. The evolution idea of Genetic Algorithm (GA) [8] can be adopted in the resampling step to solve the problem. The chromosomes and adaptability in the GA are corresponding to the particle and weight in the PF, respectively. So we can use the selection, intercross, and variation steps to deal with the particle in resampling step to increase the diversity.GA for the resampling step is shown in Algorithm3.Algorithm 3
IFNeff<NT
FOR i=1:N1
Chose randomly two particlesxak and xbk from xik,ωiki=1Ns to intercross by
x~ak=εxak+1-εxbkx~bk=εxbk+1-εxakω~ak=εωak+1-εωbkω~bk=εωbk+1-εωak
whereε~U0,1 and a,b∈1,2,…Ns
END FOR
Create a new setx~ik,ω~iki=1Ns
FORi=1:N2
Chose randomly a particlex~ck from x~ik,ω~iki=1Ns to vary by
x⌢ck=x~ck+η
whereη~NX^k,P^k
IfpZk∣x⌢ck>pZk∣x~ck, accept x⌢ck
Else ifpZk∣x⌢ck/pZk∣x~ck>u  u~U0,1, accept x⌢ck
Otherwise acceptx~ck
END IF
END FOR
Create a new setx⌢ik,ω⌢iki=1Ns
Else leap over the resampling step
END IFThe GAPF uses the GA to replace the standard resampling step in the PN and increases the diversity, which is validated by simulations.


## 4. Simulation Results

### 4.1. Hardware Condition

An Intel®Core™i3-2100@3.10 GHz CPU is used with 2.00 GB RAM and MATLAB 2010b.


### 4.2. The True Ballistic Trajectory of the Air-to-Air Missile

We assume that the aircraft attacked flies along a straight line at a constant velocityvf=200 m/s with the pitching angle θf=π/120 and the azimuth angle φf=π/6 at the original point 10000,500,10000 in north-sky-west geographic coordinate system. We take the firing point of the air-to-air missile as the coordinate origin 0,0,0 and assume that the missile moves following the Proportional Navigation law, whose proportional parameter Kθ=Kφ=4, at a constant speed vm=500 m/s with the pitching angle θm=0 and the azimuth angle φm=0. Through the missile dynamics equation set (refer to [3]), we can acquire the ballistic trajectory of the missile and the measurement set Z with a second sampling interval in Figure 2.Figure 2
The true ballistic trajectory and the measurements.


### 4.3. The Simulation Results

In Figure3, we show the missile state estimates filtered by the CV, CA, and PN models with EKF, respectively. In CV model, we assume that the original state is Xk=1=0,0,0,0,0,0; the estimate covariance matrix is Pk=1=diag⁡100,100,100,30,30,30; the state noise covariance matrix is Q=diag⁡0,0,0,50,50,50; the measurement noise covariance matrix is R=diag⁡200,0.04,0.04. In CA model, we assume that the original state and covariance matrix are Xk=1=0,0,0,0,0,0,0,0,0 and Pk=1=diag⁡100,100,100,30,30,30,10,10,10 separately; the state noise covariance matrix is Q=diag⁡0,0,0,0,0,0,30,30,30; the measurement noise covariance matrix is R=diag⁡200,0.04,0.04. In PN model, we assume that the proportional parameters is Kθk=1=Kθk=2=3 and Kφk=1=kk=2=2.5; the angles are initialized by the measurement vectors at k=1 and k=2. Figure 3 shows that the three models are all convergent, but the state estimate performances are different among them. Using 100 times Monte Carlo simulations for every model with EKF, the estimate effects are given in the form of the position Mean Square Error (RMSE) [21]. In Figure 4, we record the RMSE according to the models with EKF, for models comparison purposes. In Figure 5, we show that the average estimated values of proportional parameters Kθ and Kφ are close to the true values gradually. The plot shows that PN is the best model among them for tracking the air-to-air missile.Figure 3
The EKF track with the CV CA and PN.Figure 4
The RMSE with the CV CA and PN.Figure 5
Average proportional parameters.In Figure6, we adopt three filter algorithms, such as EKF, PF, and GAPF, to track the missile with the constant PN model. The PN-EKF algorithm is a comparison which is the same in Figures 3 and 4. The number of particles is Ns=3000 and the resampling threshold is NT=1000 in the PN-PF and PN-GAPF. The number of intercross particles is N1=0.3Ns and the number of variation particle is N2=0.005Ns in the PN-GAPF. Figure 6 shows that the GAPF performance is better than EKF and PF. The RSME curves in Figure 7 directly prove it. In Figure 8, we show that the average computation time of per step in the EKF is shortest among them, but GAPF needs more time in intercross and variation steps than the standard resampling algorithm in the PF. Further improvements in performance can be acquired by increasing the number of particles, but the computation time also becomes longer, which heavily leads to not meeting real-time tracking requirement.Figure 6
The track with PN-EKF PN-PF and PN-GAPF.Figure 7
The track RMSE with PN-EKF PN-PF and PN-GAPF.Figure 8
Average computation time of per step.


## 4.1. Hardware Condition

An Intel®Core™i3-2100@3.10 GHz CPU is used with 2.00 GB RAM and MATLAB 2010b.


## 4.2. The True Ballistic Trajectory of the Air-to-Air Missile

We assume that the aircraft attacked flies along a straight line at a constant velocityvf=200 m/s with the pitching angle θf=π/120 and the azimuth angle φf=π/6 at the original point 10000,500,10000 in north-sky-west geographic coordinate system. We take the firing point of the air-to-air missile as the coordinate origin 0,0,0 and assume that the missile moves following the Proportional Navigation law, whose proportional parameter Kθ=Kφ=4, at a constant speed vm=500 m/s with the pitching angle θm=0 and the azimuth angle φm=0. Through the missile dynamics equation set (refer to [3]), we can acquire the ballistic trajectory of the missile and the measurement set Z with a second sampling interval in Figure 2.Figure 2
The true ballistic trajectory and the measurements.


## 4.3. The Simulation Results

In Figure3, we show the missile state estimates filtered by the CV, CA, and PN models with EKF, respectively. In CV model, we assume that the original state is Xk=1=0,0,0,0,0,0; the estimate covariance matrix is Pk=1=diag⁡100,100,100,30,30,30; the state noise covariance matrix is Q=diag⁡0,0,0,50,50,50; the measurement noise covariance matrix is R=diag⁡200,0.04,0.04. In CA model, we assume that the original state and covariance matrix are Xk=1=0,0,0,0,0,0,0,0,0 and Pk=1=diag⁡100,100,100,30,30,30,10,10,10 separately; the state noise covariance matrix is Q=diag⁡0,0,0,0,0,0,30,30,30; the measurement noise covariance matrix is R=diag⁡200,0.04,0.04. In PN model, we assume that the proportional parameters is Kθk=1=Kθk=2=3 and Kφk=1=kk=2=2.5; the angles are initialized by the measurement vectors at k=1 and k=2. Figure 3 shows that the three models are all convergent, but the state estimate performances are different among them. Using 100 times Monte Carlo simulations for every model with EKF, the estimate effects are given in the form of the position Mean Square Error (RMSE) [21]. In Figure 4, we record the RMSE according to the models with EKF, for models comparison purposes. In Figure 5, we show that the average estimated values of proportional parameters Kθ and Kφ are close to the true values gradually. The plot shows that PN is the best model among them for tracking the air-to-air missile.Figure 3
The EKF track with the CV CA and PN.Figure 4
The RMSE with the CV CA and PN.Figure 5
Average proportional parameters.In Figure6, we adopt three filter algorithms, such as EKF, PF, and GAPF, to track the missile with the constant PN model. The PN-EKF algorithm is a comparison which is the same in Figures 3 and 4. The number of particles is Ns=3000 and the resampling threshold is NT=1000 in the PN-PF and PN-GAPF. The number of intercross particles is N1=0.3Ns and the number of variation particle is N2=0.005Ns in the PN-GAPF. Figure 6 shows that the GAPF performance is better than EKF and PF. The RSME curves in Figure 7 directly prove it. In Figure 8, we show that the average computation time of per step in the EKF is shortest among them, but GAPF needs more time in intercross and variation steps than the standard resampling algorithm in the PF. Further improvements in performance can be acquired by increasing the number of particles, but the computation time also becomes longer, which heavily leads to not meeting real-time tracking requirement.Figure 6
The track with PN-EKF PN-PF and PN-GAPF.Figure 7
The track RMSE with PN-EKF PN-PF and PN-GAPF.Figure 8
Average computation time of per step.


## 5. Conclusion

This paper aims to solve the problem of tracking the air-to-air missile. Firstly we establish the PN motion model based on the PN law to preferably track the air-to-air missile in 3d Cartesian coordinate system. Secondly, the nonlinear filter problem should be solved because of the nonlinear measurement equation in the state space model for tracking the missile. In the nonlinear filter algorithms, we take the standard EKF and PF to track the missile with the PN motion model. However, both filters’ results are not good because the linear approximate degree is inadequate for EKF and the diversity of particles is deficient for the standard resampling step. We introduce the intercross and variation steps in the GA to the resampling step in the standard PF for increasing the diversity. The simulation results show that the GAPN, which replaces the standard resampling step with the GA, can obtain a better performance for tracking the missile with PN. Meanwhile, the PN model proposed in this paper is a more suitable model for the air-to-air missile than CV and CA models.



---

*Source: 3921608-2016-03-20.xml*