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.

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]. 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:

Measurements of the system are observed at discrete time instant t_k expressed as: 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: 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]: 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]: where, $p\left({\mathbf{y}}_k|\mathbf{x}\left(t_k\right)\right)$ is given by:

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]:

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.

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]: 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]:

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: 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.

In this paper we shall consider Bayes’ aposteriori PDF expressed in Equation (6) as GCSM up to third order (Equation (11)): 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).