Distributed multiple model extended information filter with unbiased mixing for satellite launch vehicle tracking

Abstract

A distributed extended information filter-based interacting multiple model estimator with unbiased mixing is proposed for satellite launch vehicle tracking. In this problem, multiple heterogeneous sensors such as radars, telemetry systems receiving onboard Global Positioning System—inertial navigation system data, and electro-optical targeting systems are used. The extended information filter is used for nonlinear estimation dealing with ballistic model and spherical coordinate observation. The multiple Markov switching models comprise thrusting and coasting modes having different state vector dimensions for the launch vehicle. To effectively combine both state vectors, an unbiased mixing technique is applied and then the distributed extended information filter integrates local states and information matrix contributions. Hence, the proposed algorithm takes into account both heterogeneity of tracking sensors and multiplicity of vehicle’s dynamic model. We prove the superiority of the proposed algorithm by conducting Monte Carlo simulation with nominal trajectory data of Korea Space Launch Vehicle-1. Comparative simulation results demonstrate that the performance of the proposed method has been improved in vehicle’s position root mean square error.

Introduction

Generally, a satellite launch vehicle (SLV) flies over a very large distance for successful injection of a satellite into the designed orbit. During the flight, most SLVs experience different flight phases due to stage separation and the multiple flight modes are separated into two parts, the thrusting flight phase (TFP) and the coasting flight phase (CFP). ¹ Besides, SLV flight mission requires high fault tolerances due to safety issues and high cost. Because of the aforementioned intrinsic attributes of the SLV flight mission, SLV tracking systems (TSs) comprise multiple heterogeneous sensors such as radars, telemetry systems receiving onboard Global Positioning System-inertial navigation system (GPS-INS) data, and electro-optical targeting systems (EOTSs) installed over widespread areas. In this situation, the mission control system (MCS), ² which is a main data processing station of TSs, should concern two subjects: multiple models of SLV dynamics and distributed multiple sensor data.

Real-time SLV tracking problem has made use of single model filters such as an alpha-beta-gamma filter or Kalman filter (KF), ³ because of a reliability issue. However, since a KF only utilizes a single dynamic model, in general, the tracking performance of KF for a maneuvering launch vehicle is inferior to multiple model estimators. ⁴ Hence, the dynamic model of an SLV should be described in multiple models so that the multiple flight phases are properly accounted for. 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 short-range projectiles, sounding rockets, and SLV.^1,5–10 Short-range ballistic munitions or projectiles with multiple models such as spin-stabilized and fin-stabilized models were implemented using an interacting multiple model (IMM) estimator.^5,6 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. ⁷ The first mode was a constant axial force model for the boosting and re-entry phases. The second mode was a ballistic acceleration model that incorporated the gravitational, Coriolis, and centripetal forces for the exo-atmospheric phase. The final mode was a standard autocorrelated acceleration Singer’s model for malfunction motions of missiles such as re-entry tumbling. In Song and Han, ¹ comparison of two distinct IMM algorithms that incorporate different types of multiple models to describe the motion of SLV is presented. Here, two different multiple models of SLV are represented using linear singer model with ternary uniform mixture and nonlinear ballistic model with unbiased mixing.

Various linear and nonlinear multiple sensor fusion schemes have been developed over the past three decades. Classical linear fusion schemes are the measurement fusion (MF) ¹¹ and state vector fusion (SVF). ¹² Theoretically, the former gives optimal estimated state vector. In practical situations, however, MF is not feasible due to the communication bandwidth. ¹³ The SVF scheme had been applied in the SLV tracking problem, which is formulated using constant velocity model with multiple converted measurements of radar. ¹³ Since this algorithm integrates multiple sensor data after converting into common coordinate system, that is, earth-centered earth fixed (ECEF), it cannot be applied to different types of sensors such as angle-only or range-only measurements. Nonlinear fusion scheme based on information filtering, which is essentially a KF stated in terms of the inverse of the covariance matrix,^14,15 can be a breakthrough. Furthermore, some literatures^16,17 tried to fuse multiple model problems related to application of mobile node localization but not involved in SLV tracking. From the aforementioned literature survey, we can pay attention to the absence in those researches concerning both multiplicity of dynamic model and heterogeneity of observation simultaneously, even though the multiple model and multiple sensor data fusion of the distributed TSs help MCS have both effectiveness and robustness in the SLV flight mission.

To this end, for accurate SLV tracking, we propose a distributed extended information filter-based IMM (DEIF-IMM) estimator with unbiased mixing. The proposed algorithm is expressed via nonlinear multiple models concerning comprehensive factors of the SLV specifications to be collected in advance such as thrust, gravity, drag coefficient, Mach number, and air density.^7–10 Here, multiple motions of the SLV are separated into two parts, TFP and CFP, having different dimensions of state vector. To mix and fuse multiple state vectors and measurements, unbiased mixing ⁸ and distributed extended information filtering ¹⁴ are applied, respectively. By incorporating both approaches to the distributed SLV tracking problem, we can expect better tracking performance and stability since the proposed algorithm simultaneously concerns heterogeneity of tracking sensors and multiplicity of vehicle’s dynamic model. To show the improved performance compared to conventional tracking algorithms, we illustrate comparisons between benchmark algorithms such as MF and SVF using converted measurements and the proposed algorithm. For the practical SLV track generation, a nominal trajectory of Korea Space Launch Vehicle-1 (KSLV-1) concerning six degree of freedom (DoF) launch vehicle dynamics and stage separation events is adopted in the simulation part.

The remainder of this article is organized as follows. Section “Problem statement for SLV tracking” presents a statement of the problem for nonlinear multiple model and measurements for SLV tracking. In section “Classical linear multiple sensor fusion algorithm,” the conventional SLV tracking approach is illustrated. In section “Proposed DEIF-IMM algorithm,” the proposed DEIF-IMM algorithm is derived. Section “Simulation results” presents simulation settings and results; a comparison of the proposed algorithm and conventional algorithm is depicted as an aspect of root mean square (RMS) error of the state vector. Finally, in section “Conclusion,” the conclusions of this article are presented.

Problem statement for SLV tracking

Introduction of KSLV-1

KSLV-1 is the first South Korean SLV composed of two stages: a liquid propellant first stage developed by Khrunichev State Research and Production Space Centre and a solid propellant second stage developed in Korea by Korea Aerospace Research Institute. ¹ The objective of KSLV-1 is to inject a payload of 100 kg to low earth orbit with an inclination angle of 80 degree. KSLV-1 had launched successfully in Naro Space Center on 30 January 2013.

Multiple heterogeneous sensors in Naro Space Center

Due to wide flight range and successful tracking of SLV, TS is composed of distributed multiple heterogeneous sensors. In Naro Space Center, there are three kinds of TSs such as radar, telemetry, and EOTS. C-band pulse radars are operated on two different modes, skin mode and beacon mode. A skin mode uses returns of transmitted pulses of radar, but a beacon mode communicates with a transponder on SLV. A real SLV tracking mission applies a beacon mode because of its superior tracking performance. The onboard GPS-INS and status data of SLV can be downloaded in real-time. In addition, for the purpose of antenna pointing toward SLV, mono-pulse telemetry antenna gives azimuth and elevation angle measurements. EOTS delivers not only image data but also range measurements from Doppler radar equipped in the system. MCS gathers all tracking data for processing tracking filters and sending slaving data for preventing tracking failure. All TSs are located in remote areas as shown in Figure 1 and the distances between station #1 and stations #2 and #3 are 150 and 3100 km, respectively.

Figure 1.

Illustration of problem statement for SLV tracking.

Classical linear multiple sensor fusion algorithm

Converted measurement KF

The linear motion of the SLV is simply depicted as a discretized Wiener process acceleration model ¹⁸

$$X_{k+1}=F{X}_k+w_k$$

(1)

where the state vector $X_k\in {\Re }^9$ consists of the position, velocity, and acceleration components along the x, y, and z axes, respectively, that is, $X_k=[x_k{\stackrel{·}{x}}_k{\stackrel{··}{x}}_k{y}_k{\stackrel{·}{y}}_k{\stackrel{··}{y}}_k{z}_k{\stackrel{·}{z}}_k{\stackrel{··}{z}}_k{]}^T$ . The system matrix $F$ and covariance matrix of the system (white Gaussian) noise $w_k$ are represented as equations (2) and (3), respectively

$$f_{L=x,y,z}=\left[\begin{array}{ccc}\hfill 1\hfill & \hfill \Delta t\hfill & \hfill \Delta t^2/2\hfill \\ \multicolumn{1}{c}{0}& \hfill 1\hfill & \hfill \Delta t\hfill \\ \multicolumn{1}{c}{0}& \hfill 0\hfill & \hfill 1\hfill \end{array}\right],F=\left[\begin{array}{ccc}\hfill f_x\hfill & \hfill 0_3\hfill & \hfill 0_3\hfill \\ \multicolumn{1}{c}{0_3}& \hfill f_y\hfill & \hfill 0_3\hfill \\ \multicolumn{1}{c}{0_3}& \hfill 0_3\hfill & \hfill f_z\hfill \end{array}\right]$$

(2)

$$\begin{array}{c}E[w_k{w}_l^T]=Q{\delta }_{k-l}\hfill \\ q_{L=x,y,z}=\left[\begin{array}{ccc}\hfill \Delta t^5/20\hfill & \hfill \Delta t^4/8\hfill & \hfill \Delta t^2/6\hfill \\ \multicolumn{1}{c}{\Delta t^4/8}& \hfill \Delta t^3/3\hfill & \hfill \Delta t^2/2\hfill \\ \multicolumn{1}{c}{\Delta t^3/6}& \hfill \Delta t^2/2\hfill & \hfill \Delta t\hfill \end{array}\right],\hfill \\ Q=\left[\begin{array}{ccc}\hfill q_x\hfill & \hfill 0_3\hfill & \hfill 0_3\hfill \\ \multicolumn{1}{c}{0_3}& \hfill q_y\hfill & \hfill 0_3\hfill \\ \multicolumn{1}{c}{0_3}& \hfill 0_3\hfill & \hfill q_z\hfill \end{array}\right]\hfill \end{array}$$

(3)

A measurement model of multiple radars for the SLV gives the spherical coordinate observations such that

$$\begin{array}{c}{Z}_k^{(s)}=h(X_k)+v_k^{(s)},Z_k^{(s)}\in {\Re }^3,v_k^{(s)}\in {\Re }^3,i=1,...,n\hfill \\ Z_k^{(s)}=\left[\begin{array}{c}\hfill r_k^{(s)}\hfill \\ \multicolumn{1}{c}{{\phi }_k^{(s)}}\\ \multicolumn{1}{c}{{\theta }_k^{(s)}}\end{array}\right]=\left[\begin{array}{c}\hfill \sqrt{x_k^2+y_k^2+z_k^2}\hfill \\ \multicolumn{1}{c}{{\tan }^{-1}(y_k/x_k)}\\ \multicolumn{1}{c}{{\tan }^{-1}(y_k/\sqrt{x_k^2+y_k^2})}\end{array}\right]+\left[\begin{array}{c}\hfill v_r^{(s)}\hfill \\ \multicolumn{1}{c}{v_{\phi }^{(s)}}\\ \multicolumn{1}{c}{v_{\theta }^{(s)}}\end{array}\right]\hfill \end{array}$$

(4)

where $v_k^{(s)}~\mathbb{N}(0,R_k^{(s)})$ , $R_k^{(s)}=\mathrm{diag}({\sigma }_r^{{(s)}^2},{\sigma }_{\phi }^{{(s)}^2},{\sigma }_{\theta }^{{(s)}^2})$ , and $s$ means index of multiple sensors. Using a three-dimensional (3D) de-biased converted measurement approach, ¹⁹ we can transform the original nonlinear equation (4) into linear form as

$$\begin{array}{c}{\mathbf{z}}_k^{(s)}=H{X}_k+v_{c,k}^{(s)},v_{c,k}^{(s)}\in {\Re }^3\hfill \\ {\mathbf{z}}_k^{(s)}=\left[\begin{array}{ccc}\hfill 1\hfill & \hfill 0_{1\times 5}\hfill & \hfill 0_{1\times 3}\hfill \\ \multicolumn{1}{c}{0_{1\times 3}}& \hfill 1\hfill & \hfill 0_{1\times 5}\hfill \\ \multicolumn{1}{c}{0_{1\times 5}}& \hfill 1\hfill & \hfill 0_{1\times 2}\hfill \end{array}\right]X_k+\left[\begin{array}{c}\hfill v_{x,k}^{(s)}\hfill \\ \multicolumn{1}{c}{v_{y,k}^{(s)}}\\ \multicolumn{1}{c}{v_{z,k}^{(s)}}\end{array}\right]\hfill \end{array}$$

(5)

Here, $v_{c,k}^{(s)}$ is the converted measurement noise expressed in terms of Cartesian coordinates, that is, $v_{c,k}^{(s)}~\mathbb{N}(0,R_{c,k}^{(s)})$

$$R_{c,k}^{(s)}=\left[\begin{array}{ccc}\hfill R_{\mathrm{xx},k}^{(s)}\hfill & \hfill R_{\mathrm{xy},k}^{(s)}\hfill & \hfill R_{\mathrm{xz},k}^{(s)}\hfill \\ \multicolumn{1}{c}{R_{\mathrm{yx},k}^{(s)}}& \hfill R_{\mathrm{yy},k}^{(s)}\hfill & \hfill R_{\mathrm{yz},k}^{(s)}\hfill \\ \multicolumn{1}{c}{R_{\mathrm{zx},k}^{(s)}}& \hfill R_{\mathrm{zy},k}^{(s)}\hfill & \hfill R_{\mathrm{zz},k}^{(s)}\hfill \end{array}\right]$$

(6)

Measurement fusion and state vector fusion

Measurement fusion

MF algorithm fuses the sensor observation directly via an augmented measurement model and applies one KF as shown in Figure 2. ¹¹ The overall radar observation ${\mathbf{z}}_k$ is composed of $n$ multiple radar observation subvectors (local sensors), ${\mathbf{z}}_k^{(1)},...,{\mathbf{z}}_k^{(n)}$ , that is

$${\mathbf{z}}_k={[{{\mathbf{z}}_k^{(1)}}^T,...,{{\mathbf{z}}_k^{(n)}}^T]}^T$$

(7)

MF is the estimation of the state $X_k$ at current time measurements in mean-square sense, such that

$${\hat{X}}_{k|k}^{\mathrm{MF}}=E(X_k|{\mathbf{z}}_{[0,k]}),{\mathbf{z}}_{[0,k]}=\{{\mathbf{z}}_l,l=0,...,k\}$$

(8)

The recursion of MF is computed by equations (9) and (10). State and covariance time propagations are given as

$$\begin{array}{c}{\hat{X}}_{k+1|k}=F{\hat{X}}_{k|k}\hfill \\ P_{k+1|k}={\mathrm{FP}}_{k+1|k}{F}^T+Q\hfill \end{array}$$

(9)

The state and covariance measurement data updates are given as

$$\begin{array}{c}{K}_{k+1}=P_{k+1|k}{H}_F^T[H_F{P}_{k+1|k}{H}_F^T+R_k]\hfill \\ {\hat{X}}_{k+1|k+1}^{\mathrm{MF}}={\hat{X}}_{k+1|k}+K_{k+1}[{\mathbf{z}}_{k+1}-H_F{\hat{X}}_{k+1|k}]\hfill \\ P_{k+1|k+1}^{\mathrm{MF}}=[I_n-K_{k+1}{H}_F]P_{k+1|k}\hfill \end{array}$$

(10)

where $H_F=[H^T,...,H^T{]}^T$ and $R_k=\mathrm{diag}(R_{c,k}^{(1)},...,R_{c,k}^{(n)})$ .

Figure 2.

Measurement fusion architecture.

State vector fusion

SVF algorithm applies KF for individual local sensors and generates global fused state estimates using local estimates as shown in Figure 3. ¹² The local KF estimates are calculated by equations (9) and (10) using local sensor measurements ${\mathbf{z}}_k^{(s)},s=1,...,n$ instead of ${\mathbf{z}}_k$ . SVF of the state $X_k$ based on the overall radar observation (7) is constructed from the local KF estimates ${\hat{X}}_{k|k}^{(s)},s=1,...,n$ , ¹¹ such that

$${\hat{X}}_{k|k}^{\mathrm{SVF}}=\sum _{s=1}^n{a}_k^{(s)}{\hat{X}}_{k|k}^{(s)},\sum _{s=1}^n{a}_k^{(s)}=I_{n_l}$$

(11)

where $a_k^{(1)},...,a_k^{(n)}$ are $N\times N$ time-varying matrix weights determined from the mean-square criterion

$$J^{\mathrm{SVF}}=E\left({\Vert X_k-\sum _{s=1}^n{a}_k^{(s)}{\hat{X}}_{k|k}^{(s)}\Vert }^2\right)\to \underset{a_k^{(s)}}{min}$$

(12)

Criterion (12) to find optimal matrix weights $a_k^{(1)},...,a_k^{(n)}$ can be solved by minimizing the trace of mean square error and its derivation has been carried out various methods such as generalized Millman’s formula ¹² and Lagrange multiplier. ²⁰

Figure 3.

State vector fusion architecture.

From the linear model fusion, we can check that even though there are additional tracking data such as angle measurements from telemetry and range measurements from EOTS, the conventional linear fusion algorithm cannot handle all data due to coordinate conversion. To overcome this problem, we propose DEIF-IMM estimator with unbiased mixing.

Proposed DEIF-IMM algorithm

Multiple dynamic models of SLV

The trajectory of an SLV can be separated into two parts: thrusting and coasting. Two distinct dynamics of SLV are considerably dissimilar and characterized by different main forces in each phase. The main forces can be represented as thrust, drag, and gravity components.^1,8–10 For the thrusting mode, three forces such as thrust, drag, and gravity act on SLV simultaneously, but in the coasting mode, only gravity and drag have an effect on SLV. Therefore, the motion of SLV can be depicted as nonlinear multiple model systems. For a nonlinear SLV dynamic model, the state vector is denoted as

$${\mathbf{x}}_t=[x_t{y}_t{z}_t{\stackrel{·}{x}}_t{\stackrel{·}{y}}_t{\stackrel{·}{z}}_t{\xi }_t{\tau }_t{]}^T$$

(13)

where ${\xi }_t$ is the drag coefficient and ${\tau }_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^8,9 such that

$$\begin{array}{c}\left[\begin{array}{c}\hfill \stackrel{··}{x}\hfill \\ \multicolumn{1}{c}{\stackrel{··}{y}}\\ \multicolumn{1}{c}{\stackrel{··}{z}}\end{array}\right]=\frac{\tau }{V}\underset{\mathrm{thrust}\mathrm{term}}{\left[\begin{array}{c}\hfill \stackrel{·}{x}\hfill \\ \multicolumn{1}{c}{\stackrel{·}{y}}\\ \multicolumn{1}{c}{\stackrel{·}{z}}\end{array}\right]}+\underset{\mathrm{drag}\mathrm{term}}{\xi {\xi }_{m}D\left[\begin{array}{c}\hfill \stackrel{·}{x}\hfill \\ \multicolumn{1}{c}{\stackrel{·}{y}}\\ \multicolumn{1}{c}{\stackrel{·}{z}}\end{array}\right]}+\underset{\mathrm{gravity}\mathrm{term}}{g\left[\begin{array}{c}\hfill 0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{1}\end{array}\right]}+{\tilde{\omega }}_1\hfill \\ \stackrel{·}{\xi }={\tilde{\omega }}_2\hfill \\ \stackrel{·}{\tau }={\tilde{\omega }}_3\hfill \end{array}$$

(14)

The first term on the right-hand side of equation (14) represents the thrust $(m^2/s)$ of the SLV in the x, y, and z directions. V is the magnitude of the velocity $v=[\stackrel{·}{x}\stackrel{·}{y}\stackrel{·}{z}{]}^T$ , that is, the SLV speed $(m/s)$ . The second part of equation (14) is the drag term, which is related to velocity and altitude, that is, $D=-\rho (z)V/2$ , where $\rho (z)={\rho }_0{e}^{-\mathrm{cz}}$ is the air density (kg/m³) at an altitude z (m) and c is the air density constant $(m^{-1})$ . ⁶ $\xi$ is the drag coefficient and ${\xi }_m$ is the Mach number-dependent drag coefficient multiplier, which is approximated by the cubic spline curve shown in Figure 4. Since the SLV tested in the proposed algorithm belongs to a small-sized launch vehicle with a nose-shaped appearance, we adopt the drag coefficient of projectiles with the nose shape. ²¹ The third part of equation (14) is a gravity term. Gravity $g$ is the standard acceleration due to gravity at sea level, which is assumed constant throughout the trajectory, with a value of 9.812 m/s ² . ${\tilde{\omega }}_1$ , ${\tilde{\omega }}_2$ , and ${\tilde{\omega }}_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.^8,9 We can modify the dynamic equations (13) and (14) in a compact form such that

$${\stackrel{·}{\mathbf{x}}}_t=f[{\mathbf{x}}_t]+{\tilde{\omega }}_t$$

(15)

where

$$f[{\mathbf{x}}_t]=\left[\begin{array}{c}\hfill {\stackrel{·}{x}}_t\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{y}}_t}\\ \multicolumn{1}{c}{{\stackrel{·}{z}}_t}\\ \multicolumn{1}{c}{\tau \frac{{\stackrel{·}{x}}_t}{V_t}+{\xi }_t{D}_t{\stackrel{·}{x}}_t}\\ \multicolumn{1}{c}{\tau \frac{{\stackrel{·}{y}}_t}{V_t}+{\xi }_t{D}_t{\stackrel{·}{y}}_t}\\ \multicolumn{1}{c}{\tau \frac{{\stackrel{·}{z}}_t}{V_t}+{\xi }_t{D}_t{\stackrel{·}{z}}_t-g}\end{array}\right],{\tilde{\omega }}_t=[{\tilde{\omega }}_{1,t}{\tilde{\omega }}_{2,t}{\tilde{\omega }}_{3,t}{]}^T$$

Figure 4.

Normalized drag coefficient. ²¹

The state vector equation (15) is discretized by a second-order Taylor expansion. ²² Then, equation (15) can be written as a discretized continuous time system with white process noise such that

$${\mathbf{x}}_{k+1}={\mathbf{x}}_k+f[{\mathbf{x}}_k]\Delta t+A_{k}f[{\mathbf{x}}_k]\frac{\Delta t^2}{2}+{\omega }_k$$

(16)

where $A_k$ is the Jacobian of equation (15) evaluated at ${\mathbf{x}}_k$ ²² and ${\omega }_k$ is the discretized continuous time process noise for the sampling interval $\Delta t$ . The corresponding covariance matrix of the discretized process noise is

$$Q_k=\left[\begin{array}{ccc}\hfill Q_{\mathrm{pv}}{q}_v\hfill & \hfill 0_{6\times 1}\hfill & \hfill 0_{6\times 1}\hfill \\ \multicolumn{1}{c}{0_{1\times 6}}& \hfill \Delta t{q}_{\xi }\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0_{1\times 6}}& \hfill 0\hfill & \hfill \Delta t{q}_{\tau }\hfill \end{array}\right],Q_{\mathrm{pv}}=\left[\begin{array}{cc}\hfill \frac{\Delta t^3}{3}{I}_3\hfill & \hfill \frac{\Delta t^2}{2}{I}_3\hfill \\ \multicolumn{1}{c}{\frac{\Delta t^2}{2}{I}_3}& \hfill \Delta t{I}_3\hfill \end{array}\right]$$

(17)

where $I_3$ is the $3\times 3$ identity matrix and the continuous time process noise intensities $q_v$ , $q_{\xi }$ , and $q_{\tau }$ are the corresponding power spectral densities.

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 TFP, whereas zero thrust signifies CFP. This means we can rewrite system model (16) as Markov jump nonlinear system models, where two models of multiple model set $i=\{\mathrm{TFP},\mathrm{CFP}\}$ obey the following equation as

$${\mathbf{x}}_{i,k+1}={\mathbf{x}}_{i,k}+f[{\mathbf{x}}_{i,k}]\Delta t+A_{i,k}f[{\mathbf{x}}_{i,k}]\frac{\Delta t^2}{2}+{\omega }_{i,k}$$

(18)

where

$$\begin{array}{c}\mathrm{cov}({\omega }_{i,k})=Q_{i,k},Q_{\mathrm{TFP},k}=Q_k,Q_{\mathrm{TFP},k}=\left[\begin{array}{cc}\hfill Q_{\mathrm{pv}}{q}_v\hfill & \hfill 0_{6\times 1}\hfill \\ \multicolumn{1}{c}{0_{1\times 6}}& \hfill \Delta t{q}_{\xi }\hfill \end{array}\right]\hfill \\ \mathrm{Thrusting}\mathrm{flight}\mathrm{mode}:{\mathbf{x}}_{i=\mathrm{TFP},k}=[x_k{y}_k{z}_k{\stackrel{·}{x}}_k{\stackrel{·}{y}}_k{\stackrel{·}{z}}_k{\xi }_k{\tau }_k{]}^T\hfill \\ \mathrm{Coasting}\mathrm{flight}\mathrm{mode}:{\mathbf{x}}_{i=\mathrm{CFP},k}=[x_k{y}_k{z}_k{\stackrel{·}{x}}_k{\stackrel{·}{y}}_k{\stackrel{·}{z}}_k{\xi }_k{]}^T\hfill \end{array}$$

(19)

Multiple measurement models

Basically, the proposed algorithm uses nonlinear spherical coordinate measurement equation (4), but to apply extended information filter, we use Jacobian of equation (4) at predicted state ${\hat{\mathbf{x}}}_{k+1|k}$ for both TFP and CFP modes such that

$${\mathbf{y}}_k^{(s)}=h^{(s)}({\mathbf{x}}_{i,k})+v_k^{(s)},{\mathbf{y}}_k^{(s)}\in {\Re }^3,s=1,...,n$$

(20)

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{H}}_{\mathrm{TFP},\mathrm{rad}}^{(s)}={{[\nabla {\mathbf{x}}_k{h}^{{(\mathrm{rad})}^T}({\mathbf{x}}_k)]}^T|}_{{\mathbf{x}}_k={\hat{\mathbf{x}}}_{k+1|k}}}\\ =\left[\begin{array}{cccccccc}\hfill \frac{x_k}{\sqrt{x_k^2+y_k^2+z_k^2}}\hfill & \hfill \frac{y_k}{\sqrt{x_k^2+y_k^2+z_k^2}}\hfill & \hfill \frac{z_k}{\sqrt{x_k^2+y_k^2+z_k^2}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{\frac{-y_k}{x_k^2+y_k^2}}& \hfill \frac{x_k}{x_k^2+y_k^2}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{\frac{-x_k{z}_k}{(x_k^2+y_k^2+z_k^2)\sqrt{x_k^2+y_k^2}}}& \hfill \frac{-y_k{z}_k}{(x_k^2+y_k^2+z_k^2)\sqrt{x_k^2+y_k^2}}\hfill & \hfill \frac{\sqrt{x_k^2+y_k^2}}{x_k^2+y_k^2+z_k^2}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]\hfill \end{array}$$

(21)

$$\begin{array}{c}{\hat{H}}_{\mathrm{CFP},\mathrm{rad}}^{(s)}={{[\nabla {\mathbf{x}}_k{h}^{{(\mathrm{rad})}^T}({\mathbf{x}}_k)]}^T|}_{{\mathbf{x}}_k={\hat{\mathbf{x}}}_{k+1|k}}\hfill \\ =\left[\begin{array}{ccccccc}\hfill \frac{x_k}{\sqrt{x_k^2+y_k^2+z_k^2}}\hfill & \hfill \frac{y_k}{\sqrt{x_k^2+y_k^2+z_k^2}}\hfill & \hfill \frac{z_k}{\sqrt{x_k^2+y_k^2+z_k^2}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{\frac{-y_k}{x_k^2+y_k^2}}& \hfill \frac{x_k}{x_k^2+y_k^2}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{\frac{-x_k{z}_k}{(x_k^2+y_k^2+z_k^2)\sqrt{x_k^2+y_k^2}}}& \hfill \frac{-y_k{z}_k}{(x_k^2+y_k^2+z_k^2)\sqrt{x_k^2+y_k^2}}\hfill & \hfill \frac{\sqrt{x_k^2+y_k^2}}{x_k^2+y_k^2+z_k^2}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]\hfill \end{array}$$

(22)

Note that the aforementioned linear multiple sensor fusion algorithms cannot fuse additional heterogeneous measurements such as angle measurements from telemetry and range measurements from EOTS. The proposed algorithm, however, resolved the problem by incorporating the information filter where the fusion of multiple sensors executed simple summation of local states and information matrix contributions. Hence, we can also define Jacobian matrices for angle measurements of telemetry ${\hat{H}}_{i,\mathrm{tlm}}^{(s)}$ and range measurements of EOTS ${\hat{H}}_{i,\mathrm{eots}}^{(s)}$ as follows:

• ${\hat{H}}_{i,\mathrm{tlm}}^{(s)}$ : $2\times 8(\mathrm{or}2\times 7)$ matrix with removal of the first row of ${\hat{H}}_{i,\mathrm{rad}}^{(s)}$ where i = TFP (or i = CFP).

• ${\hat{H}}_{i,\mathrm{eot}}^{(s)}$ : $1\times 8(\mathrm{or}1\times 7)$ matrix with removal of the first and second rows of ${\hat{H}}_{i,\mathrm{rad}}^{(s)}$ where i = TFP (or i = CFP).

A DEIF-IMM algorithm for nonlinear dynamics with different sizes of the mode state vector is summarized as follows (Figure 5).

Distributed extended information filtering

• (i) Model-conditioned reinitialization (for $i=\mathrm{TFP},\mathrm{CFP}$ )

•  Predicted mode probability

$${\mu }_{i,k|k-1}\stackrel{\Delta }{=}⪻\{m_{i,k}|{\mathbf{y}}_{k-1}\}=\sum _j{\pi }_{\mathrm{ji}}{\mu }_{j,k-1}$$

(23)

• Mixing weight

$${\mu }_{j|i,k-1}\stackrel{\Delta }{=}⪻\{m_{j,k-1}|m_{i,k},{\mathbf{y}}_{k-1}\}=\sum _j{\pi }_{\mathrm{ji}}{\mu }_{j,k-1}/{\mu }_{i,k-1}$$

(24)

• Unbiased mixing estimate and covariance

$$\begin{array}{c}{\overline{\mathbf{x}}}_{i,k-1|k-1}\stackrel{\Delta }{=}E[{\mathbf{x}}_{k-1|k-1}|M_{i,k},{\mathbf{y}}_{k-1}]=\sum _j{\hat{\chi }}_{j,k-1|k-1}{\mu }_{j|i,k-1},\hfill \\ \mathrm{where}{\hat{\chi }}_{\mathrm{TFP},k-1|k-1}\stackrel{\Delta }{=}{\hat{\mathbf{x}}}_{\mathrm{TFP},k-1|k-1},{\hat{\chi }}_{\mathrm{CFP},k-1|k-1}\stackrel{\Delta }{=}[{\hat{\mathbf{x}}}_{\mathrm{CFP},k-1|k-1}^T,{\tau }_{k-1}{]}^T\hfill \\ {\overline{P}}_{i,k-1|k-1}=\sum _j[P_{i,k-1|k-1}+({\overline{\mathbf{x}}}_{i,k-1|k-1}-{\hat{\chi }}_{j,k-1|k-1}){({\overline{\mathbf{x}}}_{i,k-1|k-1}-{\hat{\chi }}_{j,k-1|k-1})}^T]{\mu }_{j|i,k-1}\hfill \end{array}$$

(25)

• (ii) Model-conditioned DEIF (for $i=\mathrm{TFP},\mathrm{CFP}$ )

• Predicted information state estimate and covariance

$$\begin{array}{c}{\mathit{y}}_{i,k|k-1}={\left(P_{i,k|k-1}\right)}^{-1}{\hat{\mathbf{x}}}_{i,k|k-1}\hfill \\ {\mathit{Y}}_{i,k|k-1}={\left(P_{i,k|k-1}\right)}^{-1}\hfill \\ ={({\Phi }_{i,k-1}{\overline{P}}_{i,k-1|k-1}{\Phi }_{i,k-1}^T+Q_{i,k-1})}^{-1},\hfill \\ {\Phi }_{i,k-1}=I+{\hat{A}}_{i,k-1}\hfill \end{array}$$

(26)

where

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{\mathbf{x}}}_{i,k|k-1}={\overline{\mathbf{x}}}_{i,k-1|k-1}+f[{\overline{\mathbf{x}}}_{i,k-1|k-1}]\Delta t}\\ +{\hat{A}}_{i,k-1}f[{\overline{\mathbf{x}}}_{i,k-1|k-1}]\frac{\Delta t^2}{2},{\hat{A}}_{i,k}={\frac{\partial f}{\partial \mathbf{x}}|}_{\mathbf{x}={\overline{\mathbf{x}}}_{i,k}}\hfill \end{array}$$

(27)

• Measurement residual

$${\tilde{\mathbf{y}}}_{i,k}^{(s)}={\mathbf{y}}_k^{(s)}-h^{(s)}({\hat{\mathbf{x}}}_{i,k|k-1})$$

(28)

• Information state contribution

$${\mathbf{i}}_{i,k}^{(s)}={\hat{H}}_i^{(s)}{TR}_k^{(s)}-1({\tilde{\mathbf{y}}}_{i,k}^{(s)}+{\hat{H}}_i^{(s)}{\hat{\mathbf{x}}}_{i,k|k-1})$$

(29)

• Information matrix contribution

$${\mathbf{I}}_{i,k}^{(s)}={\hat{H}}_i^{(s)}{TR}_k^{(s)-1}(s){\hat{H}}_i^{(s)}$$

(30)

• Update information state and information matrix

$${\mathit{y}}_{i,k|k}={\mathit{y}}_{i,k|k\_1}+\sum _{s=1}^n{\mathbf{i}}_{i,k}^{(s)},{\mathit{Y}}_{i,k|k}={\mathit{Y}}_{i,k|k-1}+\sum _{s=1}^n{\mathbf{I}}_{i,k}^{(s)}$$

(31)

• (iii) Mode probability update (for $i=\mathrm{TFP},\mathrm{CFP}$ )

• Mode likelihood

$${\Lambda }_{i,k}^{(s)}\stackrel{\Delta }{=}p[z_{i,k}^{(s)}|M_{i,k}^{(s)},{\mathbf{y}}_{k-1}^{(s)}]\stackrel{\mathbf{assume}}{=}\mathbb{N}({\mathbf{y}}_k^{(s)};0,S_{i,k}^{(s)})$$

(32)

where residual covariance is $S_{i,k}^{(s)}={\hat{H}}_i^{(s)}{P}_{i,k|k-1}^{(s)}{\hat{H}}_i^{(s)}T-R_k^{(s)}$

• Mode probability

$${\mu }_{i,k}=\frac{{\mu }_{i,k|k-1}{\Lambda }_{i,k}}{\sum _j{\mu }_{j,k|k-1}{\Lambda }_{j,k}},{\Lambda }_{i,k}=\sum _s{\Lambda }_{i,k}^{(s)}$$

(33)

• (iv) Combination (for $i=\mathrm{TFP},\mathrm{CFP}$ )

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{\mathbf{x}}}_{k|k}=\sum _i{\hat{\mathbf{x}}}_{i,k|k}{\mu }_{i,k},{\hat{P}}_{k|k}}\\ =\sum _i[P_{i,k|k}+({\hat{\mathbf{x}}}_{k|k}-{\hat{\mathbf{x}}}_{i,k|k}){({\hat{\mathbf{x}}}_{k|k}-{\hat{\mathbf{x}}}_{i,k|k})}^T]{\mu }_{i,k}\hfill \end{array}$$

(34)

Figure 5.

An architecture of the proposed DEIF-IMM algorithm for SLV tracking.

Simulation results

To demonstrate the performance of the proposed algorithm, we simulated the SLV tracking problem based on the nominal flight trajectory of KSLV-1 considering various facts such as 6-DoF launch vehicle dynamics and wind profiles in distributed sensor environment. Although we have real tracking data of KSLV-1, the real data are not applied in this simulation section. Since we cannot exactly know the realization of flight trajectory of KSLV-1, we assume that nominal flight trajectory is a realization of KSLV-1 as an alternative rather than using real data. In the simulation, only three measurements are used, that is, s = 1, 2, 3 where s = 1: radar #1, s = 2: radar #2, s = 3 telemetry #1. Measurement noise intensities of equation (4) are selected as ${\sigma }_r^{(s)}=2.5^{2}m$ , ${\sigma }_{\phi }^{(s)}=0.{01}^2\mathrm{deg}$ , and ${\sigma }_{\theta }^{(s)}=0.{01}^2\mathrm{deg}$ . When s = 3 it means telemetry’s angle-only measurements ( ${\sigma }_r^{(3)}$ does not exist). Sampling interval $\Delta t=0.1\mathrm{s}$ is used and for the multiple model setting, transition probability matrix for TFP and CFP is

$$\pi =\left[\begin{array}{cc}\hfill 0.9\hfill & \hfill 0.1\hfill \\ \multicolumn{1}{c}{0.1}& \hfill 0.9\hfill \end{array}\right]$$

with initial mode probability vector [0.9,0.1]. As shown in Table 1, the nominal flight sequence of KSLV-1 is as follows. First, the payload fairing is separated during the first stage flight at 215.0 s. After the first stage engine shutdown at 229.0 s, the upper stage is separated from the first stage and the vehicle enters the CFP mode. The second stage continues in the CFP mode until the kick motor ignition at 395.0 s, and the vehicle then enters the TFP mode. At the end of kick motor combustion (455.0 s), the upper stage enters the target orbit in CFP mode as in the previous separation. Finally, the satellite is inserted into the target orbit after it separates from the upper stage during CFP mode at 540 s ¹ and references therein. The simulation started at countdown time (CT) + 20 s, where KSLV_I flies over 20 s after lift-off and, as shown in the event timetable, the event is finished within 540 s. Hence, we focus on this regime including margin 80 s (CT + 20–620). To demonstrate the performance of DEIF-IMM, we compare RMS error of position estimation using two different conventional fusion methods, MF and SVF, with 100 Monte Carlo runs.

Table 1.

Flight sequence of KSLV-1.

Event	Flight time (s)	Altitude (km)
Lift-off	0.0	0.1
March = 1	54.0	7.2
Maximum dynamic pressure	64.0	10.5
Pairing separation	215.0	177
First stage engine cutoff	229.0	194
First stage separation	232.0	196
Second engine (kick motor) start	395.0	303
Second engine cutoff	455.0	305
Satellite separation	540.0	305

KSLV-1: Korea Space Launch Vehicle-1.

Figures 6–8 show the RMS error of Cartesian coordinate (x, y, and z axes) from three algorithms such as the proposed algorithm, MF, fusion, and SVF fusion algorithm. Although the measurement coverage of each TS is different in a real launch mission, we assume that all TSs can simultaneously cover all flight regimes in this simulation part. We can check that tracking performance of the proposed algorithm outperforms especially after CT 200 s. The reason is that in fact, tracking performance of IMM, that is, a main estimator of the proposed algorithm is better than KF when the maneuvering index ⁴ is large. In this simulation, as explained in the flight sequence of KSLV-1, the motion of the SLV is described by multiple dynamics (in our case TFP and CFP) rather than single dynamics, whereas the dynamic model of classical fusion (benchmark) algorithms used in this article is the constant acceleration (CA) model. This means that the maneuvering index of CFP is very small, that is, the system noise intensity of CFP is configured to have value very similar to benchmark algorithms, that is, KF. In the TFP, however, the maneuvering index is very large due to enormous thrust force. This means during the CFP, tracking performance of both the proposed algorithm and benchmark algorithms could be shown to be similar. On the other hand, during TFP tracking performance of the proposed algorithm should be better than benchmark algorithms. Despite this fact, we can often see the lack of performance enhancement within certain thrusting phase areas. The reason is that we adopted multistage (two stages) SLV in our simulation example. Two distinct (first and second) thrust phases have different amounts of acceleration, that is, to exactly describe multiple phases of KSLV-1 we need one more thrusting model. In this article, however, we have adopted two models for simple implementation. Furthermore, the classical algorithm cannot utilize extra measurements from TSs due to coordinate conversion. This is why the performance of the MF and SVF fusion algorithm is worse than the proposed algorithm. Nonetheless, we can see that the y-axis tracking errors shown in Figure 7 are relatively small. This result is observed because the y-axis motion of KSLV-1 in a local coordinate is relatively smaller than the other axis motion, and the CA model approximately expresses this small motion. On the other hand, we can see that only in second stage combustion, the tracking performance of MF and SVF is better than the proposed algorithm. The proposed multiple model algorithm depicts multiple modes of flight phases using different thrust and drag noise factors, but there are also quite large differences in noise factors between the first and second TFPs. To exactly describe the second TFP, one needs an additional mode with its own parameters.

Figure 6.

RMS error comparison of x-axis.

Figure 7.

RMS error comparison of y-axis.

Figure 8.

RMS error comparison of z-axis.

Hence, the contributions of the proposed DEIF-IMM algorithm fall into two categories. The first one is that a single model KF cannot exactly describe the motion of both TFP and CFP, whereas the proposed multiple model-based approaches work well. The second is that the DEIF-IMM takes advantage of auxiliary sensor data giving stability and performance enhancement. Also, Figure 9 illustrates histories of mode probabilities of the proposed algorithm described in equation (33). During 20–100 s, the proposed DEIF-IMM estimator could not find proper mode probability. The reason is that there is very small acceleration increase, especially for several seconds right after SLV’s lift-off. In this regime, due to small acceleration increase the early phase of first engine ignition pretends to be a CA motion without thrust. After this ambiguous regime, we can see that the mode probabilities are properly assigned according to the flight events. On the other hand, the proposed algorithm has its own drawback in filter robustness against system uncertainty. It originates from the independency of each multiple model. The proposed DEIF-IMM algorithm utilizes nonlinear multiple models and the main estimator, which is an IMM, determines its mode probability based on likelihood function, but it cannot be robust against sudden system or measurement uncertainties. For example, when system uncertainty occurs during CFP the uncertainty appears as a sudden acceleration (thrust). Then, likelihood function is biased to TFP, even though its status is CFP. Therefore, the proposed algorithm cannot be robust against both system and measurement uncertainties, whereas benchmark algorithms have proven to have robustness against system uncertainty. ²³ Also, Table 2 shows computational complexities of each algorithm regarding big O notation. Here n, m, and N indicate the state vector size, the measurement vector size, and the number of sensors, respectively. We can see that the proposed algorithm has reasonable complexity compared with benchmark algorithms for the simulation case $(n\ge \mathrm{Nm})$ .

Figure 9.

Histories of mode probabilities of the proposed DEIF-IMM algorithm.

Table 2.

Computational complexity $(n\ge \mathrm{Nm})$ .

Algorithm	$\mathcal{O}(·)$
CMKF-MF	$2n^3$
CMKF-SVF	$N(N^2+N+2)n^3$
DEIF-IMM	$3n^3$

CMKF-MF: converted measurement Kalman filter-measurement fusion; CMKF-SVF: converted measurement Kalman filter-state vector fusion; DEIF-IMM: distributed extended information filter-based interacting multiple model.

Conclusion

In this article, we investigated the DEIF-IMM with unbiased mixing for SLV tracking. The main idea is originated from the fact that there is no study concerning both multiplicity of dynamic model and heterogeneity of TS simultaneously in SLV tracking problem even though the multiple models and multiple sensor data fusion of the distributed TSs gives both effectiveness and robustness in SLV flight mission. By incorporating unbiased mixing ⁸ and information filter-based fusion ¹⁴ approaches to the distributed SLV tracking problem, we proposed a novel algorithm considering heterogeneity of tracking sensors and multiplicity of vehicle’s dynamic model. For verification of the proposed algorithms, tracking simulations using nominal trajectory of KSLV-1 including flight event information is used. The noisy measurements such as the radars and telemetry antenna pointing angles were generated from this nominal trajectory. To demonstrate the superiority of the proposed DEIF-IMM, in these simulations, we compared the proposed algorithm with classical multi-sensor fusion algorithms for RMS error between the nominal trajectory data and the estimated position. In the comparison of tracking performance between proposed DEIF-IMM and classical fusion algorithms, proposed DEIF-IMM algorithm shows the best performance.

Furthermore, we also noticed drawbacks of the proposed DEIF-IMM algorithm such as weakness in robustness against system uncertainty and necessity of sophisticated tuning procedures for multiple model specifications. These problems should be thoroughly concerned in the future work.