Mixture Culver Filter for Space Object Orbit Determination

1. Introduction

State estimation of orbiting objects in space such as planets, asteroids, debris, or satellites can be realized by radar tracking or optical observations using batch or sequential filtering methods. These methods improve apriori orbit determination from a set of tracking data. Batch or least square estimators provide an object’s epoch state estimates by processing complete set of observations, while sequential estimators or filters process one measurement at a time to give state vector at the time of that measurement [1]. Optimal sequential estimation for linearized orbital dynamics in Minimum Mean Square Error (MMSE) sense is provided by Linearized Kalman Filter (KF) [1,2,3]. Proficient nonlinear filtering algorithms such as Extended Kalman Filter (EKF) and improved EKF (iEKF) were derived based on Gaussian assumption of Bayes’ posterior PDF and availability of rich measurement environment [1,2,3,4,5,6]. However, majority of the aerospace systems are nonlinear which consider non-Gaussian evolution of state, for example tracking of an Exo-atmospheric Re-Entry Vehicle (ERV), space object tracking, navigation of robots or aircrafts [7,8,9,10]. Aerospace engineers and scientists need to find algorithms for real-time sequential state estimation. Methods based on Gaussian posterior PDF may be endowed with suboptimal estimates for space object state estimates due to highly nonlinear nature of orbital and measurement dynamics with multiple modes or elongated tail PDFs [11]. The Gaussian Sum Filter (GSF) [12] tackles multiple mode distributions by assuming the Bayes’ posterior PDF as Gaussian Mixture Model (GMM) [13] and can be considered as parallel banks of EKFs. Improved propagation of state uncertainty using GMM were proposed by [11,14]. Nonlinear filters for space surveillance based on such models were presented by [15,16,17]. Better alternatives to EKF include Sigma Point Filters Family (SPFF) approximations of Bayesian posterior statistics [18,19,20]. In essence, these algorithms are also based on Gaussian assumption of the nonlinear system state evolution and commonly termed as Unscented Kalman Filter (UKF). Proficient improvements to UKF based on adaptive techniques to surmount dynamic and measurement model mismatch is presented by [21]. Tightly coupled Global Positioning System (GPS) Pseudo-Range / Inertial Measurement Unit (IMU) and Precise Point Position (PPP)/IMU navigation systems [22] employed in aerospace systems show nonlinearity during large IMU misalignments and GPS outages. To keep the advantages of the nonlinear filtering methods in dealing with such nonlinear systems, a Cubature Kalman Filter (CKF) + EKF hybrid filtering method based on dual estimation framework is proposed by [23]. CKF is a nonlinear filtering method based on the spherical-radial Cubature rule [24]. Being a deterministic sampling filtering method, CKF needs 2n (n=states/parameters of the system) Cubature points to propagate the state and covariance matrix, which shows a relatively smaller computational load than the UKF, as UKF mostly needs 2n+1 sigma points for the nonlinear states’ propagation [23]. Filtering methods based on Sequential Monte Carlo (SMC) methods known as Particle Filters (PF) have also been used for aerospace systems [7,25]. PF employ ensemble of weighted samples of state variables or parameters to solve online prediction and estimation requirements in a recursive manner. There are also efficient modifications to PFs presented in [26,27,28]. Gaussian density function expanded using Hermite polynomials [29] is termed as Gram Charlier Series (GCS) [30,31]. Hermite are orthogonal polynomials with Gaussian type weighting function over $-\infty$ to $\infty$ domain. Culver used third order GCS to derive analytic solutions for nonlinear Bayesian filtering by expanding nonlinear system equations up to second order in Taylor series [32]. As an extension to use of GCS approximation for Bayes’ posterior density for filtering nonlinear systems [32,33,34], a GCS Mixture (GCSM) PDF was also proposed [35,36]. In this paper filtering technique based on high fidelity GCSM model to capture evolution of nonlinear state uncertainty is proposed by adapting the technique used by [14] for GMM. This filtering method is termed as Mixture Culver Filter (MCF). In this paper comparison of GCS filter [32] (named Culver Filter (CF)), EKF [1,2] and GSF [12] with MCF shall be presented. To authors knowledge use of GCSM for nonlinear state estimation has not been reported in filtering literature. The new filter has shown improvement/comparable performance over other methods especially under space objects’ uncertain initial conditions and sparse measurement availability.

2. Continuous Discrete Nonlinear Filtering Problem

Large number systems of concern fall under the classification of nonlinear systems. Filtering algorithm that captures characteristics of the nonlinearities is preferable than considering the nonlinear problem to that of a linear one. Consider a continuous time dynamical expressed by the nonlinear ito Stochastic Differential Equation (SDE) of the form [37].

$$d\mathbf{x}\left(t\right)=\mathbf{f}\left(\mathbf{x}\left(t\right),t\right)dt+\mathbf{G}\left(\mathbf{x}\left(t\right),t\right)d\mathsf{\beta }\left(t\right)$$

(1)

where, $\mathbf{x}\left(t\right)\in {\mathbb{R}}^d$ is state of the d-dimensional dynamical system, $\mathbf{f}\left(\mathbf{x}\left(t\right),t\right)\in {\mathbb{R}}^{d\times 1}$ is nonlinear function which describes time evolution of the dynamical system, $\mathbf{G}\left(\mathbf{x}\left(t\right),t\right)\in {\mathbb{R}}^{d\times m}$ is dispersion matrix considered as function of $\mathbf{x}\left(t\right)$ and time “t”. $\mathsf{\beta }\left(t\right)\in {\mathbb{R}}^m$ is Brownian motion of zero mean value with diffusion $\mathbf{Q}\left(t\right)$ which can be expressed as:

$$E\left[d\mathit{\beta }\left(t\right){d\mathit{\beta }\left(t\right)}^T\right]=\mathit{Q}\left(t\right)dt$$

(2)

Measurements of the system are observed at discrete time instant t_k expressed as:

$${\mathbf{y}}_k=\mathbf{h}\left(\mathbf{x}\left(t_k\right)\right)+{\mathbf{v}}_k$$

(3)

where, ${\mathbf{y}}_k\in {\mathbb{R}}^q$ is q-dimensional discrete time observation vector, $\mathbf{h}\left(\mathbf{x}\left(t_k\right)\right)\in {\mathbb{R}}^{q\times 1}$ is discrete time nonlinear measurement equation considered as a function of $\mathbf{x}\left(t_k\right)$and ${\mathbf{v}}_k\in {\mathbb{R}}^q$ is zero-mean q-dimensional Gaussian process noise.

The measurement noise covariance is given by:

$$E\left[{\mathbf{v}}_k{\mathbf{v}}_j^{\mathrm{T}}\right]={\mathsf{\delta }}_{kj}{\mathbf{R}}_k$$

(4)

where, ${\mathsf{\delta }}_{kj}$ is dirac-delta function and ${\mathbf{R}}_k$, is covariance matrix.

The requirement is to obtain conditional state estimates on availability of the measurement ${\mathbf{y}}_k$. If apriori PDF of the dynamical system of Equation (1) is available, the predictive PDF $p\left(\mathbf{x}\left(t_k\right)|{\mathbf{y}}_{k-1}\right)$conditioned on previous observation satisfies FPKE [38]:

$$\frac{\partial p}{\partial t}=-\sum _{i=1}^d\frac{\partial }{\partial x_i}\left(\mathbf{f}\left(\mathbf{x}\left(t\right),t\right)p\right)+\frac{1}{2}\sum _{i=1}^d\sum _{j=1}^d\frac{{\partial }^2}{\partial x_i\partial x_j}\left[{\left(\mathbf{G}\mathbf{Q}{\mathbf{G}}^{\mathrm{T}}\right)}_{ij}p\right]$$

(5)

where, $p=p\left(\mathbf{x}\left(t_k\right)|{\mathbf{y}}_{k-1}\right)$is termed as state transition PDF, $\mathbf{G}\mathbf{Q}{\mathbf{G}}^{\mathrm{T}}$ is diffusion matrix of SDE and subscripts $i,j\in \left\{1,\dots ,d\right\}$ are dimension subscripts.

On receipt of measurement On receipt of measurement at discrete time $t_k$ the conditional PDF known as Bayes’ aposteriori PDF is expressed as [2,19]:

$$p\left(\mathbf{x}\left(t_k\right)|{\mathbf{y}}_k\right)= \frac{p\left({\mathbf{y}}_k|\mathbf{x}\left(t_k\right)\right)p\left(\mathbf{x}\left(t_k\right)|{\mathbf{y}}_{k-1}\right)}{{\int }_{-\infty }^{+\infty }p\left({\mathbf{y}}_k|\mathbf{x}\left(t_k\right)\right)p\left(\mathbf{x}\left(t_k\right)|{\mathbf{y}}_{k-1}\right)d\mathbf{x}\left(t_k\right)}$$

(6)

where, $p\left({\mathbf{y}}_k|\mathbf{x}\left(t_k\right)\right)$ is given by:

$$p\left({\mathbf{y}}_k|\mathbf{x}\left(t_k\right)\right)= \frac{1}{{\left|2\pi {\mathbf{R}}_k\right|}^{\frac{1}{2}}}\mathrm{e}\mathrm{x}\mathrm{p}\left\{-\frac{1}{2}{\left[{\mathbf{y}}_k-\mathbf{h}\left(\mathbf{x}\left(t_k\right)\right)\right]}^{\mathrm{T}}{\mathbf{R}}_k^{-1}\left[{\mathbf{y}}_k-\mathbf{h}\left(\mathbf{x}\left(t_k\right)\right)\right]\right\}$$

(7)

Equation: 5 and 6 are commonly termed as predictor (time evolution) and corrector (measurement conditioning) equations for progression of the dynamical system’s PDF. The optimal state of the system in MMSE sense is given by computing mean of Bayes’ aposteriori PDF [19]:

$${\widehat{\mathbf{x}}}_k={\int }_{-\infty }^{+\infty }\mathbf{x}\left(t_k\right)p\left(\mathbf{x}\left(t_k\right)|{\mathbf{y}}_k\right)d\mathbf{x}\left(t_k\right)$$

(8)

In general, analytical solution of Equation (5), is achievable for linear dynamical systems [2]. Typically numerical techniques are utilized to find solutions for nonlinear systems of lower dimensions $\left(\mathrm{d}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}\le 6\right)$ [39]. Therefore, sequential state estimation of space object using numerical solution of Partial Differential Equation (PDE) expressed in Equation (5) may not be considered optimal. Numerical solution of FPKE for such a need presents excessive memory requirement due to storage and recursively progress entire Bayes’ aposteriori PDF after each time step. Therefore, a requirement for computationally feasible solution to Bayes’ aposteriori PDF is indicative to build realistic ground-based space object Orbit Determination (OD) systems.

3. Gram Charlier Series Mixture Model

Orthogonal expansion of Gaussian PDF using Hermite polynomials and higher order moments can approximate probability distributions of a nonlinear dynamical system. Third order GCS can be expressed as [31]:

$$p_{gcs}\left({\mathbf{x}}_k\right)=\mathcal{N}\left({\mathbf{x}}_k,{\mathsf{\mu }}_k,{\mathbf{P}}_k\right)\left[1+\sum _{i,j,l}\frac{{\mathrm{P}}_{ijl}^{\left(3\right)}}{3!}{h}_{ijl}\left({\mathbf{x}}_k,{\mathsf{\mu }}_k,{\mathbf{P}}_k\right)\right]$$

(9)

where, k, is subscript for discrete time ($t_k$) for the state ${\mathbf{x}}_k=\mathbf{x}\left(t_k\right))$, $\mathcal{N}\left({\mathbf{x}}_k,{\mathsf{\mu }}_k,{\mathbf{P}}_k\right)$ is Gaussian PDF with mean, ${\mathsf{\mu }}_k$ and covariance, ${\mathbf{P}}_k$, $h_{ijl}\left({\mathbf{x}}_k,{\mathsf{\mu }}_k,{\mathbf{P}}_k\right)$ is multi-dimensional third order Hermite polynomials with dimensions $i,j,l\in \left\{1,\dots d\right\}$, and ${\mathrm{P}}_{ijl}^{\left(3\right)}$ is an element of multivariate Co-skewness tensor with dimensions $i,j,l$. Time subscript $\left(k\right)$ can be omitted for Co-skewness tensor notations by considering time dependence implicitly. Rodrigues formula can be used to generate Hermite polynomials by differentiating$\mathcal{N}\left({\mathbf{x}}_k,{\mathsf{\mu }}_k,{\mathbf{P}}_k\right)$[33,34]:

$$h_{i\dots l}\left({\mathbf{x}}_k,{\mathsf{\mu }}_k,{\mathbf{P}}_k\right)={\left(-1\right)}^{i+..+l}\frac{1}{\mathcal{N}\left({\mathbf{x}}_k,{\mathsf{\mu }}_k,{\mathbf{P}}_k\right)}\frac{d^{i+..+l}}{{dx}_1^i\dots {dx}_d^l}\mathcal{N}\left({\mathbf{x}}_k,{\mathsf{\mu }}_k,{\mathbf{P}}_k\right)$$

(10)

Higher order polynomials can be used to approximate PDF of nonlinear dynamical systems. Estimation algorithms basing on single GCS [32,33,34] are specified at a particular low order moment term. Negative probability regions could be formed for lower order GCS which may not always be a suitable PDF [36]. This may result a PDF not integrated to unity. PDFs with multi-modes or centroid of such PDFs may not be captured by Single GCS, especially GCS with lower orders$\left(\mathrm{o}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{r}\le 4\right)$. Generally, GCS of the higher orders$(\mathrm{o}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{r}\ge 15)$would be required to find good approximation to such cases [35]. Computational complexity increases for multivariate systems of higher order. An increase in order of GCS adds $\left(o+d-1\right)!/(o!\left(d-1\right)!)$ moment terms where, o = order and d = multivariate dimension of the PDF. Furthermore, a point may be reached where an increase in GCS orders would not improve the approximation any further [40]. Consequently, Van Hulle suggested GCS mixture models of lower orders$\left(3\le \mathrm{o}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{r}\le 5\right)$ to overcome complexity related with single higher order series [41]. Model based on GCS Mixture (GCSM) up to third order can be expressed as:

$$p_{gcsm}\left({\mathbf{x}}_k\right)=\sum _{g=1}^G{\alpha }_k^{\left(g\right)}\mathcal{N}\left({\mathbf{x}}_k,{\mathsf{\mu }}_k^{\left(g\right)},{\mathbf{P}}_k^{\left(g\right)}\right)\left\{1+\sum _{i,j,l}\frac{{\mathrm{P}}_{ijl}^{\left(3\right)\left(g\right)}}{3!}{h}_{ijl}\left({\mathbf{x}}_k,{\mathsf{\mu }}_k^{\left(g\right)},{\mathbf{P}}_k^{\left(g\right)}\right)\right\}$$

(11)

where,${\alpha }_k^{\left(g\right)}$ = Weight of $g^{\mathrm{t}\mathrm{h}}$ component of the GCSM model, G = Total components of GCSM.

Aforesaid, one may now consider third order GCSM model as a better approximation of Bayes’ aposteriori PDF needed for state estimation of space objects.

4. Description of Mixture Culver Filter

In this paper we shall consider Bayes’ aposteriori PDF expressed in Equation (6) as GCSM up to third order (Equation (11)):

$$p\left({\mathbf{x}}_k|{\mathbf{y}}_k\right)=\frac{1}{C_k}\sum _{g=1}^G{\alpha }_{k|k-1}^{\left(g\right)}\mathcal{N}\left({\mathbf{x}}_k;{\mathsf{\mu }}_{k|k-1}^{\left(g\right)},{\mathbf{P}}_{k|k-1}^{\left(g\right)}\right)\left[1+\sum _{i,j,l}\frac{{\mathrm{P}}_{ijl}^{\left(3\right)\left(g\right)}}{3!}{h}_{ijl}\left({\mathbf{x}}_k,{\mathsf{\mu }}_{k|k-1}^{\left(g\right)},{\mathbf{P}}_{k|k-1}^{\left(g\right)}\right)\right]p\left({\mathbf{y}}_k|{\mathbf{x}}_k\right)$$

(12)

where, $k|k-1$, is subscript for state predictive PDF indicating transition from k-1 to k^th instant of time, $p\left({\mathbf{y}}_k|{\mathbf{x}}_k\right)$, is PDF of measurement conditioned on evolved state and $C_k$, is normalizing constant (expressed later).

4.1. Time Update

Predictive PDF for a nonlinear continuous time dynamical system is obtained exactly by solving FPKE (Equation (5)) between the measurements from time $t_k$ to $t_{k-1}$. However, in GCSM approximation of predictive PDF, the parameters of each GCS component, mean${\mathsf{\mu }}_{k|k-1}$, covariance, ${\mathbf{P}}_{k|k-1}$ covariance and co-skewness ${\mathrm{P}}_{ijl}^{\left(3\right)}$ tensor components shall be obtained by numerically integrating following Equations [32]:

$$\frac{d{\mathsf{\mu }}_i\left(t\right)}{dt}={\mathrm{f}}_i\left(\mathsf{\mu },t\right)+{\mathrm{A}}_{ief}\left(\mathsf{\mu },t\right){\mathrm{P}}_{ef}$$

(13)

$$\frac{d{\mathrm{P}}_{ij}\left(t\right)}{dt}=2{\left\{{{\mathrm{F}}_{ie}\left(\mathsf{\mu },t\right){\mathrm{P}}_{je}+\mathrm{A}}_{ief}\left(\mathsf{\mu },t\right){\mathrm{P}}_{jef}^{\left(3\right)}\right\}}_{\mathrm{s}}+{\mathrm{V}}_{ij}$$

(14)

$$\frac{d{\mathrm{P}}_{ijl}^{\left(3\right)}\left(t\right)}{dt}=3{\left\{{{\mathrm{F}}_{ie}\left(\mathsf{\mu },t\right){\mathrm{P}}_{jle}^{\left(3\right)}+\mathrm{A}}_{ief}\left(\mathsf{\mu },t\right)({\mathrm{P}}_{je}{\mathrm{P}}_{lf}+{\mathrm{P}}_{jf}{\mathrm{P}}_{le})\right\}}_{\mathrm{s}}$$

(15)

The Equation (13) to Equation (15) are derived using Taylor series expansion of the nonlinear function $\mathbf{f}\left[\left(\mathbf{x}\left(t\right),t\right)\right]$ (Equation (1)) up to second order and employing ito differential rule [32,37]. Therefore, ${\mathrm{F}}_{ie}\left(\widehat{\mathbf{x}},t\right)={\partial \mathrm{f}}_i\left(\widehat{\mathbf{x}},t\right)/\partial {\mathrm{x}}_e$ is a component of Jacobian matrix, ${\mathrm{A}}_{ief}\left(\widehat{\mathbf{x}},t\right)=\left(\frac{1}{2}\right){{\partial }^2\mathrm{f}}_i\left(\widehat{\mathbf{x}},t\right)/\partial {\mathrm{x}}_e\partial {\mathrm{x}}_f$ is component of Hessian matrix with tensor subscripts notations $i,e,f\in \{1\dots d\}$. $N{\left\{.\right\}}_{\mathrm{s}}$ represents number of terms in the bracket by symmetrizing with respect to all subscripts. For example, symmetric terms are expressed as; $3{\left\{{\mathrm{P}}_{ij}{\mathrm{P}}_{lm}\right\}}_{\mathrm{s}}={\mathrm{P}}_{ij}{\mathrm{P}}_{lm}+{\mathrm{P}}_{il}{\mathrm{P}}_{jm}+{\mathrm{P}}_{im}{\mathrm{P}}_{jl}$. In order to compute time update of weights ${\alpha }_{k|k-1}$ of each GCSM component of Equation: 12 methodology suggested by [14] is selected. The idea for optimal weight updates for each component of GCSM is realized by minimizing error between FPKE equation (Equation (5)) and time derivative of GCSM PDF. A continuous time notation (Equation (11)) will be used for development of time update of weights. The error in FPKE and time derivative of GCSM is expressed as:

$$e\left(\mathbf{x},t\right)=\frac{\partial p_{gcsm}\left(\mathbf{x},t\right)}{\partial t}-{\mathcal{L}}_{FP}\left(p_{gcsm}\left(\mathbf{x},t\right)\right)$$

(16)

where,${\mathcal{L}}_{FP}(.)$ = Fokker-Planck operator, is described as:

$${\mathcal{L}}_{FP}\left(p_{gcsm}\left(\mathbf{x},t\right)\right)=\sum _{g=1}^G{\alpha }_t^{\left(g\right)}{\mathcal{L}}_{FP}\left(p_{gcs}^{\left(g\right)}\right)=1_{FP}^T{\mathit{\alpha }}_t$$

(17)

where, ${\mathit{\alpha }}_t\in {\mathbb{R}}^{d\times 1}$ is the vector of weights and the elements of $1_{FP}\in {\mathbb{R}}^{d\times 1}$ are given by application of Fokker-Planck operator on individual GCS components:

$${\mathcal{L}}_{FP}\left(p_{gcs}^{\left(g\right)}\right)=-\frac{{\partial }^T{p}_{gcs}^{\left(g\right)}}{\partial t}\mathit{f}\left(\mathit{x},t\right)-p_{gcs}^{\left(g\right)}Tr\left[\frac{\partial \mathit{f}\left(\mathit{x},t\right)}{\partial \mathit{x}}\right]+\frac{1}{2}\mathrm{T}\mathrm{r}\left[\mathbb{Q}\frac{{\partial }^2{p}_{gcs}^{\left(g\right)}}{\partial \mathbf{x}\partial {\mathbf{x}}^T}\right]$$

(18)

The first term on right of Equation: 16 is obtained by taking total derivative expressed as:

$$\frac{\partial p_{gcsm}\left(\mathbf{x},t\right)}{\partial t}=\sum _{g=1}^G\left({\dot{\alpha }}_t^{\left(g\right)}{p}_{gcs}^{\left(g\right)}+{\alpha }_t^{\left(g\right)}\frac{{\partial }^T{p}_{gcs}^{\left(g\right)}}{\partial {\mathsf{\mu }}_t^{\left(g\right)}}{\dot{\mathsf{\mu }}}_t^{\left(g\right)}+{\alpha }_t^{\left(g\right)}\mathrm{T}\mathrm{r}\left[\frac{\partial p_{gcs}^{\left(g\right)}}{\partial {\mathbf{P}}_t^{\left(g\right)}}{\dot{\mathbf{P}}}_t^{\left(g\right)}\right]+{\alpha }_t^{\left(g\right)}\frac{\partial p_{gcs}^{\left(g\right)}}{\partial P_{ijl}^{\left(3\right)\left(g\right)}}{\dot{P}}_{ijl}^{\left(3\right)\left(g\right)}\right)$$

(19)

where, Tr = trace and the last term in Equation: 19 implies summation of derivatives over all indices (i,j,l) obtained as:

$${\alpha }_t^{\left(g\right)}\sum _{i,j,l}\frac{\partial p_{gcs}^{\left(g\right)}}{\partial {\mathrm{P}}_{ijl}^{\left(3\right)\left(g\right)}}{\dot{\mathrm{P}}}_{ijl}^{\left(3\right)\left(g\right)}$$

(20)

The derivative of the moments ${\dot{\mathsf{\mu }}}_t^{\left(g\right)}$, ${\dot{\mathbf{P}}}_t^{\left(g\right)}$and ${\dot{\mathrm{P}}}_{ijl}^{\left(3\right)\left(g\right)}$ for each GCS component are obtained from Equation: 13-15 and the derivative of weights, ${\dot{\alpha }}_t^{\left(g\right)}$ is obtained by time discretization using the first forward difference:

$${\dot{\alpha }}_t^{\left(g\right)}=\frac{1}{\mathsf{\Delta }t}\left({\alpha }_{t^{\text{'}}}^{\left(g\right)}-{\alpha }_t^{\left(g\right)}\right)$$

(21)

where, $t^{\mathrm{\text{'}}}=t+\mathsf{\Delta }t$

Now by substituting Equation (21) into Equation (19) one may rewrite total time derivative of GCSM as:

$$\frac{\partial p_{gcsm}\left(\mathbf{x},t\right)}{\partial t}=\sum _{g=1}^G\frac{1}{\mathsf{\Delta }t}{p}_{gcs}^{\left(g\right)}{\alpha }_{t^{\text{'}}}^{\left(g\right)}+\sum _{g=1}^G\left(\frac{{\partial }^T{p}_{gcs}^{\left(g\right)}}{\partial {\mathsf{\mu }}_t^{\left(g\right)}}{\dot{\mathsf{\mu }}}_t^{\left(g\right)}+\mathrm{T}\mathrm{r}\left[\frac{\partial p_{gcs}^{\left(g\right)}}{\partial {\mathbf{P}}_t^{\left(g\right)}}{\dot{\mathbf{P}}}_t^{\left(g\right)}\right]+\frac{\partial p_{gcs}^{\left(g\right)}}{\partial {\mathrm{P}}_{ijl}^{\left(3\right)\left(g\right)}}{\dot{\mathrm{P}}}_{ijl}^{\left(3\right)\left(g\right)}-\frac{1}{\mathsf{\Delta }t}{p}_{gcs}^{\left(g\right)}\right){\alpha }_t^{\left(g\right)}$$

(22)

$$\frac{\partial p_{gcsm}\left(\mathbf{x},t\right)}{\partial t}=\frac{1}{\mathsf{\Delta }t}{\mathbf{p}}_{gcs}^T{\mathit{\alpha }}_{t^{\text{'}}}+{\mathbf{m}}_{\mathbf{T}\mathbf{D}}^T{\mathit{\alpha }}_t$$

(23)

where, ${\mathit{\alpha }}_{t^{\text{'}}}\in {\mathbb{R}}^{d\times 1}$ is the vector of new weights which are being found out, ${\mathbf{p}}_{gcs}\in {\mathbb{R}}^{d\times 1}$ is the vector of GCS components and the elements of ${\mathbf{m}}_{\mathbf{T}\mathbf{D}}\in {\mathbb{R}}^{d\times 1}$ are expressed as second of Equation (22) (underlined). Now by substituting Equation (17) and Equation (23) into Equation (16) one would get FPKE error as:

$$e\left(t,\mathbf{x}\right)=\frac{1}{\mathsf{\Delta }t}{\mathbf{p}}_{gcs}^T{\mathit{\alpha }}_{t^{\text{'}}}+{\left({\mathbf{m}}_{\mathbf{T}\mathbf{D}}-1_{FP}\right)}^T{\mathit{\alpha }}_t$$

(24)

Furthermore, analytical expressions for different derivatives used in Equation (22) can be expressed in component wise tensor notation as:

$$\frac{\partial p_{gcs}}{\partial {\mathsf{\mu }}_a}=\left[{\mathrm{P}}_{au}^{-1}\left({\mathrm{x}}_u-{\mathsf{\mu }}_u\right)+\frac{1}{3!}{\mathrm{P}}_{ijl}^{\left(3\right)}\left({\mathrm{P}}_{au}^{-1}\left({\mathrm{x}}_u-{\mathsf{\mu }}_u\right)h_{ijl}-{\mathrm{P}}_{ia}^{-1}{h}_j{h}_l-{\mathrm{P}}_{ja}^{-1}{h}_i{h}_l-{\mathrm{P}}_{la}^{-1}{h}_i{h}_j+{\mathrm{P}}_{ia}^{-1}{\mathrm{P}}_{jl}^{-1}+{\mathrm{P}}_{ja}^{-1}{\mathrm{P}}_{il}^{-1}+{\mathrm{P}}_{la}^{-1}{\mathrm{P}}_{ij}^{-1}\right)\right]p_g$$

(25)

$$\begin{array}{l}\frac{\partial p_{gcs}}{\partial {\mathrm{P}}_{mn}}=[\frac{1}{2}\left({\mathrm{P}}_{mr}^{-1}{\mathrm{P}}_{ns}^{-1}\left({\mathrm{x}}_r-{\mathsf{\mu }}_r\right)\left({\mathrm{x}}_s-{\mathsf{\mu }}_s\right)-{\mathrm{P}}_{mn}^{-1}\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{1}{3!}{\mathrm{P}}_{ijl}^{\left(3\right)}(\frac{1}{2}\left({\mathrm{P}}_{mr}^{-1}{\mathrm{P}}_{ns}^{-1}\left({\mathrm{x}}_r-{\mathsf{\mu }}_r\right)\left({\mathrm{x}}_s-{\mathsf{\mu }}_s\right)-{\mathrm{P}}_{mn}^{-1}\right)h_{ijl}-{\mathrm{P}}_{im}^{-1}{\mathrm{P}}_{nr}^{-1}{\mathrm{P}}_{js}^{-1}{\mathrm{P}}_{lv}^{-1}\left({\mathrm{x}}_r-{\mathsf{\mu }}_r\right)\left({\mathrm{x}}_s-{\mathsf{\mu }}_s\right)\left({\mathrm{x}}_v-{\mathsf{\mu }}_v\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-{\mathrm{P}}_{jm}^{-1}{\mathrm{P}}_{ns}^{-1}{\mathrm{P}}_{ir}^{-1}{\mathrm{P}}_{lv}^{-1}\left({\mathrm{x}}_r-{\mathsf{\mu }}_r\right)\left({\mathrm{x}}_s-{\mathsf{\mu }}_s\right)\left({\mathrm{x}}_v-{\mathsf{\mu }}_v\right)-{\mathrm{P}}_{lm}^{-1}{\mathrm{P}}_{nv}^{-1}{\mathrm{P}}_{ir}^{-1}{\mathrm{P}}_{js}^{-1}\left({\mathrm{x}}_r-{\mathsf{\mu }}_r\right)\left({\mathrm{x}}_s-{\mathsf{\mu }}_s\right)\left({\mathrm{x}}_v-{\mathsf{\mu }}_v\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\mathrm{P}}_{im}^{-1}{\mathrm{P}}_{jl}^{-1}{\mathrm{P}}_{nr}^{-1}\left({\mathrm{x}}_r-{\mathsf{\mu }}_r\right)+{\mathrm{P}}_{jm}^{-1}{\mathrm{P}}_{nl}^{-1}{\mathrm{P}}_{ir}^{-1}\left({\mathrm{x}}_r-{\mathsf{\mu }}_r\right)+{\mathrm{P}}_{il}^{-1}{\mathrm{P}}_{jm}^{-1}{\mathrm{P}}_{ns}^{-1}\left({\mathrm{x}}_s-{\mathsf{\mu }}_s\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\mathrm{P}}_{im}^{-1}{\mathrm{P}}_{nl}^{-1}{\mathrm{P}}_{js}^{-1}\left({\mathrm{x}}_s-{\mathsf{\mu }}_s\right)+{\mathrm{P}}_{ij}^{-1}{\mathrm{P}}_{lm}^{-1}{\mathrm{P}}_{nt}^{-1}\left({\mathrm{x}}_v-{\mathsf{\mu }}_v\right)+{\mathrm{P}}_{im}^{-1}{\mathrm{P}}_{nj}^{-1}{\mathrm{P}}_{lv}^{-1}\left({\mathrm{x}}_v-{\mathsf{\mu }}_v\right)\left)\right]p_g\end{array}$$

(26)

$$\begin{array}{l}\frac{\partial p_{gcs}}{\partial {\mathrm{x}}_a}=[{-\mathrm{P}}_{au}^{-1}\left({\mathrm{x}}_u-{\mathsf{\mu }}_u\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{1}{3!}{\mathrm{P}}_{ijl}^{\left(3\right)}(-{\mathrm{P}}_{au}^{-1}\left({\mathrm{x}}_u-{\mathsf{\mu }}_u\right)h_{ijl}+{\mathrm{P}}_{ia}^{-1}{\mathrm{P}}_{js}^{-1}{\mathrm{P}}_{lv}^{-1}\left({\mathrm{x}}_s-{\mathsf{\mu }}_s\right)\left({\mathrm{x}}_v-{\mathsf{\mu }}_v\right){+\mathrm{P}}_{ja}^{-1}{\mathrm{P}}_{ir}^{-1}{\mathrm{P}}_{lv}^{-1}\left({\mathrm{x}}_r-{\mathsf{\mu }}_r\right)\left({\mathrm{x}}_v-{\mathsf{\mu }}_v\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\mathrm{P}}_{la}^{-1}{\mathrm{P}}_{ir}^{-1}{\mathrm{P}}_{js}^{-1}\left({\mathrm{x}}_r-{\mathsf{\mu }}_r\right)\left({\mathrm{x}}_s-{\mathsf{\mu }}_s\right)-{\mathrm{P}}_{ia}^{-1}{\mathrm{P}}_{jl}^{-1}-{\mathrm{P}}_{ja}^{-1}{\mathrm{P}}_{il}^{-1}-{\mathrm{P}}_{la}^{-1}{\mathrm{P}}_{ij}^{-1}\left)\right]p_g\end{array}$$

(27)

$$\begin{array}{l}\frac{{\partial }^2{p}_{gcs}}{\partial {\mathrm{x}}_m\partial {\mathrm{x}}_n}=[\left({\mathrm{P}}_{mr}^{-1}{\mathrm{P}}_{ns}^{-1}\left({\mathrm{x}}_r-{\mathsf{\mu }}_r\right)\left({\mathrm{x}}_s-{\mathsf{\mu }}_s\right)-{\mathrm{P}}_{mn}^{-1}\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{1}{3!}{\mathrm{P}}_{ijl}^{\left(3\right)}\left(\frac{1}{2}({\mathrm{P}}_{mr}^{-1}{\mathrm{P}}_{ns}^{-1}\left({\mathrm{x}}_r-{\mathsf{\mu }}_r\right)\left({\mathrm{x}}_s-{\mathsf{\mu }}_s\right)-{\mathrm{P}}_{mn}^{-1}\right)h_{ijl}+{\mathrm{P}}_{im}^{-1}{\mathrm{P}}_{jn}^{-1}{\mathrm{P}}_{lv}^{-1}\left({\mathrm{x}}_v-{\mathsf{\mu }}_v\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\mathrm{P}}_{im}^{-1}{\mathrm{P}}_{ln}^{-1}{\mathrm{P}}_{js}^{-1}\left({\mathrm{x}}_s-{\mathsf{\mu }}_s\right)+{\mathrm{P}}_{jm}^{-1}{\mathrm{P}}_{in}^{-1}{\mathrm{P}}_{lv}^{-1}\left({\mathrm{x}}_v-{\mathsf{\mu }}_v\right)+{\mathrm{P}}_{jm}^{-1}{\mathrm{P}}_{ln}^{-1}{\mathrm{P}}_{ir}^{-1}\left({\mathrm{x}}_r-{\mathsf{\mu }}_r\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\mathrm{P}}_{lm}^{-1}{\mathrm{P}}_{in}^{-1}{\mathrm{P}}_{js}^{-1}\left({\mathrm{x}}_s-{\mathsf{\mu }}_s\right)+{\mathrm{P}}_{lm}^{-1}{\mathrm{P}}_{jn}^{-1}{\mathrm{P}}_{ir}^{-1}\left({\mathrm{x}}_r-{\mathsf{\mu }}_r|)\right)]p_g\end{array}$$

(28)

$$\frac{\partial p_{gcs}}{\partial {\mathrm{P}}_{ijl}^{\left(3\right)}}=\left[\frac{1}{3!}{h}_{ijl}\right]p_g$$

(29)

where, component wise tensor notation used in above expressions (Equation (25) to Equation (29)) implies summation of indices, $p_g$ is a Gaussian PDF and ${\mathsf{\mu }}_u$, ${\mathrm{P}}_{au}^{-1}\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{P}}_{ijl}^{\left(3\right)}$ including similar forms, indicates individual components of mean, covariance and Coskewness tensors respectively. By propagating the mean ${\mathsf{\mu }}_t^{\left(g\right)}$, covariance ${\mathbf{P}}_t^{\left(g\right)}$and ${\mathrm{P}}_{ijl}^{\left(3\right)\left(g\right)}$ of individual GCS component using Equation (13) to Equation (15), one seeks to obtain new weights by minimizing the error in FPKE over a selected volume of state space [14]:

$${\mathrm{m}\mathrm{i}\mathrm{n}}_{{\alpha }_{t^{\text{'}}}^{\left(g\right)}}\frac{1}{2}{\int }_V^{}{e}^2\left(t,\mathbf{x}\right)d\mathbf{x}$$

(30)

$$\mathrm{s}.\mathrm{t}\sum _{g=1}^G{\alpha }_{t^{\text{'}}}^{\left(g\right)}=1$$

$${\alpha }_{t^{\text{'}}}^{\left(g\right)}\ge 0, g=1\dots G$$

The aforementioned problem can be formulated as a quadratic programming problem [14,42]:

$$\begin{array}{c}{\mathrm{m}\mathrm{i}\mathrm{n}}_{{\alpha }_{t^{\text{'}}}^{\left(g\right)}}\frac{1}{2}{\mathit{\alpha }}_{t^{\text{'}}}^T{\mathbf{M}}_{\mathrm{c}}{\mathit{\alpha }}_t+{\mathit{\alpha }}_{t^{\text{'}}}^T{\mathbf{N}}_{\mathrm{c}}{\mathit{\alpha }}_t+{\left({\mathit{\alpha }}_{t^{\text{'}}}-{\mathit{\alpha }}_t\right)}^T\left({\mathit{\alpha }}_{t^{\text{'}}}-{\mathit{\alpha }}_t\right)\\ \mathrm{s}.\mathrm{t}{1}_{d\times 1}^T{\mathit{\alpha }}_{t^{\text{'}}}=1\\ {\mathit{\alpha }}_{t^{\text{'}}}\ge 0_{d\times 1}\end{array}$$

(31)

where, $1_{d\times 1}\in {\mathbb{R}}^{d\times 1}$ is a vector of ones, $0_{d\times 1}\in {\mathbb{R}}^{d\times 1}$ is a vector of zeros and the matrices ${\mathbf{M}}_{\mathrm{c}}\in {\mathbb{R}}^{d\times d}$ and ${\mathbf{N}}_{\mathrm{c}}\in {\mathbb{R}}^{d\times d}$ are given by:

$$\begin{array}{c}{\mathbf{M}}_{\mathrm{c}}=\frac{1}{∆t^2}{\int }_V^{}{\mathbf{p}}_{gcs}{\mathbf{p}}_{gcs}^{T}d\mathbf{x}\\ {\mathbf{N}}_{\mathrm{c}}=\frac{1}{∆t}{\int }_V^{}{\mathbf{p}}_{gcs}{\left({\mathbf{m}}_{\mathbf{T}\mathbf{D}}-1_{FP}\right)}^{T}d\mathbf{x}\end{array}$$

(32)

where, “V” is domain of the ${\mathbf{p}}_{gcs}$

Analytical solutions were found out for above integrals and presented in Appendix A [35]. To author’s knowledge the adaptation of FPKE error feedback methodology for GCSM based nonlinear filter using analytical or numerical methods is new and has not appeared anywhere in space objects’ state estimation and filtering literature.

4.2. Measurement Update

Using the time updated GCSM along with new weights for each GCS component ${\alpha }_{k|k-1}^{\left(g\right)}$, one now consider treatment of Bayes’ posterior PDF (Equation (12)) for MMSE solution for nonlinear filtering problem. Firstly, the normalization constant $C_k$ in Equation (12) can be obtained as:

$$C_k\triangleq p\left({\mathbf{y}}_k|{\mathbf{y}}_{k-1}\right)=\left({\int }_{-\mathbf{\infty }}^{+\mathbf{\infty }}\sum _{g=1}^G{\alpha }_{k|k-1}^{\left(g\right)}\mathcal{N}\left({\mathbf{x}}_k,{\mathsf{\mu }}_{k|k-1}^{\left(g\right)},{\mathbf{P}}_{k|k-1}^{\left(g\right)}\right)\left[1+\sum _{i,j,l}\frac{{\mathrm{P}}_{ijl}^{\left(3\right)\left(g\right)}}{3!}{h}_{ijl}\left({\mathbf{x}}_k,{\mathsf{\mu }}_{k|k-1}^{\left(g\right)},{\mathbf{P}}_{k|k-1}^{\left(g\right)}\right)\right]p\left({\mathbf{y}}_k|{\mathbf{x}}_k\right)d{\mathbf{x}}_k\right)$$

(33a)

where, each GCS component$"g"$inside integral of Equation (33a) can be written as:

$${\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\left\{\mathcal{N}\left({\mathbf{x}}_k,{\mathsf{\mu }}_{k|k-1}^{\left(g\right)},{\mathbf{P}}_{k|k-1}^{\left(g\right)}\right)\left[1+\sum _{i,j,l}\frac{{\mathrm{P}}_{ijl}^{\left(3\right)\left(g\right)}}{3!}{h}_{ijl}\left({\mathbf{x}}_k,{\mathsf{\mu }}_{k|k-1}^{\left(g\right)},{\mathbf{P}}_{k|k-1}^{\left(g\right)}\right)\right]p\left({\mathbf{y}}_k|{\mathbf{x}}_k\right)\right\}d{\mathbf{x}}_k$$

(33b)

One may approximate $p\left({\mathbf{y}}_k|{\mathbf{x}}_k\right)$ as a multidimensional Gaussian PDF. Therefore, by linearization of the measurement function $\mathbf{h}\left(.\right)$ (Equation (3)) utilizing first order Taylor series expansion around predicted estimates ${\mathsf{\mu }}_{k|k-1}^{\left(g\right)}$, $p\left({\mathbf{y}}_k|{\mathbf{x}}_k\right)$ can be formulated as [32]:

$$p\left({\mathbf{y}}_k|{\mathbf{x}}_k\right)=\frac{1}{\sqrt{\left|2\pi \mathbf{R}\right|}}\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{1}{2}{\left({\mathbf{y}}_k-\mathbf{h}\left({\mathsf{\mu }}_{k|k-1}^{\left(g\right)}\right)-{\mathbf{H}}_k\left({\mathbf{x}}_k-{\mathsf{\mu }}_{k|k-1}^{\left(g\right)}\right)\right)}^{\mathrm{T}}{\mathbf{R}}^{-1}\left({\mathbf{y}}_k-\mathbf{h}\left({\mathsf{\mu }}_{k|k-1}^{\left(g\right)}\right){-\mathbf{H}}_k\left({\mathbf{x}}_k-{\mathsf{\mu }}_{k|k-1}^{\left(g\right)}\right)\right)\right]$$

(34)

where, ${\mathbf{H}}_k={\left.(\partial \mathbf{h}/\partial {\mathbf{x}}_k)\right|}_{{\mathbf{x}}_k={\mathsf{\mu }}_{k|k-1}^{\left(g\right)}}$

By substituting Equation (34) in Equation (33b) and computing the integral for each GCS component can be expressed as [32]:

$${\mathrm{D}}_k^{\left(g\right)}=\frac{\exp \{-\frac{1}{2}{\left({\mathbf{y}}_k-\mathbf{h}\left({\mathsf{\mu }}_{k|k-1}^{\left(g\right)}\right)\right)}^{\mathrm{T}}}{\sqrt[]{\left|\mathbf{I}+{\mathbf{H}}_{k|k-1}^{\mathrm{T}\left(g\right)}{\mathbf{R}}_k^{-1}{\mathbf{H}}_{k|k-1}^{\left(g\right)}{\mathbf{P}}_{k|k-1}^{\left(g\right)}\right|}}\times \left({\mathbf{R}}_k^{-1}-{\mathbf{R}}_k^{-1}{\mathbf{H}}_{k|k-1}^{\left(g\right)}{\mathsf{\Omega }}_k^{\left(g\right)}{\mathbf{H}}_{k|k-1}^{\mathrm{T}\left(g\right)}{\mathbf{R}}_k^{-1}\right)\left({\mathbf{y}}_k-\mathbf{h}\left({\mathsf{\mu }}_{k|k-1}^{\left(g\right)}\right)\right)\}\left(1+{\mathcal{V}}_k^{\left(g\right)}\right)$$

(35)

where, I denote identity matrix and ${\mathcal{V}}_k^{\left(g\right)}$is expressed in Table 1. This would yield the denominator $C_k$ as:

$$C_k=\sum _{g=1}^G{\alpha }_{k|k-1}^{\left(g\right)}{\mathrm{D}}_k^{\left(g\right)}$$

(36)

Now Mean, Covariance and Co-skewness of Bayes’ aposteriori PDF (Equation (12)) can be calculated using following integrals:

$${\mathbf{x}}_k^{+}=\frac{1}{C_k}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}p\left({\mathbf{x}}_k|{\mathbf{y}}_k\right){\mathbf{x}}_{k}d{\mathbf{x}}_k$$

(37a)

$${\mathbf{P}}_k^{+}=\frac{1}{C_k}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}p\left({\mathbf{x}}_k|{\mathbf{y}}_k\right)\left({\mathbf{x}}_k-{\mathbf{x}}_k^{+}\right){\left({\mathbf{x}}_k-{\mathbf{x}}_k^{+}\right)}^{\mathrm{T}}d{\mathbf{x}}_k$$

(37b)

$${\mathbf{P}}_k^{+\left(3\right)}=\frac{1}{C_k}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}p\left({\mathbf{x}}_k|{\mathbf{y}}_k\right)\left({\mathbf{x}}_k-{\mathbf{x}}_k^{+}\right){\left({\mathbf{x}}_k-{\mathbf{x}}_k^{+}\right)}^{\mathrm{T}}\otimes {\left({\mathbf{x}}_k-{\mathbf{x}}_k^{+}\right)}^{\mathrm{T}}d{\mathbf{x}}_k$$

(37c)

Firstly, one compute means ${\mathbf{x}}_k^{+}$ by rewriting Equation (37a) as:

$${\mathbf{x}}_k^{+}=\frac{1}{\sum _{g=1}^G{\alpha }_{k|k-1}^{\left(g\right)}{\mathrm{D}}_k^{\left(g\right)}}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\sum _{g=1}^G{\alpha }_{k|k-1}^{\left(g\right)}\mathcal{N}\left({\mathbf{x}}_k;{\mathsf{\mu }}_{k|k-1}^{\left(g\right)},{\mathbf{P}}_{k|k-1}^{\left(g\right)}\right)\left[1+\sum _{i,j,l}\frac{{\mathrm{P}}_{ijl}^{\left(g\right)}}{3!}{h}_{ijl}\left({\mathbf{x}}_k,{\mathsf{\mu }}_{k|k-1}^{\left(g\right)},{\mathbf{P}}_{k|k-1}^{\left(g\right)}\right)\right]p\left({\mathbf{y}}_k|{\mathbf{x}}_k\right){\mathbf{x}}_{k}d{\mathbf{x}}_k$$

(38)

Each GCS component in Equation (38) can be solved as:

$${\mathbf{x}}_k^{+\left(g\right)}=\frac{1}{{\mathrm{D}}_k^{\left(g\right)}}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{p}_{gcs}^{\left(g\right)}\left({\mathbf{x}}_{k|k-1}|{\mathbf{y}}_{k-1}\right)p\left({\mathbf{y}}_k|{\mathbf{x}}_k\right){\mathbf{x}}_{k}d{\mathbf{x}}_k$$

(39)

where,

$$p_{gcs}^{\left(g\right)}\left({\mathbf{x}}_{k|k-1}|{\mathbf{y}}_{k-1}\right)=\mathcal{N}\left({\mathbf{x}}_k;{\mathsf{\mu }}_{k|k-1}^{\left(g\right)},{\mathbf{P}}_{k|k-1}^{\left(g\right)}\right)\left[1+\sum _{i,j,l}\frac{{\mathrm{P}}_{ijl}^{\left(g\right)}}{3!}{h}_{ijl}\left({\mathbf{x}}_k,{\mathsf{\mu }}_{k|k-1}^{\left(g\right)},{\mathbf{P}}_{k|k-1}^{\left(g\right)}\right)\right]$$

(40)

Similarly, one may now compute Covariance and Co-skewness for each GCS component using following equations:

$${\mathbf{P}}_k^{+\left(g\right)}=\frac{1}{{\mathrm{D}}_k^{\left(g\right)}}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{p}_{gcs}^{\left(g\right)}\left({\mathbf{x}}_k|{\mathbf{y}}_k\right)\left({\mathbf{x}}_k-{\mathbf{x}}_k^{+\left(g\right)}\right){\left({\mathbf{x}}_k-{\mathbf{x}}_k^{+\left(g\right)}\right)}^{\mathrm{T}}d{\mathbf{x}}_k$$

(41a)

$${\mathbf{P}}_k^{+\left(3\right)\left(g\right)}=\frac{1}{{\mathrm{D}}_k^{\left(g\right)}}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{p}_{gcs}^{\left(g\right)}\left({\mathbf{x}}_k|{\mathbf{y}}_k\right)\left({\mathbf{x}}_k-{\mathbf{x}}_k^{+\left(g\right)}\right){\left({\mathbf{x}}_k-{\mathbf{x}}_k^{+\left(g\right)}\right)}^{\mathrm{T}}\otimes {\left({\mathbf{x}}_k-{\mathbf{x}}_k^{+\left(g\right)}\right)}^{\mathrm{T}}d{\mathbf{x}}_k$$

(41b)

The solution of integrals in Equation (39)-(41) are given by measurement update equations expressed as [32]:

$${\widehat{\mathrm{x}}}_i^{+\left(g\right)}={\widehat{\mathrm{x}}}_i^{-\left(g\right)}+{\mathrm{d}}_i^{\left(g\right)}+{\varphi }_i^{\left(g\right)}$$

(42a)

$${\mathrm{P}}_{ij}^{+\left(g\right)}={\mathsf{\Omega }}_{ij}^{\left(g\right)}+{\phi }_{ij}^{\left(g\right)}-{\varphi }_i^{\left(g\right)}{\varphi }_j^{\left(g\right)}$$

(42b)

$${\mathrm{P}}_{ijl}^{+\left(3\right)\left(g\right)}={\mathfrak{N}}_{ijl}^{\left(g\right)}-{\varphi }_i^{\left(g\right)}{\phi }_{jl}^{\left(g\right)}-{\varphi }_j^{\left(g\right)}{\phi }_{il}^{\left(g\right)}-{\varphi }_l^{\left(g\right)}{\phi }_{ij}^{\left(g\right)}+2{\varphi }_i^{\left(g\right)}{\varphi }_j^{\left(g\right)}{\varphi }_l^{\left(g\right)}$$

(42c)

The weight estimates in the measurement update could be found using ${\mathrm{D}}_k^{\left(g\right)}$ (Equation (35)) of each GCS component of the Bayes’ aposteriori PDF:

$${\alpha }_k^{\left(g\right)}=\frac{{\alpha }_{k|k-1}^{\left(g\right)}{\mathrm{D}}_k^{\left(g\right)}}{\sum _{g=1}^G{\alpha }_{k|k-1}^{\left(g\right)}{\mathrm{D}}_k^{\left(g\right)}}$$

(43)

MCF can be initialized using Expectation Maximization algorithm (EM) [43,44] using Matlab’s ReBEL Toolkit [45]. To simplify computation initial PDF can be assumed as GMM. An addition of Residual Resampling (RR) step [19] can be done in MCF algorithm after the weight update of Equation: 43 in order to generate children of GCS components with relatively significant weights and removal of components with insignificant weights. Thus, the effective size of weight $G_E$could be expressed as [19]:

$$G_E=\frac{1}{\sum _{g=1}^G{\left({\alpha }_k^{\left(g\right)}\right)}^2}$$

(44)

If $G_E<G_T$ where $G_T$ is required (threshold) size of weights, we would perform RR step. Time update of weights for each GCS component in MCF described earlier could become quite extensive for higher dimensional systems. Therefore, alternatively one can simplify the algorithm by keeping the weights constant between the measurements. This is essentially the same methodology used in traditional Gaussian Sum Filter (GSF). Algorithm for MCF is shown in Table 1.

5. Tracking of Space Object Using RADAR Measurements

In this section algorithms discussed for MCF would now be implemented for tracking of space object such as space debris using ground-based radars. Generally, for space debris tracking, some orbital elements, and the approximate size of the debris (radar cross-section) are available. The radar beam is pointed to a pre-determined position in space and after detection the debris is tracked, and observation vectors are collected. These observations are used to compute orbital parameters and radar signature. This mode of observation is called ‘target directed’ and is used when the uncertainty in the knowledge of a debris’ orbit is high and precise information is required for collision-avoidance maneuvers for operational spacecraft and for reentry predictions for potentially dangerous objects [46]. Equations of motion for space debris true orbital trajectory in Low Earth Orbit (LEO) used in this paper are given as:

$$\dot{\mathit{r}}=\mathbf{v}$$

(45a)

$$\dot{\mathbf{v}}=-\frac{{\mu }_E}{r^3}\mathit{r}+{\mathit{a}}_G+{\mathit{a}}_D$$

(45b)

$${\mathit{a}}_D=\frac{1}{2}\rho \frac{C_{D}A}{\mathrm{m}}{\mathrm{v}}^2\frac{\mathbf{v}}{\left|\mathbf{v}\right|}$$

(45c)

where,

$\mathit{r}={[X,Y,Z]}^{\mathrm{T}}$, $\mathbf{v}={[\dot{X,}\dot{Y},\dot{Z}]}^{\mathrm{T}}$ are position and velocity of a space debris in ECI coordinates,$\mathrm{v}= \left|\mathbf{v}\right|$, $r=\left|\mathit{r}\right|$, ${\mathit{a}}_G$ = perturbation acceleration due to zonal gravitational harmonic including up to J₄, and ${\mathit{a}}_D$ = acceleration due to atmospheric drag, $C_D$ = drag coefficient,$\rho$= atmospheric density, A = cross-sectional area of space debris and m = mass of debris [1,47].

The true trajectory shall be measured by Radar. Measurement consists of range, azimuth and elevation angles in Topocentric coordinate system [1]. Radar site is selected as Eglin US Air Force Base (AFB). Normally distributed measurement errors with following variances are selected (adapted from reference [8]):

$$\begin{array}{c}{\sigma }_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}}=25\mathrm{m},{\sigma }_{\mathrm{a}\mathrm{z}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{t}\mathrm{h}}=0.015\mathrm{d}\mathrm{e}\mathrm{g}\\ {\sigma }_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}=0.015\mathrm{d}\mathrm{e}\mathrm{g}\end{array}$$

(46)

Initial conditions of space debris to generate true trajectory are considered in LEO vicinity of International Space Station (ISS). Correct state estimation of space debris at such heights is crucial for safe operations of ISS.

5.1. MCF Results and Comparison with Other Nonlinear Filters for Shorter Durations

Space debris apriori estimates could be highly uncertain especially in case of sparsely tracked debris. Therefore, filters performance with significant uncertainty in initial position variance of ${10}^7{\mathrm{m}}^2$ and initial velocity variance of $50000{\mathrm{m}}^2.{\mathrm{s}}^{-2}$ shall be observed. Firstly, one provides the measurement data availability of 1 Hz. Due to this high frequency of measurement availability the time update using FPKE error feedback is not used in our first simulations. The time history of Instantaneous Error (IE) in ECI coordinates given Figure 1 (a & b) shows that the estimates of CF, MCF, GSF and EKF are close to each other. The convergence to lower errors of MCF is comparatively better than other three filters. The MCF and GSF are both being propagated using two GCSM and GMM components respectively. The error plots of these figures are obtained after averaging 100 Monte Carlo runs for each filter.

One of the drawbacks of filters based on mixture PDFs is the suboptimal time update of mixand weights when there are fewer or no measurements. In this situation the weights would remain constant until a measurement is received. This could possibly produce inferior estimates for filters based on mixture models. One would now incorporate optimal time update of weights in MCF algorithm using FPKE error feedback to compare filters for 0.033 Hz measurement availability (see Figure 1 (c & d)) and consider filtering performance over period once no observation is available. The time history of IE is shown in Figure 1 (c & d). The figures clearly show efficiency of MCF over CF owing to use of optimal weight updates. The error curves for position and velocity are lower for MCF. Moreover, the Root Mean Square Error (RMSE) criteria (Table: 2) and convergence to lower errors shows improvement provided by MCF over other filtering methods. Moreover, the performance of CF is slightly better than EKF. In general, the filters based on mixture PDFs (GSF and MCF) show improvement over single approximation of Bayes’ posterior PDF. These error curves are obtained by averaging 50 simulations runs for each filter. This provides a reasonable confidence over these estimation results.

One can also define instantaneous RMS error for the filters expressed as:

$${\epsilon }_i\left(k\right)=\sqrt{\frac{1}{N}\sum _{j=1}^N{\left[x_{i,j}\left(k\right)-{\widehat{x}}_{i,j}\left(k\right)\right]}^2}$$

(47)

where, $x_{i,j}\left(k\right)$ is the true state, ${\widehat{x}}_{i,j}\left(k\right)$ is the estimated, i = i^th component of state, j = j^th simulation and N = total number of simulations.

5.2. MCF Results and Comparison with Other Nonlinear Filters for Longer Durations

Now one extends the filtering performance for later orbital period i.e., once there are no observations available. Figure 2 depicts time history of IE and instantaneous RMSE (Equation (47)) over 3 orbital periods. The measurements (0.033 Hz) are only available for 4 min once the satellite is in viewing position from the radar site. These simulations are obtained from processing 50 Monte Carlo runs for each filter. The performance of CF and MCF are very close when compared for IE criteria, however, one may observe distinct improvement in RMSE results by MCF over other filtering methods.

5.3. MCF for Onboard Satellite Navigation Using GPS

LEO small satellites are equipped with a GPS receiver such as SGR-05P of Surrey Satellite Technology (SSTL), whose position and velocity accuracy are 10 m and 0.15 m/s, respectively [48]. Due to limited power resource onboard small satellites, use of GPS for OD should be done carefully to conserve the power and enhance satellite’s mission lifetime. Now we compare MCF with CF for OD of an in-orbit satellite whose position and velocity are available after 95 min (~ 1 orbital period). The time history of IE and RMSE are shown in Figure 3.

Figure 3. Comparison of CF with MCF for position/velocity errors $\left(∆\right)$ in ECI coordinates (shown in log scale) when measurement availability is once per orbital period (a to h).

Figure 3. Comparison of CF with MCF for position/velocity errors $\left(∆\right)$ in ECI coordinates (shown in log scale) when measurement availability is once per orbital period (a to h).

Figure 3 shows improvements achieved using MCF over CF. Convergence to lower errors is obtained by MCF in merely over 5 orbital periods. Therefore, by employing MCF one could increase mission lifetime of a satellite along with better OD. The error curves have also been produced by averaging using 50 Monte Carlo runs.

6. Conclusion

In this paper we have developed an improved nonlinear filter based on GCSM model approximation of Bayes’ aposteriori PDF for space object state estimation using radar measurements. Comparative analysis of our filter with other methods for tracking of space debris shows improvement/comparable performance under highly uncertain initial conditions and availability of sparse and no measurement data. MCF can also be used in other astrophysical computations such as motion and position of celestial bodies and satellites.