Geomagnetic Field Measurements based Adaptive Cubature H-infinity Filter for Microsatellite Autonomous Orbit Estimation

1. Introduction

The autonomous navigation of microsatellites pertains to the technology for real-time determination of the position and velocity of microsatellites without reliance on ground support, utilizing solely the onboard measurement equipment. With the proliferation of space programs across different nations, there has been a surge in the number of micro-satellites operating in low Earth orbit, enhancing space application capabilities while imposing significant demands on ground telemetry and control systems. To mitigate the dependence of microsatellites on ground telemetry and control and enhance their autonomous intelligent operation in orbit, as well as serving as an effective complement to traditional ground telemetry and control, research into autonomous orbit determination technology for microsatellites has emerged as a pivotal focus within the realm of microsatellite technology [1,2,3].

The autonomous orbit determination method for microsatellites can be categorized into three main approaches: astrometric navigation orbit determination based on celestial measurements, autonomous navigation orbit determination based on beacon measurements, and geodetic/geographic information matching navigation orbit determination. The first approach involves measuring the position of known natural celestial bodies relative to the microsatellite and then estimating the orbital parameters using the least squares method based on the orbital dynamics model [4]. The second approach relies on signals received from beacons located at various positions with a specific pattern, enabling estimation of the orbital parameters through measured relative distances and rate of change using an orbital dynamics model [5]. The third approach utilizes onboard instruments such as optical imaging sensors, magnetometers, or gravity gradiometers to obtain measurement information related to geophysical fields or Earth's surface features for estimating orbital parameters by matching them with a prior model [6,7]. Geomagnetic navigation is particularly suitable for medium-low precision orbit determination of near-earth micro-satellites due to abundant magnetic field resources in near-earth orbits that can be described by Gaussian spherical harmonics [8,9]. The intensity and direction of the magnetic field can be expressed as a function of satellite position, making it feasible to use magnetic measurement methods for autonomous satellite orbit determination if a magnetometer meets requirements such as light weight, low power consumption, high reliability, and low cost.

The autonomous orbit estimation of a small satellite based on geomagnetic field measurements presents a fundamental nonlinear filtering challenge. Commonly employed solutions include the extended Kalman filter (EKF) [10], the unscented Kalman filter (UKF) [11,12,13], and the cubature Kalman filter (CKF) [14,15]. Among these, the EKF's use of first-order Taylor expansion compromises filtering accuracy. The UKF approximates the distribution of the nonlinear function using a set of sigma points, thereby enhancing estimation accuracy. However, its parameter selection lacks clear theoretical basis and may yield negative weights for variable dimensions exceeding 3, potentially leading to non-positive definite covariance matrices that directly impact estimate stability. The CKF decomposes the nonlinear Gaussian weighted integral into spherical surface and radial integrals, utilizing an equal-weighted set of cubature points to achieve superior filtering accuracy and numerical stability compared to EKF and UKF, making it widely applicable in engineering [16,17,18]. During application, filters often encounter model error and imprecise statistical noise characteristics; robustness refers to a filter's ability to withstand such challenges. The H-infinity filter is designed with improved robustness in mind addresses conservative filtering due to fixed constraint levels by proposing an adaptive robust filtering algorithm that comprehensively considers filter performance [19]. Furthermore, Ref. [20] designs a robust filtering algorithm capable of improving system adaptability to changing model parameters while addressing difficulties in precise modeling and unobtainable noise statistical characteristics. Ref. [21] combines H-infintiy filtering with CKF algorithms to propose a nonlinear robust filtering algorithm aimed at enhancing system robustness against abnormal observation values.

This paper presents an adaptive cubature H-infinity filtering algorithm. Firstly, it provides the low-altitude orbital dynamics model and the Earth's magnetic field measurement model of a small satellite. Subsequently, in the context of the adaptive cubature H-infinity filtering algorithm, on one hand, the spherical integral and radial integral are computed using the spherical simplex criterion and the second-order Gaussian-Legendre quadrature criterion respectively. The spherical simplex-radial cubature quadrature rule is then introduced to enhance the approximation accuracy of nonlinear Gaussian weighted integrals. On the other hand, game theory-based extended H-infinity filtering is integrated with the cubature rule to establish an inverse relationship between constraint level and filtering innovation. The constraint level is adaptively adjusted to improve robustness against model errors within this algorithm framework. Finally, a semi-physical simulation experiment demonstrates the effectiveness of this proposed algorithm.

The remainder of this paper is organized as follows: the autonomous orbit estimation models are described in Section 2, the adaptive robust cubature H-infinity filter is deduced in Section 3, the semi-physical simulation and results are presented in Section 4, and the conclusions are drawn in Section 5.

2. Autonomous Orbit Estimation Models

2.1. Dynamics Model of Low Orbit

For satellites in low Earth orbit, the primary perturbations experienced by the satellite in the Earth-fixed coordinate system ($O-XYZ$) are the non-spherical Earth term $J_2$ and atmospheric drag. When considering these perturbations, the orbital dynamics model of the satellite is as follows:

$$\left\{\begin{array}{l}\dot{x}=v_x\\ \dot{y}=v_y\\ \dot{z}=v_z\\ {\dot{v}}_x=\frac{\mu x}{r^3}\left[J_2\left(\frac{R_e}{r}\right)\left(7.5\frac{z^2}{r^2}-1.5\right)-1\right]+{\omega }_e^{2}x+2{\omega }_e{v}_y+f_{ax}+f_x\\ {\dot{v}}_y=\frac{\mu y}{r^3}\left[J_2\left(\frac{R_e}{r}\right)\left(7.5\frac{z^2}{r^2}-1.5\right)-1\right]+{\omega }_e^{2}y-2{\omega }_e{v}_x+f_{ay}+f_y\\ {\dot{v}}_z=\frac{\mu z}{r^3}\left[J_2\left(\frac{R_e}{r}\right)\left(7.5\frac{z^2}{r^2}-4.5\right)-1\right]+f_{az}+f_z\end{array}\right.$$

(1)

where $r=\sqrt{x^2+y^2+z^2}$, ${\mathit{X}}_p={\left(x,y,z\right)}^{\mathsf{{\rm T}}}$, $X_v={\left(v_x,v_y,v_z\right)}^{\mathsf{{\rm T}}}$ represent the position vector and velocity vector of the microsatellite in the $O-XYZ$, respectively. $J_2$ denotes the coefficient of the harmonic term, $\mu$ signifies the gravitational constant of the Earth, $R_e$ stands for the Earth's radius, and ${\omega }_e$ represents the Earth's rotational angular velocity. ${\mathit{f}}_a={\left(f_{a,x},f_{a,y},f_{a,z}\right)}^{\mathsf{{\rm T}}}=-\frac{1}{2}\frac{c_{d}A}{m}{\rho }_d{v}_{rel}{\mathit{v}}_{rel}$ accounts for atmospheric drag perturbation, where $c_d$ denotes the atmospheric drag coefficient, $A/m$ represents the satellite's frontal area-to-volume ratio, and ${\rho }_d$ stands for atmospheric density. Assuming that the atmosphere rotates with the Earth, the relative velocity between a satellite and the atmosphere is expressed as ${\mathit{v}}_{rel}=\dot{\mathit{X}}-{\mathsf{\omega }}_e\times \mathit{X}$, where ${\mathsf{\omega }}_e={\omega }_e\times {\left(0,0,1\right)}^{\mathsf{{\rm T}}}$ is indicative of Earth's rotational angular velocity vector. ${\mathit{f}}_d={\left(f_x,f_y,f_z\right)}^{\mathsf{{\rm T}}}$ encompasses third-order non-spherical perturbation summing with three-body perturbation and solar radiation pressure perturbation along all three coordinate axes.

In order to facilitate the implementation of the discrete Kalman filtering orbit determination algorithm on board a spacecraft computer, it is essential to discretize the continuous-time spacecraft orbital dynamics model presented in equation (1). In this study, we employ the classical fourth-order Runge-Kutta method as the discretization technique. For the continuous equation $\dot{y}=f\left(x\right)$, its general numerical computation format is: $y_{k+1}=y_k+\frac{h}{6}\left(K_1+2K_2+2K_3+K_4\right)$,$K_1=f\left(x_k,y_k\right)$,$K_2=f\left(x_k+h/2,y_k+\left(h/2\right)K_1\right)$,$K_3=f\left(x_k+h/2,y_k+\left(h/2\right)K_2\right)$,$K_4=f\left(x_k+h,y_k+h{K}_3\right)$. Utilizing the fourth-order Runge-Kutta method enables us to express equation (32) in the following discrete form:

$${\mathit{X}}_k=\mathit{f}\left({\mathit{X}}_{k-1}\right)+{\mathit{w}}_{k-1}$$

(2)

where ${\mathit{X}}_k={\left({\mathit{X}}_{p,k}^{\mathsf{{\rm T}}},{\mathit{X}}_{v,k}^{\mathsf{{\rm T}}}\right)}^{\mathsf{{\rm T}}}={\left(x_k,y_k,z_k,v_{x,k},v_{y,k},v_{z,k}\right)}^{\mathsf{{\rm T}}}$ denotes the orbital state at time k, $f\left(\cdot \right)$ signifies the non-linear functional relationship as depicted in equation (1), and ${\mathit{w}}_k$ represents the system noise.

2.2. Geomagnetic Field Measurement Model

The magnetometer offers the advantages of compact size, lightweight design, cost-effectiveness, reliable performance, and absence of field of view restrictions. It is capable of providing real-time continuous geomagnetic field measurement information in all weather conditions and is considered a fundamental component on current low-orbit microsatellites. The theoretical expression of the spherical harmonic model for the geomagnetic field elucidates that the primary source of the geomagnetic field lies within the Earth. The International Association of Geomagnetism and Aeronomy (IAGA) periodically releases a new iteration of the International Geomagnetic Reference Field (IGRF) every five years to furnish measurement data for ongoing refinement and enhancement of our understanding about the geomagnetic field. The potential function $V$ for the geomagnetic field can be represented as:

$$V\left(r,\theta ,\lambda \right)=R_e{\displaystyle \sum _{n=1}^{\infty }{\displaystyle \sum _{m=0}^n{\left(\frac{R_e}{r}\right)}^{n+1}\left(g_n^m\cos \left(m\lambda \right)+h_n^m\sin \left(m\lambda \right)\right)P_n^m\left(\cos \theta \right)}}$$

(3)

In this equation, $R_e=6378.140\mathrm{km}$ represents the reference equatorial radius of the Earth, R denotes the geocentric distance, λ stands for geographic longitude, $\theta$ signifies geocentric colatitude, $\theta =\pi /2-\phi$, and φ indicates geocentric latitude. $P_n^m$ refers to an nth-order mth-degree quasi-associated Legendre function, while $g_n^m$ and $h_n^m$ are denoted as the Gaussian coefficient. The IGRF model provides a table of corresponding Gaussian coefficients.

$$\left\{\begin{array}{l}{P}_n^m(\cos \theta )=\frac{1}{2^{n}n!}\sqrt{\frac{C_m(n-m)!{(1-{\cos }^2\theta )}^m}{(n+m)!}}\frac{d^{n+m}{({\cos }^2\theta -1)}^n}{d\cos {\theta }^{n+m}}\\ C_m=\left\{\begin{array}{l}1(m=0)\\ 2(m\ge 1)\end{array}\right.\end{array}\right.$$

(4)

The geomagnetic field vector $B_m$ can be expressed as the negative gradient of the geomagnetic potential function $V$ as follows.

$$B_m={\left(B_r,B_{\theta },B_{\lambda }\right)}^{\mathsf{{\rm T}}}=-\nabla V\left(r,\theta ,\lambda \right)$$

(5)

The relationship between the geomagnetic field vector and satellite position can be derived from equations (3) and (5). Specifically, the components of the geomagnetic field vector $B_m$ in the Earth-fixed coordinate system are as follows:

$$\left\{\begin{array}{l}{B}_r=-\frac{\partial V}{\partial r}={\displaystyle \sum _{n=1}^k{\left(\frac{R_e}{r}\right)}^{n+2}\left(n+1\right){\displaystyle \sum _{m=0}^n\left(g_n^m\cos \left(m\lambda \right)+h_n^m\sin \left(m\lambda \right)\right)P_n^m\left(\cos \theta \right)}}\\ B_{\theta }=-\frac{1}{r}\frac{\partial V}{\partial \theta }=-{\displaystyle \sum _{n=1}^k{\left(\frac{R_e}{r}\right)}^{n+2}{\displaystyle \sum _{m=0}^n\left(g_n^m\cos \left(m\lambda \right)+h_n^m\sin \left(m\lambda \right)\right)\frac{\partial P_n^m\left(\cos \theta \right)}{\partial \theta }}}\\ B_{\lambda }=-\frac{1}{r\sin \theta }\frac{\partial V}{\partial \lambda }=-\frac{1}{\sin \theta }{\displaystyle \sum _{n=1}^k{\left(\frac{R_e}{r}\right)}^{n+2}{\displaystyle \sum _{m=0}^n\left(-g_n^m\cos \left(m\lambda \right)+h_n^m\sin \left(m\lambda \right)\right)P_n^m\left(\cos \theta \right)}}\end{array}\right.$$

(6)

The higher the order of expansion for equation (6), the greater the model accuracy, and correspondingly, the greater the computational complexity. Let $\mathit{Z}={[B_r,B_{\theta },B_{\lambda }]}^{\mathsf{{\rm T}}}$, then equation (6) can be expressed in the form of the following discrete measurement equation.

$${\mathit{Z}}_k=h\left(X_k\right)+v_k$$

(7)

where ${\mathit{Z}}_k$ represents the three-axis geomagnetic field measured at time k by the magnetometer, $h\left(\cdot \right)$ denotes the nonlinear model function describing the three-axis geomagnetic field on the right-hand side, and $v_k$ signifies the measurement noise of the magnetometer at time k.

3. Adaptive Robust Cubature H-Infinity Filter

3.1. Suboptimal Solution of the Extended H-Infinity Filter

Consider the discrete nonlinear dynamic system described by equations (2) and (7).

$$\left\{\begin{array}{l}{\mathit{X}}_k=\mathit{f}\left({\mathit{X}}_{k-1}\right)+{\mathit{w}}_{k-1}\\ {\mathit{Z}}_k=\mathit{h}\left({\mathit{X}}_k\right)+{\mathit{v}}_k\end{array}\right.$$

(8)

where we make the assumption that the system noise ${\mathit{w}}_k$ and the measurement noise ${\mathit{v}}_k$ are bounded energy signals with unknown statistical properties. Specifically, we impose the following constraints: ${\displaystyle {\sum }_{k=0}^{\infty }{\mathit{w}}_k^{\mathsf{{\rm T}}}{\mathit{w}}_k}<\infty$, ${\displaystyle {\sum }_{k=0}^{\infty }{\mathit{v}}_k^{\mathsf{{\rm T}}}{\mathit{v}}_k}<\infty$.

Let $T_k\left(\mathit{F}\right)$ denote the mapping operation that associates the unknown perturbation, ${\mathit{X}}_0-{\widehat{\mathit{X}}}_0^{+}$,${\left\{{\mathit{w}}_i\right\}}_{i=0}^k$,${\left\{{\mathit{v}}_i\right\}}_{i=0}^k$ , with the estimated error ${\displaystyle {\sum }_{k=0}^{\infty }{\mathit{w}}_k^{\mathsf{{\rm T}}}{\mathit{w}}_k}<\infty$. Here, ${\mathit{X}}_0$represents the initial state of the system, ${\widehat{\mathit{X}}}_0$denotes its estimated value, and its variance ${\mathit{P}}_0$reflects the approximation accuracy of ${\widehat{\mathit{X}}}_0^{+}$to${\mathit{X}}_0$. Moreover, the objective of the H-infinity filter is to determine an estimation strategy $\mathit{F}$ such that minimizes the H-infinity norm of this operation $T_k\left(\mathit{F}\right)$. The H-infinity norm of this transfer operation is defined as:

$${‖T‖}_{\infty }=\underset{u\in l_2,u\ne 0}{\sup }\frac{{‖Tu‖}_2}{{‖u‖}_2}$$

(9)

In the equation, "sup" represents the supremum, and $u\in l_2\left[0,+\infty \right)$ is a bounded signal, that is ${\displaystyle {\sum }_{k=0}^{\infty }{u}_k^{\mathsf{{\rm T}}}{u}_k}<\infty$.

Nevertheless, in certain exceptional scenarios, closed-form solutions for the optimal H-infinity filtering can be derived, wherein a cost function is stipulated [22].

$$\mathit{J}=\frac{{\displaystyle \sum _{k=1}^N{‖{\mathit{X}}_k-{\widehat{\mathit{X}}}_k‖}^2}}{{‖{\mathit{X}}_0-{\widehat{\mathit{X}}}_0‖}_{{\mathit{P}}_0^{-1}}^2+{\displaystyle \sum _{k=1}^N\left({‖{\mathit{w}}_k‖}_{{\mathit{Q}}_k^{-1}}^2+{‖{\mathit{v}}_k‖}_{{\mathit{R}}_k^{-1}}^2\right)}}$$

(10)

where ${\mathit{Q}}_k$and ${\mathit{R}}_k$are weight matrices that can be adjusted based on performance requirements. In practical applications, ${\mathit{Q}}_k$and${\mathit{R}}_k$ can be regarded as estimations of the corresponding noise covariance matrices, with a norm ${‖\mathit{e}‖}_{{\mathit{S}}_k}^2={\mathit{e}}^{\mathsf{{\rm T}}}{\mathit{S}}_k\mathit{e}$.

The H-infinity filtering problem aims to minimize the normalized energy of the estimation error under various conditions, while also minimizing the energy of the input noise and initial error, thus satisfying:

$${‖\mathit{J}‖}_{\infty }=\sup \frac{{\displaystyle \sum _{k=1}^N{‖{\mathit{X}}_k-{\widehat{\mathit{X}}}_k‖}^2}}{{‖{\mathit{X}}_0-{\widehat{\mathit{X}}}_0‖}_{{\mathit{P}}_0^{-1}}^2+{\displaystyle \sum _{k=1}^N\left({‖{\mathit{w}}_k‖}_{{\mathit{Q}}_k^{-1}}^2+{‖{\mathit{v}}_k‖}_{{\mathit{R}}_k^{-1}}^2\right)}}\le {\gamma }^2$$

(11)

where $\gamma$ represents a predetermined threshold value known as the constraint level.

An analytical solution for suboptimal H-infinity filtering based on Krein space Kalman filtering, which is a suboptimal filter compared to the standard Kalman filter, differs from the latter primarily in the computation of the posterior covariance matrix. In suboptimal H-infinity filtering, the method for computing the posterior covariance matrix is as follows:"

$${\mathit{P}}_k^{+}={\left({\left({\mathit{P}}_k^{-}\right)}^{-1}+{\mathit{H}}_k^{\mathsf{{\rm T}}}{\mathit{R}}_k^{-1}{\mathit{H}}_k-{\gamma }^{-2}{\mathit{I}}_n\right)}^{-1}$$

(12)

It is evident that the covariance matrix of the extended H-infinity filter (EH∞F) incorporates an additional regulating term ${\gamma }^{-2}{\mathit{I}}_n$ in comparison to the covariance matrix of the standard Kalman filter. The iterative calculation steps for the EH∞F algorithm are as follows:

Compute the prior state estimation${\widehat{\mathit{X}}}_k^{-}=\mathit{f}\left({\widehat{\mathit{X}}}_{k-1}^{+}\right)$;

Compute the prior error covariance matrix${\mathit{P}}_k^{-}={\mathit{F}}_k{\mathit{P}}_{k-1}^{+}{\mathit{F}}_k^{\mathsf{{\rm T}}}+{\mathit{Q}}_k$;

Calculate the inverse of the posterior error covariance matrix${\left({\mathit{P}}_k^{+}\right)}^{-1}={\left({\mathit{P}}_k^{-}\right)}^{-1}+{\mathit{H}}_k^{\mathsf{{\rm T}}}{\mathit{R}}_k^{-1}{\mathit{H}}_k-{\gamma }^{-2}{\mathit{I}}_n$;

Compute the posterior state estimation${\widehat{\mathit{X}}}_k^{+}={\widehat{\mathit{X}}}_k^{-}+{\mathit{K}}_k\left({\mathit{Z}}_k-\mathit{h}\left({\widehat{\mathit{X}}}_k^{-}\right)\right)$;

Determine the Kalman filter gain matrix${\mathit{K}}_k={\mathit{P}}_k^{-}{\mathit{H}}_k^{\mathsf{{\rm T}}}{\left({\mathit{H}}_k{\mathit{P}}_k^{-}{\mathit{H}}_k^{\mathsf{{\rm T}}}+{\mathit{R}}_k\right)}^{-1}$.

where ${\mathit{F}}_k$and${\mathit{H}}_k$ denote the Jacobian matrices evaluated at ${\widehat{\mathit{X}}}_{k-1}^{+}$and ${\widehat{\mathit{X}}}_k^{-}$, respectively.

It is evident that the EH∞F exhibits a computational structure akin to that of the EKF, particularly when $\gamma \to \infty$. As such, the EH∞F simplifies to the EKF. Consequently, the EH∞F can be considered a more general algorithm, with $\gamma$ serving as a tuning parameter that balances robustness and precision in control algorithms.

3.2. Adaptive Determination of Constraint Level Values

The degree of constraint significantly influences the accuracy and robustness of the filtering process. Traditionally, it is assigned a fixed value based on engineering experience, leading to conservative filter performance. To strike a balance between filtering accuracy and robustness, an adaptive optimization approach should be employed.

To incorporate the EH∞F into the CKF computation framework, it is essential to adopt the following linear error propagation principle for approximating the measurement error covariance matrix and cross-error covariance matrix, ${\mathit{P}}_{xz,k}\approx {\mathit{P}}_k^{-}{\mathit{H}}_k^{\mathsf{{\rm T}}}$. Subsequently, the measurement update step of the cubature H∞ filter is computed as follows:

$${\mathit{H}}_k^{\mathsf{{\rm T}}}{\mathit{R}}_k^{-1}{\mathit{H}}_k={\left({\mathit{P}}_k^{-}\right)}^{-1}{\mathit{P}}_{xz,k}{\mathit{R}}_k^{-1}{\mathit{P}}_{xz,k}^{\mathsf{{\rm T}}}{\left({\mathit{P}}_k^{-}\right)}^{-1}$$

(13)

The inverse matrix of the posterior error covariance matrix for EH∞F is derived as follows.

$${\left({\mathit{P}}_k^{+}\right)}^{-1}={\left({\mathit{P}}_k^{-}\right)}^{-1}+{\left({\mathit{P}}_k^{-}\right)}^{-1}{\mathit{P}}_{xz,k}{\mathit{R}}_k^{-1}{\mathit{P}}_{xz,k}^{\mathsf{{\rm T}}}{\left({\mathit{P}}_k^{-}\right)}^{-1}-{\gamma }^{-2}{\mathit{I}}_n$$

(14)

The constraint level $\gamma$directly impacts the robust filter's performance. A smaller value enhances system robustness but reduces filtering accuracy, while a larger value compromises system stability and may lead to divergence. The selection of the constraint level $\gamma$ must satisfy the following Riccati inequality.

$${\left({\mathit{P}}_k^{-}\right)}^{-1}+{\left({\mathit{P}}_k^{-}\right)}^{-1}{\mathit{P}}_{xz,k}{\mathit{R}}_k^{-1}{\mathit{P}}_{xz,k}^{\mathsf{{\rm T}}}{\left({\mathit{P}}_k^{-}\right)}^{-1}-{\gamma }^{-2}{\mathit{I}}_n>0$$

(15)

To achieve the adaptive selection of the constraint level $\gamma$, we initially introduce the following lemma.

Lemma 1. [23] Suppose A and B are two n-order Hermitian matrices, $\mathit{A}>0$,$\mathit{B}\ge 0$, and then $\mathit{A}>\mathit{B}\iff {\lambda }_{\max }\left(\mathit{B}{\mathit{A}}^{-1}\right)<1$ ,where ${\lambda }_{\max }\left(\mathit{A}\right)$ represents the maximum eigenvalue.

The square of the filter innovation ${\mathit{e}}_k^{\mathsf{{\rm T}}}{\mathit{e}}_k$ can be utilized to denote the actual estimation error of the filter. A large value of ${\mathit{e}}_k^{\mathsf{{\rm T}}}{\mathit{e}}_k$indicates poor system stability, necessitating a reduction of $\gamma$ to enhance filter robustness. Conversely, a small value signifies high estimation accuracy, allowing for an increase in $\gamma$ to reduce the need for filter robustness. Thus, there exists an inverse relationship between the constraint level $\gamma$and ${\mathit{e}}_k^{\mathsf{{\rm T}}}{\mathit{e}}_k$.

The formula (15) is written in the following equivalent form:

$${\gamma }^2\mathit{I}>{\left[{\left({\mathit{P}}_k^{-}\right)}^{-1}+{\left({\mathit{P}}_k^{-}\right)}^{-1}{\mathit{P}}_{xz,k}{\mathit{R}}^{-1}{\mathit{P}}_{xz,k}^{\mathsf{{\rm T}}}{\left({\mathit{P}}_k^{-}\right)}^{-1}\right]}^{-1}$$

(16)

After defining the matrix $\mathit{A}={\left[{\left({\mathit{P}}_k^{-}\right)}^{-1}+{\left({\mathit{P}}_k^{-}\right)}^{-1}{\mathit{P}}_{xz,k}{\mathit{R}}^{-1}{\mathit{P}}_{xz,k}^{\mathsf{{\rm T}}}{\left({\mathit{P}}_k^{-}\right)}^{-1}\right]}^{-1}$, we obtain ${\gamma }^2\mathit{I}>\mathit{A}$, and according to Lemma 1, we have ${\lambda }_{\max }\left(\mathit{A}{\left({\gamma }^2\mathit{I}\right)}^{-1}\right)<1\Rightarrow {\lambda }_{\max }\left(\mathit{A}\right)<{\gamma }^2$, thus enabling us to establish:

$$\gamma =\alpha \sqrt{{\lambda }_{\max }\left(\mathit{A}\right)}$$

(17)

where $\alpha >1$, as there exists an inverse proportional relationship between $\gamma$and ${\mathit{e}}_k^{\mathsf{{\rm T}}}{\mathit{e}}_k$, let

$$a=1+\frac{\beta }{\sqrt{{\mathit{e}}_k^{\mathsf{{\rm T}}}{\mathit{e}}_k}}$$

(18)

where $\beta >0$ represents the phase relationship coefficient, typically determined through experimental means. Once $\beta$ is established, the coefficient $\alpha$ becomes solely dependent on the new information filter, effectively quantifying the relationship between $\gamma$ and ${\mathit{e}}_k^{\mathsf{{\rm T}}}{\mathit{e}}_k$.

3.3. Adaptive H-Infinity Filtering Algorithm with Adjustable Constraint Level

3.3.1. Spherical Simplex-Radial Cubature Quadrature Rule

Given the integral $I\left(\mathit{f}\right)={\displaystyle {\int }_{{\mathrm{R}}^n}\mathit{f}\left(\mathit{x}\right)e^{-{\mathit{x}}^{\mathsf{{\rm T}}}\mathit{x}}d\mathit{x}}$, let $\mathit{x}=r\mathit{y}$, where , $\mathit{y}$ satisfies ${\mathit{y}}^{\mathsf{{\rm T}}}\mathit{y}=1$ represents the surface of the unit sphere, and $r\ge 0$is the radius of the sphere. Then $I\left(\mathit{f}\right)$can be decomposed into a spherical surface integral $S\left(r\right)={\displaystyle {\int }_{U_n}\mathit{f}\left(r\mathit{y}\right)d\sigma \left(\mathit{y}\right)}$and a radial integral$R={\displaystyle {\int }_0^{\infty }S\left(r\right)r^{n-1}{e}^{-r^2}dr}$; however, obtaining analytical solutions for these integrals is generally challenging. Therefore, numerical approximation methods are considered. According to [24], a three-order spherical simplex rule consisting of $2n+2$ points can be used to approximate the spherical surface integral as follows:

$$S\left(r\right)=\frac{A_n}{2\left(n+1\right)}{\displaystyle \sum _{l=1}^{n+1}\left[\mathit{f}\left(r{\mathit{a}}_l\right)+\mathit{f}\left(-r{\mathit{a}}_l\right)\right]}$$

(19)

where $A_n=2\sqrt{{\pi }^n}/\mathsf{\Gamma }\left(n/2\right)$ represents the surface area of the n-dimensional unit sphere, while ${\mathit{a}}_l={\left[a_{l,1},a_{l,2},\cdots ,a_{l,n}\right]}^{\mathsf{{\rm T}}}$,$l=1,2,\cdots ,n+1$ denote the set of vertices of the n-dimensional simplex, with their elements defined as follows:

$$a_{l,m}=\left\{\begin{array}{l}-\sqrt{\frac{n+1}{n\left(n-m+2\right)\left(n-m+1\right)}},m<l\\ \sqrt{\frac{\left(n+1\right)\left(n-l+1\right)}{n\left(n-l+2\right)}},m=l\\ 0,m>l\end{array}\right.$$

(20)

For the radial integral$R={\displaystyle {\int }_0^{\infty }S\left(r\right)r^{n-1}{e}^{-r^2}dr}$, let $r^2=t$, resulting in $r=\sqrt{t}$ . Subsequently, we have $R=\frac{1}{2}{\displaystyle {\int }_0^{\infty }S\left(\sqrt{t}\right)t^{\frac{n-2}{2}}{e}^{-t}dt}$. Set $g\left(t\right)=S\left(\sqrt{t}\right)$,$\beta =\frac{n-2}{2}$and we obtain $R=\frac{1}{2}{\displaystyle {\int }_0^{\infty }g\left(t\right)t^{\beta }{e}^{-t}dt}$. This approach enables us to approximate the integral term using the Gauss-Laguerre quadrature rule [25] as follows.

$${\displaystyle {\int }_0^{\infty }g\left(t\right)t^{\beta }{e}^{-t}dt}\approx {\displaystyle \sum _{j=1}^p{A}_{j}g\left(t_j\right)}$$

(21)

where $t_j$ denotes the integration points, while $A_j$ represents the corresponding weights. The integration points can be derived from the solutions of p-th order Chebyshev-Laguerre polynomials as follows.

$$L_p^{\beta }\left(t\right)={\left(-1\right)}^p{t}^{-\beta }{e}^t\frac{d^p}{d{t}^p}\left(t^{\beta +p}{e}^{-t}\right)=0$$

(22)

The corresponding weights can be derived by solving the following equation.

$$A_j=\frac{p!\mathsf{\Gamma }\left(\beta +p+1\right)}{t_j{\left[{\dot{L}}_p^{\beta }\left(t_j\right)\right]}^2}$$

(23)

where the Gamma function $\mathsf{\Gamma }$ is utilized. The approximate precision of this rule relies on the number of points employed for integration, and for $p=2$, the subsequent approximation can be derived.

$$R=\frac{1}{2}\left[A_{1}S\left(\sqrt{t_1}\right)+A_{2}S\left(\sqrt{t_2}\right)\right]$$

(24)

Solving equations (22) and (23) yields the values of $t_1$，$t_2$，$A_1$ and $A_2$ as follows:

$$\left\{\begin{array}{l}{t}_1=\frac{n}{2}+1+\sqrt{\frac{n}{2}+1},A_1=\frac{n\mathsf{\Gamma }\left(\frac{n}{2}\right)}{2n+4+2\sqrt{2n+4}}\\ t_2=\frac{n}{2}+1-\sqrt{\frac{n}{2}+1},A_2=\frac{n\mathsf{\Gamma }\left(\frac{n}{2}\right)}{2n+4-2\sqrt{2n+4}}\end{array}\right.$$

(25)

The essence of the nonlinear Kalman filter lies in the computation of the Gaussian weighted integral $I_N={\displaystyle {\int }_{{\mathrm{R}}^n}\mathit{f}\left(\mathit{x}\right)N\left(\mathit{x};\overline{\mathit{x}},{\mathit{P}}_x\right)d\mathit{x}}$. Here, $N\left(\cdot \right)$ represents the Gaussian probability density function. The integral takes on the following identity form:

$$I_N=\frac{1}{\sqrt{{\pi }^n}}{\displaystyle {\int }_{R^n}\mathit{f}\left(\sqrt{2{\mathit{P}}_x}\mathit{x}+\overline{\mathit{x}}\right)e^{-{\mathit{x}}^{\mathsf{{\rm T}}}\mathit{x}}d\mathit{x}}$$

(26)

By combining equation (24) with $t_1$,$t_2$,$A_1$,$A_2$ we can derive the following spherical simplex-radial cubature quadrature rule for calculating $I_N$ as follows :

$$\begin{array}{l}{I}_N=\frac{n}{4\left(n+1\right)\left(n+2+\sqrt{2n+4}\right)}{\displaystyle \sum _{i=1}^{n+1}\left[\mathit{f}\left(\overline{\mathit{x}}+\sqrt{\left(n+2+\sqrt{2n+4}\right){\mathit{P}}_x}{\mathit{a}}_i\right)+\right.}\left.\mathit{f}\left(\overline{\mathit{x}}-\sqrt{\left(n+2+\sqrt{2n+4}\right){\mathit{P}}_x}{\mathit{a}}_i\right)\right]+\\ \frac{n}{4\left(n+1\right)\left(n+2-\sqrt{2n+4}\right)}{\displaystyle \sum _{i=1}^{n+1}\left[\mathit{f}\left(\overline{\mathit{x}}+\sqrt{\left(n+2-\sqrt{2n+4}\right){\mathit{P}}_x}{\mathit{a}}_i\right)+\right.}\left.\mathit{f}\left(\overline{\mathit{x}}-\sqrt{\left(n+2-\sqrt{2n+4}\right){\mathit{P}}_x}{\mathit{a}}_i\right)\right]\end{array}$$

(27)

In particular, when the condition in Equation (21) is met $p=1$, we can solve for $t_1=n/2$,$A_1=\mathsf{\Gamma }\left(n/2\right)$ , and subsequently derive the spherical simplex-radial cubature quadrature rule as presented in the literature [24]. This demonstrates that the proposed spherical simplex-radial cubature quadrature rule in this paper exhibits higher accuracy than the previously established spherical simplex-radial cubature rule.

The $\mathit{a}=\left[{\mathit{a}}_1,{\mathit{a}}_2,\cdots ,{\mathit{a}}_{n+1}\right]$ represents the matrix composed of ${\mathit{a}}_l$, ${\left[\cdot \right]}_i$ is used to represent the i-th column of the matrix, using formula (27) to construct the cubature points ${\widehat{\mathit{x}}}^{\left(i\right)}$and weights ${\omega }_i$as follows:

$${\widehat{\mathit{x}}}^{\left(i\right)}=\left\{\begin{array}{l}\overline{\mathit{x}}+\sqrt{\left(n+2+\sqrt{2n+4}\right){\mathit{P}}_x}{\left[\mathit{a},-\mathit{a}\right]}_i,i=1,2,\cdots ,2n+2\\ \overline{\mathit{x}}+\sqrt{\left(n+2-\sqrt{2n+4}\right){\mathit{P}}_x}{\left[\mathit{a},-\mathit{a}\right]}_{i-2n-2},i=2n+3,\cdots ,4n+4\end{array}\right.$$

(28)

$${\omega }_i=\left\{\begin{array}{l}\frac{n}{4\left(n+1\right)\left(n+2+\sqrt{2n+4}\right)},i=1,2,\cdots ,2n+2\\ \frac{n}{4\left(n+1\right)\left(n+2-\sqrt{2n+4}\right)},i=2n+3,\cdots ,4n+4\end{array}\right.$$

(29)

3.3.2. Adaptive H-Infinity Filter with Constrained Level Adjustment

In this paper, based on the system model equation (2), it is evident that n equals 6. Substituting this value into the equations (28) and (29) yields the cubature points and corresponding weights as follows:

$${\widehat{\mathit{x}}}^{\left(i\right)}=\left\{\begin{array}{l}\overline{\mathit{x}}+2\sqrt{3{\mathit{P}}_x}{\left[\mathit{a},-\mathit{a}\right]}_i,i=1,2,\cdots ,14\\ \overline{\mathit{x}}+2\sqrt{{\mathit{P}}_x}{\left[\mathit{a},-\mathit{a}\right]}_{i-14},i=15,\cdots ,28\end{array}\right.$$

(30)

$${\omega }_i=\left\{\begin{array}{l}\frac{1}{56},i=1,2,\cdots ,14\\ \frac{3}{56},i=15,\cdots ,28\end{array}\right.$$

(31)

The specific calculation steps for cubature H-infinity filtering with adaptive constraint level adjustment are provided, incorporating an adaptive adjustment algorithm for the constraint level.

Step 1: Filter initialization.

Given the initial values of the filter ${\widehat{\mathit{X}}}_0^{+}$，${\mathit{P}}_0^{+}$，$\gamma$;

In the context of recursion, $k=1,2,\cdots$, complete the subsequent recursive update

Step 2: Time update.

Compute the cubature points by utilizing the preceding orbital state${\widehat{\mathit{X}}}_{k-1}^{+}$ at time k-1 and the prior error covariance${\mathit{P}}_{k-1}^{+}$ as ${\widehat{\mathit{X}}}_{k-1}^{\left(i\right)}=\left\{\begin{array}{l}{\widehat{\mathit{X}}}_{k-1}^{+}+2\sqrt{3{\mathit{P}}_{k-1}^{+}}{\left[\mathit{a},-\mathit{a}\right]}_i,i=1,2,\cdots ,14\\ {\widehat{\mathit{X}}}_{k-1}^{+}+2\sqrt{{\mathit{P}}_{k-1}^{+}}{\left[\mathit{a},-\mathit{a}\right]}_{i-14},i=15,\cdots ,28\end{array}\right.$.

Determine the non-linear transfer points of orbital states${\mathit{X}}_k^{\left(i\right)}=\mathit{f}\left({\widehat{\mathit{X}}}_{k-1}^{\left(i\right)}\right)$.

Compute the prior estimation of the orbital state at time k as${\widehat{\mathit{X}}}_k^{-}=\frac{1}{56}{\displaystyle \sum _{i=1}^{14}{\mathit{X}}_k^{\left(i\right)}}+\frac{3}{56}{\displaystyle \sum _{i=15}^{28}{\mathit{X}}_k^{\left(i\right)}}$.

Calculate the prior error covariance matrix for orbital state as${\mathit{P}}_k^{-}=\frac{1}{56}{\displaystyle \sum _{i=1}^{14}{\omega }_i\left({\mathit{X}}_k^{\left(i\right)}-{\widehat{\mathit{X}}}_k^{-}\right){\left({\mathit{X}}_k^{\left(i\right)}-{\widehat{\mathit{X}}}_k^{-}\right)}^{\mathsf{{\rm T}}}}+\frac{3}{56}{\displaystyle \sum _{i=15}^{28}{\omega }_i\left({\mathit{X}}_k^{\left(i\right)}-{\widehat{\mathit{X}}}_k^{-}\right){\left({\mathit{X}}_k^{\left(i\right)}-{\widehat{\mathit{X}}}_k^{-}\right)}^{\mathsf{{\rm T}}}}+{\mathit{Q}}_{k-1}$;

Step 3: Measurement Update

Individually, utilize ${\widehat{\mathit{X}}}_k^{-}$and${\mathit{P}}_k^{-}$ to calculate the cubature points in accordance with the specified cubature criteria as ${\widehat{\mathit{X}}}_k^{\left(i\right)}=\left\{\begin{array}{l}{\widehat{\mathit{X}}}_k^{-}+2\sqrt{3{\mathit{P}}_k^{-}}{\left[\mathit{a},-\mathit{a}\right]}_i,i=1,2,\cdots ,14\\ {\widehat{\mathit{X}}}_k^{-}+2\sqrt{{\mathit{P}}_k^{-}}{\left[\mathit{a},-\mathit{a}\right]}_{i-14},i=15,\cdots ,28\end{array}\right.$.

Determine the non-linear transfer points of orbital states as${\mathit{Z}}_k^{\left(i\right)}=\mathit{h}\left({\widehat{\mathit{X}}}_k^{\left(i\right)}\right)$.

Calculate the predicted measurement of magnetometer as${\widehat{\mathit{Z}}}_k=\frac{1}{56}{\displaystyle \sum _{i=1}^{14}{\omega }_i{\mathit{Z}}_k^{\left(i\right)}}+\frac{3}{56}{\displaystyle \sum _{i=15}^{28}{\omega }_i{\mathit{Z}}_k^{\left(i\right)}}$.

Compute the matrix of measurement error covariance${\mathit{P}}_{z,k}=\frac{1}{56}{\displaystyle \sum _{i=1}^{14}{\omega }_i\left({\mathit{Z}}_k^{\left(i\right)}-{\widehat{\mathit{Z}}}_k\right){\left({\mathit{Z}}_k^{\left(i\right)}-{\widehat{\mathit{Z}}}_k\right)}^{\mathsf{{\rm T}}}}+\frac{3}{56}{\displaystyle \sum _{i=15}^{28}{\omega }_i\left({\mathit{Z}}_k^{\left(i\right)}-{\widehat{\mathit{Z}}}_k\right){\left({\mathit{Z}}_k^{\left(i\right)}-{\widehat{\mathit{Z}}}_k\right)}^{\mathsf{{\rm T}}}}+{\mathit{R}}_k$.

Compute the cross-covariance matrix. ${\mathit{P}}_{xz,k}=\frac{1}{56}{\displaystyle \sum _{i=1}^{14}{\omega }_i\left({\widehat{\mathit{X}}}_k^{\left(i\right)}-{\widehat{\mathit{X}}}_k^{-}\right){\left({\mathit{Z}}_k^{\left(i\right)}-{\widehat{\mathit{z}}}_k\right)}^{\mathsf{{\rm T}}}}+\frac{3}{56}{\displaystyle \sum _{i=15}^{28}{\omega }_i\left({\widehat{\mathit{X}}}_k^{\left(i\right)}-{\widehat{\mathit{X}}}_k^{-}\right){\left({\mathit{Z}}_k^{\left(i\right)}-{\widehat{\mathit{z}}}_k\right)}^{\mathsf{{\rm T}}}}$.

Compute the matrix of Kalman filter gain ${\mathit{K}}_k={\mathit{P}}_{xz,k}{\mathit{P}}_{z,k}^{-1}$;

Utilize the magnetic field measurements ${\mathit{Z}}_k$ obtained from the magnetometer at time k to refine the posterior estimation of the computed state${\widehat{\mathit{X}}}_k^{+}={\widehat{\mathit{X}}}_k^{-}+{\mathit{K}}_k\left({\mathit{Z}}_k-{\widehat{\mathit{Z}}}_k\right)$.

Compute the posterior error covariance matrix for the orbital state${\mathit{P}}_k^{+}={\left({\left({\mathit{P}}_k^{-}\right)}^{-1}+{\left({\mathit{P}}_k^{-}\right)}^{-1}{\mathit{P}}_{xz,k}{\mathit{R}}_k^{-1}{\mathit{P}}_{xz,k}^{\mathsf{{\rm T}}}{\left({\mathit{P}}_k^{-}\right)}^{-1}-{\gamma }^{-2}{\mathit{I}}_n\right)}^{-1}$.

Step 4: Adaptive calculation of the level value constraint

Conduct iterative recursive computations of the value $\gamma$ based on equations (17) and (18).

4. Simulation and Results

A semi-physical simulation experiment system for autonomous orbit determination of micro-satellite magnetic field has been established, as depicted in the figure. The experimental setup primarily comprises a micro-satellite simulator, a magnetic field coil, an air bearing turntable, a high-pressure air compressor, and relevant upper-level computer software. The magnetic field coil is capable of generating a desired strength magnetic field, and various disturbances can be introduced to simulate the in-orbit conditions of the micro-satellite as realistically as possible.

Figure 1. Microsatellite Autonomous Orbit Estimation Semi-physical Simulation Experiment System; (a) Micro-satellite Simulator and Earth Magnetic Field Coil; (b) Complementary Parts.

Figure 1. Microsatellite Autonomous Orbit Estimation Semi-physical Simulation Experiment System; (a) Micro-satellite Simulator and Earth Magnetic Field Coil; (b) Complementary Parts.

The micro-satellite simulator is equipped with a magnetometer to measure the magnetic field strength generated by the magnetic field coil, as depicted in the figure. Additionally, it incorporates a star sensor and gyroscope for attitude determination, along with a flywheel and a magnetic torque device for attitude control. The simulator utilizes the HPOP (High Precision Orbit Propagation) orbit prediction algorithm to run tests, accounting for perturbations such as atmospheric drag, solar pressure, and three-body gravity. The reference orbit's epoch is 1 Jun 2024 16:32:00 (UTCG), with a semi-major axis of 6878.137 km, eccentricity of 0, inclination of 97.035°, right ascension of the ascending node at 279.135°, argument of perigee at 0°, and true anomaly at 0°. The magnetometer boasts a measurement accuracy of 167 nT.

Figure 2. Magnetometer, gyroscope and flywheel.

Figure 2. Magnetometer, gyroscope and flywheel.

The root mean square error (RMSE) of position and velocity is employed for the assessment of the autonomous orbit estimation. The RMSE is defined as follows:

$$E_{position}(k)={\left[\frac{1}{N}{\displaystyle \sum _{n=1}^N\left({(r_x^n-{\widehat{r}}_x^n)}^2+{(r_y^n-{\widehat{r}}_y^n)}^2+\left.{(r_z^n-{\widehat{r}}_z^n)}^2\right)\right.}\right]}^{\frac{1}{2}}$$

(32)

where N denotes the number of Monte Carlo simulations, and velocity RMSE is defined analogously to equation (32).

The simulation results are presented in Figure 3 and Figure 4, demonstrating that the estimated RMSE of ACH∞F is smaller than that of CKF, indicating the superior resilience of ACH∞F against unmodeled errors and its enhanced filtering performance under such conditions. To quantitatively characterize the autonomous orbit estimation accuracy of the microsatellites, the average RMSE of orbit determination is statistically summarized in Table 1 over the time interval of 5000~18,000s. The comparison reveals a notable improvement in positioning accuracy and velocity accuracy by 7.49% and 3.1%, respectively, for ACH∞F as compared to the CKF algorithm.

5. Conclusions

In order to address the issue of reduced orbit estimation accuracy resulting from disparities between in-orbit measured geomagnetic field values and model values in autonomous orbit estimation of microsatellites based on geomagnetic field measurements, this paper proposes an adaptive cubature H-infinity filtering algorithm. Firstly, it presents the low-altitude orbital dynamics model and geomagnetic field measurement model for microsatellites. Subsequently, within the adaptive cubature H-infinity filtering algorithm, a spherical polyhedron criterion and a second-order Gaussian-Legendre quadrature criterion are utilized to compute the spherical integral and radial integral. Additionally, a spherical simplex-radial quadrature integration criterion is introduced to enhance the approximation accuracy of nonlinear Gaussian weighted integrals. Furthermore, by integrating extended H-infinity filtering based on game theory with the cubature criterion, an inverse proportional relationship between constraint level and innovation in filtering is established. This allows for adaptive adjustment of constraint level to enhance algorithm robustness against model errors. Finally, through semi-physical simulation experiments, it is demonstrated that compared with traditional algorithms, the proposed algorithm improves autonomous orbit estimation accuracy of position and velocity of microsatellites by 7.49% and 3.1%, thus validating its effectiveness.