# Square-Root Higher-Order Unscented Estimators for Robust Orbit Determination

Yang Yang Dr Y. Yang is with the School of Mechanical and Manufacturing
Engineering, University of New South Wales, Sydney , NSW 2052, Australia.
E-mail<EMAIL_ADDRESS>
###### Abstract

Orbit determination (OD) is a fundamental problem in space surveillance and
tracking, crucial for ensuring the safety of space assets. Real-world ground-
based optical tracking scenarios often involve challenges such as limited
measurement time, short visible arcs, and the presence of outliers, leading to
sparse and non-Gaussian observational data. Additionally, the highly
perturbative and nonlinear orbit dynamics of resident space objects (RSOs) in
low Earth orbit (LEO) add further complexity to the OD problem. This paper
introduces a novel variant of the higher-order unscented Kalman estimator
(HOUSE) called $w$-HOUSE, which employs a square-root formulation and
addresses the challenges posed by nonlinear and non-Gaussian OD problems. The
effectiveness of $w$-HOUSE was demonstrated through synthetic and real-world
measurements, specifically outlier-contaminated angle-only measurements
collected for the Sentinel 6A satellite flying in LEO. Comparative analyses
are conducted with the original HOUSE (referred to as $\delta$-HOUSE),
unscented Kalman filters (UKF), conjugate unscented transformation (CUT)
filters, and precise orbit determination solutions estimated via onboard
global navigation satellite systems measurements. The results reveal that the
proposed $w$-HOUSE filter exhibits greater robustness when dealing with
varying values of the dependent parameter compared to the original
$\delta$-HOUSE. Moreover, it surpasses all other filters in terms of
positioning accuracy, achieving three-dimensional root-mean-square errors of
less than $60\text{\,}\mathrm{m}$ in a three-day scenario. This research
suggests that the new $w$-HOUSE filter represents a viable alternative to UKF
and CUT filters, offering improved positioning performance in handling the
nonlinear and non-Gaussian OD problem associated with LEO RSOs.

##  I Introduction

The escalating presence of space debris in Earth’s orbit poses a significant
and ever-growing threat to the safety of space assets. In order to mitigate
the risk of collisions with these hazardous objects, it is imperative to
obtain highly accurate and reliable orbital information about them. Ground-
based optical tracking has emerged as a crucial component of space
surveillance, given the proliferation of resident space objects (RSOs).
However, the operation of traditional optical sensors is contingent on
specific visibility conditions, where RSOs are sunlit while the background is
dark. Consequently, these facilities primarily operate during nighttime hours.
Moreover, various practical challenges complicate the process of gathering
orbital data through ground-based optical tracking. These challenges include
constraints on the field of view of telescopes, the rapid motion of RSOs,
particularly those in low Earth orbit (LEO), and the limited scheduling
opportunities for tracking each object. As a result, the number of visible
arcs, or tracks of measurements, is often limited, and each track contains
only short intervals of effective measurements, often spanning just a few
minutes. The complexity of orbit determination (OD) is further amplified when
there is a substantial data gap, such as a day or more, between consecutive
tracks of measurements. Additionally, optical sensor measurements are
typically subject to various types of noises, including background noise, shot
noise, readout noise, dark current noise, and quantisation noise [1], which
can introduce non-Gaussian characteristics into angular measurements like
right ascension and declination. Furthermore, the imperfect nonlinear orbital
dynamics degrades the prediction accuracy and contributes to filter
divergence.

The combination of these challenges transforms the task of OD for LEO RSOs
into a formidable problem, characterised by nonlinearity and non-Gaussianity.
To address such complex estimation problems, various filtering techniques have
been developed over the years. The historical roots of online, real-time state
estimation based on dynamic models can be traced back to Rudolf E. Kalman’s
pioneering work in 1960 [2], where he introduced the Kalman filter (KF) for
linear stochastic systems with Gaussian assumptions. This enabled first
publicly known application in the Apollo lunar mission for the estimation of
the position of the spacecraft [3]. Since then, numerous KF variants have been
devised and widely employed in space navigation and OD problems, some capable
of handling nonlinear and non-Gaussian estimation scenarios.

Nonlinearities in the system dynamics can be handled by approximating them via
first-order Taylor expansions [4, 5], leading to the development of extended
Kalman filters (EKFs). However EKFs may fail when the simple linearisation
causes considerably inaccurate approximation. To address these limitations,
alternative approaches, such as Gaussian higher-order Kalman filters (GHKFs)
that take second- or even higher-order approximation of the nonlinearity, are
proposed [6]. Another set of alternatives are sigma point Kalman filters
(SPKF), including the unscented Kalman filter (UKF) [7], the central
difference filter [8], the divided difference filter [9], the cubature Kalman
filter [10], the conjugate unscented Kalman filter [11], etc. SPKFs
essentially are sampling-based filters and can yield more accurate results
than EKFs or GHKFs which are both derivative dependent, especially when
accurate initial condition states are poorly known [12]. The Gaussian sum
Kalman filter (GSKF) is also proposed to approximate a posterior non-Gaussian
probability density function (PDF) of the state using a weighted sum of
Gaussian PDFs [8, 13]. For the general OD problem, if measurements are subject
to gross errors or outliers, which is the case of optical tracking usually,
non-Gaussian behaviours due to both the orbital state and measurements have to
be handled in the sequential filtering framework. The characteristics of these
non-Gaussian PDFs can be accurately captured by introducing higher-order
moments such as skewness and kurtosis, as seen in the high-order unscented
estimator (HOUSE) [14]. Beyond addressing non-Gaussinity, the nonlinearity
inherent in orbit propagation can also be mitigated by using orbital elements,
such as the modified equinoctial orbital elements (MEEs) [15], instead of
Cartesian coordinates. This format has demonstrated the capability to reduce
the nonlinearity of orbital dynamics over more extended time periods [16, 17],
which, in turn, curtails the growth of non-Gaussianity due to the nonlinearity
of the orbit propagator.

This work aims to tackle the nonlinear and non-Gaussian OD problem by
introducing a novel version of HOUSE, dubbed $w$-HOUSE. $w$-HOUSE leverages a
square-root formulation for covariance evolution, eliminating the constraint
of high kurtosis values for the state to be estimated. To assess the efficacy
of $w$-HOUSE, comparisons are made with other filtering methods, including UKF
and square-root UKF (SRUKF), CUT-4 and CUT-6, and $\delta$-HOUSE. The
evaluation encompasses three numerical examples: a projectile example
inherited from previous work [14], a ground-based radar tracking scenario, and
a real-world ground-based optical tracking scenario.

This paper makes several significant contributions to the OD field:

  * •

It introduces a new $w$-HOUSE algorithm based on a square-root formulation,
removing the requirement for high kurtosis in the state being estimated
required by the orignial $\delta$-HOUSE. This leads to enhanced robustness
concerning the varying control parameter, i.e., $w$ for $w$-HOUSE compared to
$\delta$ for $\delta$-HOUSE.

  * •

It presents a ground-based angle-only OD algorithm that achieves state-of-the-
art performance in terms of accuracy and efficiency, even with outlier-
contaminated measurements. An orbit propagator based on MEEs is employed
across all filters. The proposed $w$-HOUSE outperforms the UKF, CUT-4, and
CUT-6 filters in terms of positioning errors, achieving better than
$60\text{\,}\mathrm{m}$ of three-dimensional (3D) root mean square errors
(RMSE) in a three-day scenario. Its runtime is comparable to UKFs and
significantly more efficient than CUT filters.

  * •

It provides an OD toolkit written in C++ and will be released to public along
the publication of this paper. This toolkit offers the community valuable
insights into advanced filtering techniques for the OD problem, a fundamental
challenge for numerous satellite missions and space situational awareness
applications.

The rest of the paper is organised as follows. Section II introduces the
fundamental Bayesian theorem as applied to the specific OD problem, where
orbit dynamics and measurement models are presented. Section III revisits the
sigma point sampling based linear least mean squares (LLMS) filters, including
UKF, CUT, and HOUSE, and particular emphasis is placed on the orignial
$\delta$-HOUSE. The novel $w$-HOUSE method is introduced by leveraging the
square-root formulation in Section IV. Section V presents numerical examples,
including a projectile example inherited from prior work, as well as two OD
scenarios featuring synthetic and real-world measurements. Results and
analyses of the OD outcomes are provided in Section VI, with concluding
remarks to follow in Section VII.

##  II Bayesian Filtering for Orbit Determination

###  II-A Systematic Equations with Initial Conditions

Consider the discrete-time state space models [18]:

$\displaystyle\vec{\boldsymbol{x}}(k)=\boldsymbol{F}(\vec{\boldsymbol{x}}(k-1),\vec{\boldsymbol{\omega}}(k-1)),$
(1)
$\displaystyle\vec{\boldsymbol{z}}(k)=\boldsymbol{h}(\vec{\boldsymbol{x}}(k),\vec{\boldsymbol{\nu}}(k)),$
(2)

where $\vec{\boldsymbol{x}}(k-1)$ indicates the orbital state at time
$t_{k-1}$; $\vec{\boldsymbol{z}}(k)$ is the measurement at time $t_{k}$;
$\vec{\boldsymbol{\omega}}$ and $\vec{\boldsymbol{\nu}}$ represent the noise
associated with the orbital state and measurement variables, respectively,
both of which follow a specific PDF such as a Gaussian PDF and a non-Gaussion
Pearson Type PDF [14]; $\boldsymbol{F}$ represents the nonlinear dynamics
driving the discrete-time evolution of the states; and $\boldsymbol{h}$ is the
nonlinear measurement function. The dimension of the state and measurement are
represented by $n_{x}$ and $n_{z}$, respectively. It is assumed that
$\vec{\boldsymbol{\omega}}$ and $\vec{\boldsymbol{\nu}}$ are independent of
$\vec{\boldsymbol{x}}$.

Given the known PDFs $p(\vec{\boldsymbol{\omega}})$,
$p(\vec{\boldsymbol{\nu}})$, $p_{0|0}(\vec{\boldsymbol{x}})$ for the
stochastic processes $\vec{\boldsymbol{\omega}}(k)$,
$\vec{\boldsymbol{\nu}}(k)$ and the initial state $\vec{\boldsymbol{x}}(0)$,
respectively, the objective of the sequential Bayesian filtering for OD is to
find the PDF $p_{k|k}(\vec{\boldsymbol{x}})$ from
$p_{k|k-1}(\vec{\boldsymbol{x}})$ given
$p_{k}(\vec{\boldsymbol{z}}|\vec{\boldsymbol{x}})$ in the following Eq. 9, so
that the estimate of the orbital state can be obtained. The operation of the
filter can be divided into two recursive steps: time update and measurement
update. During the time update step, the objective is to determine the PDF of
$\vec{\boldsymbol{x}}(k|k-1)$ given the PDFs of $\vec{\boldsymbol{x}}(k-1)$
and $\vec{\boldsymbol{\omega}}(k-1)$. During the measurement update step, the
objective is to determine the PDF of $\vec{\boldsymbol{x}}(k)$ given a
measurement $\vec{\boldsymbol{z}}(k)$, the PDF of $\vec{\boldsymbol{\nu}}(k)$,
and the prior predicted PDF of $\vec{\boldsymbol{x}}(k|k-1)$. This step
refines the state PDF by incorporating the measurement information and
adjusting the predicted PDF accordingly. By alternating between these two
steps, the filter recursively updates and refines the state estimate as new
measurements become available, improving the accuracy of the estimated state
over time.

###  II-B Time Update for $t\in[t_{k-1},t_{k}],\ k\geq 1$

Following the most general form by the Chapman-Kolmogorov equation, the time
update step involves propagating the state PDF forward in time, so the prior
PDF at $t_{k}$ is given by:

$\displaystyle p_{k|k-1}(\vec{\boldsymbol{x}})=\int
p_{k|k-1}(\vec{\boldsymbol{x}}|\vec{\boldsymbol{x}^{\prime}})p_{k-1}(\vec{\boldsymbol{x}^{\prime}})\text{d}\vec{\boldsymbol{x}^{\prime}},$
(3)

where $p_{k|k-1}(\vec{\boldsymbol{x}})$ is the predicted or prior PDF for the
oribal state $\vec{\boldsymbol{x}}(k)$, and
$p_{k|k-1}(\vec{\boldsymbol{x}}|\vec{\boldsymbol{x}^{\prime}})$ is the
conditional PDF for $\vec{\boldsymbol{x}}(k)$ given
$\vec{\boldsymbol{x}}(k-1)=\vec{\boldsymbol{x}^{\prime}}$. The PDF
$p_{k-1}(\vec{\boldsymbol{x}^{\prime}})$ is the prior PDF for
$\vec{\boldsymbol{x}}(k-1)$, which is obtained from the previous measurement
update step or the initial value for the first time epoch. The PDF
$p_{k}(\vec{\boldsymbol{x}})$ is then used in the measurement update step.
According to Eq. 1, the transition from $\vec{\boldsymbol{x}}(k-1)$ to
$\vec{\boldsymbol{x}}(k)$ is solely affected by the process noise
$\vec{\boldsymbol{\omega}}(k-1)$. Hence, the PDF
$p_{k|k-1}(\vec{\boldsymbol{x}}|\vec{\boldsymbol{x}^{\prime}})$ can be written
based on a Dirac delta function ($\delta_{\text{Dirac}}$):

$\displaystyle
p_{k|k-1}(\vec{\boldsymbol{x}}|\vec{\boldsymbol{x}^{\prime}})=\int\delta_{\text{Dirac}}(\vec{\boldsymbol{x}}-\boldsymbol{F}(\vec{\boldsymbol{x}^{\prime}},\vec{\boldsymbol{\omega}}))p_{k-1}(\vec{\boldsymbol{\omega}})\text{d}\vec{\boldsymbol{\omega}}.$
(4)

There are various ways of arriving at the discrete form of the orbit dynamics
equation above. For example, Eq. 1 can be derived from a continuous-time model
of orbit dynamics in the Cartesian coordinates:

$\displaystyle\frac{\text{d}}{\text{d}t}\vec{\boldsymbol{x}}(t)=\boldsymbol{f}(\vec{\boldsymbol{x}}(t),\vec{\boldsymbol{\omega}}(t)),$
(5)

where ${\vec{\boldsymbol{x}}}$ comprises the position ${\vec{\boldsymbol{r}}}$
and ${\vec{\boldsymbol{v}}}$ vectors. Usually these orbital states are
represented in the Earth centred inertial (ECI) system. For satellites in the
near earth orbit, gravitational and non-gravitational perturbations have to be
accounted for. In general the accelerations acting on satellites consist of
terms for the Earth gravitational force, the 3rd-body gravitational attraction
of the Sun and Moon, the solar radiation pressure and atmosphere drag on the
satellite, if no active orbital manoeuvre is performed. The exact formulations
for each term included in $\boldsymbol{f}(\vec{\boldsymbol{x}})$ can be
obtained from classical books, for example, [19] or author’s previous work
[20]. Given a sampling interval $T$, it is assumed the variation of the
unmeasured disturbances can be adequately approximated using piecewise
constant functions of the form
$\vec{\boldsymbol{\omega}}(t)=\vec{\boldsymbol{\omega}}(k),\ \text{for}\
t\in[t_{k-1},t_{k}]$. As a result, the discrete-time model in Eq. 1 can be
obtained by integrating Eq. 5 over the interval $[t_{k-1},t_{k}]$
($t_{k}=t_{k-1}+T$):

$\displaystyle\boldsymbol{F}(\vec{\boldsymbol{x}}(k-1),\vec{\boldsymbol{\omega}}(k-1))=\vec{\boldsymbol{x}}(k-1)+\int_{(k-1)T}^{kT}\boldsymbol{f}(\vec{\boldsymbol{x}}(\tau),\vec{\boldsymbol{\omega}}(k))\text{d}\tau.$
(6)

Instead of formulating the orbital dynamics in the Cartesian coordinate
system, orbital elements are used to describe the orbital dynamics and present
advantegous features in eliminating non-linearity. It is favourable for
handling orbital uncertainty propagation within a filtering framework. The MEE
is introduced addtionally to represent the motion of RSOs with respect to the
ECI frame [15, 21]. The state vector is
$\vec{\boldsymbol{x}}_{\text{mee}}=(p,f,g,h,k,L)^{T}$, where the first five
elements are real numbers while the last one is a circular variable. The MEE
is defined via the classical orbital elements (COE)
$\vec{\boldsymbol{x}}_{\text{coe}}=(a,e,i,\omega,\Omega,\nu)^{T}$:

$\displaystyle p$ $\displaystyle=a(1-e^{2}),$ (7) $\displaystyle f$
$\displaystyle=e\cos(\omega+\Omega),$ $\displaystyle g$
$\displaystyle=e\sin(\omega+\Omega),$ $\displaystyle h$
$\displaystyle=\tan(i/2)\cos(\Omega),$ $\displaystyle k$
$\displaystyle=\tan(i/2)sin(\Omega),$ $\displaystyle L$
$\displaystyle=\Omega+\omega+\nu,$

where $p$ is semiparameter, $a$ is semi-major axis, $e$ is orbital
eccentricity, $i$ is the orbital inclination, $\omega$ is argument of perigee,
$\Omega$ is right ascension of the ascending node, $\nu$ is true anomaly, and
$L$ is true longitude. There are no similarities with zero eccentricity and
inclination by using MEE, except when an orbit has an inclination of $\pm
180^{\circ}$. These extreme retrograde orbits are rare and can be ignored by
this work. As a result, the following ordinary differential equations are used
to describe the orbital dynamics:

$\displaystyle\dot{p}$ $\displaystyle=\frac{2p}{q}\sqrt{\frac{p}{\mu}}a_{t},$
(8) $\displaystyle\dot{f}$
$\displaystyle=\sqrt{\frac{p}{\mu}},\Bigg{\\{}\frac{[(1+q)\cos(L)+f]a_{t}-g(h\sin(L)-k\cos(L))a_{n}}{q}+a_{r}\sin(L)\Bigg{\\}},$
$\displaystyle\dot{g}$
$\displaystyle=\sqrt{\frac{p}{\mu}},\Bigg{\\{}\frac{[(1+q)\sin(L)+g]a_{t}+f(h\sin(L)-k\cos(L))a_{n}}{q}-a_{r}\cos(L)\Bigg{\\}},$
$\displaystyle\dot{h}$
$\displaystyle=\sqrt{\frac{p}{\mu}}\frac{s^{2}a_{n}}{2q}\cos(L),$
$\displaystyle\dot{k}$
$\displaystyle=\sqrt{\frac{p}{\mu}}\frac{s^{2}a_{n}}{2q}\sin(L),$
$\displaystyle\dot{L}$ $\displaystyle=\sqrt{\mu
p}\frac{z^{2}}{p^{2}}+\sqrt{\frac{p}{\mu}}\frac{(h\sin(L)-k\cos(L))a_{n}}{q},$

where $s^{2}=1+h^{2}+k^{2}$, $q=1+f\cos(L)+g\sin(L)$, and $a_{r}$, $a_{t}$ and
$a_{n}$ are non-two-body acceleration components due to perturbations in the
radial, tangential and normal directions, respectively. The radial direction
corresponds to the geocentric radius vector of the RSO, measured positively in
an outward direction from the Earth centre. The tangential direction is
perpendicular to this radius vector and measured positively in the direction
of orbital motion. The normal direction is positive along the angular momentum
vector of the RSO’s orbit.

###  II-C Measurement Update for $t=t_{k}$

The prior PDF is then updated via Bayesian theorem once a measurement is
available:

$\displaystyle
p_{k|k}(\vec{\boldsymbol{x}})=\frac{p_{k}(\vec{\boldsymbol{z}}|\vec{\boldsymbol{x}})p_{k|k-1}(\vec{\boldsymbol{x}})}{\int
p_{k}(\vec{\boldsymbol{z}}|\vec{\boldsymbol{x}})p_{k|k-1}(\vec{\boldsymbol{x}})},$
(9)

where $p_{k}(\vec{\boldsymbol{z}}|\vec{\boldsymbol{x}})$ is the measurement
likelihood funciton given by

$\displaystyle
p_{k}(\vec{\boldsymbol{z}}|\vec{\boldsymbol{x}})=\int\delta_{\text{Dirac}}(\vec{\boldsymbol{z}}-\boldsymbol{h}(\vec{\boldsymbol{x}}(k),\vec{\boldsymbol{\nu}}(k)))p_{\nu}\text{d}\vec{\boldsymbol{\nu}}.$
(10)

Normally there is no analytical solution to the integral in Eq. 9. Moments are
calulated to approximate the prior and posterior PDFs of the state, which
results in various filters.

Based on a traditional astronaumical telescope, measurements of topocentric
right ascension $\alpha_{\text{RA}}$ and declination $\delta_{\text{Dec}}$ can
be obtained. They are defined in the topocentric coordinate system and are
calculated from Cartesian coordinates of the RSO and tracking sensor as
follows:

$\displaystyle\alpha_{\text{RA}}=\tan^{-1}\Big{(}\frac{y_{\text{rso}}-y_{\text{stn}}}{x_{\text{rso}}-x_{\text{stn}}}\Big{)},$
(11)
$\displaystyle\delta_{\text{Dec}}=\sin^{-1}\Big{(}\frac{z_{\text{rso}}-z_{\text{stn}}}{\rho}\Big{)},$
(12)

where $\rho$ is the range between the RSO (indicated by the subscript ”rso”)
and the sensor (indicated by the subscript ”stn”). A Radar sensor can measure
the range and the range rate as additioanl two observables shown as below:

$\displaystyle\rho=\sqrt{(x_{\text{rso}}-x_{\text{stn}})^{2}+(y_{\text{rso}}-y_{\text{stn}})^{2}+(z_{\text{rso}}-z_{\text{stn}})^{2}},$
(13)
$\displaystyle\dot{\rho}=\frac{\vec{\boldsymbol{r}}_{\text{rso}}\cdot\vec{\boldsymbol{v}}_{\text{rso}}}{\rho},$
(14)

where $\vec{\boldsymbol{r}}_{\text{rso}}$ and
$\vec{\boldsymbol{v}}_{\text{rso}}$ are vectors of position and velocity,
respectively. When the orbit is propagated via MEE, a transformation is
necessary betweeen MEE to the Cartesian coordinate. Hence, the aforementioned
measurement equations can be rewritten w.r.t. MEE easily so are skipped here.

##  III Sigma Point Sampling Based Linear Least Mean Squares Filters

To obtain analytical solutions to the sequential Bayesian filtering for OD
introduced in Section II, one type of methods are to approximate the posterior
PDF given by Eq. 9 and apply one optimisation criteria to obtain the optimal
estimates from the posterior PDF [18]. Minimising the mean square error is a
widely used criteria which was first develped by Kalman [2] for linear
systems. This has been extended to nonlinear systems by using various ways,
for example, approximation based on Taylor series expansion and LLMS
approximation. LLMS approximation is a technique to approximate the
distribution of a nonlinear function of a random variable by linearising it
around its expected value [6, 22], which are shown as a better alternative for
approximating the nonlinear function than employing the Taylor expansion based
linearisation. LLMS estimators can work with Eqns. 1-2 directly and can also
be applied for non-Gaussian PDFs if these PDFs are known.

###  III-A Sigma Point Sampling Based LLMS Estimators

Sigma point sampling based LLMS estimators, also known as sigma point filters,
are a class of LLMS estimators that use sigma points to approximate moments of
a PDF, such as mean, variance, skewness, and kurtosis, which are essential for
characterising the shape and properties of the PDF. Different choices of sigma
points lead to different sigma point filters as introduced in Section I.
However, the general concept in these filters is sets of predicted state
values and measurement values are used to calculate the estimated covariance
matrices, which are then used to calculate the Kalman gain matrix in a linear
manner. Hence, the optimal state estimate is constructed recursively by
combining the predicted estimate with measurements. These steps are summarised
in the Appendix.

The UKF, proposed by Julier and Uhlman [23], uses an unscented transformation
(UT) to approximate the state PDF. It is based on the assumption that the
state follows the Gaussian distribution. The nonlinear UT employs sigma points
generated by using Eq. 15 to approximate the first two orders of moments,
i.e., mean and covariance, of the state:

$\displaystyle{\vec{\boldsymbol{x}}}^{j}=\mu+\sqrt{n+\kappa}\boldsymbol{c}^{j}\quad
w^{j}=\frac{1}{2(n+\kappa)},$
$\displaystyle{\vec{\boldsymbol{x}}}^{j+n}=\mu-\sqrt{n+\kappa}\boldsymbol{c}^{j}\quad
w^{j+n}=\frac{1}{2(n+\kappa)},$ (15)
$\displaystyle{\vec{\boldsymbol{x}}}^{2n+1}=\mu\quad
w^{2n+1}=\kappa/(n+\kappa),$

where $\mu$ is the mean of the $n$-dimensional state variable
${\vec{\boldsymbol{x}}}$, the $\boldsymbol{c}^{j}$ is the $j^{th}$ column of
original Cholesky factor of ${\boldsymbol{P}_{x}}$, $\kappa$ is a tuning
parameter, and $w$ is the corresponding weight. A tutorial of the original UKF
algorithm is summarised in [24], and its square-root form is presented in
[25].

Besides first and second orders of moments, Adurthi and et al. [11] proposed a
unified approach, labelled as CUT, to generate sigma points that approximate
higher-order moments for Gaussian and uniform PDFs. These filters are called
CUT filters. CUT introduces additional axes, labeled as conjugate axes, such
that the resultant points are more ”spread out” symmetrically in space. In
contrast, the UT selects sigma points that lie on carefully selected axes,
such as the principal axes or the eigenvectors of the covariance matrix. The
weights for the CUT and UT are calculated using different formulas, which are
derived based on the moment matching criteria. Therefore, the CUT and UT
differ in the selection of sigma points and corresponding weights, which can
affect the accuracy and efficiency of the approximation.

Still considering high-order moments, such as the state skewness and kurtosis
besides the mean and covariance, Stojanovski and Savransky [14] extended the
UKF based on the method by Ponomareva et al [26]. The generation of sigma
points is based on a more general form shown in Eq. 16:

$\displaystyle{\vec{\boldsymbol{x}}}^{j}=\mu+\vec{\boldsymbol{\alpha}}\circ\boldsymbol{c}^{j}\quad
w^{j}=\frac{1}{\vec{\boldsymbol{\alpha}}\cdot(\vec{\boldsymbol{\alpha}}+\vec{\boldsymbol{\beta}})n},$
$\displaystyle{\vec{\boldsymbol{x}}}^{j+n}=\mu-\vec{\boldsymbol{\beta}}\circ\boldsymbol{c}^{j}\quad
w^{j+n}=\frac{1}{\vec{\boldsymbol{\beta}}\cdot(\vec{\boldsymbol{\alpha}}+\vec{\boldsymbol{\beta}})n},$
(16) $\displaystyle{\vec{\boldsymbol{x}}}^{2n+1}=\mu\quad
w^{2n+1}=1-\sum_{j=1}^{2n}w^{j}.$

where $\vec{\boldsymbol{\alpha}}$ and $\vec{\boldsymbol{\beta}}$ are
coefficient vectors chosen to preserve the marginal third-and fourth-order
moments of $\vec{\boldsymbol{x}}$ and can be calculated via the following Eq.
48, and the operator $\circ$ indicates element-wise multiplication between two
vectors. The resulted filter is called the Higher-Order Unscented Estimator
[14]. It is designed to be more robust in cases where the initial
distribution, process noise and measurement noise have high kurtosis. It is
also indicated in [14] that HOUSE is computationally more efficient than the
fourth and higher-order CUT filters for three given simulation examples.

###  III-B High-Order Unscented Estimater

HOUSE accounts for non-Gaussian and nonuniform distributions by propagating
the state skewness and kurtosis in addition to the mean and covariance.
Skewness represents the asymmetry of the PDF, indicating whether the
probability mass is concentrated more to the left or right of the mean.
Kurtosis, on the other hand, measures the heaviness of the tails of the
distribution, highlighting the presence of extreme values, e.g., outliers
associated with measurements. Explicit solutions for sample points and weights
can be obtained via Eq. 16, making evaluations via HOUSE easier than CUT
filters. However, the original HOUSE filter is suffering from possible
negative value of the last weight given by Eq. 16. The solution proposed in
[14] is to introduce a $\delta$ parameter to intentionally increase the small
values of kurtosis. The working procedures of the original $\delta$-HOUSE are
introduced as below.

####  III-B1 Time Update

The time update process begins with the generation of sigma points. This is
achieved by considering the distribution of the augmented state vector, which
is represented by its mean:

$\displaystyle{\vec{\boldsymbol{x}}_{a}}(k-1)=\begin{bmatrix}{\vec{\boldsymbol{x}}}_{m}(k-1)\\\
{\vec{\boldsymbol{0}}_{\omega}}(k-1)\end{bmatrix},$ (17)

its covariace:

$\displaystyle{\boldsymbol{P}}_{a}(k-1)=\begin{bmatrix}{\boldsymbol{P}}_{x}(k-1)&\boldsymbol{0}\\\
\boldsymbol{0}&{\boldsymbol{P}}_{\omega}(k-1)\end{bmatrix},$ (18)

its skewness:

$\displaystyle{\vec{\boldsymbol{\gamma}}_{a}}(k-1)=\begin{bmatrix}{\vec{\boldsymbol{\gamma}}}_{x}(k-1)\\\
{\vec{\boldsymbol{\gamma}}}_{\omega}(k-1)\end{bmatrix},$ (19)

and its kurtosis:

$\displaystyle{\vec{\boldsymbol{\kappa}}_{a}}(k-1)=\begin{bmatrix}{\vec{\boldsymbol{\kappa}}}_{x}(k-1)\\\
{\vec{\boldsymbol{\kappa}}}_{\omega}(k-1)\end{bmatrix},$ (20)

where $k$ indicates the $k^{th}$ time epoch and $k\geq 1$, and
$\vec{\boldsymbol{0}}_{\omega}$ is an $n_{\omega}$-dimensional zero vector.
This results in $2(n_{x}+n_{\omega})+1$ sigma points, which are represented
as:

$\displaystyle{\vec{\boldsymbol{x}}}^{(j)}_{a}(k-1)=\begin{bmatrix}{\vec{\boldsymbol{x}}}^{(j)}(k-1)\\\
{\vec{\boldsymbol{\omega}}}^{(j)}(k-1)\end{bmatrix},$ (21)

where $1\leq j\leq 2(n_{x}+n_{\omega})+1$. The calculation of the vector of
weights ${\vec{\boldsymbol{w}}}$111Note the weight vector
$\vec{\boldsymbol{w}}$ is different from the process noise vector
${\vec{\boldsymbol{\omega}}}(k)$. for all sigma points can be referred to the
Section III.B of [14] and is also presented in the following Eq. 48. To
guarantee the nonnegativeness of the last weight $w^{2n+1}$, one way is to
increase the small values of kurtosis. For any component $\kappa_{i}$, if it
is smaller than $\frac{n}{1-\delta}+\gamma_{i}^{2}$, let

$\displaystyle\kappa_{i}=\frac{n}{1-\delta}+\gamma_{i}^{2},$ (22)

where $\delta$ is a predefined small positive number, $n$ is the dimention of
the random variable, i.e., the state ${\vec{\boldsymbol{x}}}$ or the process
noise ${\vec{\boldsymbol{\omega}}}$, and $\gamma_{i}$ is the corresponding
skewness component. As a result, the last weight $w^{2n+1}=\delta$ and its
positiveness is guaranteed.

The next phase involves the propagation of each sigma point using the orbital
dynamics models:

$\displaystyle{\vec{\boldsymbol{x}}}_{a}^{(j)}(k|k-1)=\boldsymbol{F}\big{(}\vec{\boldsymbol{x}}^{(j)}_{a}(k-1)\big{)}.$
(23)

Subsequently, the predicted mean and covariance are computed as follows:

$\displaystyle{\vec{\boldsymbol{x}}}_{a,m}(k|k-1)=\sum_{j=1}^{2(n_{x}+n_{\omega})+1}w^{(j)}\vec{\boldsymbol{x}}_{a}^{(j)}(k|k-1),$
(24)
$\displaystyle\vec{\delta\boldsymbol{x}}_{a}^{(j)}(k|k-1)=\vec{\boldsymbol{x}}_{a}^{(j)}(k|k-1)-{\vec{\boldsymbol{x}}}_{a,m}(k|k-1),$
(25)
$\displaystyle{\boldsymbol{P}}_{x}(k|k-1)=\sum_{j=1}^{2(n_{x}+n_{\omega})+1}w^{(j)}\vec{\delta\boldsymbol{x}}^{(j)}(k|k-1)\Big{[}\vec{\delta\boldsymbol{x}}^{(j)}(k|k-1)\Big{]}^{T},$
(26)

where $\vec{\delta\boldsymbol{x}}^{(j)}(k|k-1)$ is the first $n_{x}$ rows of
$\vec{\delta\boldsymbol{x}}_{a}^{(j)}(k|k-1)$.

Following this, the normalised states for the sigma points are computed:

$\displaystyle\vec{\tilde{\boldsymbol{x}}}^{(j)}(k|k-1)=\Big{(}\sqrt{{\boldsymbol{P}}_{x}(k|k-1)}\Big{)}^{-1}\vec{\delta\boldsymbol{x}}^{(j)}(k|k-1).$
(27)

Finally, the skewness and kurtosis of the normalised state are computed:

$\displaystyle\vec{\boldsymbol{\gamma}}_{x}(k|k-1)=\sum_{j=1}^{2(n_{x}+n_{\omega})+1}w^{(j)}\Big{[}\vec{\tilde{\boldsymbol{x}}}^{(j)}(k|k-1)\Big{]}^{3},$
(28)
$\displaystyle\vec{\boldsymbol{\kappa}}_{x}(k|k-1)=\sum_{j=1}^{2(n_{x}+n_{\omega})+1}w^{(j)}\Big{[}\vec{\tilde{\boldsymbol{x}}}^{(j)}(k|k-1)\Big{]}^{4}.$
(29)

####  III-B2 Measurement Update

The measurement update process begins by generating sigma points. This is
achieved by considering the distribution of the augmented state vector, which
is represented by its mean:

$\displaystyle{\vec{\boldsymbol{x}}}_{a}(k)=\begin{bmatrix}{\vec{\boldsymbol{x}}_{m}}(k|k-1)\\\
{\vec{\boldsymbol{0}}_{\nu}}(k)\end{bmatrix},$ (30)

its covariace:

$\displaystyle{{\boldsymbol{P}}_{a}}(k)=\begin{bmatrix}{\boldsymbol{P}}_{x}(k|k-1)&\boldsymbol{0}\\\
\boldsymbol{0}&{\boldsymbol{P}}_{\nu}(k)\end{bmatrix},$ (31)

its skewness:

$\displaystyle{\vec{\boldsymbol{\gamma}}_{a}}(k)=\begin{bmatrix}{\vec{\boldsymbol{\gamma}}}_{x}(k|k-1)\\\
{\vec{\boldsymbol{\gamma}}}_{\nu}(k)\end{bmatrix},$ (32)

and its kurtosis:

$\displaystyle{\vec{\boldsymbol{\kappa}}_{a}}(k)=\begin{bmatrix}{\vec{\boldsymbol{\kappa}}}_{x}(k|k-1)\\\
{\vec{\boldsymbol{\kappa}}}_{\nu}(k)\end{bmatrix},$ (33)

where ${\vec{\boldsymbol{x}}_{m}}(k|k-1)$, ${\boldsymbol{P}_{x}}(k|k-1)$,
${\vec{\boldsymbol{\gamma}}_{x}}(k|k-1)$ and
${\vec{\boldsymbol{\kappa}}_{x}}(k|k-1)$ can be extracted from the
corresponding augmented vectors and matrix in the time update, and
$\vec{\boldsymbol{0}}_{\nu}$ is an $n_{\nu}$-dimensional zero vector. This
results in $2(n_{x}+n_{\nu})+1$ sigma points, which are represented as:

$\displaystyle{\vec{\boldsymbol{x}}}^{(j)}_{a}(k)=\begin{bmatrix}{\vec{\boldsymbol{x}}}^{(j)}(k)\\\
{\vec{\boldsymbol{\nu}}}^{(j)}(k)\end{bmatrix},$ (34)

where $1\leq j\leq 2(n_{x}+n_{\nu})+1$. The next phase involves the
calculation of the corresponding measurement for each sigma point using the
measurement model:

$\displaystyle{\vec{\boldsymbol{z}}}^{(j)}(k)=\boldsymbol{h}\big{(}\vec{\boldsymbol{x}}^{(j)}_{a}(k)\big{)}.$
(35)

Subsequently, the mean of sigma points of the augmented state and the
deviations of sigma points from the mean are calculated as:

$\displaystyle{\vec{\boldsymbol{x}}}_{a,m}(k)=\sum_{j=1}^{2(n_{x}+n_{\nu})+1}w^{(j)}\vec{\boldsymbol{x}}_{a}^{(j)}(k),$
(36)
$\displaystyle\vec{\delta\boldsymbol{x}}_{a}^{(j)}(k)=\vec{\boldsymbol{x}}_{a}^{(j)}(k)-{\vec{\boldsymbol{x}}}_{a,m}(k).$
(37)

Following this, the measurement mean and covariance are computed:

$\displaystyle{\vec{\boldsymbol{z}}}_{m}(k)=\sum_{j=1}^{2(n_{x}+n_{\nu})+1}w^{(j)}\vec{\boldsymbol{z}}^{(j)}(k),$
(38)
$\displaystyle\vec{\delta\boldsymbol{z}}^{(j)}(k)=\vec{\boldsymbol{z}}^{j}(k)-{\vec{\boldsymbol{z}}}_{m}(k),$
(39)
$\displaystyle{\boldsymbol{P}}_{z}(k)=\sum_{j=1}^{2(n_{x}+n_{\nu})+1}w^{(j)}\vec{\delta\boldsymbol{z}}^{(j)}(k)\Big{[}\vec{\delta\boldsymbol{z}}^{(j)}(k)\Big{]}^{T},$
(40)
$\displaystyle{\boldsymbol{P}}_{xz}(k)=\sum_{j=1}^{2(n_{x}+n_{\nu})+1}w^{(j)}\vec{\delta\boldsymbol{x}}^{(j)}(k)\Big{[}\vec{\delta\boldsymbol{z}}^{(j)}(k)\Big{]}^{T},$
(41)

where $\vec{\delta\boldsymbol{x}}^{(j)}(k)$ is the first $n_{x}$ rows of
$\vec{\delta\boldsymbol{x}}_{a}^{(j)}(k)$. The Kalman gain is then calculated
as:

$\displaystyle{\boldsymbol{K}}(k)={\boldsymbol{P}}_{xz}(k){\boldsymbol{P}}_{z}^{-1}(k).$
(42)

The mean and covariance are updated using the following linear forms:

$\displaystyle{\vec{\boldsymbol{x}}}_{m}(k)={\vec{\boldsymbol{x}}}_{m}(k|k-1)+{\boldsymbol{K}}(k)(\vec{\boldsymbol{z}}(k)-{\vec{\boldsymbol{z}}}_{m}(k)),$
(43)
$\displaystyle{\boldsymbol{P}}_{x}(k)={\boldsymbol{P}}_{x}(k|k-1)-{\boldsymbol{K}}(k){\boldsymbol{P}}_{xz}^{T}(k).$
(44)

Based on the deviations in Eq. 37, normalised sigma points are calculated:

$\displaystyle\vec{\tilde{\boldsymbol{x}}}^{(j)}(k)=\Big{(}\sqrt{{\boldsymbol{P}}_{x}(k)}\Big{)}^{-1}\vec{\delta\boldsymbol{x}}^{(j)}(k).$
(45)

Accordingly, the skewness and kurtorsis vectors are computed as:

$\displaystyle\vec{\boldsymbol{\gamma}}_{x}(k)=\sum_{j=1}^{2(n_{x}+n_{\nu})+1}w^{(j)}\Big{[}\vec{\tilde{\boldsymbol{x}}}^{(j)}(k)\Big{]}^{3},$
(46)

and

$\displaystyle\vec{\boldsymbol{\kappa}}_{x}(k)=\sum_{j=1}^{2(n_{x}+n_{\nu})+1}w^{(j)}\Big{[}\vec{\tilde{\boldsymbol{x}}}^{(j)}(k)\Big{]}^{4}.$
(47)

The time update and measurement update are running recursively from the first
time epoch till the last one. The $\delta$-HOUSE algorithm is summarised in
Algorithm 1.

Data: Initial time epoch, initial state PDF represented by the mean
$\vec{\boldsymbol{x}}_{m}(0)$, covariance
$\textbf{cholesky}({\boldsymbol{P}}_{x}(0))$, skewness
$\vec{\boldsymbol{\gamma}}_{x}(0)$ and kurtoiss
$\vec{\boldsymbol{\kappa}}_{x}(0)$, the process noise PDF represented by the
covariance ${\boldsymbol{Q}}$, skewness $\vec{\boldsymbol{\gamma}}_{Q}$ and
kurtoiss $\vec{\boldsymbol{\kappa}}_{Q}$, the measurement noise PDF
represented by the covariance ${\boldsymbol{R}}$, skewness
$\vec{\boldsymbol{\gamma}}_{R}$ and kurtoiss $\vec{\boldsymbol{\kappa}}_{R}$,
and the pre-defined $\delta$ value.

Result: Updated state PDF for $k\geq 1$, represented by the mean
${\vec{\boldsymbol{x}}}_{m}(k)$, covariance ${\boldsymbol{P}}_{x}(k)$,
skewness $\vec{\boldsymbol{\gamma}}_{x}(k)$, and kurtosis
$\vec{\boldsymbol{\kappa}}_{x}(k)$.

Initialisation;

while  _Time epochs not finished_ do

Time Update:

  * •

Update the kurtosis component according to Eq. 22 if necessary before
generating

sigma points for the augmented state in Eq. 17 according to Eq. 48;

* •

Propagate sigma points through the nonlinear state transition function in Eq.
6; • Generate the predicted PDF of the state via Eqs. 17-29. Measurement
Update: • Update the kurtosis component according to Eq. 22 if necessary,
before generating

sigma points for the augmented state in Eq. 30 according to Eq. 48;

* •

Propagate sigma points through the nonlinear measurement function in Eq. 2; •
Generate the predicted PDF of the measurement via Eqs. 30-47.

end while

Algorithm 1 $\delta$-HOUSE Algorithm

##  IV New $w$-HOUSE in a Square-Root Form

One drawback of using the orignial $\delta$-HOUSE is that the ${(2n+1)}^{th}$
weight of sigma points may become negative according to Eq. 16, which will
damage the positive semi-definiteness of covariance update, leading to failure
of $\delta$-HOUSE. Meanwhile, it is proved the square-root formulation can
improve the numerical stability and guarantee the positive semi-definiteness
of the state covariances [25]. To mitigate the aforementioned issue, this work
develops a new variant of HOUSE, dubbed $w$-HOUSE, by utilising the recursive
update of square-root form of the covariance matrix rather than the covariance
matrix itself. A rank-one Cholesky update or downdate is performed
subsequently, depending on the sign of the ${(2n+1)}^{th}$ weight, which lifts
the constraint of having a nonnegative ${(2n+1)}^{th}$ weight so that avoids
distorted higher moments by defining minmum kurtorsis threshold values for the
orignial $\delta$-HOUSE in [14]. This threshold value of kurtorsis is still
useful when kurtorsis of state tends to grow along with time. Like the
original $\delta$-HOUSE, a new parameter $w$ is introduced to prevent the
${(2n+1)}^{th}$ weight of sigma points from growing wildly negative. More
specifically, if $w^{2n+1}$ exceeds the threshold value $w$, it will be reset
as zero and other $2n$ weights will also be recalculated. Revisit the
calculation of weights derived from the general form of Eq. 16:

$\displaystyle
w^{j}=\frac{1}{\vec{\boldsymbol{\alpha}}^{j}(\vec{\boldsymbol{\alpha}}^{j}+\vec{\boldsymbol{\beta}}^{j})},$
$\displaystyle
w^{j+n}=\frac{1}{\vec{\boldsymbol{\beta}}^{j}(\vec{\boldsymbol{\alpha}}^{j}+\vec{\boldsymbol{\beta}}^{j})},$
(48) $\displaystyle w^{2n+1}=1-\sum_{j=1}^{2n}w^{j},$

where

$\displaystyle\vec{\boldsymbol{\alpha}}^{j}=\frac{\vec{\boldsymbol{\gamma}}^{j}+\sqrt{4\vec{\boldsymbol{\kappa}}^{j}-3\times(\vec{\boldsymbol{\gamma}}^{j}\circ\vec{\boldsymbol{\gamma}}^{j})}}{2}$

and

$\displaystyle\vec{\boldsymbol{\beta}}^{j}=\frac{-\vec{\boldsymbol{\gamma}}^{j}+\sqrt{4\vec{\boldsymbol{\kappa}}^{j}-3\times(\vec{\boldsymbol{\gamma}}^{j}\circ\vec{\boldsymbol{\gamma}}^{j})}}{2}$

are dependent on the skewness $\vec{\boldsymbol{\gamma}}^{j}$ and kurtosis
$\vec{\boldsymbol{\kappa}}^{j}$ of the normalised state
$\vec{\tilde{\boldsymbol{x}}}^{j}$ in Eqs. 27 and 45 in time update and
measurement update, respectively. Note weights $w^{j}$ and $w^{n+j}$ ($1\leq
j\leq n$) are always positive while $w^{2n+1}$ is not. Instead, the proposed
$w$-HOUSE does not require $w^{2n+1}$ to be positive as always in the square-
root form by using the cholupdate algorithm introduced as below.

If ${\boldsymbol{S}}$ is the original Cholesky factor of ${\boldsymbol{P}}$,
then the Cholesky factor of the rank-one update or downdate
${\boldsymbol{P}}\pm\sqrt{\omega}{\boldsymbol{u}}{\boldsymbol{u}}^{T}$ can be
obtained by
${\boldsymbol{S}}=\textbf{cholupdate}{({\boldsymbol{S}},{\boldsymbol{u}},\pm\omega})$.
The operation will be executed for $m$ consecutive updates of Choleskey factor
if ${\boldsymbol{u}}$ is an $m-$column matrix [25]. An efficient
implementation of the cholupdate algorithm is given in [27]. More
specifically, instead of using Eq. 18 in the time update, its original
Cholesky factor is used throughout. It is then predicted via the QR
decomposition given the deviations of sigma points from the mean and
associated weights:

$\displaystyle{\boldsymbol{S}}_{a}(k|k-1)=\textbf{qr}\Bigg{(}\sum_{j=2}^{2(n_{x}+n_{\omega})+1}\sqrt{w^{(j)}}\vec{\delta\boldsymbol{x}}_{a}^{(j)}(k|k-1)\Bigg{)},$
(49)

followed by the rank-one update or downdate by the
${(2(n_{x}+n_{\omega})+1)}^{th}$ sigma point:

$\displaystyle{\boldsymbol{S}}_{a}(k|k-1)=\textbf{cholupdate}\Big{(}{\boldsymbol{S}}_{a}(k|k-1),\vec{\delta\boldsymbol{x}}_{a}^{(1)}(k|k-1),w^{({2(n_{x}+n_{\omega})+1})}\Big{)}.$
(50)

The normalised sigma points are calculated as:

$\displaystyle\vec{\tilde{\boldsymbol{x}}}^{(j)}(k|k-1)=[{\boldsymbol{S}}_{x}(k|k-1)]^{-1}\vec{\delta\boldsymbol{x}}^{(j)}(k|k-1),$
(51)

where $\vec{\delta\boldsymbol{x}}^{(j)}(k|k-1)$ is the first $n_{x}$ rows of
$\vec{\delta\boldsymbol{x}}_{a}^{(j)}(k|k-1)$ and
${\boldsymbol{S}}_{x}(k|k-1)$ is the first $n_{x}\times n_{x}$ block of
${\boldsymbol{S}}_{a}(k|k-1)$. Based on thee normalised sigma points, the
predicted skewness and kurtorsis of the state are calculated.

In the measurement update, instead of using the state covirance in Eq. 31, its
orginal Cholesky factor is used throughout. Likewise, the Choleskey factor of
the measurement covariance is calculated via the QR decomposition:

$\displaystyle{\boldsymbol{S}}_{z}(k)=\textbf{qr}\Bigg{(}\sum_{j=2}^{2(n_{x}+n_{\nu})+1}\sqrt{w^{j}}\vec{\delta\boldsymbol{z}}^{(j)}(k)\Bigg{)},$
(52)

followed by the rank-one update or downdate by the
${(2(n_{x}+n_{\nu})+1)}^{th}$ sigma point:

$\displaystyle{\boldsymbol{S}}_{z}(k)=\textbf{cholupdate}({\boldsymbol{S}}_{z}(k),\vec{\delta\boldsymbol{z}}^{(1)}(k),w^{(2(n_{x}+n_{\nu})+1)}).$
(53)

Then the Kalman gain is calculated as:

$\displaystyle\boldsymbol{K}={\boldsymbol{P}}_{xz}(k)/{\boldsymbol{S}}_{z}^{T}(k)/{\boldsymbol{S}}_{z}(k).$
(54)

The Choleskey factor of the covariance of the state is updated as:

$\displaystyle{\boldsymbol{S}}_{x}(k)=\textbf{cholupdate}({\boldsymbol{S}}_{x}(k|k-1),\boldsymbol{K}{\boldsymbol{S}}_{z}(k),-1).$
(55)

Accordingly, the normalised sigma points are calculated as:

$\displaystyle\vec{\hat{\boldsymbol{x}}}^{(j)}(k)=[{\boldsymbol{S}}_{x}(k)]^{-1}\vec{\delta\boldsymbol{x}}^{(j)}(k),$
(56)

followed by the updated skewness and kurtorsis of the state. Accordingly, the
new $w$-HOUSE algorithm is summarised in Algorithm 2.

Data: Initial time epoch, initial state PDF represented by the mean
$\vec{\boldsymbol{x}}_{m}(0)$, Cholesky factorisation of the covariance
${\boldsymbol{S}}_{x}(0)=\textbf{cholesky}({\boldsymbol{P}}_{x}(0))$, the
skewness $\vec{\boldsymbol{\gamma}}_{x}(0)$ and kurtoiss
$\vec{\boldsymbol{\kappa}}_{x}(0)$, the process noise PDF represented by the
Cholesky factorisations of the covariance
${\boldsymbol{S}}_{Q}=\textbf{cholesky}({\boldsymbol{Q}})$, skewness
$\vec{\boldsymbol{\gamma}}_{Q}$ and kurtoiss $\vec{\boldsymbol{\kappa}}_{Q}$,
the measurement noise PDF represented by the Cholesky factorisations of the
covariance ${\boldsymbol{S}}_{R}=\textbf{cholesky}({\boldsymbol{R}})$,
skewness $\vec{\boldsymbol{\gamma}}_{R}$ and kurtoiss
$\vec{\boldsymbol{\kappa}}_{R}$, and the pre-defined $w$ value.

Result: Updated state PDF for $k\geq 1$, represented by the mean
${\vec{\boldsymbol{x}}}_{m}(k)$, Cholesky factorisation of the covariance
${\boldsymbol{P}}_{x}(k)$, skewness $\vec{\boldsymbol{\gamma}}_{x}(k)$, and
kurtosis $\vec{\boldsymbol{\kappa}}_{x}(k)$.

Initialisation;

while  _Time epochs not finished_ do

Time Update:

  * •

Update the kurtosis component according to Eq. 22 if necessary before
generating

sigma points for the augmented state in Eq. 17 according to Eq. 48;

* •

Propagate sigma points through the nonlinear state transition function in Eq.
6; • Generate the predicted PDF similar to $\delta$-HOUSE. Eqs. 17, 19-25,
49-51, and 28-29. Measurement Update: • Update the kurtosis component
according to Eq. 22 if necessary before generating

sigma points for the augmented state in Eq. 30 according to Eq. 48;

* •

Propagate sigma points through the nonlinear measurement function in Eq. 2; •
Generate the updated PDF of the state via Eqs. 30, 32-39, 52-53, 41, 54, 43,
55-56,

and 46-47.

end while

Algorithm 2 $w$-HOUSE Algorithm

##  V Numerical Simulations and Experiements

To demosntrate the effectiveness of the proposed $w$-HOUSE, numerical
simulations and experiments were conducted in this section. The first example
was inherited from the projectile example presented in [14]. The second
example was a synthetic OD scenario with radar measurements. The third example
was a real-world OD scenario with optical angles.

###  V-A Projectile Example

As the first example, this work implemented the proposed $w$-HOUSE into the
Projectile Example presented in [14]. The equations of motion are expressed
as:

$\displaystyle\dot{x}$ $\displaystyle=v_{x},$ (57) $\displaystyle\dot{y}$
$\displaystyle=v_{y},$ $\displaystyle\dot{z}$ $\displaystyle=v_{z},$
$\displaystyle\dot{v}_{x}$ $\displaystyle=-bv\dot{v}_{x}+f_{x},$
$\displaystyle\dot{v}_{y}$ $\displaystyle=-bv\dot{v}_{y}+f_{y},$
$\displaystyle\dot{v}_{z}$ $\displaystyle=-bv\dot{v}_{z}+f_{z}-g,$

where $b$ is a paramter to calculate the atmospehric drag, which is set as a
constant of $0.001\text{\,}{\mathrm{m}}^{-1}$,
$g=$9.807\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-2}$$ is the Earth
gravitational acceleration, $v$ is the speed, and $f_{x}$, $f_{y}$ and $f_{z}$
are other specific accelerations. These specific accelerations are taken as
independent variables with mean zeros and standard deviations of
$0.01\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-2}$. The initial state is
$[1000,1000,0,500,0,500]^{T}$ with the unit of $\mathrm{m}$ for position
components and $\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ for velocity components.
The associated standard deviation values are $250\text{\,}\mathrm{m}$ for
position components and $100\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ for
velocity components. Azimuth and elevation angles were simulated at
$5\text{\,}\mathrm{Hz}$ as measurements by the following formulas:

$\displaystyle\alpha_{\text{Az}}=\arctan(y,-x)+\nu_{\alpha},$ (58)
$\displaystyle\epsilon_{\text{El}}=\arctan(z,\sqrt{x^{2}+y^{2}})+\nu_{\epsilon},$

where $\nu_{\alpha}$ and $\nu_{\epsilon}$ indicate additive measurement
noises, both with mean zeros and standard deviations of
$1\text{\,}\mathrm{\SIUnitSymbolArcminute}$. Gaussian and non-Gaussian Pearson
type IV distributions were considered in two cases, respectively. For the
latter case, all distributions have a kurtosis of 30.0; the initial state and
process noise have a skewness of 1.0 and the measurement noise has a skewness
of -1.0. For $\delta$-HOUSE, $\delta$ was set to 0. For $w$-HOUSE, $w$ was set
to -0.1.

###  V-B Synthetic OD Example with Radar Measurements

A numberical example similiar to the one presented in [17] was simulated, in
which a RSO in LEO is observed by a ground station at the Earth’s north pole.
The only difference is the measurements were generated according to the
elevation threshlod instead of randomly choosing the observational time
window. The initial orbital state at
$0\text{\,}\mathrm{h}$$0\text{\,}\mathrm{min}$$0\text{\,}\mathrm{s}$ of
January 4, 2010 is given in the ECI coordinate system as below:

$\displaystyle\vec{\boldsymbol{x}}_{0}=[7007.2175\times
10^{3},0,0,0,0.6606\times 10^{3},7.5509\times 10^{3}],$

with the units of $\mathrm{m}$ and $\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ for
position and velocity components, respectively. A non-Gaussian Pearson Type IV
distribution is considered for the initial state. Its covariance matrix is
given as:

$\centering\mathbf{P}_{0}=\begin{bmatrix}1.481\times 10^{8}&0&0&0&-9.237\times
10^{4}&-5.333\times 10^{4}\\\ 0&2.885\times 10^{7}&9.994\times
10^{6}&-3.121\times 10^{4}&0&0\\\ 0&9.994\times 10^{6}&5.770\times
10^{6}&-1.242\times 10^{4}&0&0\\\ 0&-3.121\times 10^{4}&-1.242\times
10^{4}&3.687\times 10^{1}&0&0\\\ -9.237\times 10^{4}&0&0&0&6.798\times
10^{1}&3.145\times 10^{1}\\\ -5.333\times 10^{4}&0&0&0&3.145\times
10^{1}&3.166\times 10^{1}\end{bmatrix}.\@add@centering$

The initial skewness and kurtosis are set as -1.6 and 10.0, respectively, for
all state variables.

As regards the orbital dynamics models for the ground truth, a $40\times 40$
Earth gravitational potential field was considered using the GGM03S model
[28], together with third-body attractions by Sun and Moon. No surface forces
were considered. These orbital dynamics models were used to generate synthetic
measurements at an interval of $30\text{\,}\mathrm{s}$. The ground truth of
the trajectory is shown in Fig. 1.

Figure 1: Scatter Plot of the 3D Reference Orbit for the Synthetic OD Example

A radar sensor was considered in this scenario so a set of range, range rate
and two angles were simulated when the satellite is above the horizon. A
Pearson Type IV distribution was considered, with standard deviations of
$100\text{\,}\mathrm{m}$, $0.3\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$
and $100\text{\,}\mathrm{\SIUnitSymbolArcsecond}$, and skewness of -1.0 and
kurtosis of 30.0 for all measurement variables. The uncorrupted measurements
are shown in Fig. 2.

Figure 2: Uncorrupted Measurements for the Synthetic OD Example

The filter only used a $10\times 10$ Earth gravity model for orbit prediction.
A Gaussian process noise was considered and its covariance is given as a
diagnal matrix with the standard deviations of $3.162\text{\,}\mathrm{m}$ for
all position components and
$0.01\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ for all velocity
components. For $\delta$-HOUSE, $\delta$ was set to 0. For $w$-HOUSE, $w$ was
set to -0.1.

###  V-C Real-World OD Example with Optical Angles

The Chinese Academy of Sciences has developed and installed the a couple of
electro-optical arrays at Changchun Observatory that monitor RSOs cross
different orbital domains from ground, ranging from LEO to higher orbits.
Their tracking facilities are operating in different modes, either to track a
mount of specific targets or to survey a complete region [29]. In general, the
all-sky survey yields short sequences of angular measurements, covering only a
small portion of the whole orbit. Hence, usually it is impossible to retrieve
the accurate and reliable orbital state information from a single track. After
track-to-track association for different time epochs [29], multiple tracks of
right ascension and declination were obtained for the Sentinel 6A satellite in
LEO and tested by six filters in the OD process.

The precise orbit determination (POD) solutions of this satetellite are
accessible via the European Space Agency Copernicus Open Access Hub [30]. They
were solved via processing measurements from onboard global navigation
satellite systems and have reached decimetre-level of accuracy for the three
dimensional (3D) position. These POD solutions were interpolated to the
optical measurement epochs and used as a benchmark for all filters to compare
to.

The MEE propgator aligned with the intervals of two adjcent measurements for
state prediction. The satellite mass is $961.831\text{\,}\mathrm{kg}$. Values
of 1.23 and 2.0 were used for the coefficients of solar radiation pressure
(SRP) and atmospheric drag calculations, respectively. They are the two main
surface forces acting on the LEO satellites. The effective areas are
$8.73\text{\,}{\mathrm{m}}^{2}$ for SRP and $4.09\text{\,}{\mathrm{m}}^{2}$
for drag. As for field forces, the $100\times 100$ GRACE Gravity Model GGM03S
was considered, together with the third-body attractions from the Sun and Moon
whose locations are obtained via the Jet Propulsion Laboratory ephemeris DE440
[31]. The solid Earth tide, ocean tide loading and relativistic effect were
also considered for as field forces [20]. The Earth rotational parameters were
taken from the Center for Orbit Determination in Europe [32] to facilitate the
transformations between the ECI frame and the Earth centred Earth fixed (ECEF)
frame. It is indciated one-day orbit propagation via high-fidelity models can
reach less than $20\text{\,}\mathrm{m}$ in 3D position error and less than
$5\text{\,}\mathrm{mm}\text{\,}{\mathrm{s}}^{-1}$ in 3D velocity error for the
Sentinel 6A satellite.

During the measurement outage, the maximum time step was fixed at
$180\text{\,}\mathrm{s}$ for the propagator to work with. Other time steps can
also be chosen, but the process noise covariance matrix need to be adjusted
accordingly. The initial condition for the satellite orbit can be generated
via the method in [29] or using the two-line element (TLE). In this work, a
set of TLE on the first day of the tracking window has been downloaded from
www.spacetrack.org and propagated to the first observational epoch. The
generated initial orbital state in the ECI frame is:

$\displaystyle\vec{\boldsymbol{x}}_{0}=[-3669575.624,1040418.514,6705990.196,-5051.328,-4686.718,-2035.774]^{T},$
(59)

where the unit is $\mathrm{m}$ for position components and
$\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ for velocity components. Without
knowing the exact statistical information of the TLE data, the initial
covariance matrix is assumed as a diagnoal matrix, given by:

$\displaystyle\mathbf{P}_{0}=\text{diag}(1.0\times 10^{6},1.0\times
10^{4},4.0\times 10^{4},0.01,0.25,0.04),$ (60)

with the unit of $$\mathrm{m}$^{2}$ for position components and
${\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-2}$ for velocity components. The
oscculating MEE for the initial state is obtained by converting the Cartesian
coordinate in Eq. 59 into oscculating COE first, followed by the conversion of
COE to MEE via Eq. 7:

$\displaystyle\vec{\boldsymbol{x}}_{\text{MEE},0}=\begin{bmatrix}$7.708\,68\text{\times}{10}^{6}$,$7.28204e{-4}$,$-3.56713e{-4}$,$0.530\,838$,$0.376\,335$,$8.785\,92$\end{bmatrix}^{T},$
(61)

where the unit for the first element is $\mathrm{m}$, the unit for the last
element is $\mathrm{rad}$, and the rest elements are unitless. The
corresponding covariance matrices of the initial state
$\mathbf{P}_{\text{MEE},0}$ and covariance of the process noise are obtained
via the unscented transformation [33]:

$\displaystyle\mathbf{P}_{\text{MEE},0}=\begin{bmatrix}$1.5519\text{\times}{10}^{6}$&$-0.063\,262\,5$&$0.140\,824$&$-0.045\,545\,6$&$-0.020\,913\,8$&$0.092\,296\,2$\\\
$-0.063\,262\,5$&$7.587\,18\text{\times}{10}^{-9}$&$-3.176\,09\text{\times}{10}^{-9}$&$1.885\,19\text{\times}{10}^{-9}$&$-1.851\,11\text{\times}{10}^{-9}$&$-1.712\,62\text{\times}{10}^{-9}$\\\
$0.140\,824$&$-3.176\,09\text{\times}{10}^{-9}$&$1.472\,01\text{\times}{10}^{-8}$&$-4.514\,09\text{\times}{10}^{-9}$&$-3.203\,95\text{\times}{10}^{-9}$&$1.000\,16\text{\times}{10}^{-8}$\\\
$-0.045\,545\,6$&$1.885\,19\text{\times}{10}^{-9}$&$-4.514\,09\text{\times}{10}^{-9}$&$1.821\,68\text{\times}{10}^{-9}$&$4.914\text{\times}{10}^{-10}$&$-3.245\,68\text{\times}{10}^{-9}$\\\
$-0.020\,913\,8$&$-1.851\,11\text{\times}{10}^{-9}$&$-3.203\,95\text{\times}{10}^{-9}$&$4.914\text{\times}{10}^{-10}$&$2.117\,49\text{\times}{10}^{-9}$&$-2.245\,64\text{\times}{10}^{-9}$\\\
$0.092\,296\,2$&$-1.712\,62\text{\times}{10}^{-9}$&$1.000\,16\text{\times}{10}^{-8}$&$-3.245\,68\text{\times}{10}^{-9}$&$-2.245\,64\text{\times}{10}^{-9}$&$7.210\,86\text{\times}{10}^{-9}$\\\
\end{bmatrix},$ (62)

$\displaystyle\mathbf{Q}_{\text{MEE}}=\begin{bmatrix}$461.406$&$-4.795\,89\text{\times}{10}^{-5}$&$3.562\,82\text{\times}{10}^{-5}$&$1.862\,65\text{\times}{10}^{-9}$&$-9.313\,23\text{\times}{10}^{-10}$&$-2.980\,23\text{\times}{10}^{-8}$\\\
$-4.795\,89\text{\times}{10}^{-5}$&$5.676\,55\text{\times}{10}^{-12}$&$-2.777\,37\text{\times}{10}^{-12}$&$-2.440\,54\text{\times}{10}^{-16}$&$1.814\,95\text{\times}{10}^{-16}$&$2.769\,49\text{\times}{10}^{-15}$\\\
$3.562\,82\text{\times}{10}^{-5}$&$-2.777\,37\text{\times}{10}^{-12}$&$4.001\,53\text{\times}{10}^{-12}$&$-4.984\,89\text{\times}{10}^{-16}$&$3.7058\text{\times}{10}^{-16}$&$3.908\,77\text{\times}{10}^{-15}$\\\
$1.862\,65\text{\times}{10}^{-9}$&$-2.440\,54\text{\times}{10}^{-16}$&$-4.984\,89\text{\times}{10}^{-16}$&$6.332\,16\text{\times}{10}^{-13}$&$-4.685\,97\text{\times}{10}^{-13}$&$-6.821\,21\text{\times}{10}^{-13}$\\\
$-9.313\,23\text{\times}{10}^{-10}$&$1.814\,95\text{\times}{10}^{-16}$&$3.706\,35\text{\times}{10}^{-16}$&$-4.686\,25\text{\times}{10}^{-13}$&$3.502\,48\text{\times}{10}^{-13}$&$5.098\,14\text{\times}{10}^{-13}$\\\
$1.490\,12\text{\times}{10}^{-8}$&$2.769\,49\text{\times}{10}^{-15}$&$3.9092\text{\times}{10}^{-15}$&$-6.821\,21\text{\times}{10}^{-13}$&$5.0937\text{\times}{10}^{-13}$&$7.673\,86\text{\times}{10}^{-13}$\\\
\end{bmatrix}$ (63)

An initial skewness of -1.6 and an initial kurtosis of 15.0 are given for each
of six MEEs.

The coordidate of the ground station in the ECEF frame is:

$\displaystyle\vec{\boldsymbol{x}}_{\text{stn}}=[-2730.0,3714.0,4394.0,0,0,0]^{T},$
(64)

where the unit is $\mathrm{km}$ for position components and
$\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ for velocity components. During the
three-day tracking campaign starting on 15/05/2021, the Sentinel 6A satellite
was only observed by three short arcs, each spanning approximately
$2\text{\,}\mathrm{min}$. The pre-residuals of the right ascension and
declination angles were calculated based on POD solutions of Sentinel 6A
mentioned above, shown in Fig. 3. Note that the declination residuals were
shiftted by $30\text{\,}\mathrm{min}$ intentionally for better visualisation.
It is indicated that both angles are subject to some errors and large
outliers, which lead to non-Gaussian pheanomenon of these measurements shown
by the distribution plots in the bottom subfigures of Fig. 3. The statistics
of two angular residuals is summarised in Table I.

Figure 3: Pre-residuals of Right Ascension and Declination Measurements Based on Precise Orbit Ephemeris of Sentinel 6A Satellite TABLE I: Statistics of Right Ascension and Declination Pre-residuals | Mean ($\mathrm{\SIUnitSymbolArcsecond}$) | Std. Dev. ($\mathrm{\SIUnitSymbolArcsecond}$) | RMS ($\mathrm{\SIUnitSymbolArcsecond}$) | Skewness | Kurtosis
---|---|---|---|---|---
Right Ascension | $4.650\text{\,}$ | $20.141\text{\,}$ | $20.671\text{\,}$ | -6.528 | 84.349
Declination | $1.460\text{\,}$ | $10.319\text{\,}$ | $10.421\text{\,}$ | -1.196 | 36.684

##  VI Results and Discussions

###  VI-A Projectile Example

####  VI-A1 Accuracy

100 Monte Carlo simulations (MCS) have been implemented by six filters,
including $w$-HOUSE, $\delta$-HOUSE, SRUKF, UKF, CUT-4 and CUT-6. Estimation
errors are shown in Fig. 4. The upper subfigures indicate the 3D position
errors while the lower subfigures indicate the 3D velocity errors. Note the
$96^{th}$ trial has been excluded for generating these plots due to large
errors for the first four filters. For the Gaussian case in the left column,
CUT-4 and CUT-6 result in both smaller position and velocity errors compared
to other four filters including $w$-HOUSE, generating two groups of curves.
For the non-Gaussion Pearson type IV case, two HOUSE filters stay close to
each other and generate more accurate estimates compared to other filters. In
particular, the errors quickly decrease by both HOUSE filters. Overall, these
results are a duplicate of those in Reference [14]. In two cases, SRUKF
overlaps with UKF for all types of errors, and $w$-HOUSE yields very similiar
results to $\delta$-HOUSE.

Figure 4: Estimation Errors for the Projectile Example

Fig. 5 shows the distribution of estimation errors for different filters. The
plot reproudces the results in Reference [14]. $w$-HOUSE generates similar
distributions to $\delta$-HOUSE and both are robust to outliers in the Pearson
type IV case.

Figure 5: Distribution of Estimation Error for the Projectile Example

####  VI-A2 Consumption Time

The average runtime over 100 MCS for all six filters are compared in Fig. 6.
Two HOUSE filters take slightly more computational time than two UKFs. The
runtime for CUT-4 and CUT6 are significantly larger. The $w$-HOUSE is even
slightly more efficient than $\delta$-HOUSE.

Figure 6: Runtime for the Projectile Example

###  VI-B Synthetic OD Example with Radar Measurements

####  VI-B1 Accuracy

Averaged over 100 Monte Carlo trials, the 3D position and velocity RMSE are
depicted in Fig. 7 for the synthetic OD example. All filters consistently
produce converged position and velocity estimates, with 3D RMSE decreasing
from the first track to the third track. Notably, the error curves for all
filters substantially overlap, with no discernible differences in both 3D
position and velocity errors. Table II presents the average RMSE across all
time epochs. The average 3D position RMSE for all filters falls within the
range of $4421\text{\,}\mathrm{m}4422\text{\,}\mathrm{m}$, while the average
3D velocity RMSE ranges from
$4.219\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}4.223\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$.
It is worth mentioning that the CUT filters exhibit the largest average 3D
position RMSE, whereas $\delta$-HOUSE yields the smallest value. The
$\delta$-HOUSE filter does exhibit slightly larger errors compared to the
SRUKF and UKF. These results can be attributed to the radar measurements used
in this numerical example, which encompass range, range rate, and two angles.
The information provided by these measurements effectively mitigates the non-
Gaussian phenomena, resulting in similar outcomes for all filters. The next
subsection presents a more challenging scenario where only two angles are
available.

Figure 7: 3D Position and Velocity RMSE for All Filters TABLE II: Average 3D Position and 3D Velocity Errors for All Filters (Units: $\mathrm{m}$ and $\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$) Filter Type | Average 3D Position Error | Average 3D Velocity Error
---|---|---
SRUKF | 4422.117 | 4.219
UKF | 4422.117 | 4.219
CUT4 | 4422.916 | 4.220
CUT6 | 4422.812 | 4.220
$\delta$-HOUSE | 4421.885 | 4.222
$w$-HOUSE | 4422.313 | 4.223

####  VI-B2 Time Consumption

Fig. 8 depicts the average run times for all filters across 100 MCS. Notably,
the SRUKF and UKF filters demonstrate the highest efficiency, while the CUT
filters are the most time-consuming. The $w$-HOUSE and $\delta$-HOUSE filters
exhibit comparable time consumption, with their runtime being approximately
twice that of two UKF filters.

Figure 8: Runtime of Synthetic OD for All Filters

###  VI-C Real-World OD Example with Optical Angles

The OD accuracy and precision via six filters are presented and analysed in
the following section. filters were run with MEE.

####  VI-C1 Absolute Position and Velocity Errors

The absolute position and velocity errors for three coordinates are shown in
Fig. 9 and in Fig. 10, respectively, each row of subfigures illustrating each
component. Three tracking arcs are shown separately in three columns. In each
subfigure, absolute errors generated by six filters are illustrated together.
Position and velocity errors decrease from the first track to the last. For
the first track, the absolute position and velocity errors by all filters stay
very close. For the last track, different filters yield larger discrepencies,
particularly for velocity estimates. At the begining phases of tracks 2 and 3
when measurements become available agin, all filters are subject to large
fluctuations. $\delta$-HOUSE yeilds some obvious and steep decreases on all
position components in the last track. The velocity estimates generated by
CUT-4 and CUT-6 scatter around. Given the uncertainty associated with the
orbital state does not match well with a Gaussian or uniform PDF, CUT filters
present even worse results compared to lower-order UKFs, i.e., the normal UKF
and SRUKF. Overall the proposed $w$-HOUSE method presents more smooth position
and velocity estimates than all other filters.

Figure 9: Absolute Position Errors for All Filters Figure 10: Absolute
Velocity Errors for All Filters

####  VI-C2 RMSE Statistics

The RMSE for three positon components and three velocity components are shown
in Table III and Table IV, respectively. All filters yield converged position
and velocity estimates in terms of each component and 3D values. The 3D
position RMSE decreases from the first track at around
$232.5\text{\,}\mathrm{m}$ to the last track betwen
$58.8\text{\,}\mathrm{m}70.1\text{\,}\mathrm{m}$. The 3D velocity RMSE
decreases from the first track at around
$0.5\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ to the last track at around
$0.1\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. Although the velocity RMSE
values are very close across all filters, the $w$-HOUSE filter yields the
smallest RMSE values for position estimates, followed by SRUKF and
$\delta$-HOUSE. CUT-4 and CUT-6 yield the largest RMSE values.

TABLE III: Position RMSE for All Filters (Unit: $\mathrm{m}$) Filter Type | Track | $\mathbf{x}$ | $\mathbf{y}$ | $\mathbf{z}$ | 3D Pos
---|---|---|---|---|---
$w$-HOUSE | 1 | $120.730$ | $114.857$ | $162.229$ | $232.564$
2 | $83.648$ | $91.689$ | $98.958$ | $158.734$
3 | $21.633$ | $23.765$ | $49.282$ | $58.834$
$\delta$-HOUSE | 1 | $120.708$ | $114.861$ | $162.342$ | $232.634$
2 | $82.335$ | $91.120$ | $99.745$ | $158.212$
3 | $19.492$ | $25.461$ | $50.657$ | $59.953$
SRUKF | 1 | $120.745$ | $114.853$ | $162.238$ | $232.576$
2 | $85.164$ | $92.946$ | $99.336$ | $160.498$
3 | $21.763$ | $24.071$ | $49.897$ | $59.521$
UKF | 1 | $120.732$ | $114.840$ | $162.223$ | $232.552$
2 | $83.209$ | $91.185$ | $98.620$ | $158.001$
3 | $25.404$ | $27.170$ | $51.579$ | $63.592$
CUT4 | 1 | $120.714$ | $114.821$ | $162.271$ | $232.568$
2 | $79.955$ | $88.787$ | $98.543$ | $154.876$
3 | $22.397$ | $29.474$ | $53.930$ | $65.412$
CUT6 | 1 | $120.695$ | $114.791$ | $162.160$ | $232.465$
2 | $85.800$ | $93.868$ | $99.644$ | $161.560$
3 | $26.641$ | $32.489$ | $56.165$ | $70.141$
TABLE IV: Velocity RMSE for All Filters (Unit: $\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$) Filter Type | Track | $\mathbf{vx}$ | $\mathbf{vy}$ | $\mathbf{vz}$ | $\mathbf{3DVel}$
---|---|---|---|---|---
$w$-HOUSE | 1 | $0.085$ | $0.495$ | $0.216$ | $0.547$
2 | $0.091$ | $0.085$ | $0.140$ | $0.188$
3 | $0.069$ | $0.054$ | $0.047$ | $0.099$
$\delta$-HOUSE | 1 | $0.085$ | $0.495$ | $0.216$ | $0.546$
2 | $0.099$ | $0.076$ | $0.127$ | $0.178$
3 | $0.090$ | $0.037$ | $0.028$ | $0.101$
SRUKF | 1 | $0.085$ | $0.495$ | $0.216$ | $0.547$
2 | $0.092$ | $0.085$ | $0.144$ | $0.191$
3 | $0.068$ | $0.055$ | $0.050$ | $0.101$
UKF | 1 | $0.085$ | $0.495$ | $0.216$ | $0.547$
2 | $0.095$ | $0.084$ | $0.140$ | $0.189$
3 | $0.064$ | $0.055$ | $0.053$ | $0.100$
CUT4 | 1 | $0.085$ | $0.495$ | $0.216$ | $0.547$
2 | $0.114$ | $0.071$ | $0.123$ | $0.182$
3 | $0.097$ | $0.030$ | $0.040$ | $0.109$
CUT6 | 1 | $0.085$ | $0.495$ | $0.216$ | $0.546$
2 | $0.097$ | $0.081$ | $0.139$ | $0.188$
3 | $0.083$ | $0.042$ | $0.023$ | $0.096$

####  VI-C3 Post-residuals

The post-residuals of two angles, i.e., right ascension and delination, based
on the position estimates were also calculated and they are presented in Fig.
11. Interestingly, the post-residuals for the second track are larger than
those for the first and third tracks, which is also indicated by the root mean
square (RMS) residuals of two angels shown in Table V. This is due to the fact
that right ascension angles of the second track are more contaminated by the
systematic biases indicated by the pre-residual analysis in Fig. 3. Smallest
post-residuals are obtained for the third track, the RMS values of the right
ascension are still larger than $9\text{\,}\mathrm{\SIUnitSymbolArcsecond}$
for all filters, while the RMS values of the declination drop down to less
than $4\text{\,}\mathrm{\SIUnitSymbolArcsecond}$. The $w$-HOUSE,
$\delta$-HOUSE, SRUKF, and UKF show consistent performance in terms of RMS
residuals across all tracks for both right ascension and declination angles.
The CUT4 filter demonstrates slightly higher RMS values while the CUT6 filter
stands out with lower RMS values compared to the other filters, indicating its
highest precision. Overall the differences are not significant among all
filters.

Figure 11: Post-residuals of Two Angles for All Filters TABLE V: RMS Residuals for All Filters (Unit: $\mathrm{\SIUnitSymbolArcsecond}$) Filter Type | Track | RA | Dec
---|---|---|---
$w$-HOUSE | 1 | 10.038 | 11.895
2 | 29.605 | 10.631
3 | 9.738 | 3.868
$\delta$-HOUSE | 1 | 10.038 | 11.896
2 | 29.458 | 10.613
3 | 9.523 | 3.865
SRUKF | 1 | 10.039 | 11.895
2 | 29.615 | 10.624
3 | 9.726 | 3.870
UKF | 1 | 10.038 | 11.894
2 | 29.671 | 10.656
3 | 9.461 | 3.876
CUT4 | 1 | 10.039 | 11.897
2 | 29.499 | 10.601
3 | 9.642 | 3.856
CUT6 | 1 | 10.037 | 11.893
2 | 29.345 | 10.539
3 | 9.294 | 3.867

####  VI-C4 Paramertic Analysis of HOUSE Filters

Sigma points generation in both $\delta$-HOUSE and $w$-HOUSE is in relation to
a parameter, respectively. A total of 100 MCS are conducted to analyse the
performances of two filters given different values of these parameters. Given
the requirement of nonnegativeness of $\delta$ values in $\delta$-HOUSE, 100
values are chosen uniformly from the range of $[0,0.2]$. In contrast, 100
values are chosen uniformly from the range of $[-0.1,0.1]$ for $w$ in
$w$-HOUSE, where the nonnegativeness of the last weight is not strictly
necessary. For $\delta$-HOUSE, only 13 sets of $\delta$ values out of 100
generate valid OD solutions, among which 2 values even lead to divergent
solutions. Shown in Fig. 12, trial numbers 73 and 95 yield very large 3D
position RMSE values since the second track. The rest trial numbers yield
converged OD solutions and it is indicated the 3D position RMSE values
decrease from the first track to the last one, reaching
$232.525\text{\,}\mathrm{m}$, $158.96\text{\,}\mathrm{m}$,
$65.089\text{\,}\mathrm{m}$ for each time window. On the contrary, all
selected 100 $w$ values lead to the same converged OD solution at
$232.564\text{\,}\mathrm{m}$, $158.734\text{\,}\mathrm{m}$,
$58.834\text{\,}\mathrm{m}$ for three tracks, respectively. It is explained by
the fact that weights of sigma points at each single step of $w$-HOUSE do not
exceed the threshold value set up by the $w$ value so the reset of the last
weight $w_{2n+1}$ is not triggered at all. Therefore, $w$-HOUSE shows superior
robustness and stability in terms of parameter dependency in contract to
$\delta$-HOUSE. In other words, it is easier for the newly developed $w$-HOUSE
to choose a parameter of $w$ for converged estimates than the original
$\delta$-HOUSE.

Figure 12: 3D Position RMSE Values w.r.t. Various $\delta$ Values in
$\delta$-HOUSE

####  VI-C5 Time Consumption

The runtimes for all filters are shown in Fig. 13. Likewise, SRUKF and UKF are
the most efficient, while CUT-4 and CUT-6 are the most time expensive. The
$w$-HOUSE and $\delta$-HOUSE filters consume similar amounts of time, less
than two times of that of two UKFs. The $w$-HOUSE is just slightly more time
consuming than $\delta$-HOUSE.

Figure 13: Runtime of Real-world OD for All Filters

##  VII Conclusion

This paper introduces a novel variant of the high-order unscented estimator
called $w$-HOUSE. Building upon the original $\delta$-HOUSE, the author
proposes the adoption of the square-root form of the covariance matrix. This
modification ensures a more numerically stable evolution of the probability
density function and alleviates the constraint imposed by the large kurtosis
of the state distribution on the moving weight of one sigma point.

The performance of the newly proposed $w$-HOUSE wss evaluated and compared
with several existing filters, including square-root unscented Kalman filter,
unscented Kalman filter, the conjugate unscented transformation family (CUT4
and CUT6), and the original $\delta$-HOUSE. The evaluation was conducted in
the context of three numerical examples, including a real-world angle-only
orbit determination problem of resident space objects, by considering
estimation accuracy & percision, robustness, and computational complexity.
Additionally, a comparative analysis was carried out between $\delta$-HOUSE
and $w$-HOUSE regarding their parametric dependency.

To facilitate orbit prediction within the OD problem, the author develops a
high-fidelity orbital propagator based on modified equinoctial elements. Real-
world tracking data in the format of two angles collected for the Sentinel 6A
satellite in low Earth orbit are utilised in the scenarios. The processing
results demonstrate that the proposed $w$-HOUSE method outperforms other
filters in terms of positioning accuracy. Remarkably, it achieves a root mean
square error of less than $60\text{\,}\mathrm{m}$ for 3D position estimation,
even with imperfect angle-only measurements contaminated with frequent
outliers, in a three-day OD scenario.

Both $w$-HOUSE and $\delta$-HOUSE exhibit significantly reduced computational
complexity comparable to CUT filters. Furthermore, the paper reveals that
$w$-HOUSE demonstrates greater robustness than $\delta$-HOUSE concerning the
selection of the dependent parameter. Consequently, the proposed $w$-HOUSE
emerges as a favourable choice for nonlinear and non-Gaussian estimation
problems, particularly in vital OD.

The availability of an OD toolbox equipped with a collection of aforementioned
derivative-free filters will provide the astrodynamics community with valuable
insights into emerging OD techniques for RSOs in LEO and other orbits.

##  VIII Acknowledgement

The author would acknowledge Dr Zhenwei Li from Changchun Observatory of
National Astronomical Observatory, Chinese Academy of Sciences for providing
valuable optical space tracking measurements and discussion on the measurement
quality. The preliminary implementation of the orbit determination toolbox by
the undergraduate student James Appleton supported by the UNSW Engineering
Taste of Research grant is also appreciated. The author would also acknowledge
the use of ChatGPT 3.5 for polishing the writing of this paper.

## Appendix: Sigma Point Sampling Based Linear Least Mean Squares Estimators

The process of sigma point sampling based linear least mean squares estimation
involves several steps. The first step is the generation of samples for the
augmented vector of state variables and measurement noise. Here, $N$ sigma
points are generated based on the moments of the PDF of the augmented state,
along with the associated weights for each sigma point. The set of sigma
points and weights is represented as
$\\{(\vec{\boldsymbol{x}}^{(j)}(k-1),\vec{\boldsymbol{\omega}}^{(j)}(k-1),w^{(j)}):j=1,2,...N\\}$,
with the weights normalised such that $\sum_{j=1}^{N}w^{j}=1$.

Following the generation of samples, the next step is the time update. In this
step, the sigma points are propagated through the nonlinear dynamics model to
obtain the predicted sigma points $\vec{\boldsymbol{x}}^{(j)}(k|k-1)$. This is
mathematically represented as:

$\displaystyle\vec{\boldsymbol{x}}^{(j)}(k|k-1)=\boldsymbol{F}\big{(}\vec{\boldsymbol{x}}^{(j)}(k-1),\vec{\boldsymbol{\omega}}^{(j)}(k-1)\big{)}.$
(65)

The predicted mean and covariance of these sigma points are then calculated as
follows:

$\displaystyle{\vec{\boldsymbol{x}}}_{m}(k|k-1)=\sum_{j=1}^{N}w^{(j)}\vec{\boldsymbol{x}}^{(j)}(k|k-1),$
(66)
$\displaystyle{\boldsymbol{P}}_{x}(k|k-1)=\sum_{j=1}^{N}w^{(j)}\vec{\delta\boldsymbol{x}}^{(j)}(k|k-1)\Big{[}\vec{\delta\boldsymbol{x}}^{(j)}(k|k-1)\Big{]}^{T}.$
(67)

The process then moves to the second round of sample generation for the
augmented vector of state variables and measurement noise. The set of sigma
points and weights is represented as
$\\{(\vec{\boldsymbol{x}}^{(j)}(k),\vec{\boldsymbol{\nu}}^{(j)}(k),w^{(j)}):j=1,2,...N\\}$,
with the weights normalised such that $\sum_{j=1}^{N}w^{j}=1$.

The final step is the measurement update. The deviations of the sigma points
are calculated as follows:

$\displaystyle\vec{\delta\boldsymbol{x}}^{(j)}(k)=\vec{\boldsymbol{x}}^{(j)}(k)-{\vec{\boldsymbol{x}}}_{m}(k|k-1).$
(68)

The predicted sigma points are propagated through the measurement model to
obtain the predicted measurements $\vec{\boldsymbol{z}}^{(j)}(k)$:

$\displaystyle{\vec{\boldsymbol{z}}}^{(j)}(k)=\boldsymbol{h}\big{(}\vec{\boldsymbol{x}}^{(j)}(k),\vec{\boldsymbol{\nu}}^{(j)}(k)\big{)},$
(69)

followed by the calculation of the predicted mean and covariance of
measurements as follows:

$\displaystyle{\vec{\boldsymbol{z}}}_{m}(k)=\sum_{j=1}^{N}w^{(j)}\vec{\boldsymbol{z}}^{(j)}(k),$
(70)
$\displaystyle\vec{\delta\boldsymbol{z}}^{(j)}(k)=\vec{\boldsymbol{z}}^{j}(k)-{\vec{\boldsymbol{z}}}_{m}(k),$
(71)
$\displaystyle{\boldsymbol{P}}_{z}(k)=\sum_{j=1}^{N}w^{(j)}\vec{\delta\boldsymbol{z}}^{(j)}(k)\Big{[}\vec{\delta\boldsymbol{z}}^{(j)}(k)\Big{]}^{T}.$
(72)

Subsequently, the Kalman gain is calculated as follows:

$\displaystyle{\boldsymbol{K}}(k)={\boldsymbol{P}}_{x}(k|k-1){\boldsymbol{P}}_{z}^{-1}(k),$
(73)

based on which the state estimate is updated as follows:

$\displaystyle\vec{\boldsymbol{x}}(k)=\vec{\boldsymbol{x}}(k|k-1)+{\boldsymbol{K}}(k)(\vec{\boldsymbol{z}}(k)-{\vec{\boldsymbol{z}}}_{m}(k)).$
(74)

In addition to the above, other orders of moments can also be updated. For
instance, the update of the second-order moment, i.e., the covariance, is also
necessary by using the Kalman gain as follows:

$\displaystyle{\boldsymbol{P}}_{x}(k)={\boldsymbol{P}}_{x}(k-1)-{\boldsymbol{K}}(k){\boldsymbol{P}}_{z}^{-1}(k){\boldsymbol{K}}^{T}(k).$
(75)

## References

  * [1] Z. Tan, Y. Yang, Y. Zhang, X. Zhang, and B. Lin, “Quantifying uncertainties of space objects centroid position based on optical observation,” in _Proc. 8th European Conference on Space Debris (virtual), Darmstadt, Germany, 20–23 April 2021_ , S. L. . F. S. T. Flohrer, Ed., no. 114. The European Space Agency Space Debris Office, 2021.
  * [2] R. E. Kalman, “A new approach to linear filtering and prediction problems,” _Journal of Basic Engineering_ , vol. 82, no. 1, pp. 35–45, 1960.
  * [3] S. F. Schdmit, “The Kalman filter - its recognition and development for aerospace applications,” _Journal of Guidance and Control_ , vol. 4, no. 1, pp. 4–7, 1981. [Online]. Available: https://doi.org/10.2514/3.19713
  * [4] G. Welch, G. Bishop _et al._ , “An introduction to the kalman filter,” 1995.
  * [5] Q. Lam, D. Junker, D. Anhalt, and D. Vallado, _Analysis of an extended Kalman filter based orbit determination system_ , 2012. [Online]. Available: https://arc.aiaa.org/doi/abs/10.2514/6.2010-7600
  * [6] A. Gelb, _Applied optimal estimation_. MIT press, 1974.
  * [7] S. J. Julier and J. K. Uhlmann, “Unscented filtering and nonlinear estimation,” _Proceedings of the IEEE_ , vol. 92, no. 3, pp. 401–422, 2004.
  * [8] K. Ito and K. Xiong, “Gaussian filters for nonlinear filtering problems,” _Automatic Control, IEEE Transactions on_ , vol. 45, no. 5, pp. 910–927, 2000.
  * [9] M. NøRgaard, N. K. Poulsen, and O. Ravn, “New developments in state estimation for nonlinear systems,” _Automatica_ , vol. 36, no. 11, pp. 1627–1638, 2000.
  * [10] I. Arasaratnam and S. Haykin, “Cubature Kalman filters,” _IEEE Transactions on automatic control_ , vol. 54, no. 6, pp. 1254–1269, 2009.
  * [11] N. Adurthi, P. Singla, and T. Singh, “Conjugate unscented transformation: applications to estimation and control,” _Journal of Dynamic Systems, Measurement, and Control_ , vol. 140, no. 3, 11 2017, 030907. [Online]. Available: https://doi.org/10.1115/1.4037783
  * [12] J. L. Crassidis, “Sigma-point Kalman filtering for integrated GPS and inertial navigation,” _Aerospace and Electronic Systems, IEEE Transactions on_ , vol. 42, no. 2, pp. 750–756, 2006.
  * [13] Y. Yang and H. Cai, “Analysis of adaptive Gaussian mixture unscented Kalman filter using sparse optical observations for initial orbit determination.” Proceedings of the 68th International Astronautical Conference, Adelaide, Australia, 2017.
  * [14] Z. Stojanovski and D. Savransky, “Higher-order unscented estimator,” _Journal of Guidance, Control, and Dynamics_ , vol. 44, no. 12, pp. 2186–2198, 2021. [Online]. Available: https://doi.org/10.2514/1.G006109
  * [15] M. J. H. Walker, B. Ireland, and J. Owens, “A set modified equinoctial orbit elements,” _Celestial mechanics_ , vol. 36, no. 4, pp. 409–419, 1985. [Online]. Available: https://doi.org/10.1007/BF01227493
  * [16] M. L. Psiaki, “Gaussian mixture filter for angles-only orbit determination in modified equinoctial elements,” _Journal of Guidance, Control, and Dynamics_ , vol. 45, no. 1, pp. 73–83, 2022. [Online]. Available: https://doi.org/10.2514/1.G005891
  * [17] S. Yun, R. Zanetti, and B. A. Jones, “Kernel-based ensemble Gaussian mixture filtering for orbit determination with sparse data,” _Advances in Space Research_ , vol. 69, no. 12, pp. 4179–4197, 2022.
  * [18] S. C. Patwardhan, S. Narasimhan, P. Jagadeesan, B. Gopaluni, and S. L. Shah, “Nonlinear Bayesian state estimation: A review of recent developments,” _Control Engineering Practice_ , vol. 20, no. 10, pp. 933–953, 2012, 4th Symposium on Advanced Control of Industrial Processes (ADCONIP).
  * [19] O. Montenbruck and E. Gill, _Satellite Orbits: Models, Methods and Applications_. Springer Science & Business Media, 2012.
  * [20] Y. Yang, E. Klein, C. He, H. Cai, and Y. Zhang, “PHiFA – a tool for numerical propagation of high-fidelity astrodynamics,” _Astrodynamics_ , vol. 6, no. 2, pp. 189–204, 2022. [Online]. Available: https://doi.org/10.1007/s42064-021-0113-1
  * [21] N. Das, V. Deshpande, and R. Bhattacharya, “Optimal-transport-based tracking of space objects using range data from a single ranging station,” _Journal of Guidance, Control, and Dynamics_ , vol. 42, no. 6, pp. 1237–1249, 2019.
  * [22] T. Söderström, _Discrete-time stochastic systems: estimation and control_. Springer Science & Business Media, 2002.
  * [23] S. J. Julier and J. K. Uhlmann, “New extension of the Kalman filter to nonlinear systems,” in _Signal processing, sensor fusion, and target recognition VI_ , vol. 3068. Spie, 1997, pp. 182–193.
  * [24] G. A. Terejanu, “Unscented Kalman filter tutorial,” _University at Buffalo, Buffalo_ , 2011.
  * [25] R. Van der Merwe and E. Wan, “The square-root unscented Kalman filter for state and parameter-estimation,” in _2001 IEEE International Conference on Acoustics, Speech, and Signal Processing. Proceedings (Cat. No.01CH37221)_ , vol. 6, 2001, pp. 3461–3464 vol.6.
  * [26] K. Ponomareva, P. Date, and Z. Wang, “A new unscented Kalman filter with higher order moment-matching,” in _Proceedings of Mathematical Theory of Networks and Systems (MTNS 2010), Budapest_ , 2010.
  * [27] O. Krause and C. Igel, “A more efficient rank-one covariance matrix update for evolution strategies,” in _Proceedings of the 2015 ACM Conference on Foundations of Genetic Algorithms XIII_ , ser. FOGA ’15. New York, NY, USA: Association for Computing Machinery, 2015, p. 129–136. [Online]. Available: https://doi.org/10.1145/2725494.2725496
  * [28] B. Tapley, J. Ries, S. Bettadpur, D. Chambers, M. Cheng, F. Condi, and S. Poole, “The GGM03 mean earth gravity model from GRACE,” in _AGU Fall Meeting Abstracts_ , vol. 2007, 2007, pp. G42A–03.
  * [29] Y. Yang, Z. Li, and H. Cai, “Single-/multiple-revolution Lambert’s problem constraints for optical track-to-track association.” Advances in the Astronautical Sciences AAS/AIAA Astrodynamics 2022, volume 177, 2021.
  * [30] J. F. Sanchez, M. Fernandez, P. Femenias, Y. Andres _et al._ , “Sentinel-6 orbit determination at the Copernicus POD service,” Copernicus Meetings, Tech. Rep., 2022.
  * [31] R. S. Park, W. M. Folkner, J. G. Williams, and D. H. Boggs, “The JPL planetary and lunar ephemerides DE440 and DE441,” _The Astronomical Journal_ , vol. 161, no. 3, p. 105, feb 2021. [Online]. Available: https://dx.doi.org/10.3847/1538-3881/abd414
  * [32] J. Hefty, M. Rothacher, T. Springer, R. Weber, and G. Beutler, “Analysis of the first year of earth rotation parameters with a sub-daily resolution gained at the code processing center of the igs,” _Journal of Geodesy_ , vol. 74, pp. 479–487, 2000.
  * [33] C. R. Binz and L. M. Healy, “Direct uncertainty estimates from angles-only initial orbit determination,” _Journal of Guidance, Control, and Dynamics_ , vol. 41, no. 1, pp. 34–46, 2018. [Online]. Available: https://doi.org/10.2514/1.G002227