# Multiple Model-Based Synchronization Approaches for Time Delayed Slaving Data in a Space Launch Vehicle Tracking System

Authors: Haryong Song; Yongtae Choi
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Due to the inherent characteristics of the flight mission of a space launch vehicle (SLV), which is required to fly over very large distances and have very high fault tolerances, in general, SLV tracking systems (TSs) comprise multiple heterogeneous sensors such as radars, GPS, INS, and electrooptical targeting systems installed over widespread areas. To track an SLV without interruption and to hand over the measurement coverage between TSs properly, the mission control system (MCS) transfers slaving data to each TS through mission networks. When serious network delays occur, however, the slaving data from the MCS can lead to the failure of the TS. To address this problem, in this paper, we propose multiple model-based synchronization (MMS) approaches, which take advantage of the multiple motion models of an SLV. Cubic spline extrapolation, prediction through anα-β-γ filter, and a single model Kalman filter are presented as benchmark approaches. We demonstrate the synchronization accuracy and effectiveness of the proposed MMS approaches using the Monte Carlo simulation with the nominal trajectory data of Korea Space Launch Vehicle-I.


---

## Body

## 1. Introduction

The range safety system (RSS) [1] for the flight mission of a space launch vehicle (SLV) consists of multiple heterogeneous tracking systems (TSs) with mission control systems (MCSs). Critical launch mission details such as time-space-position information (TSPI), launch mission status data, that is, quick look message (QLM), and flight safety information can be acquired from the RSS. Tracking an extensive mission trajectory of an SLV requires widespread multiple TSs so that the RSS covers the entire trajectory of the SLV flight mission. Generally, multiple TSs are spread out over different sites and they automatically hand over observation coverage according to the flight of the SLV. In this circumstance, one of the most important roles of the RSS is to distribute slaving data to each TS for continuous tracking of the SLV. If a critical network delay results in time delayed slaving data to be sent to the TSs, the MCS will not receive accurate SLV tracking data. This problem can lead to significant difficulties for the SLV mission progress and analysis. Therefore, the motivation of this research is to enhance slaving data accuracy by compensating for possible network time delays. The basic solution to this problem is simple (linear or nonlinear) extrapolation of the filtered data. In this case, extrapolation methods simply propagate the TSPI data without regard to the system dynamics of the SLV. On the other hand, a Kalman filter (KF) and its prediction capability [2, 3] can reflect the system dynamics of the SLV, which results in better synchronization performance. However, since a KF only utilizes a single dynamic model, in general, the tracking performance of a KF for a maneuvering target is inferior to multiple model estimators [4]. In addition, due to stage separation, the flight phase of an SLV is separated into two parts, the propelled flight phase (PFP) and the coasting flight phase (CFP). Hence, the dynamic model of an SLV can be described using multiple models so that the multiple flight phases are properly accounted for. To adaptively select one of the multiple dynamic models according to the flight phase of the SLV, multiple model estimators (MME) such as an interacting multiple model (IMM) [5] and a multiple model adaptive estimator (MMAE) [6] could significantly improve SLV tracking when a network delay results in delayed slaving data transmission.In the past several decades, considerable research has been undertaken in the field of launch vehicle tracking based on multiple dynamic models; researchers have shown interest in various applications such as the tracking of reentry vehicles, short-range projectiles, and sounding rockets [7–13]. A reentry vehicle tracking problem known as highly nonlinear dynamics was conducted using a modified IMM, with a different algorithm cycle compared with an ordinary IMM with multiple modes of diverse ballistic coefficients [7]. Short-range ballistic munitions or projectiles with multiple models such as spin-stabilized and fin-stabilized models were implemented using an IMM [8, 9]. Both research alternatives can be applied in the impact point prediction of projectiles. Tactical ballistic missile tracking was carried out using an IMM estimator with three modes [10]. The first mode was a constant axial force model for the boosting and reentry phases. The second mode was a ballistic acceleration model that incorporated the gravitational, Coriolis, and centripetal forces for the exoatmospheric phase. The final mode was a standard autocorrelated acceleration Singer’s model for malfunction motions of missiles such as reentry tumbling. A multiple IMM algorithm with an unbiased mixing approach for multiple modes of thrusting or for ballistic projectiles was presented [11, 12]. A sounding rocket with multiple modes of propelled flight or free fall flight was tracked using a multiple model adaptive estimation approach [14].In this paper, multiple model-based synchronization (MMS) approaches are proposed to synchronize the time delayed slaving data of the RSS. The proposed approaches can be expressed via two distinct multiple models, a nonlinear model and a linear model. The nonlinear model considers comprehensive factors such as thrust, gravity, drag coefficient, Mach number, and air density [10–13]. Although the nonlinear model precisely describes the motion of the SLV, it requires complex information concerning the SLV specifications to be collected in advance. In contrast, in the case of the linear dynamic model, a simple constant velocity (CV) or constant acceleration (CA) model with multiple hypotheses, which takes advantage of Singer’s model [16, 17], can be utilized [14]. To describe the motion of the SLV, the motion modes of both models are separated into two parts, PFP and CFP. We propose a slaving data synchronization approach for the RSS based on MME so that the MCS can adaptively find an appropriate dynamic model at an arbitrary time index, where time delay has occurred. The performance of slaving data synchronization is compared to various benchmark methods such as cubic spline extrapolation, prediction through an α-β-γ filter, and a single model KF to demonstrate the effectiveness of the proposed algorithm.The remainder of the paper is organized as follows. Section2 presents a statement of the problem for delayed slaving data in RSS. In Section 3, conventional synchronization approaches are illustrated. In Section 4, two proposed MME-based synchronization approaches are derived. Section 5 presents simulation settings and results; the comparison between different types of synchronization approaches is depicted as an aspect of RMS error of the state vector. Finally, in Section 6, the conclusions of this paper are presented.


## 2. Problem Statement for Delayed Slaving Data in RSS

The transmission of slaving data from the MCS to multiple TSs facilitates seamless tracking of the SLV in a sparsely located multiple TS environment. If a data transmission delay problem occurs, it can cause an SLV tracking failure. As depicted in Figure1, the antennas of TSs are pointing at the SLV by controlling their attitude according to slaving data from the MCS. In this situation, a slight time delay in slaving data can give rise to large differences between the antenna beam and the SLV due to the fast movement of the SLV. To solve this problem, we propose MMS approaches. Hence, the goal of this paper is to find a synchronized state X
^
k
+
s
s
y
n
c at time k
+
s (wheres is a known delay) based on delayed slaving data X
^
k
d
e
l
a
y at timek such that(1)
X
^
k
+
s
s
y
n
c
=
f
p
s
,
X
^
k
d
e
l
a
y
,
X
=
x
y
z
T
,where f
p is a linear or nonlinear propagation function and X
=
x
y
z
T is a slaving state vector that is composed ofx-axis,y-axis, andz-axis positions but is not limited to the position components.Figure 1
Illustration of problem statement for delayed slaving data in RSS.


## 3. Conventional Synchronization Approaches for Slaving Data

### 3.1. Synchronization Using Cubic Spline Extrapolation [18, 19]

Let us assume a synchronized slaving state vector to be an unknown function of the delayed slaving state vector whose values are known only until timek. We then define a cubic spline extrapolation function f
n
e
p
(
k
+
s
), where n
=
x
,
y
,
z at a specific synchronization time k
+
s that is extrapolated based on the known function f
n
k that has a real value, with N
+
1 points, where k
0
≤
k
≤
k
N.We approximatef
n
e
p
k
+
s as a three-order polynomial based on the interval k
i
,
k
i
+
1, where i
=
0
,
…
,
N
-
1. Then f
n
e
p
k
+
s can be defined as follows:(2)
f
n
e
p
k
+
s
=
f
n
k
,
k
∈
k
i
,
k
i
+
1
,
i
=
0
,
…
,
N
-
1
,
f
n
k
=
a
n
,
i
k
3
+
b
n
,
i
k
2
+
c
n
,
i
k
+
d
n
,
i
,where a
n
,
i, b
n
,
i, c
n
,
i, and d
n
,
i are unknown coefficients. To find the unknown coefficients, we make the assumption that f
n
k should be continuous in k
0
,
k
N. Therefore, the following equation containing unknown coefficients a
n
,
i, b
n
,
i, c
n
,
i, and d
n
,
i is obtained: (3)
f
n
,
i
k
i
=
f
n
,
i
+
1
k
i
f
n
,
i
′
k
i
=
f
n
,
i
+
1
′
k
i
,
f
n
,
i
′
′
k
i
=
f
n
,
i
+
1
′
′
k
i
n
=
x
,
y
,
z
,
i
=
1
,
…
,
N
-
1
.On combining (2) and (3), the cubic polynomials f
n
e
p
(
k
+
s
) are reconstructed by solving the linear equations obtained. Once we find the coefficients a
n
,
i, b
n
,
i, c
n
,
i, and d
n
,
i, where i
=
0
,
…
,
N
-
1, we can evaluate f
n
e
p
k
+
s, where s is an arbitrary lead-time for future points in x
0
,
x
N
+
s.


### 3.2. Synchronization Using anα-β-γ Filter

When the state estimation covariance for a time invariant system converges under suitable conditions to a steady-state value, explicit expressions of the steady-state covariance and filter gain can be obtained. The resulting steady-state filters for noisy kinematic models are known asα-β and α-β-γ filters [20]. In this paper, a combined α-β-γ filter using both an expanding memory polynomial filter (EMF) and a fading memory polynomial filter (FMF) was used [21, 22]. At first, the EMF is represented as follows:(4)
p
n
,
k
E
=
p
n
,
k
-
1
E
+
Δ
t
v
n
,
k
-
1
E
+
Δ
t
2
2
a
n
,
k
-
1
E
+
3
3
k
2
+
3
k
+
2
k
+
3
k
+
2
k
+
1
ε
n
,
k
E
,
v
n
,
k
E
=
v
n
,
k
-
1
E
+
Δ
a
v
n
,
k
-
1
E
+
1
Δ
t
18
2
k
+
1
k
+
3
k
+
2
k
+
1
ε
n
,
k
E
,
a
n
,
k
E
=
a
n
,
k
-
1
E
+
1
Δ
t
2
60
k
+
3
k
+
2
k
+
1
ε
n
,
k
E
,
ε
n
,
k
E
=
y
n
,
k
-
p
n
,
k
-
1
E
-
Δ
t
v
n
,
k
-
1
E
-
Δ
t
2
2
a
n
,
k
-
1
E
,where Δ
t is the sampling time and p
n
,
k
E, v
n
,
k
E, and a
n
,
k
E are the position, velocity, and acceleration estimates of the EMF, respectively. In addition, the variance reduction factor (VRF) [22] of the EMF, that is, V
R
F
E, can be represented as(5)
V
R
F
E
=
9
k
2
+
27
k
+
24
k
k
+
1
k
-
1
.On the other hand, the FMF and its VRF (V
R
F
F) can be written as follows: (6)
p
n
,
k
F
=
p
n
,
k
-
1
F
+
Δ
t
v
n
,
k
-
1
F
+
Δ
t
2
2
a
n
,
k
-
1
F
+
1
-
λ
3
ε
n
,
k
F
,
v
n
,
k
F
=
v
n
,
k
-
1
F
+
Δ
a
v
n
,
k
-
1
F
+
3
2
Δ
t
1
-
λ
2
1
+
λ
ε
n
,
k
F
,
a
n
,
k
F
=
a
n
,
k
-
1
F
+
1
Δ
t
2
1
-
λ
3
ε
n
,
k
F
,
ε
n
,
k
F
=
y
n
,
k
-
p
n
,
k
-
1
F
-
Δ
t
v
n
,
k
-
1
F
-
Δ
t
2
2
a
n
,
k
-
1
F
,
V
R
F
F
=
1
-
λ
1
+
λ
5
19
+
24
λ
+
16
λ
2
+
6
λ
3
+
λ
4
,where 0
<
λ
<
1. Both filters are conducted in parallel but only one of them is selected as a final estimate by comparing the VRFs. During the early tracking phase, the EMF is selected as the final estimate; however, at a certain time, where V
R
F
F is larger than V
R
F
E, FMF is selected as the final estimate: (7)
X
^
n
,
k
α
β
γ
=
f
n
X
^
n
,
k
-
1
α
β
γ
+
K
y
k
-
H
α
β
γ
X
^
n
,
k
-
1
α
β
γ
,
p
n
,
k
v
n
,
k
a
n
,
k
=
1
Δ
t
Δ
t
2
2
0
1
Δ
t
0
0
1
p
n
,
k
-
1
v
n
,
k
-
1
a
n
,
k
-
1
+
α
β
Δ
t
γ
Δ
t
2
y
k
-
H
α
β
γ
p
n
,
k
-
1
v
n
,
k
-
1
a
n
,
k
-
1
,
where
X
n
,
k
α
β
γ
=
p
n
,
k
v
n
,
k
a
n
,
k
T
,
H
α
β
γ
=
1
0
0
,
V
R
F
E
≥
V
R
F
F
⟶
p
n
,
k
=
p
n
,
k
E
,
v
n
,
k
=
v
n
,
k
E
,
a
n
,
k
=
a
n
,
k
E
,
V
R
F
F
>
V
R
F
E
⟶
p
n
,
k
=
p
n
,
k
F
,
v
n
,
k
=
v
n
,
k
F
,
a
n
,
k
=
a
n
,
k
F
.Finally, synchronization using an α-β-γ filter is completed by linear propagation using the system matrix f
n such that(8)
X
^
n
,
k
+
s
α
β
γ
=
∏
L
=
1
s
f
n
·
X
^
n
,
k
α
β
γ
.


### 3.3. Kalman Predictor

The motion of the SLV is simply depicted as a discretized Wiener process acceleration model [20]: (9)
X
k
+
1
=
F
X
k
+
w
k
,where the state vector X
k
∈
R
9 consists of the position, velocity, and acceleration components along x-axis, y-axis, and z-axis, respectively; that is, X
k
=
x
p
,
k
x
v
,
k
x
a
,
k
y
p
,
k
y
v
,
k
y
a
,
k
z
p
,
k
z
v
,
k
z
a
,
k
T. The system matrix and covariance matrix of the system noise w
k are represented as (10) and (11), respectively:(10)
f
n
=
x
,
y
,
z
=
1
Δ
t
Δ
t
2
2
0
1
Δ
t
0
0
1
,
F
=
f
x
0
3
0
3
0
3
f
y
0
3
0
3
0
3
f
z
,
(11)
E
w
k
w
l
T
=
Q
δ
k
-
l
,
q
L
=
x
,
y
,
z
=
Δ
t
5
20
Δ
t
4
8
Δ
t
2
6
Δ
t
4
8
Δ
t
3
3
Δ
t
2
2
Δ
t
3
6
Δ
t
2
2
Δ
t
,
Q
=
q
x
0
3
0
3
0
3
q
y
0
3
0
3
0
3
q
z
.A radar measurement for the SLV gives the spherical coordinate observations such that (12)
Z
k
=
h
X
k
+
v
s
,
k
,
Z
k
∈
R
3
,
v
s
,
k
∈
R
3
,
Z
k
=
r
k
φ
k
θ
k
=
x
p
,
k
2
+
y
p
,
k
2
+
z
p
,
k
2
tan
-
1
⁡
y
p
,
k
x
p
,
k
tan
-
1
⁡
y
p
,
k
x
p
,
k
2
+
y
p
,
k
2
+
v
r
v
φ
v
θ
,where v
s
,
k
~
N
0
,
R
s
,
k and R
s
,
k
=
d
i
a
g
σ
r
2
,
σ
φ
2
,
σ
θ
2. Using a 3D debiased converted measurement [23], we can transform the original nonlinear equations (12) into their linear form as (13)
z
k
=
H
X
k
+
v
k
,
v
k
∈
R
3
,
z
k
=
x
k
y
k
z
k
=
1
0
1
×
5
0
1
×
3
0
1
×
3
1
0
1
×
5
0
1
×
5
1
0
1
×
2
X
k
+
v
x
,
k
v
y
,
k
v
z
,
k
.Here, v
c
,
k is the converted measurement noise expressed in terms of Cartesian coordinates; that is, v
k
~
N
0
,
R
k:(14)
R
k
=
R
x
x
,
k
R
x
y
,
k
R
x
z
,
k
R
y
x
,
k
R
y
y
,
k
R
y
z
,
k
R
z
x
,
k
R
z
y
,
k
R
z
z
,
k
.Prediction or fixed-lead prediction in mean square means the synchronization of the slaving data is the estimation of the state at a future time k
+
s, where s
>
0 beyond the observation interval; that is, based on data up to an earlier time [20, 24], (15)
X
^
k
+
s
∣
k
=
E
X
k
+
s
∣
k
∣
Z
k
,
Z
k
=
z
0
,
…
,
z
k
.The optimal predictor or synchronized state X
^
k
+
s
∣
k
≜
X
^
k
+
s
∣
k
K
P and its error covariance P
k
+
s
∣
k
≜
P
k
+
s
∣
k
K
P are given by the Kalman predictor equations [2, 24]:(16)
X
^
k
+
s
∣
k
K
P
=
F
k
+
s
-
1
X
^
k
+
s
-
1
∣
k
K
P
=
⋯
=
Φ
k
+
s
,
k
X
^
k
∣
k
K
P
,
X
^
k
∣
k
K
P
≜
X
^
k
∣
k
K
F
,
P
k
+
s
∣
k
K
P
=
F
k
+
s
-
1
P
k
+
s
-
1
∣
k
K
P
F
k
+
s
-
1
T
+
Q
k
+
s
-
1
=
Φ
k
+
s
,
k
P
k
∣
k
K
P
Φ
k
+
s
,
k
T
+
∑
j
=
0
s
-
1
Φ
k
+
s
,
k
+
j
+
1
Q
k
+
j
Φ
k
+
s
,
k
+
j
+
1
T
,
P
k
∣
k
K
P
≜
P
k
∣
k
K
F
,
Φ
k
+
s
,
k
=
F
k
+
s
-
1
F
k
+
s
-
2
⋯
F
k
,
Φ
k
,
k
=
I
n
,
k
>
1
,where X
^
k
∣
k
K
F and P
k
∣
k
K
F are the KF estimate and covariance, respectively.


## 3.1. Synchronization Using Cubic Spline Extrapolation [18, 19]

Let us assume a synchronized slaving state vector to be an unknown function of the delayed slaving state vector whose values are known only until timek. We then define a cubic spline extrapolation function f
n
e
p
(
k
+
s
), where n
=
x
,
y
,
z at a specific synchronization time k
+
s that is extrapolated based on the known function f
n
k that has a real value, with N
+
1 points, where k
0
≤
k
≤
k
N.We approximatef
n
e
p
k
+
s as a three-order polynomial based on the interval k
i
,
k
i
+
1, where i
=
0
,
…
,
N
-
1. Then f
n
e
p
k
+
s can be defined as follows:(2)
f
n
e
p
k
+
s
=
f
n
k
,
k
∈
k
i
,
k
i
+
1
,
i
=
0
,
…
,
N
-
1
,
f
n
k
=
a
n
,
i
k
3
+
b
n
,
i
k
2
+
c
n
,
i
k
+
d
n
,
i
,where a
n
,
i, b
n
,
i, c
n
,
i, and d
n
,
i are unknown coefficients. To find the unknown coefficients, we make the assumption that f
n
k should be continuous in k
0
,
k
N. Therefore, the following equation containing unknown coefficients a
n
,
i, b
n
,
i, c
n
,
i, and d
n
,
i is obtained: (3)
f
n
,
i
k
i
=
f
n
,
i
+
1
k
i
f
n
,
i
′
k
i
=
f
n
,
i
+
1
′
k
i
,
f
n
,
i
′
′
k
i
=
f
n
,
i
+
1
′
′
k
i
n
=
x
,
y
,
z
,
i
=
1
,
…
,
N
-
1
.On combining (2) and (3), the cubic polynomials f
n
e
p
(
k
+
s
) are reconstructed by solving the linear equations obtained. Once we find the coefficients a
n
,
i, b
n
,
i, c
n
,
i, and d
n
,
i, where i
=
0
,
…
,
N
-
1, we can evaluate f
n
e
p
k
+
s, where s is an arbitrary lead-time for future points in x
0
,
x
N
+
s.


## 3.2. Synchronization Using anα-β-γ Filter

When the state estimation covariance for a time invariant system converges under suitable conditions to a steady-state value, explicit expressions of the steady-state covariance and filter gain can be obtained. The resulting steady-state filters for noisy kinematic models are known asα-β and α-β-γ filters [20]. In this paper, a combined α-β-γ filter using both an expanding memory polynomial filter (EMF) and a fading memory polynomial filter (FMF) was used [21, 22]. At first, the EMF is represented as follows:(4)
p
n
,
k
E
=
p
n
,
k
-
1
E
+
Δ
t
v
n
,
k
-
1
E
+
Δ
t
2
2
a
n
,
k
-
1
E
+
3
3
k
2
+
3
k
+
2
k
+
3
k
+
2
k
+
1
ε
n
,
k
E
,
v
n
,
k
E
=
v
n
,
k
-
1
E
+
Δ
a
v
n
,
k
-
1
E
+
1
Δ
t
18
2
k
+
1
k
+
3
k
+
2
k
+
1
ε
n
,
k
E
,
a
n
,
k
E
=
a
n
,
k
-
1
E
+
1
Δ
t
2
60
k
+
3
k
+
2
k
+
1
ε
n
,
k
E
,
ε
n
,
k
E
=
y
n
,
k
-
p
n
,
k
-
1
E
-
Δ
t
v
n
,
k
-
1
E
-
Δ
t
2
2
a
n
,
k
-
1
E
,where Δ
t is the sampling time and p
n
,
k
E, v
n
,
k
E, and a
n
,
k
E are the position, velocity, and acceleration estimates of the EMF, respectively. In addition, the variance reduction factor (VRF) [22] of the EMF, that is, V
R
F
E, can be represented as(5)
V
R
F
E
=
9
k
2
+
27
k
+
24
k
k
+
1
k
-
1
.On the other hand, the FMF and its VRF (V
R
F
F) can be written as follows: (6)
p
n
,
k
F
=
p
n
,
k
-
1
F
+
Δ
t
v
n
,
k
-
1
F
+
Δ
t
2
2
a
n
,
k
-
1
F
+
1
-
λ
3
ε
n
,
k
F
,
v
n
,
k
F
=
v
n
,
k
-
1
F
+
Δ
a
v
n
,
k
-
1
F
+
3
2
Δ
t
1
-
λ
2
1
+
λ
ε
n
,
k
F
,
a
n
,
k
F
=
a
n
,
k
-
1
F
+
1
Δ
t
2
1
-
λ
3
ε
n
,
k
F
,
ε
n
,
k
F
=
y
n
,
k
-
p
n
,
k
-
1
F
-
Δ
t
v
n
,
k
-
1
F
-
Δ
t
2
2
a
n
,
k
-
1
F
,
V
R
F
F
=
1
-
λ
1
+
λ
5
19
+
24
λ
+
16
λ
2
+
6
λ
3
+
λ
4
,where 0
<
λ
<
1. Both filters are conducted in parallel but only one of them is selected as a final estimate by comparing the VRFs. During the early tracking phase, the EMF is selected as the final estimate; however, at a certain time, where V
R
F
F is larger than V
R
F
E, FMF is selected as the final estimate: (7)
X
^
n
,
k
α
β
γ
=
f
n
X
^
n
,
k
-
1
α
β
γ
+
K
y
k
-
H
α
β
γ
X
^
n
,
k
-
1
α
β
γ
,
p
n
,
k
v
n
,
k
a
n
,
k
=
1
Δ
t
Δ
t
2
2
0
1
Δ
t
0
0
1
p
n
,
k
-
1
v
n
,
k
-
1
a
n
,
k
-
1
+
α
β
Δ
t
γ
Δ
t
2
y
k
-
H
α
β
γ
p
n
,
k
-
1
v
n
,
k
-
1
a
n
,
k
-
1
,
where
X
n
,
k
α
β
γ
=
p
n
,
k
v
n
,
k
a
n
,
k
T
,
H
α
β
γ
=
1
0
0
,
V
R
F
E
≥
V
R
F
F
⟶
p
n
,
k
=
p
n
,
k
E
,
v
n
,
k
=
v
n
,
k
E
,
a
n
,
k
=
a
n
,
k
E
,
V
R
F
F
>
V
R
F
E
⟶
p
n
,
k
=
p
n
,
k
F
,
v
n
,
k
=
v
n
,
k
F
,
a
n
,
k
=
a
n
,
k
F
.Finally, synchronization using an α-β-γ filter is completed by linear propagation using the system matrix f
n such that(8)
X
^
n
,
k
+
s
α
β
γ
=
∏
L
=
1
s
f
n
·
X
^
n
,
k
α
β
γ
.


## 3.3. Kalman Predictor

The motion of the SLV is simply depicted as a discretized Wiener process acceleration model [20]: (9)
X
k
+
1
=
F
X
k
+
w
k
,where the state vector X
k
∈
R
9 consists of the position, velocity, and acceleration components along x-axis, y-axis, and z-axis, respectively; that is, X
k
=
x
p
,
k
x
v
,
k
x
a
,
k
y
p
,
k
y
v
,
k
y
a
,
k
z
p
,
k
z
v
,
k
z
a
,
k
T. The system matrix and covariance matrix of the system noise w
k are represented as (10) and (11), respectively:(10)
f
n
=
x
,
y
,
z
=
1
Δ
t
Δ
t
2
2
0
1
Δ
t
0
0
1
,
F
=
f
x
0
3
0
3
0
3
f
y
0
3
0
3
0
3
f
z
,
(11)
E
w
k
w
l
T
=
Q
δ
k
-
l
,
q
L
=
x
,
y
,
z
=
Δ
t
5
20
Δ
t
4
8
Δ
t
2
6
Δ
t
4
8
Δ
t
3
3
Δ
t
2
2
Δ
t
3
6
Δ
t
2
2
Δ
t
,
Q
=
q
x
0
3
0
3
0
3
q
y
0
3
0
3
0
3
q
z
.A radar measurement for the SLV gives the spherical coordinate observations such that (12)
Z
k
=
h
X
k
+
v
s
,
k
,
Z
k
∈
R
3
,
v
s
,
k
∈
R
3
,
Z
k
=
r
k
φ
k
θ
k
=
x
p
,
k
2
+
y
p
,
k
2
+
z
p
,
k
2
tan
-
1
⁡
y
p
,
k
x
p
,
k
tan
-
1
⁡
y
p
,
k
x
p
,
k
2
+
y
p
,
k
2
+
v
r
v
φ
v
θ
,where v
s
,
k
~
N
0
,
R
s
,
k and R
s
,
k
=
d
i
a
g
σ
r
2
,
σ
φ
2
,
σ
θ
2. Using a 3D debiased converted measurement [23], we can transform the original nonlinear equations (12) into their linear form as (13)
z
k
=
H
X
k
+
v
k
,
v
k
∈
R
3
,
z
k
=
x
k
y
k
z
k
=
1
0
1
×
5
0
1
×
3
0
1
×
3
1
0
1
×
5
0
1
×
5
1
0
1
×
2
X
k
+
v
x
,
k
v
y
,
k
v
z
,
k
.Here, v
c
,
k is the converted measurement noise expressed in terms of Cartesian coordinates; that is, v
k
~
N
0
,
R
k:(14)
R
k
=
R
x
x
,
k
R
x
y
,
k
R
x
z
,
k
R
y
x
,
k
R
y
y
,
k
R
y
z
,
k
R
z
x
,
k
R
z
y
,
k
R
z
z
,
k
.Prediction or fixed-lead prediction in mean square means the synchronization of the slaving data is the estimation of the state at a future time k
+
s, where s
>
0 beyond the observation interval; that is, based on data up to an earlier time [20, 24], (15)
X
^
k
+
s
∣
k
=
E
X
k
+
s
∣
k
∣
Z
k
,
Z
k
=
z
0
,
…
,
z
k
.The optimal predictor or synchronized state X
^
k
+
s
∣
k
≜
X
^
k
+
s
∣
k
K
P and its error covariance P
k
+
s
∣
k
≜
P
k
+
s
∣
k
K
P are given by the Kalman predictor equations [2, 24]:(16)
X
^
k
+
s
∣
k
K
P
=
F
k
+
s
-
1
X
^
k
+
s
-
1
∣
k
K
P
=
⋯
=
Φ
k
+
s
,
k
X
^
k
∣
k
K
P
,
X
^
k
∣
k
K
P
≜
X
^
k
∣
k
K
F
,
P
k
+
s
∣
k
K
P
=
F
k
+
s
-
1
P
k
+
s
-
1
∣
k
K
P
F
k
+
s
-
1
T
+
Q
k
+
s
-
1
=
Φ
k
+
s
,
k
P
k
∣
k
K
P
Φ
k
+
s
,
k
T
+
∑
j
=
0
s
-
1
Φ
k
+
s
,
k
+
j
+
1
Q
k
+
j
Φ
k
+
s
,
k
+
j
+
1
T
,
P
k
∣
k
K
P
≜
P
k
∣
k
K
F
,
Φ
k
+
s
,
k
=
F
k
+
s
-
1
F
k
+
s
-
2
⋯
F
k
,
Φ
k
,
k
=
I
n
,
k
>
1
,where X
^
k
∣
k
K
F and P
k
∣
k
K
F are the KF estimate and covariance, respectively.


## 4. Proposed MME-Based Synchronization Approaches

### 4.1. Synchronization Using an IMM with Singer’s Linear Model

Singer [16] described 2D manned maneuvering targets in range-bearing coordinates. This model can be adapted to the SLV kinematics in 3D Cartesian coordinates with multiple flight phases: (17)
p
˙
n
v
˙
n
a
˙
n
=
F
p
n
v
n
a
n
+
G
w
n
,
F
=
0
1
0
0
0
1
0
0
-
α
n
,
G
=
0
0
0
,where w
n
∈
R
3
×
1 is a white noise process along the Cartesian axis n
=
x
,
y
,
z
. The parameter β
=
1
/
α
n is the maneuver correlation time constant, and σ
a
n
2 is the acceleration variance describing maneuver intensity. In a steady state,(18)
σ
w
n
2
=
2
α
σ
a
n
2
.To describe SLV kinematics, the model must cope with crucial nonzero mean acceleration maneuvers during the propelled phase. In addition, after each stage’s engine burns out, a multiple model approach is applied to the coasting flight; that is, one model describes PFP, whereas the other depicts the CFP. Empirically tuned, independent probability density functions (PDFs) represented by TUM describe the accelerations of the SLV in the local coordinate frame. Figure 2 shows the means and variances of the acceleration processes of the SLV in this paper. The discrete-time model with state transition matrix Ψ
n
α
n
,
Δ
t is as follows:(19)
p
n
,
k
+
1
v
n
,
k
+
1
a
n
,
k
+
1
=
Ψ
n
α
n
,
Δ
t
p
n
,
k
v
n
,
k
a
n
,
k
+
w
n
,
k
,
Ψ
n
α
n
,
Δ
t
=
e
F
Δ
t
1
Δ
t
α
n
Δ
t
-
1
+
e
-
α
n
Δ
t
0
1
1
-
e
-
α
n
Δ
t
α
n
0
0
e
-
α
n
Δ
t
,
w
n
,
k
=
∫
k
Δ
t
k
+
1
Δ
t
e
F
k
+
1
Δ
t
-
r
G
w
n
d
r
.In the case of the PFP, w
n
,
k is a nonzero mean white noise sequence caused by the nonzero mean acceleration μ
a
n seen in Figure 2. A nonzero mean white noise sequence for the PFP should be considered in the target kinematics when implementing the state propagation stage in the KF. Thus, the deterministic input u
n
,
k caused by w
n
,
k along x-axis, y-axis, and z-axis is derived as follows: (20)
u
x
,
k
=
E
w
x
,
k
=
E
∫
k
Δ
t
k
+
1
Δ
t
e
F
k
+
1
Δ
t
-
r
G
w
x
d
r
≈
E
∫
k
Δ
t
k
+
1
Δ
t
e
-
α
x
k
+
1
Δ
t
-
r
0
0
w
x
d
r
,
u
x
,
k
=
μ
a
x
α
n
1
-
e
-
α
x
Δ
t
0
0
.
y-axis and z-axis can be derived in the same manner as shown in (20). The maneuver excitation covariance [10], which represents the uncertainty of the SLV kinematics model, is (21)
Q
n
,
k
α
n
,
Δ
t
=
E
w
n
,
k
-
E
w
n
,
k
w
n
,
k
-
E
w
n
,
k
T
=
2
α
n
σ
a
n
2
q
11
q
12
q
13
q
12
q
22
q
23
q
13
q
23
q
33
.The specific components of q
11
,
…
,
q
33 are illustrated in [25]. In addition, the measurement matrix H
L for Singer’s model can be depicted as (22)
z
L
,
k
=
H
L
x
L
,
k
+
ν
k
,
H
L
=
1
0
0
0
0
3
×
2
1
0
3
×
2
0
0
3
×
2
0
0
1
,
x
L
,
k
=
x
k
x
˙
k
x
¨
k
y
k
y
˙
k
y
¨
k
z
k
z
˙
k
z
¨
k
T
,where ν
k is measurement noise (as shown in (14)) with error covariance.Figure 2
Nonzero mean acceleration PDF in the PFP model (a) alongx-axis, (b) along y-axis, and (c) along z-axis. Zero mean acceleration PDF in the CFP Model (d) along x-axis, y-axis, and z-axis.
(a)
(b)
(c)
(d)


#### 4.1.1. IMM with Singer’s Linear Model

From (19)–(22) in Section 4.1, we can rewrite the Markov jump linear systems, where the ith model of the finite multiple model set M
=
{
m
1
,
…
,
m
M
} obeys the following equations:(23)
x
L
,
k
+
1
=
Ψ
k
i
x
L
,
k
+
w
k
i
,
z
L
,
k
=
H
L
x
L
,
k
+
ν
k
,where(24)
c
o
v
w
k
i
=
Q
k
i
,
c
o
v
ν
k
=
R
k
,
Ψ
k
i
=
d
i
a
g
Ψ
x
,
k
α
x
i
,
Δ
t
,
Ψ
y
,
k
α
y
i
,
Δ
t
,
Ψ
z
,
k
α
z
i
,
Δ
t
,
Q
k
i
=
d
i
a
g
Q
x
,
k
α
x
i
,
Δ
t
,
Q
y
,
k
α
y
i
,
Δ
t
,
Q
z
,
k
α
z
i
,
Δ
t
.Here, σ
a
n
2 for the SLV can be represented by TUM as in Figure 2. The PDFs of Figures 2(a)–2(d) are experimentally sampled from the nominal acceleration profile of the SLV. The superscript (i) denotes quantities pertinent to model m
(
i
) in M, and the jumps of the system mode are assumed to have transition probabilities:(25)
Pr
⁡
m
k
+
1
j
∣
m
k
i
≜
π
i
j
,where m
k
i denotes the event in which model m
i matches the system mode in effect at time k. In our application, M
=
2, and m
1 and m
2 are the propelled and coasting modes, respectively.Finally, complete recursion of the IMM with mode matched KF for the SLV tracking is summarized as follows:(i)
Model-conditioned reinitialization (fori
=
1,2
,
…
,
M):
(a)
predicted mode probability:μ
L
,
k
∣
k
-
1
i
≜
Pr
⁡
{
m
k
i
∣
z
L
,
k
-
1
}
=
∑
j
π
j
i
μ
L
,
k
-
1
j,
(b)
mixing weight:μ
L
,
k
-
1
j
∣
i
≜
Pr
⁡
{
m
k
-
1
j
∣
m
k
i
,
z
L
,
k
-
1
}
=
∑
j
π
j
i
μ
L
,
k
-
1
j
/
μ
L
,
k
-
1
i,
(c)
mixing estimate and covariance:(26)
x
-
L
,
k
-
1
∣
k
-
1
i
≜
E
x
L
,
k
-
1
∣
k
-
1
∣
m
k
i
,
z
L
,
k
-
1
=
∑
j
x
^
L
,
k
-
1
∣
k
-
1
j
μ
L
,
k
-
1
j
∣
i
,
P
-
L
,
k
-
1
∣
k
-
1
i
=
∑
j
P
L
,
k
-
1
∣
k
-
1
i
+
x
-
L
,
k
-
1
∣
k
-
1
i
-
x
^
L
,
k
-
1
∣
k
-
1
j
x
-
L
,
k
-
1
∣
k
-
1
i
-
x
^
L
,
k
-
1
∣
k
-
1
j
T
μ
L
,
k
-
1
j
∣
i
.
(ii)
Model-conditioned filtering (fori
=
1,2
,
…
,
M):
(a)
predicted estimate and covariance:(27)
x
^
L
,
k
∣
k
-
1
i
=
Ψ
k
-
1
i
x
-
L
,
k
-
1
∣
k
-
1
i
+
u
k
-
1
i
,
where
u
k
-
1
1
=
u
x
,
k
-
1
T
u
y
,
k
-
1
T
u
z
,
k
-
1
T
T
,
u
k
-
1
2
=
0
3
×
1
,
P
L
,
k
∣
k
-
1
i
=
Ψ
k
-
1
i
P
-
L
,
k
-
1
∣
k
-
1
i
Ψ
k
-
1
i
T
+
Q
k
-
1
i
,
(b)
measurement residual:z
~
L
,
k
(
i
)
=
z
L
,
k
-
H
L
x
^
L
,
k
∣
k
-
1
(
i
)
-
ν
k,
(c)
residual covariance:S
L
,
k
i
=
H
L
P
L
,
k
∣
k
-
1
i
H
L
T
-
R
k,
(d)
filter gain:K
L
,
k
i
=
P
L
,
k
∣
k
-
1
i
H
L
T
S
L
,
k
i
-
1,
(e)
update of state and covariance:x
^
L
,
k
∣
k
i
=
x
^
L
,
k
∣
k
-
1
i
+
K
L
,
k
i
z
~
L
,
k
i
,
P
L
,
k
∣
k
i
=
P
L
,
k
∣
k
-
1
i
-
K
L
,
k
i
S
L
,
k
i
K
L
,
k
i
T.
(iii)
Mode probability update (fori
=
1,2
,
…
,
M):
(a)
mode likelihood:Λ
L
,
k
i
≜
p
z
L
,
k
i
∣
m
k
i
,
Z
L
,
k
-
1
=
a
s
s
u
m
e
N
z
L
,
k
i
;
0
,
S
L
,
k
i,
(b)
mode probability:μ
L
,
k
i
=
μ
L
,
k
∣
k
-
1
i
Λ
L
,
k
i
/
∑
j
μ
L
,
k
∣
k
-
1
j
Λ
L
,
k
j.
(iv)
Combination (fori
=
1,2
,
…
,
M):(28)
x
^
L
,
k
∣
k
=
∑
i
x
^
L
,
k
∣
k
i
μ
L
,
k
i
,
P
^
L
,
k
∣
k
i
=
∑
i
P
L
,
k
∣
k
i
+
x
^
L
,
k
∣
k
-
x
^
L
,
k
∣
k
i
x
^
L
,
k
∣
k
-
x
^
L
,
k
∣
k
i
T
μ
L
,
k
i
.Singer’s MMS of time delayed slaving data is completed by propagating the combined estimate based on the current mode’s dynamic model such that (29)
x
^
L
,
k
+
s
∣
k
s
y
n
c
=
Φ
k
+
s
,
k
s
y
n
c
x
^
L
,
k
∣
k
,
Φ
k
+
s
,
k
s
y
n
c
=
F
k
+
s
-
1
s
y
n
c
F
k
+
s
-
2
s
y
n
c
⋯
F
k
s
y
n
c
,
Φ
k
,
k
s
y
n
c
=
I
n
,
k
>
1
,where F
k
+
s
-
1
s
y
n
c
,
F
k
+
s
-
2
s
y
n
c
,
…
,
F
k
s
y
n
c are system matrices depending on a flight phase mode at current time k.


### 4.2. Synchronization Using an IMM with a Nonlinear Ballistic Model

For a nonlinear ballistic model, the state vector for the propelled mode is denoted as(30)
x
t
=
x
t
y
t
z
t
x
˙
t
y
˙
t
z
˙
t
ξ
t
τ
t
T
,where ξ
t is the drag coefficient and τ
t is the thrust. Generally, the drag coefficient varies significantly with the Mach number regime: subsonic, transonic, and supersonic. Therefore, we take advantage of the dynamic model considering a Mach number-dependent multiplier [11, 12] such that(31)
x
¨
y
¨
z
¨
=
τ
V
x
˙
y
˙
z
˙
thrust  term
+
ξ
ξ
m
D
x
˙
y
˙
z
˙
drag  term
+
g
0
0
1
g
r
a
v
i
t
y
term
+
ω
~
1
,
ξ
˙
=
ω
~
2
,
τ
˙
=
ω
~
3
.The first term on the right side of (31) represents the thrust (m
2
/
s) of the SLV in x, y, and z directions. Two distinct multiple modes of a nonlinear SLV model can be separated by the existence of thrust. In other words, the existence of thrust signifies propelled flight mode, whereas zero thrust signifies coasting flight mode. V is the magnitude of the velocity v
=
x
˙
y
˙
z
˙
T, that is, the SLV speed (m
/
s). The second part of (31) is the drag term, which is related to velocity and altitude; that is, D
=
-
ρ
z
V
/
2, where ρ
z
=
ρ
0
e
-
c
z is the air density (k
g
/
m
3) at an altitude z
(
m
) and c is the air density constant (m
-
1) [9]. ξ is the drag coefficient and ξ
m is the Mach number-dependent drag coefficient multiplier, which is approximated by the cubic spline curve shown in Figure 3. In this paper, for the drag characteristics of the SLV, which are applied in the subsequent simulation section, drag coefficients of the Saturn V launch vehicle [15] are used. The third part of (31) is a gravity term. Gravity g is the standard acceleration due to gravity at sea level, which is assumed to be constant throughout the trajectory, with a value of 9.812 m
/
s
2. ω
~
1, ω
~
2, and ω
~
3 are assumed to be continuous time zero mean white Gaussian noises. The drag coefficient and thrust acceleration are represented as Wiener processes with a slow variation [11, 12]. We can modify the dynamic equations (30) and (31) as a compact form such that (32)
x
˙
t
=
f
x
t
+
ω
~
t
,where (33)
f
x
t
=
x
˙
t
y
˙
t
z
˙
t
τ
x
˙
t
V
t
+
ξ
t
D
t
x
˙
t
τ
y
˙
t
V
t
+
ξ
t
D
t
y
˙
t
τ
z
˙
t
V
t
+
ξ
t
D
t
z
˙
t
-
g
,
ω
~
t
=
ω
~
1
,
t
ω
~
2
,
t
ω
~
3
,
t
T
.The state vector equation (33) is discretized by a second-order Taylor expansion [26]. Then, (33) can be written as a discretized continuous time system with white process noise such that(34)
x
k
+
1
=
x
k
+
f
x
k
Δ
t
+
A
k
f
x
k
Δ
t
2
2
+
ω
k
,where A
k is the Jacobian of (33) evaluated at x
k [26] and ω
k is the discretized continuous time process noise for the sampling interval Δ
t. The corresponding covariance matrix of the discretized process noise is(35)
Q
=
Q
1
q
v
0
6
×
1
0
6
×
1
0
1
×
6
Δ
t
q
ξ
0
0
1
×
6
0
Δ
t
q
τ
,
Q
1
=
Δ
t
3
3
I
3
Δ
t
2
2
I
3
Δ
t
2
2
I
3
Δ
t
I
3
,where I
3 is the 3
×
3 identity matrix and the continuous time process noise intensities q
v, q
ξ, and q
τ are the corresponding power spectral densities.Figure 3
Normalized drag coefficient [15].The measurement matrixH
N
L for the nonlinear ballistic multiple model can be depicted as(36)
z
N
L
,
k
=
H
N
L
x
N
L
,
k
+
ν
k
,
H
N
L
=
I
3
0
,
x
N
L
,
k
=
x
k
y
k
z
k
x
˙
k
y
˙
k
z
˙
k
ξ
k
τ
k
T
,where ν
k is measurement noise with error covariance R
k
δ
k
-
j
=
E
ν
k
ν
j
T, that is, (14).An IMM algorithm for nonlinear dynamics with different sizes of the mode state vector is summarized as follows [25]:(v)
Model-conditioned reinitialization (fori
=
1,2
,
…
,
M):
(a)
predicted mode probability:μ
N
L
,
k
∣
k
-
1
i
≜
Pr
⁡
m
k
i
∣
z
N
L
,
k
-
1
=
∑
j
π
j
i
μ
N
L
,
k
-
1
j,
(b)
mixing weight:(37)
μ
N
L
,
k
-
1
j
∣
i
≜
P
r
m
k
-
1
j
∣
m
k
i
,
z
N
L
,
k
-
1
=
∑
j
π
j
i
μ
N
L
,
k
-
1
j
μ
N
L
,
k
-
1
i
,
(c)
unbiased mixing estimate and covariance:(38)
x
-
N
L
,
k
-
1
∣
k
-
1
i
≜
E
x
N
L
,
k
-
1
∣
k
-
1
∣
m
k
i
,
z
N
L
,
k
-
1
=
∑
j
χ
^
N
L
,
k
-
1
∣
k
-
1
j
μ
N
L
,
k
-
1
j
∣
i
,
where
χ
^
N
L
,
k
-
1
∣
k
-
1
1
≜
x
^
N
L
,
k
-
1
∣
k
-
1
1
,
χ
^
N
L
,
k
-
1
∣
k
-
1
2
≜
x
^
N
L
,
k
-
1
∣
k
-
1
2
T
,
τ
k
-
1
T
,
P
-
N
L
,
k
-
1
∣
k
-
1
i
=
∑
j
P
N
L
,
k
-
1
∣
k
-
1
i
+
x
-
N
L
,
k
-
1
∣
k
-
1
i
-
χ
^
N
L
,
k
-
1
∣
k
-
1
j
x
-
N
L
,
k
-
1
∣
k
-
1
i
-
χ
^
N
L
,
k
-
1
∣
k
-
1
j
T
μ
N
L
,
k
-
1
j
∣
i
.
(vi)
Model-conditioned filtering (fori
=
1,2
,
…
,
M):
(a)
predicted estimate and covariance:(39)
x
^
N
L
,
k
∣
k
-
1
i
=
x
^
N
L
,
k
-
1
∣
k
-
1
i
+
f
x
^
N
L
,
k
-
1
∣
k
-
1
i
Δ
t
+
A
^
N
L
,
k
-
1
i
f
x
^
N
L
,
k
-
1
∣
k
-
1
i
Δ
t
2
2
,
A
^
k
i
=
∂
f
∂
x
x
=
x
^
k
i
,
P
N
L
,
k
∣
k
-
1
i
=
Ω
k
-
1
i
P
-
N
L
,
k
-
1
∣
k
-
1
i
Ω
k
-
1
i
T
+
Q
k
-
1
,
Ω
k
-
1
i
=
I
+
A
^
N
L
,
k
-
1
i
,
(b)
measurement residual:z
~
N
L
,
k
(
i
)
=
z
N
L
,
k
-
H
N
L
x
^
N
L
,
k
∣
k
-
1
(
i
)
-
ν
k,
(c)
residual covariance:S
N
L
,
k
i
=
H
N
L
P
N
L
,
k
∣
k
-
1
i
H
N
L
T
-
R
k,
(d)
filter gain:K
N
L
,
k
i
=
P
N
L
,
k
∣
k
-
1
i
H
N
L
T
S
N
L
,
k
i
-
1,
(e)
update of state and covariance:(40)
x
^
N
L
,
k
∣
k
i
=
x
^
N
L
,
k
∣
k
-
1
i
+
K
N
L
,
k
i
z
~
N
L
,
k
i
,
P
N
L
,
k
∣
k
i
=
P
N
L
,
k
∣
k
-
1
i
-
K
N
L
,
k
i
S
N
L
,
k
i
K
N
L
,
k
i
T
.
(vii)
Mode probability update (fori
=
1,2
,
…
,
M):
(a)
mode likelihood:(41)
Λ
N
L
,
k
i
≜
p
z
N
L
,
k
i
∣
m
k
i
,
Z
N
L
,
k
-
1
=
a
s
s
u
m
e
N
z
N
L
,
k
i
;
0
,
S
N
L
,
k
i
,
(b)
mode probability:μ
N
L
,
k
i
=
μ
N
L
,
k
∣
k
-
1
i
Λ
N
L
,
k
i
/
∑
j
μ
N
L
,
k
∣
k
-
1
j
Λ
N
L
,
k
j.
(viii)
Combination (fori
=
1,2
,
…
,
M): (42)
x
^
N
L
,
k
∣
k
=
∑
i
x
^
N
L
,
k
∣
k
i
μ
N
L
,
k
i
,
P
^
N
L
,
k
∣
k
i
=
∑
i
P
N
L
,
k
∣
k
i
+
x
^
N
L
,
k
∣
k
-
x
^
N
L
,
k
∣
k
i
x
^
N
L
,
k
∣
k
-
x
^
N
L
,
k
∣
k
i
T
μ
N
L
,
k
i
.The nonlinear MMS of time delayed slaving data is completed by propagating the combined estimate based on the current mode’s estimated vector such that (43)
x
^
N
L
,
k
+
s
∣
k
=
x
^
N
L
,
k
∣
k
+
f
x
^
N
L
,
k
∣
k
·
s
+
A
^
N
L
,
k
f
x
^
N
L
,
k
∣
k
·
s
2
2
,where s is a lead-time for synchronization (s
=
n
·
Δ
t). If the SLV is in the PFP (μ
N
L
,
k
1
>
μ
N
L
,
k
2), the current state estimate includes the thrust term; that is, the state vector x
N
L
,
k
=
x
k
y
k
z
k
x
˙
k
y
˙
k
z
˙
k
ξ
k
τ
k
T. If the SLV is not in the PFP, the state vector does not include the thrust term (τ
k
=
0).


## 4.1. Synchronization Using an IMM with Singer’s Linear Model

Singer [16] described 2D manned maneuvering targets in range-bearing coordinates. This model can be adapted to the SLV kinematics in 3D Cartesian coordinates with multiple flight phases: (17)
p
˙
n
v
˙
n
a
˙
n
=
F
p
n
v
n
a
n
+
G
w
n
,
F
=
0
1
0
0
0
1
0
0
-
α
n
,
G
=
0
0
0
,where w
n
∈
R
3
×
1 is a white noise process along the Cartesian axis n
=
x
,
y
,
z
. The parameter β
=
1
/
α
n is the maneuver correlation time constant, and σ
a
n
2 is the acceleration variance describing maneuver intensity. In a steady state,(18)
σ
w
n
2
=
2
α
σ
a
n
2
.To describe SLV kinematics, the model must cope with crucial nonzero mean acceleration maneuvers during the propelled phase. In addition, after each stage’s engine burns out, a multiple model approach is applied to the coasting flight; that is, one model describes PFP, whereas the other depicts the CFP. Empirically tuned, independent probability density functions (PDFs) represented by TUM describe the accelerations of the SLV in the local coordinate frame. Figure 2 shows the means and variances of the acceleration processes of the SLV in this paper. The discrete-time model with state transition matrix Ψ
n
α
n
,
Δ
t is as follows:(19)
p
n
,
k
+
1
v
n
,
k
+
1
a
n
,
k
+
1
=
Ψ
n
α
n
,
Δ
t
p
n
,
k
v
n
,
k
a
n
,
k
+
w
n
,
k
,
Ψ
n
α
n
,
Δ
t
=
e
F
Δ
t
1
Δ
t
α
n
Δ
t
-
1
+
e
-
α
n
Δ
t
0
1
1
-
e
-
α
n
Δ
t
α
n
0
0
e
-
α
n
Δ
t
,
w
n
,
k
=
∫
k
Δ
t
k
+
1
Δ
t
e
F
k
+
1
Δ
t
-
r
G
w
n
d
r
.In the case of the PFP, w
n
,
k is a nonzero mean white noise sequence caused by the nonzero mean acceleration μ
a
n seen in Figure 2. A nonzero mean white noise sequence for the PFP should be considered in the target kinematics when implementing the state propagation stage in the KF. Thus, the deterministic input u
n
,
k caused by w
n
,
k along x-axis, y-axis, and z-axis is derived as follows: (20)
u
x
,
k
=
E
w
x
,
k
=
E
∫
k
Δ
t
k
+
1
Δ
t
e
F
k
+
1
Δ
t
-
r
G
w
x
d
r
≈
E
∫
k
Δ
t
k
+
1
Δ
t
e
-
α
x
k
+
1
Δ
t
-
r
0
0
w
x
d
r
,
u
x
,
k
=
μ
a
x
α
n
1
-
e
-
α
x
Δ
t
0
0
.
y-axis and z-axis can be derived in the same manner as shown in (20). The maneuver excitation covariance [10], which represents the uncertainty of the SLV kinematics model, is (21)
Q
n
,
k
α
n
,
Δ
t
=
E
w
n
,
k
-
E
w
n
,
k
w
n
,
k
-
E
w
n
,
k
T
=
2
α
n
σ
a
n
2
q
11
q
12
q
13
q
12
q
22
q
23
q
13
q
23
q
33
.The specific components of q
11
,
…
,
q
33 are illustrated in [25]. In addition, the measurement matrix H
L for Singer’s model can be depicted as (22)
z
L
,
k
=
H
L
x
L
,
k
+
ν
k
,
H
L
=
1
0
0
0
0
3
×
2
1
0
3
×
2
0
0
3
×
2
0
0
1
,
x
L
,
k
=
x
k
x
˙
k
x
¨
k
y
k
y
˙
k
y
¨
k
z
k
z
˙
k
z
¨
k
T
,where ν
k is measurement noise (as shown in (14)) with error covariance.Figure 2
Nonzero mean acceleration PDF in the PFP model (a) alongx-axis, (b) along y-axis, and (c) along z-axis. Zero mean acceleration PDF in the CFP Model (d) along x-axis, y-axis, and z-axis.
(a)
(b)
(c)
(d)


### 4.1.1. IMM with Singer’s Linear Model

From (19)–(22) in Section 4.1, we can rewrite the Markov jump linear systems, where the ith model of the finite multiple model set M
=
{
m
1
,
…
,
m
M
} obeys the following equations:(23)
x
L
,
k
+
1
=
Ψ
k
i
x
L
,
k
+
w
k
i
,
z
L
,
k
=
H
L
x
L
,
k
+
ν
k
,where(24)
c
o
v
w
k
i
=
Q
k
i
,
c
o
v
ν
k
=
R
k
,
Ψ
k
i
=
d
i
a
g
Ψ
x
,
k
α
x
i
,
Δ
t
,
Ψ
y
,
k
α
y
i
,
Δ
t
,
Ψ
z
,
k
α
z
i
,
Δ
t
,
Q
k
i
=
d
i
a
g
Q
x
,
k
α
x
i
,
Δ
t
,
Q
y
,
k
α
y
i
,
Δ
t
,
Q
z
,
k
α
z
i
,
Δ
t
.Here, σ
a
n
2 for the SLV can be represented by TUM as in Figure 2. The PDFs of Figures 2(a)–2(d) are experimentally sampled from the nominal acceleration profile of the SLV. The superscript (i) denotes quantities pertinent to model m
(
i
) in M, and the jumps of the system mode are assumed to have transition probabilities:(25)
Pr
⁡
m
k
+
1
j
∣
m
k
i
≜
π
i
j
,where m
k
i denotes the event in which model m
i matches the system mode in effect at time k. In our application, M
=
2, and m
1 and m
2 are the propelled and coasting modes, respectively.Finally, complete recursion of the IMM with mode matched KF for the SLV tracking is summarized as follows:(i)
Model-conditioned reinitialization (fori
=
1,2
,
…
,
M):
(a)
predicted mode probability:μ
L
,
k
∣
k
-
1
i
≜
Pr
⁡
{
m
k
i
∣
z
L
,
k
-
1
}
=
∑
j
π
j
i
μ
L
,
k
-
1
j,
(b)
mixing weight:μ
L
,
k
-
1
j
∣
i
≜
Pr
⁡
{
m
k
-
1
j
∣
m
k
i
,
z
L
,
k
-
1
}
=
∑
j
π
j
i
μ
L
,
k
-
1
j
/
μ
L
,
k
-
1
i,
(c)
mixing estimate and covariance:(26)
x
-
L
,
k
-
1
∣
k
-
1
i
≜
E
x
L
,
k
-
1
∣
k
-
1
∣
m
k
i
,
z
L
,
k
-
1
=
∑
j
x
^
L
,
k
-
1
∣
k
-
1
j
μ
L
,
k
-
1
j
∣
i
,
P
-
L
,
k
-
1
∣
k
-
1
i
=
∑
j
P
L
,
k
-
1
∣
k
-
1
i
+
x
-
L
,
k
-
1
∣
k
-
1
i
-
x
^
L
,
k
-
1
∣
k
-
1
j
x
-
L
,
k
-
1
∣
k
-
1
i
-
x
^
L
,
k
-
1
∣
k
-
1
j
T
μ
L
,
k
-
1
j
∣
i
.
(ii)
Model-conditioned filtering (fori
=
1,2
,
…
,
M):
(a)
predicted estimate and covariance:(27)
x
^
L
,
k
∣
k
-
1
i
=
Ψ
k
-
1
i
x
-
L
,
k
-
1
∣
k
-
1
i
+
u
k
-
1
i
,
where
u
k
-
1
1
=
u
x
,
k
-
1
T
u
y
,
k
-
1
T
u
z
,
k
-
1
T
T
,
u
k
-
1
2
=
0
3
×
1
,
P
L
,
k
∣
k
-
1
i
=
Ψ
k
-
1
i
P
-
L
,
k
-
1
∣
k
-
1
i
Ψ
k
-
1
i
T
+
Q
k
-
1
i
,
(b)
measurement residual:z
~
L
,
k
(
i
)
=
z
L
,
k
-
H
L
x
^
L
,
k
∣
k
-
1
(
i
)
-
ν
k,
(c)
residual covariance:S
L
,
k
i
=
H
L
P
L
,
k
∣
k
-
1
i
H
L
T
-
R
k,
(d)
filter gain:K
L
,
k
i
=
P
L
,
k
∣
k
-
1
i
H
L
T
S
L
,
k
i
-
1,
(e)
update of state and covariance:x
^
L
,
k
∣
k
i
=
x
^
L
,
k
∣
k
-
1
i
+
K
L
,
k
i
z
~
L
,
k
i
,
P
L
,
k
∣
k
i
=
P
L
,
k
∣
k
-
1
i
-
K
L
,
k
i
S
L
,
k
i
K
L
,
k
i
T.
(iii)
Mode probability update (fori
=
1,2
,
…
,
M):
(a)
mode likelihood:Λ
L
,
k
i
≜
p
z
L
,
k
i
∣
m
k
i
,
Z
L
,
k
-
1
=
a
s
s
u
m
e
N
z
L
,
k
i
;
0
,
S
L
,
k
i,
(b)
mode probability:μ
L
,
k
i
=
μ
L
,
k
∣
k
-
1
i
Λ
L
,
k
i
/
∑
j
μ
L
,
k
∣
k
-
1
j
Λ
L
,
k
j.
(iv)
Combination (fori
=
1,2
,
…
,
M):(28)
x
^
L
,
k
∣
k
=
∑
i
x
^
L
,
k
∣
k
i
μ
L
,
k
i
,
P
^
L
,
k
∣
k
i
=
∑
i
P
L
,
k
∣
k
i
+
x
^
L
,
k
∣
k
-
x
^
L
,
k
∣
k
i
x
^
L
,
k
∣
k
-
x
^
L
,
k
∣
k
i
T
μ
L
,
k
i
.Singer’s MMS of time delayed slaving data is completed by propagating the combined estimate based on the current mode’s dynamic model such that (29)
x
^
L
,
k
+
s
∣
k
s
y
n
c
=
Φ
k
+
s
,
k
s
y
n
c
x
^
L
,
k
∣
k
,
Φ
k
+
s
,
k
s
y
n
c
=
F
k
+
s
-
1
s
y
n
c
F
k
+
s
-
2
s
y
n
c
⋯
F
k
s
y
n
c
,
Φ
k
,
k
s
y
n
c
=
I
n
,
k
>
1
,where F
k
+
s
-
1
s
y
n
c
,
F
k
+
s
-
2
s
y
n
c
,
…
,
F
k
s
y
n
c are system matrices depending on a flight phase mode at current time k.


## 4.1.1. IMM with Singer’s Linear Model

From (19)–(22) in Section 4.1, we can rewrite the Markov jump linear systems, where the ith model of the finite multiple model set M
=
{
m
1
,
…
,
m
M
} obeys the following equations:(23)
x
L
,
k
+
1
=
Ψ
k
i
x
L
,
k
+
w
k
i
,
z
L
,
k
=
H
L
x
L
,
k
+
ν
k
,where(24)
c
o
v
w
k
i
=
Q
k
i
,
c
o
v
ν
k
=
R
k
,
Ψ
k
i
=
d
i
a
g
Ψ
x
,
k
α
x
i
,
Δ
t
,
Ψ
y
,
k
α
y
i
,
Δ
t
,
Ψ
z
,
k
α
z
i
,
Δ
t
,
Q
k
i
=
d
i
a
g
Q
x
,
k
α
x
i
,
Δ
t
,
Q
y
,
k
α
y
i
,
Δ
t
,
Q
z
,
k
α
z
i
,
Δ
t
.Here, σ
a
n
2 for the SLV can be represented by TUM as in Figure 2. The PDFs of Figures 2(a)–2(d) are experimentally sampled from the nominal acceleration profile of the SLV. The superscript (i) denotes quantities pertinent to model m
(
i
) in M, and the jumps of the system mode are assumed to have transition probabilities:(25)
Pr
⁡
m
k
+
1
j
∣
m
k
i
≜
π
i
j
,where m
k
i denotes the event in which model m
i matches the system mode in effect at time k. In our application, M
=
2, and m
1 and m
2 are the propelled and coasting modes, respectively.Finally, complete recursion of the IMM with mode matched KF for the SLV tracking is summarized as follows:(i)
Model-conditioned reinitialization (fori
=
1,2
,
…
,
M):
(a)
predicted mode probability:μ
L
,
k
∣
k
-
1
i
≜
Pr
⁡
{
m
k
i
∣
z
L
,
k
-
1
}
=
∑
j
π
j
i
μ
L
,
k
-
1
j,
(b)
mixing weight:μ
L
,
k
-
1
j
∣
i
≜
Pr
⁡
{
m
k
-
1
j
∣
m
k
i
,
z
L
,
k
-
1
}
=
∑
j
π
j
i
μ
L
,
k
-
1
j
/
μ
L
,
k
-
1
i,
(c)
mixing estimate and covariance:(26)
x
-
L
,
k
-
1
∣
k
-
1
i
≜
E
x
L
,
k
-
1
∣
k
-
1
∣
m
k
i
,
z
L
,
k
-
1
=
∑
j
x
^
L
,
k
-
1
∣
k
-
1
j
μ
L
,
k
-
1
j
∣
i
,
P
-
L
,
k
-
1
∣
k
-
1
i
=
∑
j
P
L
,
k
-
1
∣
k
-
1
i
+
x
-
L
,
k
-
1
∣
k
-
1
i
-
x
^
L
,
k
-
1
∣
k
-
1
j
x
-
L
,
k
-
1
∣
k
-
1
i
-
x
^
L
,
k
-
1
∣
k
-
1
j
T
μ
L
,
k
-
1
j
∣
i
.
(ii)
Model-conditioned filtering (fori
=
1,2
,
…
,
M):
(a)
predicted estimate and covariance:(27)
x
^
L
,
k
∣
k
-
1
i
=
Ψ
k
-
1
i
x
-
L
,
k
-
1
∣
k
-
1
i
+
u
k
-
1
i
,
where
u
k
-
1
1
=
u
x
,
k
-
1
T
u
y
,
k
-
1
T
u
z
,
k
-
1
T
T
,
u
k
-
1
2
=
0
3
×
1
,
P
L
,
k
∣
k
-
1
i
=
Ψ
k
-
1
i
P
-
L
,
k
-
1
∣
k
-
1
i
Ψ
k
-
1
i
T
+
Q
k
-
1
i
,
(b)
measurement residual:z
~
L
,
k
(
i
)
=
z
L
,
k
-
H
L
x
^
L
,
k
∣
k
-
1
(
i
)
-
ν
k,
(c)
residual covariance:S
L
,
k
i
=
H
L
P
L
,
k
∣
k
-
1
i
H
L
T
-
R
k,
(d)
filter gain:K
L
,
k
i
=
P
L
,
k
∣
k
-
1
i
H
L
T
S
L
,
k
i
-
1,
(e)
update of state and covariance:x
^
L
,
k
∣
k
i
=
x
^
L
,
k
∣
k
-
1
i
+
K
L
,
k
i
z
~
L
,
k
i
,
P
L
,
k
∣
k
i
=
P
L
,
k
∣
k
-
1
i
-
K
L
,
k
i
S
L
,
k
i
K
L
,
k
i
T.
(iii)
Mode probability update (fori
=
1,2
,
…
,
M):
(a)
mode likelihood:Λ
L
,
k
i
≜
p
z
L
,
k
i
∣
m
k
i
,
Z
L
,
k
-
1
=
a
s
s
u
m
e
N
z
L
,
k
i
;
0
,
S
L
,
k
i,
(b)
mode probability:μ
L
,
k
i
=
μ
L
,
k
∣
k
-
1
i
Λ
L
,
k
i
/
∑
j
μ
L
,
k
∣
k
-
1
j
Λ
L
,
k
j.
(iv)
Combination (fori
=
1,2
,
…
,
M):(28)
x
^
L
,
k
∣
k
=
∑
i
x
^
L
,
k
∣
k
i
μ
L
,
k
i
,
P
^
L
,
k
∣
k
i
=
∑
i
P
L
,
k
∣
k
i
+
x
^
L
,
k
∣
k
-
x
^
L
,
k
∣
k
i
x
^
L
,
k
∣
k
-
x
^
L
,
k
∣
k
i
T
μ
L
,
k
i
.Singer’s MMS of time delayed slaving data is completed by propagating the combined estimate based on the current mode’s dynamic model such that (29)
x
^
L
,
k
+
s
∣
k
s
y
n
c
=
Φ
k
+
s
,
k
s
y
n
c
x
^
L
,
k
∣
k
,
Φ
k
+
s
,
k
s
y
n
c
=
F
k
+
s
-
1
s
y
n
c
F
k
+
s
-
2
s
y
n
c
⋯
F
k
s
y
n
c
,
Φ
k
,
k
s
y
n
c
=
I
n
,
k
>
1
,where F
k
+
s
-
1
s
y
n
c
,
F
k
+
s
-
2
s
y
n
c
,
…
,
F
k
s
y
n
c are system matrices depending on a flight phase mode at current time k.


## 4.2. Synchronization Using an IMM with a Nonlinear Ballistic Model

For a nonlinear ballistic model, the state vector for the propelled mode is denoted as(30)
x
t
=
x
t
y
t
z
t
x
˙
t
y
˙
t
z
˙
t
ξ
t
τ
t
T
,where ξ
t is the drag coefficient and τ
t is the thrust. Generally, the drag coefficient varies significantly with the Mach number regime: subsonic, transonic, and supersonic. Therefore, we take advantage of the dynamic model considering a Mach number-dependent multiplier [11, 12] such that(31)
x
¨
y
¨
z
¨
=
τ
V
x
˙
y
˙
z
˙
thrust  term
+
ξ
ξ
m
D
x
˙
y
˙
z
˙
drag  term
+
g
0
0
1
g
r
a
v
i
t
y
term
+
ω
~
1
,
ξ
˙
=
ω
~
2
,
τ
˙
=
ω
~
3
.The first term on the right side of (31) represents the thrust (m
2
/
s) of the SLV in x, y, and z directions. Two distinct multiple modes of a nonlinear SLV model can be separated by the existence of thrust. In other words, the existence of thrust signifies propelled flight mode, whereas zero thrust signifies coasting flight mode. V is the magnitude of the velocity v
=
x
˙
y
˙
z
˙
T, that is, the SLV speed (m
/
s). The second part of (31) is the drag term, which is related to velocity and altitude; that is, D
=
-
ρ
z
V
/
2, where ρ
z
=
ρ
0
e
-
c
z is the air density (k
g
/
m
3) at an altitude z
(
m
) and c is the air density constant (m
-
1) [9]. ξ is the drag coefficient and ξ
m is the Mach number-dependent drag coefficient multiplier, which is approximated by the cubic spline curve shown in Figure 3. In this paper, for the drag characteristics of the SLV, which are applied in the subsequent simulation section, drag coefficients of the Saturn V launch vehicle [15] are used. The third part of (31) is a gravity term. Gravity g is the standard acceleration due to gravity at sea level, which is assumed to be constant throughout the trajectory, with a value of 9.812 m
/
s
2. ω
~
1, ω
~
2, and ω
~
3 are assumed to be continuous time zero mean white Gaussian noises. The drag coefficient and thrust acceleration are represented as Wiener processes with a slow variation [11, 12]. We can modify the dynamic equations (30) and (31) as a compact form such that (32)
x
˙
t
=
f
x
t
+
ω
~
t
,where (33)
f
x
t
=
x
˙
t
y
˙
t
z
˙
t
τ
x
˙
t
V
t
+
ξ
t
D
t
x
˙
t
τ
y
˙
t
V
t
+
ξ
t
D
t
y
˙
t
τ
z
˙
t
V
t
+
ξ
t
D
t
z
˙
t
-
g
,
ω
~
t
=
ω
~
1
,
t
ω
~
2
,
t
ω
~
3
,
t
T
.The state vector equation (33) is discretized by a second-order Taylor expansion [26]. Then, (33) can be written as a discretized continuous time system with white process noise such that(34)
x
k
+
1
=
x
k
+
f
x
k
Δ
t
+
A
k
f
x
k
Δ
t
2
2
+
ω
k
,where A
k is the Jacobian of (33) evaluated at x
k [26] and ω
k is the discretized continuous time process noise for the sampling interval Δ
t. The corresponding covariance matrix of the discretized process noise is(35)
Q
=
Q
1
q
v
0
6
×
1
0
6
×
1
0
1
×
6
Δ
t
q
ξ
0
0
1
×
6
0
Δ
t
q
τ
,
Q
1
=
Δ
t
3
3
I
3
Δ
t
2
2
I
3
Δ
t
2
2
I
3
Δ
t
I
3
,where I
3 is the 3
×
3 identity matrix and the continuous time process noise intensities q
v, q
ξ, and q
τ are the corresponding power spectral densities.Figure 3
Normalized drag coefficient [15].The measurement matrixH
N
L for the nonlinear ballistic multiple model can be depicted as(36)
z
N
L
,
k
=
H
N
L
x
N
L
,
k
+
ν
k
,
H
N
L
=
I
3
0
,
x
N
L
,
k
=
x
k
y
k
z
k
x
˙
k
y
˙
k
z
˙
k
ξ
k
τ
k
T
,where ν
k is measurement noise with error covariance R
k
δ
k
-
j
=
E
ν
k
ν
j
T, that is, (14).An IMM algorithm for nonlinear dynamics with different sizes of the mode state vector is summarized as follows [25]:(v)
Model-conditioned reinitialization (fori
=
1,2
,
…
,
M):
(a)
predicted mode probability:μ
N
L
,
k
∣
k
-
1
i
≜
Pr
⁡
m
k
i
∣
z
N
L
,
k
-
1
=
∑
j
π
j
i
μ
N
L
,
k
-
1
j,
(b)
mixing weight:(37)
μ
N
L
,
k
-
1
j
∣
i
≜
P
r
m
k
-
1
j
∣
m
k
i
,
z
N
L
,
k
-
1
=
∑
j
π
j
i
μ
N
L
,
k
-
1
j
μ
N
L
,
k
-
1
i
,
(c)
unbiased mixing estimate and covariance:(38)
x
-
N
L
,
k
-
1
∣
k
-
1
i
≜
E
x
N
L
,
k
-
1
∣
k
-
1
∣
m
k
i
,
z
N
L
,
k
-
1
=
∑
j
χ
^
N
L
,
k
-
1
∣
k
-
1
j
μ
N
L
,
k
-
1
j
∣
i
,
where
χ
^
N
L
,
k
-
1
∣
k
-
1
1
≜
x
^
N
L
,
k
-
1
∣
k
-
1
1
,
χ
^
N
L
,
k
-
1
∣
k
-
1
2
≜
x
^
N
L
,
k
-
1
∣
k
-
1
2
T
,
τ
k
-
1
T
,
P
-
N
L
,
k
-
1
∣
k
-
1
i
=
∑
j
P
N
L
,
k
-
1
∣
k
-
1
i
+
x
-
N
L
,
k
-
1
∣
k
-
1
i
-
χ
^
N
L
,
k
-
1
∣
k
-
1
j
x
-
N
L
,
k
-
1
∣
k
-
1
i
-
χ
^
N
L
,
k
-
1
∣
k
-
1
j
T
μ
N
L
,
k
-
1
j
∣
i
.
(vi)
Model-conditioned filtering (fori
=
1,2
,
…
,
M):
(a)
predicted estimate and covariance:(39)
x
^
N
L
,
k
∣
k
-
1
i
=
x
^
N
L
,
k
-
1
∣
k
-
1
i
+
f
x
^
N
L
,
k
-
1
∣
k
-
1
i
Δ
t
+
A
^
N
L
,
k
-
1
i
f
x
^
N
L
,
k
-
1
∣
k
-
1
i
Δ
t
2
2
,
A
^
k
i
=
∂
f
∂
x
x
=
x
^
k
i
,
P
N
L
,
k
∣
k
-
1
i
=
Ω
k
-
1
i
P
-
N
L
,
k
-
1
∣
k
-
1
i
Ω
k
-
1
i
T
+
Q
k
-
1
,
Ω
k
-
1
i
=
I
+
A
^
N
L
,
k
-
1
i
,
(b)
measurement residual:z
~
N
L
,
k
(
i
)
=
z
N
L
,
k
-
H
N
L
x
^
N
L
,
k
∣
k
-
1
(
i
)
-
ν
k,
(c)
residual covariance:S
N
L
,
k
i
=
H
N
L
P
N
L
,
k
∣
k
-
1
i
H
N
L
T
-
R
k,
(d)
filter gain:K
N
L
,
k
i
=
P
N
L
,
k
∣
k
-
1
i
H
N
L
T
S
N
L
,
k
i
-
1,
(e)
update of state and covariance:(40)
x
^
N
L
,
k
∣
k
i
=
x
^
N
L
,
k
∣
k
-
1
i
+
K
N
L
,
k
i
z
~
N
L
,
k
i
,
P
N
L
,
k
∣
k
i
=
P
N
L
,
k
∣
k
-
1
i
-
K
N
L
,
k
i
S
N
L
,
k
i
K
N
L
,
k
i
T
.
(vii)
Mode probability update (fori
=
1,2
,
…
,
M):
(a)
mode likelihood:(41)
Λ
N
L
,
k
i
≜
p
z
N
L
,
k
i
∣
m
k
i
,
Z
N
L
,
k
-
1
=
a
s
s
u
m
e
N
z
N
L
,
k
i
;
0
,
S
N
L
,
k
i
,
(b)
mode probability:μ
N
L
,
k
i
=
μ
N
L
,
k
∣
k
-
1
i
Λ
N
L
,
k
i
/
∑
j
μ
N
L
,
k
∣
k
-
1
j
Λ
N
L
,
k
j.
(viii)
Combination (fori
=
1,2
,
…
,
M): (42)
x
^
N
L
,
k
∣
k
=
∑
i
x
^
N
L
,
k
∣
k
i
μ
N
L
,
k
i
,
P
^
N
L
,
k
∣
k
i
=
∑
i
P
N
L
,
k
∣
k
i
+
x
^
N
L
,
k
∣
k
-
x
^
N
L
,
k
∣
k
i
x
^
N
L
,
k
∣
k
-
x
^
N
L
,
k
∣
k
i
T
μ
N
L
,
k
i
.The nonlinear MMS of time delayed slaving data is completed by propagating the combined estimate based on the current mode’s estimated vector such that (43)
x
^
N
L
,
k
+
s
∣
k
=
x
^
N
L
,
k
∣
k
+
f
x
^
N
L
,
k
∣
k
·
s
+
A
^
N
L
,
k
f
x
^
N
L
,
k
∣
k
·
s
2
2
,where s is a lead-time for synchronization (s
=
n
·
Δ
t). If the SLV is in the PFP (μ
N
L
,
k
1
>
μ
N
L
,
k
2), the current state estimate includes the thrust term; that is, the state vector x
N
L
,
k
=
x
k
y
k
z
k
x
˙
k
y
˙
k
z
˙
k
ξ
k
τ
k
T. If the SLV is not in the PFP, the state vector does not include the thrust term (τ
k
=
0).


## 5. Simulation Results

To demonstrate the performance of the proposed synchronization approaches for delayed slaving data, we simulated the SLV tracking problem based on the nominal flight trajectory of the Korea Space Launch Vehicle-I (KSLV-I). In the simulation, the radar measurement noise intensities of (12) are selected as σ
r
=
15 m, σ
φ
=
0.01 deg, and σ
θ
=
0.01 deg. The nominal flight sequence of KSLV-I is as follows. First, the payload fairing separates during the first stage flight at 215.4 s. After the first stage engine shutdown at 228.7 s, the upper stage separates from the first stage and enters the CFP. The second stage continues in the CFP until the kick motor ignition at 395 s, and the vehicle then enters the PFP. At the end of kick motor combustion (452.7 s), the upper stage enters the target orbit in CFP as in the previous separation. Finally, the satellite is inserted into the target orbit after it separates from the upper stage during CFP at 540 s. Since delayed slaving data from the MCS to the TSs can occur in both the PFP and CFP, synchronization simulations are conducted at arbitrary points of the flight phases. Figure 4 shows the synchronization error of the delayed slaving data with respect to the benchmark and the proposed approaches for delays of s
=
0.1,0.2
,
…
,
1 s. Figures 4(a)–4(c) illustrate synchronization errors, where delayed slaving data occurred at 155 s, which is the first stage of PFP, and Figures 4(d)–4(f) show synchronization errors, where delayed slaving data occurred at 280 s, which is the first stage of CFP. As shown in Figures 4(a)–4(c), the extrapolation and α-β-γ filter-based approach cannot find the proper synchronized slaving data. The synchronization errors of these approaches exponentially increase according to increasing delay s; in particular, the error of the extrapolation approaches an out-of-figure bound (2500 m). On the other hand, the Kalman prediction and the proposed approaches have stable synchronization errors even if the delay time s is larger. Furthermore, we can observe that the errors in the synchronization approach of the Kalman prediction are relatively larger than the proposed approaches, especially in the PFP. As explained in the flight sequence of KSLV-I, in general, the motion of the SLV is described by multiple dynamics (in our case PFP and CFP) rather than single dynamics. However, the dynamic model of the KF and prediction capability used in this paper is the constant acceleration (CA) model. This means that the filter is optimized with respect to CA motion, which is actually not occurring in our application. This is why the performance of the Kalman prediction-based synchronization is worse than the proposed approaches. Nonetheless, we can see that the y-axis synchronization errors shown in Figures 4(b) and 4(e) are small. This result is observed because the y-axis motion of KSLV-I in a local coordinate is relatively smaller than the other axis motion, and the CA model approximately expresses this small motion. A single KF prediction model cannot exactly describe the motion of both PFP and CFP, whereas the proposed multiple model-based approaches work well. As shown in Figure 4, regardless of the delay s
, the proposed multiple model-based approaches find synchronized position vectors. In the comparison of synchronization performance between linear IMM and nonlinear IMM synchronization, the nonlinear IMM-based synchronization approach shows the best performance, except for y-axis, where CA motion is dominant. In addition, the difference between the proposed multiple model-based approaches is very small as shown in Figure 4, but the complexities of the algorithms for real-time applications are quite dissimilar. Hence, the operator may adaptively select one of the proposed approaches according to one’s environment.Figure 4
Synchronization errors at PFP and CFP.


## 6. Conclusions

In this paper, we investigated the time synchronization approaches of delayed slaving data in the RSS for SLV tracking. One of the most important roles of the RSS is to distribute slaving data to each TS for continuous tracking of the SLV. If there is a critical network delay resulting in time delayed slaving data being sent to each TS, the MCS will not receive accurate SLV tracking data. This problem can give rise to significant difficulties for the SLV mission progress and analysis. To overcome this problem, we proposed MMS approaches which take advantage of the multiple motion models of an SLV. The linear IMM-based synchronization approach was developed using Singer’s model with ternary uniform mixtures and the nonlinear IMM-based synchronization approach was derived from a nonlinear ballistic model with a drag coefficient. For verification of the proposed algorithms, SLV tracking simulations using KSLV-I and the radar measurement data generated from nominal trajectory were conducted. To demonstrate the superiority of time synchronization performance in these simulations, we compared the proposed algorithm with benchmark approaches for absolute error between the nominal trajectory data and the synchronized slaving data; the simulation results demonstrated that the proposed MMS approaches performed competitively.



---

*Source: 9319282-2016-08-23.xml*