Autonomous State Estimation and Observability Analysis for the Taiji Formation Using High-Precision Optical Sensors

Abstract

In certain observation periods of navigation missions for the Taiji formation, ground observation stations are unable to observe the spacecraft, while the state of the spacecraft can be estimated through the utilization of dynamic equations simulated on prior knowledge. However, this method cannot accurately track the spacecraft. In this paper, we focus on appropriately selecting the available onboard measurement to estimate the state of the spacecraft of the Taiji formation. We design two schemes to explore the performance of the state estimation based on the interspacecraft interferometry measurements and the measurements obtained from the Sun sensor and the radial velocity sensor. The observability of the system is numerically analyzed using the singular value decomposition method. Furthermore, we analyze error covariance propagation using the cubature Kalman filter. The results show that using high-precision interspacecraft angle measurement can improve significantly the observability of the system. The absolute position and velocity of the spacecraft can be estimated respectively with an accuracy of about 3.1 km and 0.14 m/s in the first scheme, where the prior information of the precision of the position and velocity is respectively 100 km and 1 m/s. When the measurement from the radial velocity sensor is used in the second scheme, the estimation accuracy of the velocity can be improved about 18 times better than that in the first scheme.

1. Introduction

The Taiji mission [1] is a type of space-based gravitational wave (GW) detection program that aims to detect GWs within a frequency range spanning from 0.1 mHz to 1.0 Hz. The Taiji mission consists of three spacecraft (SC), forming an SC formation known as the Taiji formation. Each SC follows a heliocentric orbit, forming a massive equilateral triangle with sides spanning approximately three million kilometers. Generally, navigation methods for such formations, which are located in heliocentric orbits, rely on a combination of deep space network (DSN) and other ground observations to estimate the absolute state of the SC [2,3]. However, due to the limitation of the observation arc, neither of these methods can provide continuous coverage for the entire observation period of the formation. As the Taiji formation is like that of LISA [4], we can borrow from the coverage analysis of observation of LISA. In the case of the LISA formation, there may be only a few hours a week of ground-based observation time, and the state of the SC is usually estimated based on the evolution of dynamical equations and a priori knowledge in the absence of ground observations [5]. These methods cannot track the SC well, especially when the SC performs unpredictable maneuvers. However, in the Taiji mission, each SC is equipped with a high-precision laser interspacecraft link measurement system [6] and other sensors; thus, these measurements may be used to perform the state estimation. These measurements include a laser measurement system to measure the relative ranges and range rates between SCs, and a differential wavefront sensing (DWS) system to measure the relative angles of SCs. Further, each SC is also equipped with a digital Sun sensor to determine the relative pointing angle between the SC and Sun. It is also possible to carry a spectrometer sensor to measure the radial velocity of the SC relative to Sun. Therefore, in this paper, we will focus on how to combine the available onboard measurements to autonomously estimate the state of the SC when the ground observations are insufficient.

The autonomous state estimation problem refers to estimating the SC state using only onboard measurement information. For the Taiji formation, there is currently no existing literature addressing the autonomous state estimation of the formation based on onboard sensors of the SC. Research approaches from another space-based GW detection formation with a design similar to the Taiji formation include the LISA formation [7], which can serve as valuable references. For both of these formations, the current approaches are to perform state estimation or orbit determination of the SC typically involving utilizing ground-based measurements (such as DSN and very long baseline interferometry (VLBI) measurement) in combination with interspacecraft measurements. For the LISA formation, Chung [2] investigated the absolute orbit determination and relative orbit determination of the SC utilizing ground observations and measurements of relative range and relative range rate between the SCs. This paper shows that when attempting to estimate the state of the SC solely based on interspacecraft range and range rate measurements, the filtering results would diverge, and the system is unobservable. Wang et al. [8] discovered that relying solely on interspacecraft measurements yields poor estimation of the absolute position state, while relative distance estimation remains accurate. For other formations with similar configurations, Hu et al. [9] found that the observability of the system using angles-only measurements is enhanced when there is a significant difference in eccentricity or inclination among SC orbits. On the other hand, studies [10,11,12] focusing on other interspacecraft measurements, such as relative range or relative range rate, indicate that using these measurements alone can determine the initial absolute orbit of the SC. However, when the SC orbits are coplanar, the state of the system using range or range rate measurement is usually unobservable.

Based on the aforementioned study on autonomous state estimation, we find that the observability of the system only using interspacecraft measurements becomes poor when the difference in eccentricity or inclination of SC orbits is small. For the Taiji or LISA formation, their nominal Keplerian orbit elements also exhibit these features, such as having the same inclination and eccentricity [13,14]. Consequently, the autonomous state estimation system may become unobservable or weakly observable in such cases. Studies [15,16,17] show that the addition of interspacecraft measurements can improve the accuracy of the state estimation. In addition to interspacecraft measurements, other onboard sensors can also be utilized. In these papers, they also show that, when using the measurements from the Sun sensor and the spectrometer sensor, the estimation accuracy of the absolute state of the SC can also be improved. Thus, like the study [18] that proposes a sensor selection algorithm to assess and compare the effects of different sensors on the accuracy of state estimation, for the Taiji mission, choosing the right sensors to model observation schemes for estimating the state of the SCs is an important study for the mission design.

In this paper, we investigate the problem of autonomous state estimation in the Taiji formation by utilizing all available onboard sensors. We specifically address state estimation in a multibody perturbation environment. Considering the payload configuration and sensor costs, for instance, the Taiji formation may opt not to equip a spectrometer for measuring relative velocity. We establish two observation schemes: the first case contains all measurements except the radial velocity of the spacecraft relative to Sun, and the second case contains all measurements. Furthermore, we explore the impact of high-precision optical angle sensors on the observability of the system using the singular value decomposition (SVD) method. Our findings demonstrate that the inclusion of high-precision angle sensors can enhance the observability of the system. Additionally, we analyze the propagation of error covariance of the state using the square cubature Kalman filter, particularly in a scenario with poor prior knowledge. Through our research, we achieve an estimated accuracy of approximately 3 km for the absolute position of the SC and approximately 0.1 m/s for the absolute velocity, where the prior information of the precision of the position and velocity is respectively 100 km and 1 m/s. These results contribute valuable a priori knowledge for the tasks in the Taiji mission, which need precision state information of SCs.

This paper is structured as follows: In Section 2, we introduce the dynamics and observation models of the Taiji mission. Section 3 focuses on the analysis of system observability using different schemes and measurement accuracies. In Section 4, we design a square root format of the cubature Kalman filter and then perform the simulation experiments of state estimation for the Taiji mission. Finally, Section 5 provides a summary of the main results of this paper.

2. System Model

The Taiji formation consists of three identical spacecraft (SC1, SC2, SC3). These three SCs form an approximately equilateral triangle formation with a side length of $3\times {10}^6$ km, orbiting Sun [19], as shown in Figure 1. The reference coordinate frame for the formation is selected from the solar center of a mass inertial coordinate system. In this coordinate system, the origin is positioned at the center of the mass of Sun, the fundamental plane aligns with the plane of the ecliptic, the x-axis points towards the equinox, the y-axis lies within the fundamental plane and is perpendicular to the x-axis, and the z-axis extends perpendicular to the origin, following the right-hand rule.

The nonlinear dynamical model characterizing the Taiji formation is expressed as follows:

$$\left\{\begin{array}{cc}\hfill \dot{\mathit{x}}\left(t\right)& =\mathit{f}\left(\mathit{x}\right(t\left)\right)+\mathit{w}\left(t\right)\hfill \\ \hfill \mathit{y}\left(k\right)& =\mathit{h}\left(\mathit{x}\right(k\left)\right)+\mathit{v}\left(k\right)\hfill \end{array}\right..$$

(1)

where $\mathit{x}$ represents the state vector of the system, $\mathit{f}\left(\mathit{x}\right(t\left)\right)$ corresponds to the nonlinear continuous dynamical model of the system, $\mathit{w}\left(t\right)$ represents the process noise of the system, $\mathit{h}\left(\mathit{x}\right(k\left)\right)$ denotes the discrete observation model of the system, and $\mathit{v}\left(k\right)$ represents the measurement noise of the state. A comprehensive explanation of these variables will be provided in the subsequent subsections.

2.1. Dynamics Model

The state vector of the system $\mathit{x}$ is defined as follows:

$$\begin{array}{cc}\hfill {\mathit{r}}_{\mathit{i}}& ={[x_i,y_i,z_i]}^T(i=1,2,3),\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \dot{{\mathit{r}}_{\mathit{i}}}& ={[\dot{x_i},\dot{y_i},\dot{z_i}]}^T(i=1,2,3),\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \mathit{x}& ={[{\mathit{r}}_{\mathbf{1}},{\mathit{r}}_{\mathbf{2}},{\mathit{r}}_{\mathbf{3}},\dot{{\mathit{r}}_{\mathbf{1}}},\dot{{\mathit{r}}_{\mathbf{2}}},\dot{{\mathit{r}}_{\mathbf{3}}}]}^T\in {\mathbb{R}}^{18},\hfill \end{array}$$

(4)

where i represents the index of the SC in the formation, ${\mathit{r}}_{\mathit{i}}$ represents the position state of the SC in Cartesian form, and $\dot{{\mathit{r}}_{\mathit{i}}}$ represents the velocity state of the SC in the coordinate frame.

$\mathit{f}\left(\mathit{x}\right(t\left)\right)$ represents the nonlinear dynamical equations of the system:

$$\begin{array}{cc}\hfill \mathit{f}\left(\mathit{x}\right(t\left)\right)& =\left[\begin{array}{c}\dot{{\mathit{r}}_{\mathbf{1}}}\hfill \\ \dot{{\mathit{r}}_{\mathbf{2}}}\hfill \\ \dot{{\mathit{r}}_{\mathbf{3}}}\hfill \\ \ddot{{\mathit{r}}_{\mathbf{1}}}\hfill \\ \ddot{{\mathit{r}}_{\mathbf{2}}}\hfill \\ \ddot{{\mathit{r}}_{\mathbf{3}}}\hfill \end{array}\right]\in {\mathbb{R}}^{18}.\hfill \end{array}$$

(5)

where

$$\ddot{{\mathit{r}}_{\mathit{i}}}=-\frac{{\mu }_s}{r_i^3}{\mathit{r}}_{\mathit{i}}+\sum _{j=1}^{n_p}{\mu }_j\left[\frac{{\mathit{r}}_{pj}-{\mathit{r}}_{\mathit{i}}}{\parallel {\mathit{r}}_{pj}-{\mathit{r}}_{\mathit{i}}{\parallel }^3}-\frac{{\mathit{r}}_{pj}}{r_{pj}^3}\right](i=1,2,3,j=1,2,3,\dots ,6),$$

(6)

where ${\mu }_s$ represents the gravitational constant of Sun, $r_i=\sqrt{x_i^2+y_i^2+z_i^2}$ represents the range from the i-th SC to the center of Sun, ${\mu }_j$ represents the gravitational constant of the j-th perturbing planet, ${\mathit{r}}_{pj}$ represents the heliocentric coordinate system position vector of the j-th perturbing planet, ${\mathit{r}}_{pj}-{\mathit{r}}_{\mathit{i}}$ represents the position vector of the j-th perturbing planet relative to the i-th SC, and $n_p=6$ is the number of perturbing planets. For each SC, we have considered the gravitational influence of six celestial bodies, including Mercury, Venus, Earth, Mars, Jupiter, and Saturn.

$\mathit{w}\left(t\right)$ and $\mathit{v}\left(t\right)$ are modeled as the zero-mean Gaussian white noise and are not related to each other:

$$\left\{\begin{array}{cc}\hfill \mathit{w}\left(t\right)& \sim (0,Q_c)\hfill \\ \hfill \mathit{v}\left(k\right)& \sim (0,R_k)\hfill \\ \hfill E\left(\mathit{w}\left(t\right)\mathit{w}{(t-\tau )}^T\right)& =Q_c\delta (t-\tau )\hfill \\ \hfill E\left(\mathit{v}\left(k\right)\mathit{v}{(k-\tau )}^T\right)& =R_k\delta (k-\tau )\hfill \\ \hfill E\left(\mathit{w}\left(t\right)\mathit{v}{\left(k\right)}^T\right)& =0\hfill \end{array}\right..$$

(7)

where $Q_c$ is the covariance of the process noise distribution, and $R_k$ is the covariance of the measurement noise distribution.

2.2. Observational Model

In our study, we employ a high-precision interspacecraft measurement system along with other onboard sensors on the SC to gather measurement information for the purpose of state estimation. The interspacecraft measurements for the SC primarily encompass relative range measurements, relative range rate measurements, and relative angle measurements. In addition, we make use of the digital Sun sensor and spectrometer installed on the SC as supplementary measurement sensors. It is important to note that the measurement information considered in our study undergoes preprocessing and time synchronization. Consequently, we provide detailed descriptions of the specific instrumentation used and present simplified mathematical models for interpreting these measurement data.

2.2.1. Interspacecraft Measurements

Space-based GW detection formations employ a laser interferometry instrument [20] with extremely high measurement accuracy to capture the information between spacecraft. For the sake of convenience, we initially denote the relative position and relative velocity vectors between the spacecraft:

$$\begin{array}{cc}\hfill {\mathit{r}}_{\mathit{ij}}& =[x_i-x_j,y_i-y_j,z_i-z_j],\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \dot{{\mathit{r}}_{\mathit{ij}}}& =[\dot{x_i}-\dot{x_j},\dot{y_i}-\dot{y_j},\dot{z_i}-\dot{z_j}],\hfill \end{array}$$

(9)

where i and j represent the indices of the SC. The relative range between the i-th SC and the j-th SC is measured by the laser ranging system, which can be modeled as

$$\begin{array}{c}\hfill L_{ij}=\left|\right|{\mathit{r}}_{\mathit{ij}}\left|\right|,\end{array}$$

(10)

The relative range rate measured by the Doppler shift measurement system can be modeled as

$$\dot{L_{ij}}={\left(\frac{{\mathit{r}}_{\mathit{ij}}}{\left|\right|{\mathit{r}}_{\mathit{ij}}\left|\right|}\right)}^T\dot{{\mathit{r}}_{\mathit{ij}}}.$$

(11)

The relative angles between spacecraft are measured using the DWS system, a widely recognized technique for precisely quantifying the relative wavefront misalignment between two beams with high sensitivity. In this context, the measured relative angles include the azimuth angle denoted by $\lambda$ and the elevation angle denoted by $\varphi$, which are generally measured using CCD and quadrant photodiode (QPD) sensors. In the Taiji program, relative angles can be measured with a very high accuracy of 1 nrad [21,22]. The DWS measurement of these angles can be modeled as follows:

$$\begin{array}{cc}\hfill {\lambda }_{ij}& =arctan\left(\frac{y_{ij}}{x_{ij}}\right),\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill {\varphi }_{ij}& =arcsin\left(\frac{z_{ij}}{\left|\right|{\mathit{r}}_{\mathit{ij}}\left|\right|}\right),\hfill \end{array}$$

(13)

where $x_{ij}=x_i-x_j$ and $y_{ij}$, $z_{ij}$ are similar to it. Moreover, the relative angle measurements can also be alternatively modeled using the light-of-sight (LOS) vector, which provides an equivalent representation:

$$A_{ij}=\frac{{\mathit{r}}_{\mathit{ij}}}{\left|\right|{\mathit{r}}_{\mathit{ij}}\left|\right|}.$$

(14)

2.2.2. Other Measurements Onboard the Spacecraft

The measurement of the relative Sun-pointing angles between Sun and the i-th SC is obtained using the digital Sun sensor, which can be modeled as follows:

$$A_{si}=\frac{{\mathit{r}}_{\mathit{i}}}{\left|\right|{\mathit{r}}_{\mathit{i}}\left|\right|},$$

(15)

The radial velocity can be measured using the spectrometer on the SC, which detects the Doppler shift caused by the relative motion between Sun and the i-th SC. This measurement can be modeled as follows:

$$V_i={\left(\frac{{\mathit{r}}_{\mathit{i}}}{\left|\right|{\mathit{r}}_{\mathit{i}}\left|\right|}\right)}^T\dot{{\mathit{r}}_{\mathit{i}}}.$$

(16)

2.2.3. Measurement Equations

Considering the payload constraints of the SC, we have developed two observation model schemes as demonstrated in

Table 1. Measurement cases.

	Interspacecraft Measurements			Sun Angle Measurements	Radial Velocity Measurements
	Range	Range Rate	Angle
case 1	✓	✓	✓	✓
case 2	✓	✓	✓	✓	✓

. The first scheme, referred to as case 1, does not include the radial velocity information provided by the spectrometer. On the other hand, the second scheme, known as case 2, incorporates this additional measurement.

These two observation models of the system can be written as

$${\mathit{h}}_{\mathbf{1}}\left(\mathit{x}\right)={\left[L_{12},L_{13},L_{23},\dot{L_{12}},\dot{L_{13}},\dot{L_{23}},A_{12},A_{13},A_{23},A_{s1},A_{s2},A_{s3}\right]}^T\in {\mathbb{R}}^{24},$$

(17)

$${\mathit{h}}_{\mathbf{2}}\left(\mathit{x}\right)={\left[{\mathit{h}}_{\mathbf{1}}\left(\mathit{x}\right),V_1,V_2,V_3\right]}^T\in {\mathbb{R}}^{27}.$$

(18)

3. Observability Analysis

The observability analysis of a system provides a crucial theoretical foundation for state estimation. It helps determine whether the system’s state can be accurately estimated based on available measurement information. If the system is observable, the convergence of the state error occurs gradually as the measurement information is continually updated. The observability of nonlinear systems is typically assessed using the rank criterion approach, where the rank of the observability matrix serves as an indicator of the system’s observability level. In this study, we employ the observability Gramian matrix [23] to compute the system’s observability matrix. Specifically, for discrete systems, the observability Gramian matrix without considering measurement noise can be expressed as follows:

$$\begin{array}{c}\hfill {\mathit{G}}_{0,n}=\sum _{k=0}^n{\mathrm{\Phi }}_{0,k}^T{\mathit{H}}_k^T{\mathit{H}}_k{\mathrm{\Phi }}_{0,k},\end{array}$$

(19)

where n represents the total number of measurement instances starting from $t_0$. ${\mathrm{\Phi }}_{0,k}$ denotes the system transfer matrix from $t_0$ to $t_k$, while ${\mathit{H}}_k$ represents the Jacobian derivative matrix derived from the observation model at time $t_k$.

$$\begin{array}{cc}\hfill {\mathit{A}}_k& =\frac{\partial \mathit{f}}{\partial \mathit{x}}{|}_{{\mathit{x}}_k},\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill {\mathit{H}}_k& =\frac{\partial \mathit{h}}{\partial \mathit{x}}{|}_{{\mathit{x}}_k},\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_k& =e^{{\mathit{A}}_k\mathrm{\Delta }t},\hfill \end{array}$$

(22)

In the provided equation, for small sampling times $\mathrm{\Delta }t$, the system transfer matrix ${\mathrm{\Phi }}_k$ can be approximated as ${\mathrm{\Phi }}_k=\mathit{I}+\mathrm{\Delta }t{\mathit{A}}_k$. To numerically quantify the initial observability of the system at time $t_k$, we employ the singular value decomposition (SVD) method to analyze the observability Gramian matrix ${\mathit{G}}_{0,n}$ as follows:

$$\begin{array}{c}\hfill {\mathit{G}}_{0,n}=U\Sigma V^T.\end{array}$$

(23)

In the last equation, the matrix $\Sigma$ represents the diagonal matrix composed of the singular values of the observability Gramian matrix ${\mathit{G}}_{0,n}$, while U and V correspond to the left and right unitary matrices, respectively. We employ two metrics to measure the degree of observability of the system:

$$\begin{array}{cc}\hfill OI& ={\sigma }_{min},\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill CN& =\frac{{\sigma }_{max}}{{\sigma }_{min}}.\hfill \end{array}$$

(25)

The variable ${\sigma }_{min}$ represents the smallest individual singular value obtained from the matrix $\Sigma$, and ${\sigma }_{max}$ represents the largest individual singular value obtained from that. The observability index (OI) signifies the measure of observability for the observability matrix, while the condition number (CN) denotes the numerical condition of the observability matrix. A lower OI value or a higher CN value signifies reduced system observability. These two indicators serve as evaluative metrics to assess the degree of observability of the system. For our study, we numerically calculate these two metrics, OI and CN, for both case 1 and case 2, using the initial state of the system. The initial state vector for each SC within the Taiji formation in the inertial frame can be derived from the initial Keplerian orbital elements [24].

One initial Keplerian orbital parameter for the Taiji formation is present in

Table 2. Initial Keplerian orbital elements for the Taiji formation.

	a (m)	e	i (rad)	$\mathbf{\Omega }$ (rad)	$\mathit{\omega }$ (rad)	f (rad)
SC1	149596345656	0.0057765	0.41885699 04957726	6.2789100 17676982	4.5455776 28265889	3.14159265 3580112
SC2	149596558992	0.0057787	0.40563848 15646019	0.0236941 87073916	0.3314021 192709651	1.05695599 21295202
SC3	149596971939	0.0057800	0.40266269 23819615	6.2634647 12324839	2.4647099 56430478	−1.0565641 717284677

, encompassing the semimajor axis a, eccentricity e, inclination i, right ascension of the ascending node $\mathrm{\Omega }$, argument of the perihelion $\omega$, and true anomaly f. To determine the true value of the measurement information, we calculate it indirectly based on the state value at the initial time $t_0$. For both case 1 and case 2, we compute the condition number (CN) and observability index (OI) of the initial observability Gramian matrix ${\mathit{G}}_{0,0}$ as demonstrated in

Table 3. The observability metrics for two observation schemes.

	OI	CN
case 1	$3.251\times {10}^{-32}$	$9.227\times {10}^{31}$
case 2	$2.404\times {10}^{-31}$	$1.248\times {10}^{31}$

.

Our findings indicate that the observability of the system employing an observation scheme similar to case 1 is slightly weaker compared with the system employing the case 2 scheme. This disparity may be attributed to the principle that the inclusion of a greater variety of observation types enhances the effectiveness of state estimation.

The observability of a system is also influenced by the accuracy of the measurement information. For instance, Geller et al. [25] discovered that when the angle measurement accuracy is on the order of 1 mrad, certain initial orbit determination challenges become difficult to resolve. Additionally, Yim et al. [26] found that there is a linear relationship between the observability of a system utilizing angles-only measurements and the angle measurement accuracy. In the Taiji formation, each SC is equipped with high-precision angle measurement sensors, such as CCD and QPD. In this study, our focus lies in investigating the extent to which the accuracy of the angle measurement sensor impacts the observability of the system.

We also use the observability Gramian matrix with measurement noise to investigate this issue:

$$\begin{array}{c}\hfill {\mathit{L}}_{0,n}=\sum _{k=0}^n{\mathrm{\Phi }}_{0,k}^T{\mathit{H}}_k^T{\mathit{R}}_k^{-1}{\mathit{H}}_k{\mathrm{\Phi }}_{0,k}\end{array}$$

(26)

Here, ${\mathit{R}}_k$ represents the measurement covariance of the system at time k. To quantify the impact of measurement accuracy on the observability of the system, we employ a ratio metric proposed by Hu et al. [9] to measure the effect of measurement accuracy on the observability of the system, where a smaller value indicates weaker observability of the system.

$$O{I}_r=\frac{O{I}_L}{O{I}_G}=\frac{min\sigma \left({\mathit{L}}_{0,k}\right)}{min\sigma \left({\mathit{G}}_{0,k}\right)}.$$

(27)

In our research, we focus on the case 2 observation scheme, which involves the utilization of almost all available measurement sensors on the SC. To assess the observability of the system, we compare three levels of angle measurement accuracy. These accuracy levels represent different types of sensors: low-precision sensors, reflecting current sensor capabilities; common sensors, planned to be installed on the SC; and high-precision sensors that may be used in the future. These accuracy levels are separated by a factor of ${10}^3$. We assume that the measurement noise of the sensors follows a Gaussian distribution, and the standard deviation of the sensors for each level is presented in

Table 4. The measurement standard deviation of measurement sensors.

	Interspacecraft Measurements			Sun Angle Measurements	Radial Velocity Measurements
	Range (m)	Range Rate (m/s)	Angle (rad)	Sun-Pointing Angle (rad)	Range Rate to Sun (m/s)
level 1 (low-precision)	1	$1\times {10}^{-6}$	$1\times {10}^{-3}$	$1\times {10}^{-5}$	$1\times {10}^{-2}$
level 2 (common)			$1\times {10}^{-6}$
level 3 (high-precision)			$1\times {10}^{-9}$

:

To assess and compare the observability of the system, we calculate the observability index ratio ($O{I}_r$) using 1000 data points within each interval. The measurement standard deviation step for each level is set as $1\times {10}^{-9}$ from level 3 to level 2. From level 2 to level 1, the measurement standard deviation step for a sensor is set as $1\times {10}^{-6}$. The $O{I}_r$ values for the system with varying measurement accuracy are presented below:

Table 5. OI and CN for different measurement accuracy levels.

	OI	CN
level 1	$2.090\times {10}^{-26}$	$1.436\times {10}^{26}$
level 2	$1.542\times {10}^{-23}$	$1.945\times {10}^{23}$
level 3	$1.244\times {10}^{-20}$	$2.412\times {10}^{20}$

below presents OI and CN of the observability matrix of the system at different levels of measurement accuracy.

From Figure 2, we observe that the degree of observability of the system demonstrates a nearly linear growth pattern as the accuracy of angle measurement increases, particularly from level 2 to level 1. However, we find that the observability of the system shows minimal improvement from level 3 to level 2, despite the increase in angle measurement accuracy from approximately $3\times {10}^{-9}$ rad to $1\times {10}^{-9}$ rad. There are several factors contributing to these observations. First, in the absence of ground-based observations, the system tends to exhibit weak observability. Additionally, the unique orbit configurations of the Taiji formation, such as having orbits with similar inclinations or eccentricities, also contribute to the reduced observability of the system.

Through the above analysis, we find that high-precision angle measurement information can significantly improve the observability of the system. In the Taiji formation, each SC is equipped with high-precision angle sensors like CCD and QPD. To further investigate the impact of angle measurement accuracy on system observability, numerical simulation experiments were conducted using the nonlinear Kalman filter as described in Section 4. Additionally, a comparative evaluation of the state estimation accuracy was performed for two observation schemes.

4. Numerical Simulation

4.1. Simulation Initialization

Since the Taiji formation has not yet been launched, there is a lack of real data available for research. We employ a numerical simulation method to investigate the problem of autonomous state estimation. The initial state of the formation is computed using the Keplerian orbits, as illustrated in Table 2. Subsequently, we employ the fourth-order Runge–Kutta method to solve differential Equation (5) and simulate the true state of the SC from the initial time. For this simulation, we use a discrete step length of $1/3$ s. The true measurement values for case 1 and case 2 at each time instant k are computed using Equations (17) and (18), respectively. It is assumed that the measurement information has been preprocessed, including time synchronization. In the state estimation process, the measurement information employed comprises the true measurement value with added noise. The process noise $\mathit{w}\left(t\right)$ and the measurement noise $\mathit{v}\left(k\right)$ are assumed to adhere to Equation (7). Specific noise deviations and simulation parameters are provided in

Table 6. The simulation parameters.

Parameter	Value
process noise standard deviation	$1\times {10}^{-12}$
relative range noise standard deviation	1 m
relative range rate noise standard deviation	$1\times {10}^{-6}$ m/s
relative angle noise standard deviation	$1\times {10}^{-6}$ rad
Sun-pointing angle noise standard deviation	$1\times {10}^{-5}$ rad
radial velocity noise standard deviation	$1\times {10}^{-2}$ m/s
simulation time step	$1/3$ s
simulation duration	$2^{16}/3$ s
measurement frequency	once per $1/3$ s

:

The covariance matrix of the process noise, denoted as Q, is set as $Q=\mathrm{diag}\left([I_{1\times 18}\times {10}^{-12}]\right)$. Similarly, the measurement noise matrix, denoted as R, is set as $R=\mathrm{diag}\left([I_{1\times 3}·1,I_{1\times 3}·{10}^{-12},I_{1\times 9}·{10}^{-12},I_{1\times 9}·{10}^{-10},I_{1\times 3}·{10}^{-4}]\right)$. In general, the a priori estimation error of the spacecraft’s position state is assumed to be 10 km, while the a priori knowledge error of the velocity state is 0.1 m/s. However, to consider a worst-case scenario with less accurate prior state knowledge, we increase the uncertainty by a factor of 10. Consequently, the covariance matrix of the state uncertainty is set as $P_0=\mathrm{diag}\left([I_{1\times 9}·{10}^{10},I_{1\times 9}·1]\right)$.

We employ the cubature Kalman filter (CKF) [27,28] to estimate the state, leveraging its robustness against divergence in highly nonlinear systems. Considering numerical stability, we adopt the square root form of the filter to perform covariance analysis in the case of weakly observable systems. The algorithmic process of the square root cubature Kalman filter (SR-CKF) for a hybrid system, consisting of continuous dynamics equations and discrete observation equations, can be described as follows:

Given the initial state estimation ${\widehat{\mathit{x}}}_{0|0}$ and the covariance matrix of the uncertainty estimation ${\mathit{P}}_{0|0}$ obtained from prior knowledge, the initial square root of the error covariance ${\mathit{S}}_{0|0}$ can be computed using the Cholesky decomposition of ${\mathit{P}}_{0|0}$.

$$\begin{array}{c}\hfill {\mathit{P}}_{0|0}={\mathit{S}}_{0|0}{\mathit{S}}_{0|0}^T\end{array}$$

(28)

The time update process can be described as follows:

1.

Calculate the cubature points ${\mathit{X}}_{i,k-1\mid k-1}$ (where $i=1,2,\dots ,m$) using the square root factor ${\mathit{S}}_{k-1\mid k-1}$ and the common cubature points ${\mathit{\zeta }}_i$:

$${\mathit{X}}_{i,k-1\mid k-1}={\mathit{S}}_{k-1\mid k-1}{\mathit{\zeta }}_i+{\widehat{\mathit{x}}}_{k-1\mid k-1}$$

(29)

Here, m is twice the dimension $n_x$ of the state variable $\mathit{x}$, and ${\mathit{\zeta }}_i$ can be defined as

$${\mathit{\zeta }}_i=\left\{\begin{array}{ccc}\sqrt{\frac{m}{2}}·{\mathit{ϵ}}_i\hfill & \hfill & i=1,2\cdots n_x\hfill \\ -\sqrt{\frac{m}{2}}·{\mathit{ϵ}}_{i-n_x}\hfill & \hfill & i=n_x+1,n_x+2\cdots m\hfill \end{array}\right.$$

(30)

where ${\mathit{ϵ}}_i$ denotes the i-th column of the unit matrix ${\mathit{I}}_{{\mathit{n}}_{\mathit{x}}}$.

2.

Evaluate the propagated cubature points:

$${\mathit{X}}_{i,k\mid k-1}^{\ast }=\mathit{f}\left({\mathit{X}}_{i,k-1\mid k-1}\right)$$

(31)

3.

Estimate the predicted state:

$${\widehat{\mathit{x}}}_{k\mid k-1}=\frac{1}{m}\sum _{i=1}^m{\mathit{X}}_{i,k\mid k-1}^{\ast }$$

(32)

4.

Estimate the square root factor of the predicted error covariance:

$$\begin{array}{cc}\hfill {\mathit{S}}_{k\mid k-1}& =\mathit{qr}\left(\left[{\mathcal{X}}_{k|k-1}^{\ast }{\mathit{S}}_{Q,k-1}\right]\right)\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill {\mathit{Q}}_{k-1}& ={\mathit{S}}_{Q,k-1}{\mathit{S}}_{Q,k-1}^T\hfill \end{array}$$

(34)

Here, $\mathit{qr}$ denotes the QR decomposition of the matrix, ${\mathit{S}}_{Q,k-1}$ represents the square root factor of the process noise matrix ${\mathit{Q}}_{k-1}$ at time k, and ${\mathcal{X}}_{k|k-1}^{}$ is the centered weighted matrix, which can be calculated as

$$\begin{array}{c}\hfill {\mathcal{X}}_{k\mid k-1}^{\ast }=\frac{1}{\sqrt{m}}\left[{\mathit{X}}_{1,k\mid k-1}^{\ast }-{\widehat{\mathit{x}}}_{k\mid k-1}{\mathit{X}}_{2,k\mid k-1}^{\ast }-{\widehat{\mathit{x}}}_{k\mid k-1}\cdots {\mathit{X}}_{m,k\mid k-1}^{\ast }-{\widehat{\mathit{x}}}_{k\mid k-1}\right]\end{array}$$

(35)

The measurement update process is as follows:

1.

Calculate the cubature points again:

$${\mathit{X}}_{i,k\mid k-1}={\mathit{S}}_{k\mid k-1}{\mathit{\zeta }}_i+{\widehat{\mathit{x}}}_{k\mid k-1}$$

(36)

2.

Evaluate the propagated cubature points:

$${\mathit{Y}}_{i,k\mid k-1}=\mathit{h}\left({\mathit{X}}_{i,k\mid k-1}\right)$$

(37)

3.

Estimate the predicted measurement:

$${\widehat{\mathit{y}}}_{k\mid k-1}=\frac{1}{m}\sum _{i=1}^m{\mathit{Y}}_{i,k\mid k-1}^{\ast }$$

(38)

4.

Estimate the square root factor of the innovation covariance matrix:

$$\begin{array}{cc}\hfill {\mathit{S}}_{yy,k\mid k-1}& =\mathit{qr}\left(\left[{\mathcal{Y}}_{k|k-1}{\mathit{S}}_{R,k-1}\right]\right)\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill {\mathit{R}}_{k-1}& ={\mathit{S}}_{R,k-1}{\mathit{S}}_{R,k-1}^T\hfill \end{array}$$

(40)

Here, ${\mathit{S}}_{R,k-1}$ is a square root factor of the measurement noise matrix ${\mathit{R}}_{k-1}$ at time k, and ${\mathcal{Y}}_{k|k-1}$ is the centered weighted matrix, calculated as

$$\begin{array}{c}\hfill {\mathcal{Y}}_{k\mid k-1}=\frac{1}{\sqrt{m}}\left[{\mathit{Y}}_{1,k\mid k-1}-{\widehat{\mathit{y}}}_{k\mid k-1}{\mathit{Y}}_{2,k\mid k-1}-{\widehat{\mathit{y}}}_{k\mid k-1}\cdots {\mathit{Y}}_{m,k\mid k-1}-{\widehat{\mathit{y}}}_{k\mid k-1}\right]\end{array}$$

(41)

Additionally, the centered weighted matrix ${\mathcal{X}}_{k\mid k-1}$ can be calculated as

$$\begin{array}{c}\hfill {\mathcal{X}}_{k\mid k-1}=\frac{1}{\sqrt{m}}\left[{\mathit{X}}_{1,k\mid k-1}-{\widehat{\mathit{x}}}_{k\mid k-1}{\mathit{X}}_{2,k\mid k-1}-{\widehat{\mathit{x}}}_{k\mid k-1}\cdots {\mathit{X}}_{m,k\mid k-1}-{\widehat{\mathit{x}}}_{k\mid k-1}\right]\end{array}$$

(42)

5.

Calculate the cross-covariance matrix:

$${\mathit{P}}_{xy,k\mid k-1}={\mathcal{X}}_{k\mid k-1}{\mathcal{Y}}_{k\mid k-1}^T$$

(43)

6.

Estimate the Kalman gain:

$${\mathit{K}}_k=({\mathit{P}}_{xy,k\mid k-1}/{\mathit{S}}_{yy,k\mid k-1})/{\mathit{S}}_{yy,k\mid k-1}$$

(44)

7.

Update the state estimation:

$${\widehat{\mathit{x}}}_{k\mid k}={\widehat{\mathit{x}}}_{k\mid k-1}+{\mathit{K}}_k({\mathit{y}}_k-{\widehat{\mathit{y}}}_{k\mid k-1})$$

(45)

Then the square root factor of the corresponding error covariance can be estimated as

$${\mathit{S}}_{k\mid k}=\mathit{qr}\left([{\mathcal{X}}_{k\mid k-1}-{\mathit{K}}_k{\mathcal{Y}}_{k\mid k-1}{\mathit{K}}_k{\mathit{S}}_{R,k}]\right)$$

(46)

This state estimation algorithm is executed recursively for each sampling time, denoted as k. The simulation results are presented in the subsequent subsection, as demonstrated below.

4.2. Simulation Results

As depicted in Table 4, three levels of accuracy are considered for relative angle measurement. Initially, we conduct state estimation experiments on case 2, varying the accuracy of angle measurements, to validate the results obtained from the observability analysis of the system presented in Section 3. Subsequently, we compare the influence of case 1 and case 2 on the accuracy of state estimation to investigate the significance of the radial velocity sensor.

As illustrated in Figure 3, Figure 4 and Figure 5 and

Table 7. The state estimation errors of the spacecraft in case 2.

	Measurement Level	x (km)	y (km)	z (km)	$\dot{\mathit{x}}$ (m/s)	$\dot{\mathit{y}}$ (m/s)	$\dot{\mathit{z}}$ (m/s)
SC1	level 1	$61.089$	$202.740$	$100.296$	$0.077$	$-0.065$	$0.124$
	level 2	$12.657$	$55.100$	$9.132$	$4.197\times {10}^{-3}$	$0.987\times {10}^{-3}$	$-2.164\times {10}^{-3}$
	level 3	$1.987$	$1.415$	$-0.746$	$1.733\times {10}^{-3}$	$3.475\times {10}^{-3}$	$-7.778\times {10}^{-3}$
SC2	level 1	$41.785$	$212.048$	$113.512$	$0.076$	$-0.066$	$0.126$
	level 2	$12.650$	$55.104$	$9.136$	$4.248\times {10}^{-3}$	$1.174\times {10}^{-3}$	$-2.203\times {10}^{-3}$
	level 3	$1.988$	$1.416$	$-0.747$	$1.732\times {10}^{-3}$	$3.480\times {10}^{-3}$	$-7.767\times {10}^{-3}$
SC3	level 1	$39.591$	$209.159$	$88.955$	$0.073$	$-0.068$	$0.126$
	level 2	$12.647$	$55.088$	$9.130$	$4.200\times {10}^{-3}$	$0.818\times {10}^{-3}$	$-2.131\times {10}^{-3}$
	level 3	$1.985$	$1.414$	$-0.745$	$1.734\times {10}^{-3}$	$3.478\times {10}^{-3}$	$-7.768\times {10}^{-3}$

, we can see that the worst-case estimation errors for position and velocity states are approximately 2.0 km and 0.0078 m/s, respectively; under the highest relative angular measurement accuracy of 1 nrad, the state estimation outcomes for the three SCs in the Taiji formation exhibit similar trends. That may be due to the fact that the three spacecraft have similar orbital configurations and the same sensor configuration. We find that the estimation error of the position state of the SC gradually increases at level 1 but converges progressively at level 2 and level 3. The velocity state error remains stable across all three levels of measurement accuracy. This stability can be attributed to the direct availability of velocity measurement information from the radial velocity sensor, whereas the position state cannot be directly measured. Notably, as the relative angular measurement accuracy improves, the accuracy of SC state estimation also increases.

Considering the SC payload limitations and other practical considerations, it is possible that the SC in the Taiji formation may not be equipped with an extensive array of sensors specifically designed for solar measurement purposes, such as a spectrometer capable of measuring the spacecraft’s velocity with respect to Sun. To assess the impact of the radial velocity sensor on state estimation, we conduct identical simulation experiments in case 1. This allows for a direct comparison and evaluation of the effect of the radial velocity sensor on the accuracy of state estimation.

As illustrated in Figure 5 and Figure 6 and Table 7 and

Table 8. The state estimation errors of the spacecraft with the accuracy of level 3 in case 1.

	x (km)	y (km)	z (km)	$\dot{\mathit{x}}$ (m/s)	$\dot{\mathit{y}}$ (m/s)	$\dot{\mathit{z}}$ (m/s)
SC1	$1.290$	$-3.120$	$-0.885$	$-0.103$	$-0.143$	$-0.024$
SC2	$1.291$	$-3.118$	$-0.883$	$-0.101$	$-0.142$	$-0.027$
SC3	$1.289$	$-3.121$	$-0.883$	$-0.103$	$-0.146$	$-0.022$

, we can see that the worst-case estimation errors for position and velocity states are approximately 3.1 km and 0.14 m/s, respectively, under the highest relative angular measurement accuracy of 1 nrad. When comparing the state estimation results between case 1 and case 2 with the same relative angle measurement accuracy, there is little difference in the accuracy of position state estimation of the SC. However, there is a significant disparity in the estimation of velocity of the SC. The velocity estimation results in case 2 are approximately 18 times more precise than those in case 1. This is likely due to the inclusion of radial velocity measurements in case 2, which directly provides velocity information of the SC, resulting in significant improvement of velocity estimation. Further, the filter is more stable in case 2 compared with that in case 1.

5. Conclusions

This paper presents an analysis of the autonomous state estimation problem of the Taiji formation. To assess the importance of a radial velocity sensor and considering the sensor configuration of the formation, we propose two observation schemes: case 1 contains all measurements except the radial velocity measurement, and case 2 contains all measurements. Considering that there are high-precision optical relative angle sensors in the Taiji formation, we study the influence of three distinct angular measurement accuracies on state estimation precision. By the observability analysis, we find that the utilization of high-precision angle measurement sensors can significantly enhance the observability of the system. By numerical simulation experiments, we find that, in case 2 with a relative angle measurement accuracy of 1 nrad, even with poor prior state knowledge, in which the prior information of the precision of the position and velocity is respectively 100 km and 1 m/s, the position of the SC can be estimated with an accuracy of approximately 2.0 km, while the velocity can reach approximately 0.0078 m/s. In the absence of radial velocity measurement in case 1, with a relative angle measurement accuracy of 1 nrad, the estimated accuracy of the position is approximately 3.1 km. However, the estimation accuracy of velocity is approximately 0.14 m/s. This could be attributed to the absence of measured velocity state information of the SC. It is noteworthy that the Taiji formation may not carry the spectrometer to measure the radial velocity of the SC relative to Sun. For certain missions, such as the alignment of laser between two distant SCs, of the Taiji formation, the sensor configuration of case 1 may also meet the requirements for SC state estimation accuracy, especially when the SC is equipped with a high-precision relative angle measurement sensor. The study proposed in this paper can be applied to address the autonomous state estimation problem when ground measurement information is unavailable and can assist in selecting the appropriate sensor configuration for the state estimation of the Taiji formation.