Angle-Only Cooperative Orbit Determination Considering Attitude Uncertainty

Abstract

In this paper, a novel concept for cooperative orbit determination (OD) using inter-spacecraft angle-only measurements is proposed. Different from the conventional cooperative OD that only estimates orbit states, the attitude of the observer spacecraft is considered by incorporating the attitude into the estimated vector. The observability of a two-spacecraft system is analyzed based on the observability matrix. Observability analysis reveals that inter-spacecraft angle-only measurements are inadequate to estimate both the attitude and the orbit states in two-body dynamics. The observability of the two-spacecraft system can be improved by considering high-order gravitational perturbation or executing an attitude maneuver on the observer spacecraft. This is the first time that we present the observability analysis and orbit estimation results for a two-spacecraft system considering attitude uncertainty for the observer. Finally, simulation results demonstrate the effectiveness of the proposed method. The results in this paper can be potentially useful for autonomous managements of a spacecraft constellation and formation.

1. Introduction

Autonomous orbit determination (OD), due to its considerable value in space systems engineering, has been of widespread interest in the last several decades [1,2,3,4]. The ability of spacecraft to determine their own states, without the help of ground-based tracking equipment, can improve their ‘intelligence’ and survivability and may also reduce operational management costs [5,6,7,8]. With the current plans and development of future multispacecraft constellations, the latter issue, i.e., efficient constellation management, becomes particularly important, and autonomous OD for a spacecraft constellation is highly desirable [9,10,11].

Filter algorithms have been widely used in many OD applications since R.E. Kalman proposed his famous recursive method Kalman filter (KF) to solve discrete linear filtering problems in 1960 [12,13,14]. Those filter algorithms could utilize the provided measurements coming from different sensors to obtain the required estimates of the orbit states [15,16,17,18]. Thus, current filter algorithms rely on the accurate measurement information measured by different sensors equipped on the spacecraft. In this situation, many filter-based methods have been proposed to solve the autonomous OD problem for a single spacecraft using GPS, radar, guidepost and magnetic field vector measurements [19,20,21,22]. Although the above methods have been proven to be efficient in providing considerable estimation, these measurements rely on complicated measurement equipment, which is impossible for constellation spacecraft.

Recently, a variety of methods have been proposed to solve the autonomous OD problem for a satellite constellation by using only interspacecraft relative measurements [23,24,25,26]. However, they all face the rank deficiency, meaning that the multispacecraft system with interspacecraft measurements alone is unobservable [27]. Interspacecraft ranging systems are widely used to obtain range-only measurements between two spacecraft, but they lack observability in some specific situations [23]. It has been concluded that in the multibody system, due to the nonsymmetric gravity from the noncentral gravitational body, all the orbit elements of two spacecraft can be obtained using only intersatellite range measurements [28,29]. However, intersatellite range measurements are not adequate to obtain all orbit elements in two-body dynamics because the rotation of the orbits with respect to the inertial reference system cannot be determined, and only the shape, size, and relative orientation of two orbits can be obtained using only intersatellite range measurements.

Another well-known method employs onboard cameras to obtain the relative angle measurements, which only include the relative direction information [24,25]. The angle-only navigation is simple, robust, and well proven in many applications [30]. However, the inertial states of the spacecraft cannot be observed without additional measurements. At least one beacon (usually noted as attitude-reference spacecraft) must be added to the system to ensure that there is a reference in inertial space or that other sources of measurements should be added [31]. The basic requirement for orbit determination using angle-only measurements is that the line-of-sight vector must be measured with respect to the inertial frame, which requires several types of sensors, such as a laser sensor to measure the relative range, an optical camera to measure the relative direction, and a star sensor to measure the inertial attitude [24,25,32]. However, due to the size, mass, and power consumption problem of active sensors such as star sensors, they might be unable to be on-boarded for small and microspacecraft to estimate the inertial attitude. Therefore, the inertial attitude of observing spacecraft should also be added into the OD problem when using intersatellite line-of-sight angle measurements. Fortunately, the conclusion of OD using angle-only measurements in two-body dynamics does not hold in other dynamical cases. Specifically, the obtainable orbit elements deviate from the prior conclusion under a more complex dynamic. For example, perturbations affect the orbits and contribute to deviations of real orbits from ideal two-body dynamics, which may improve the observability of a system. Moreover, with a calibrated thrust maneuver, observability can be guaranteed for possible situations [33].

Observability analysis is necessary to determine whether a particular measurement system is sufficient to solve the state estimation problem [34,35,36], and a number of researchers have carried out this analysis in the context of the OD problem [37,38,39]. One unwieldy approach to establish observability is to use classical nonlinear techniques that rely on Lie-derivative criterion methods [24,25]. In general, it is difficult to provide further results on the system observability based on the Lie derivative criterion because the analytical expressions of higher-order Lie derivatives become progressively more complicated. Hence, discussing observability for nonlinear systems using Lie algebra and differential geometry methods becomes quite difficult. Another type of approach for establishing observability is the observability matrix (OM) [1,11]. The OM method can be used to determine whether a system can be observed using sequential measurements, which are easy to obtain through numerical propagation.

In this paper, a novel concept for autonomous OD is proposed based on a two-spacecraft system with angle-only measurements. Differently from the traditional autonomous OD using intersatellite angle measurements, the inertial attitude of the observer is also added to be estimated. In this manuscript, we present for the first time the observability analysis and the orbit estimation results for a two-spacecraft system with different additional conditions. Five conditions are considered and for each condition, the observability of the OD system is investigated using the OM.

The remainder of this paper is organized as follows. The state model and the measurement model used is firstly presented in Section 2, followed by brief reviews of the OM and the fifth-order CKF. Five conditions, as well as the corresponding OMs, are discussed in Section 3. The observability analysis and filter estimation results are shown in Section 4. Finally, conclusions are given in Section 5.

2. Attitude and Orbit Determination Model

In this section, the basic mathematical model, namely, the state model and the observation model, for the autonomous attitude and orbit determination problem of a two-spacecraft system using angle-only measurements is first presented. In addition, the observability matrix is introduced to analyze the observability of the system. Finally, the widely used fifth-order cubature Kalman filter is reviewed for later estimation.

2.1. State Model

Two spacecraft, defined as $S_1$ and $S_2$, are considered and shown in Figure 1. Assume that $S_1$ is able to observe the line-of-sight angle between $S_1$ and $S_2$. Note that $S_1$ is not considered as the attitude-reference spacecraft, which is different from recent works.

The attitude-orbit determination developed in this paper aims to determine the absolute states (i.e., the absolute orbit of both $S_1$ and $S_2$, together with the absolute attitude of $S_2$) of the spacecraft system, which for Earth orbiting bodies are usually described in the Earth-centered inertial (ECI) frame (coordinate system $O_e-X_e{Y}_e{Z}_e$ in Figure 1). The analysis developed in this paper builds primarily on the two-body dynamics of the spacecraft with Earth as the primary body. Define the state of the i-th spacecraft in the ECI frame as:

(1)${\mathit{X}}_i={\left[{\mathit{r}}_i^{\mathrm{T}},{\mathit{v}}_i^{\mathrm{T}}\right]}^{\mathrm{T}}={\left[x_i,y_i,z_i,{\dot{x}}_i,{\dot{y}}_i,{\dot{z}}_i\right]}^{\mathrm{T}}$

The state equation of the spacecraft orbit under the two-body dynamics can be expressed in general form as:

(2)${\dot{\mathit{X}}}_i=\mathit{f}\left({\mathit{X}}_i\right)=\left[\begin{array}{c}{\dot{\mathit{r}}}_i\\ {\ddot{\mathit{r}}}_i\end{array}\right]$

In the case of ideal two-body dynamics (i.e., particle dynamics model), Equation (2) is given by:

(3)${\ddot{\mathit{r}}}_i=-\frac{{\mu }_e}{∥{\mathit{r}}_i∥}{\mathit{r}}_i$

In the case of considering $J_2$, $J_3$, and $J_4$ perturbations, Equation (3) is rewritten as:

(4)${\ddot{\mathit{r}}}_i=-\frac{{\mu }_e}{∥{\mathit{r}}_i∥}{\mathit{r}}_i+{\mathit{a}}_{J_2}+{\mathit{a}}_{J_3}+{\mathit{a}}_{J_4}$

where ${\mu }_e$ is the Earth’s gravitational parameter, and ${\mathit{a}}_{J_2}$, ${\mathit{a}}_{J_3}$ and ${\mathit{a}}_{J_4}$ are the perturbation acceleration of $J_2$, $J_3$, and $J_4$, which can be obtained by derivative of potential function with respect to position. The potential function is given by:

(5)$U=\frac{\mu }{r}\left[1-J_2{\left(\frac{R_e}{r}\right)}^2{P}_2(sin\phi )-J_3{\left(\frac{R_e}{r}\right)}^3{P}_3(sin\phi )-J_4{\left(\frac{R_e}{r}\right)}^4{P}_4(sin\phi )\right]$

where $\phi$ is the latitude of the spacecraft’s ground trace; $P_2(sin\phi )=\frac{3}{2}{(sin\phi )}^2-\frac{1}{2}$, $P_3(sin\phi )=\frac{5}{2}{(sin\phi )}^3-\frac{1}{2}(sin\phi )$ and $P_4(sin\phi )=\frac{35}{8}{(sin\phi )}^4-\frac{30}{8}{(sin\phi )}^2+\frac{3}{8}$.

The attitude of spacecraft $S_1$ is represented through the quaternion, defined by:

(6)$\mathit{q}=\left[\begin{array}{c}\hfill q_{\mathbf{0}}\hfill \\ \hfill q_{13}\hfill \end{array}\right]$

where:

(7)$q_{13}=\left[\begin{array}{c}\hfill q_1\hfill \\ \hfill q_2\hfill \\ \hfill q_3\hfill \end{array}\right]=\widehat{\mathit{n}}sin(\alpha /2)=\widehat{\mathsf{n}}sin(\alpha /2)$

(8)$q_0=cos(\alpha /2)$

where $\widehat{\mathit{n}}$ is a unit vector corresponding to the principal axis of rotation and $\alpha$ is the angle of rotation. The quaternion kinematics are derived through the spacecraft’s angular velocity as follows:

(9)$\dot{\mathit{q}}={\mathit{g}}_1\left(\mathit{q}\right)=\frac{1}{2}\mathsf{\Omega }\left(\omega \right)\mathit{q}$

(10)$\dot{\mathit{\omega }}={\mathit{g}}_2\left(\mathit{\omega }\right)={\mathit{J}}^{-1}\left({\mathit{\omega }}^{\times }\mathit{J}\mathit{\omega }+\mathit{M}\right)$

where $\mathit{\omega }={\left[{\omega }_x{\omega }_y{\omega }_z\right]}^{\mathrm{T}}$ denotes the angular velocity, $\mathit{J}$ represents the moment of inertia of the spacecraft, $\mathit{M}$ is the sum of external moments, and matrix $\mathsf{\Omega }\left(\mathit{\omega }\right)$ is defined as:

(11)$\mathsf{\Omega }\left(\mathit{\omega }\right)=\left[\begin{array}{cc}0& -{\mathit{\omega }}^T\\ \mathit{\omega }& -{\mathit{\omega }}^{\times }\end{array}\right],{\omega }^{\times }=\left[\begin{array}{ccc}0& -{\omega }_z& {\omega }_y\\ {\omega }_z& 0& -{\omega }_x\\ -{\omega }_y& {\omega }_x& 0\end{array}\right]$

The four quaternion elements satisfy the following normalization constraint:

(12)${\mathit{q}}^{\mathrm{T}}\mathit{q}={\mathit{q}}_{13}^{\mathrm{T}}{\mathit{q}}_{13}+q_0^2=1$

Combining the elements for the above orbits and attitude, the state vector to be estimated is:

(13)$\begin{array}{c}\hfill \mathit{X}={\left[{\mathit{X}}_1^{\mathrm{T}},{\mathit{X}}_2^{\mathrm{T}},{\mathit{q}}^{\mathrm{T}},{\omega }^{\mathrm{T}}\right]}^{\mathrm{T}}\hfill \\ \hfill ={\left[x_1,y_1,z_1,{\dot{x}}_1,{\dot{y}}_1,{\dot{z}}_1,x_2,y_2,z_2,{\dot{x}}_2,{\dot{y}}_2,{\dot{z}}_2,q_0,q_1,q_2,q_3,{\omega }_x,{\omega }_y,{\omega }_x\right]}^{\mathrm{T}}\hfill \end{array}$

The state model is given by:

(14)$\mathit{X}=\mathit{F}\left(\mathit{X}\right)={\left[{\dot{\mathit{X}}}_1^{\mathrm{T}},{\dot{{\mathit{X}}_2}}^{\mathrm{T}},{\dot{\mathit{q}}}^{\mathrm{T}},{\dot{\mathit{\omega }}}^{\mathrm{T}}\right]}^{\mathrm{T}}={\left[\mathit{f}{\left({\mathit{X}}_1\right)}^{\mathrm{T}},\mathit{f}{\left({\mathit{X}}_2\right)}^{\mathrm{T}},{\mathit{g}}_1{\left(\dot{\mathit{q}}\right)}^{\mathrm{T}},{\mathit{g}}_2\left(\mathit{\omega }\right){}^{\mathrm{T}}\right]}^{\mathrm{T}}$

2.2. Observation Model

As shown in Figure 1, the spacecraft $S_2$ is assumed to be observed by the spacecraft $S_1$. The inertial inter-spacecraft angle measurements can be obtained if the observer $S_1$ is an attitude-reference spacecraft. The inertial inter-spacecraft angle measurements can be represented by two angulars $\alpha$ and $\beta$, given as:

(15)$\left\{\begin{array}{c}\hfill \alpha ={tan}^{-1}(\Delta y/\Delta x)+{\epsilon }_{\alpha }\hfill \\ \hfill \beta ={sin}^{-1}\left(\Delta z/\sqrt{\Delta x^2+\Delta y^2+\Delta z^2}\right)+{\epsilon }_{\beta }\hfill \end{array}\right.$

where $\Delta x=x_2-x_1$, $\Delta y=y_2-y_1$ and $\Delta z=z_2-z_1$; ${\epsilon }_{\alpha }$ and ${\epsilon }_{\beta }$, respectively, denote the randomly distributed noise for the two angulars. To be convenient, the angular measurement equations in Equation (15) can be further replaced by the line-of-sight model, given by:

(16)${\mathit{y}}^{\prime }=\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-\mathit{r}∥}+{\epsilon }^{\prime }$

where ${\mathit{r}}_1={[x_1,y_1,z_1]}^T$ and ${\mathit{r}}_2={[x_2,y_2,z_2]}^T$; ${\mathbf{\epsilon }}^{\prime }$ denotes the randomly distributed noise for vector measurement ${\mathit{y}}^{\prime }$.

Note that in this paper, the attitude of spacecraft $S_1$ also needs to be estimated. Therefore, $S_1$ can only measure the line-of-sight angle in the spacecraft body coordinate system (i.e., coordinate system $O_B-X_B{Y}_B{Z}_B$), and the real measurement $\mathit{y}$ is given by the following:

(17)$\mathit{y}={\mathbf{R}}_E^B{\mathit{y}}^{\prime }+\epsilon ={\mathbf{R}}_E^B\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥}+\epsilon$

where $\epsilon$ is the corresponding measurement noise vector characterized by a normal distribution with zero mean and covariance $\mathbf{R}\in {\mathbb{R}}^{3\times 3}$ and ${\mathbf{R}}_E^B$ is the rotation matrix from ECI (coordinate system $O_E-X_E{Y}_E{Z}_E$) to the spacecraft body coordinate system $O_B-X_B{Y}_B{Z}_B$:

(18)${\mathbf{R}}_E^B=\left[\begin{array}{ccc}\hfill q_0^2+q_1^2-q_2^2-q_3^2\hfill & \hfill 2\left(q_0{q}_3+q_1{q}_2\right)\hfill & \hfill 2\left(q_1{q}_3-q_0{q}_2\right)\hfill \\ \hfill 2\left(-q_0{q}_3+q_1{q}_2\right)\hfill & \hfill q_1^2-q_1^2+q_2^2-q_3^2\hfill & \hfill 2\left(q_2{q}_3+q_{q_1}{q}_1\right)\hfill \\ \hfill 2\left(q_1{q}_2+q_1{q}_3\right)\hfill & \hfill 2\left(q_2{q}_3-q_0{q}_1\right)\hfill & \hfill q_0^2-q_1^2-q_2^2+q_3^2\hfill \end{array}\right]$

Hence, the navigation problem is given by:

(19)$\left\{\begin{array}{c}\hfill \dot{\mathit{X}}={\left[{\dot{\mathit{X}}}_1^{\mathrm{T}},{\dot{\mathit{X}}}_2^{\mathrm{T}},{\dot{\mathit{q}}}^{\mathrm{T}}\right]}^{\mathrm{T}}={\left[\mathit{f}{\left({\mathit{X}}_1\right)}^{\mathrm{T}},\mathit{f}{\left({\mathit{X}}_1\right)}^{\mathrm{T}},\mathit{g}{\left(\mathit{q}\right)}^{\mathrm{T}}\right]}^{\mathrm{T}}\hfill \\ \hfill \mathit{y}={\mathbf{R}}_E^B\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥}+\epsilon \hfill \end{array}\right.$

2.3. Observability Matrix

In this paper, the observability matrix (OM) is taken as a metric to evaluate the feasibility of the two-spacecraft system. With measurements collected k times sequentially, denoted as $\{{\mathit{y}}_0,{\mathit{y}}_1,\cdots ,{\mathit{y}}_{k-1}\}$ from time epoch $t_0$ to $t_{k-1}$, the OM is represented as:

(20)$\mathbf{N}=\left[\begin{array}{c}{\tilde{\mathbf{H}}}_0\\ ⋮\\ {\tilde{\mathbf{H}}}_{k-1}\end{array}\right]$

(21)${\tilde{\mathbf{H}}}_k={\mathbf{H}}_k\mathsf{\Phi }\left(t_k,t_0\right)$

where $\mathsf{\Phi }(t_k,t_0)$ is the state transformation matrix (STM) from $t_0$ to $t_k$ and

(22)${\mathbf{H}}_k={\left.\frac{\partial \mathit{y}}{\partial \mathit{X}}\right|}_{t=t_k}$

The differential equation of the STM is as follows:

(23)$\dot{\mathsf{\Phi }}\left(t_k,t_0\right)={\left.\frac{\partial \mathit{F}\left(\mathit{X}\right)}{\partial \mathit{X}}\right|}_{t=t_k}\mathsf{\Phi }\left(t_k,t_0\right)$

and is initialized by:

(24)$\mathsf{\Phi }\left(t_0,t_0\right)={\mathbf{I}}_{n_{\times }n}$

where n is the dimension of the state vector to be estimated and ${\mathbf{I}}_{n\times n}$ is an n-dimensional unit matrix. Note that the differential term $\partial \mathit{F}\left(\mathit{X}\right)/\partial \mathit{X}$ is also known as the Jacobi matrix of state model (14).

An OM with a full rank (i.e., $\mathrm{rank}\left(\mathbf{N}\right)=n$) indicating that the two-spacecraft system is observable using the given measurements [1]. Moreover, the observability degree of the two-spacecraft system can be described by the condition number (CN) of OM, represented by $\mathrm{cond}\left(\mathbf{N}\right)=∥\mathbf{N}∥·∥{\mathbf{N}}^{-1}∥$. The smaller the CN, the better the observability [24].

2.4. Review of Fifth-Order CKF

In this section, the fifth-order CKF algorithm is briefly summarized. First, consider a discrete nonlinear system as:

(25)$\begin{array}{c}\hfill {\mathit{x}}_k=\mathit{f}\left({\mathit{x}}_{k-1}\right)+{\mathit{w}}_{k-1}\hfill \\ \hfill {\mathit{z}}_k=\mathit{h}\left({\mathit{x}}_k\right)+{\mathit{v}}_k\hfill \end{array}$

where ${\mathit{x}}_k\in {\mathbb{R}}^n$ is the state vector at time epoch k and ${\mathit{z}}_k\in {\mathbb{R}}^m$ is the measurement. ${\mathit{w}}_{k-1}\in {\mathbb{R}}^n$ and ${\mathit{v}}_k\in {\mathbb{R}}^m$ denote the independent system and measurement noise, respectively, and are both considered independent and white Gaussian distributions with covariances ${\mathbf{Q}}_{k-1}$ and ${\mathbf{R}}_k$, respectively.

The optimal Bayesian filters contains two steps: the prediction step and the update step. Both of the two steps require us to calculate the Gaussian weighted integration ${\int }_{{\mathbb{R}}^n}\mathit{g}\left(\mathit{x}\right)\mathit{N}(\mathit{x};\overline{\mathit{x}},\mathbf{P})d\mathit{x}$, where $\mathit{g}\left(\mathit{x}\right)$ is a nonlinear function. The integral with respect to the general Gaussian distribution $\mathit{N}(\mathit{x};\overline{\mathit{x}},\mathbf{P})$ can be further approximated by:

(26)$\begin{array}{cc}\hfill {\int }_{{\mathbb{R}}^n}\mathit{g}\left(x\right)\mathit{N}(\mathit{x};\overline{\mathit{x}},\mathbf{P})d\mathit{x}& \hfill ={\int }_{{\mathbb{R}}^n}g(\mathbf{A}\mathit{x}+\overline{\mathit{x}})\mathit{N}(\mathit{x};\overline{\mathit{x}},\mathbf{P})d\mathit{x}\hfill \\ \hfill & \hfill \approx \sum _{i=1}^{N_p}{W}_i\mathit{g}\left(\mathit{A}{\gamma }_i+\overline{\mathit{x}}\right)\hfill \end{array}$

where $\mathbf{P}=\mathbf{S}{\mathbf{S}}^{\mathrm{T}}$, $N_P$ is the total number of points, and ${\gamma }_i$ and $W_i$ are the quadrature points and weights, respectively, corresponding to the Gaussian distribution $\mathit{N}(\mathit{x};\overline{\mathit{x}},\mathbf{P})$. Specifically, the integral with respect to $\mathit{N}(\mathit{x};\mathbf{0},\mathbf{I})$ can be approximated by the following quadrature rule:

(27)${\int }_{{\mathbb{R}}^n}\mathit{g}\left(\mathit{x}\right)N(\mathit{x};\mathbf{0},\mathbf{I})d\mathit{x}\approx \sum _{i=1}^{N_p}{W}_i\mathit{g}\left({\gamma }_i\right)$

According to the cubature rule and Equation (26), in fifth-order CKF, the cubature points ${\mathbf{\zeta }}_i$ are given by [21]:

(28)$\left\{\begin{array}{c}{\mathbf{\zeta }}_0={\mathbf{0}}_{n\times 1}\hfill \\ {\mathbf{\zeta }}_i={\mathbf{A}}_k\cdot \sqrt{n+2}{\mathit{e}}_i\hfill \\ {\mathbf{\zeta }}_{i+n}=-{\mathbf{A}}_k\cdot \sqrt{n+2}{\mathit{e}}_i\hfill \\ {\mathbf{\zeta }}_{j+2n}={\mathbf{A}}_k\cdot \sqrt{n+2}{\mathsf{s}}_j^{+}\hfill \\ {\mathbf{\zeta }}_{j+2n+n\left(n-1\right)/2}=-{\mathbf{A}}_k\cdot \sqrt{n+2}{\mathit{s}}_j^{+}\hfill \\ {\mathbf{\zeta }}_{j+2n+n\left(n-1\right)}={\mathbf{A}}_k\cdot \sqrt{n+2}{\mathit{s}}_j^{-}\hfill \\ {\mathbf{\zeta }}_{j+2n+3n\left(n-1\right)/2}=-{\mathbf{A}}_k\cdot \sqrt{n+2}{\mathit{s}}_j^{-}\hfill \\ i=1,2,\dots ,n;j=1,2,\dots ,n\left(n-1\right)/2\hfill \end{array}\right.$

where n is the state dimension of the system to be estimated, ${\mathbf{A}}_k$ is the Cholesky decomposition of covariance matrix ${\mathbf{P}}_k$ at epoch $t_k$ and $P_k=A_k{A}_k^T$, and ${\mathit{e}}_i$ is the i-th column of the n-th identity matrix ${\mathbf{I}}_n$. The point sets ${\mathit{s}}_j^{+}$ and ${\mathit{s}}_j^{-}$ are given as follows.

(29)$\begin{array}{cc}\hfill & \hfill \left\{{\mathit{s}}_j^{+}\right\}=\left\{({\mathit{e}}_p+{\mathit{e}}_q)/\sqrt{2}\left|p<q,p,q=1,2,\dots ,n\right.\right\}\hfill \\ \hfill & \hfill \left\{{\mathit{s}}_j^{-}\right\}=\left\{({\mathit{e}}_p-{\mathit{e}}_q)/\sqrt{2}\left|p<q,p,q=1,2,\dots ,n\right.\right\}\hfill \end{array}$

The corresponding weight $w_i$ of each Cubature point ${\zeta }_i$ is given by

(30)$w_i=\left\{\begin{array}{cc}\hfill \frac{2}{n+2}\hfill & \hfill i=0\hfill \\ \hfill \frac{1}{{(n+2)}^2}\hfill & \hfill i=1,\dots ,2n(n-1)\hfill \\ \hfill \frac{4-n}{{(n+2)}^2}\hfill & \hfill i=2n(n-1)+1,\dots ,2n^2\hfill \end{array}\right.$

Then, the fifth-order CKF algorithm is summarized as follows:

i. Time updating:

(31)$\begin{array}{cc}\hfill & \hfill {\mathrm{A}}_k=Cholesky\left({\mathbf{P}}_{k-1|k-1}\right)\hfill \\ \hfill & \hfill {\chi }_{i,k-1|k-1}={\mathbf{\zeta }}_i+{\mathit{x}}_{k-1|k-1}\hfill \\ \hfill & \hfill {\chi }_{i,k|k-1}^{*}=\mathit{f}\left({\chi }_{i,k-1|k-1}\right)\hfill \\ \hfill & \hfill {\mathit{x}}_{k|k-1}=\sum _{i=0}^{2n^2}{w}_i{\chi }_{i,k|k-1}^{*}\hfill \\ \hfill & \hfill {\mathbf{P}}_{k|k-1}=\sum _{i=0}^{2n^2}{w}_i({\chi }_{i,k|k-1}^{*}-{\mathit{x}}_{k|k-1}){({\chi }_{i,k|k-1}^{*}-{\mathit{x}}_{k|k-1})}^{\mathrm{T}}+{\mathbf{Q}}_k\hfill \end{array}$

ii. Measurement updating

(32)$\begin{array}{cc}\hfill & \hfill {\mathbf{A}}_k=Cholesky\left({\mathbf{P}}_{k|k-1}\right)\hfill \\ \hfill & \hfill {\chi }_{i,k|k-1}={\mathbf{\zeta }}_i+{\mathit{x}}_{k|k-1}\hfill \\ \hfill & \hfill {\mathit{z}}_{i,k|k-1}^{*}=\mathit{h}\left({\chi }_{i,k|k-1}\right)\hfill \\ \hfill & \hfill {\mathit{z}}_{k|k-1}=\sum _{i=0}^{2n^2}{w}_i{\mathit{z}}_{i,k|k-1}^{*}\hfill \\ \hfill & \hfill {\mathbf{P}}_{\mathit{xz},k|k-1}=\sum _{i=0}^{2n^2}{w}_i({\chi }_{i,k|k-1}^{*}-{\mathit{x}}_{k|k-1}){({\mathit{z}}_{i,k|k-1}^{*}-{\mathit{z}}_{k|k-1})}^{\mathrm{T}}\hfill \\ \hfill & \hfill {\mathbf{P}}_{\mathit{zz},k|k-1}=\sum _{i=0}^{2n^2}{w}_i({\mathit{z}}_{i,k|k-1}^{*}-{\mathit{z}}_{k|k-1}){({\mathit{z}}_{i,k|k-1}^{*}-{\mathit{z}}_{k|k-1})}^{\mathrm{T}}+{\mathbf{R}}_k\hfill \\ \hfill & \hfill {\mathbf{K}}_k={\mathbf{P}}_{\mathit{xz},k|k-1}{\mathbf{P}}_{\mathit{zz},k|k-1}^{-1}\hfill \\ \hfill & \hfill {\mathit{x}}_{k|k}={\mathit{x}}_{k|k-1}+{\mathbf{K}}_k({\mathit{z}}_k-{\mathit{z}}_{k|k-1})\hfill \\ \hfill & \hfill {\mathbf{P}}_{k|k}={\mathbf{P}}_{k|k-1}-{\mathbf{K}}_k{\mathbf{P}}_{\mathit{zz},k|k-1}{\mathbf{K}}_k^{\mathrm{T}}\hfill \end{array}$

where Cholesky(·) represents the Cholesky decomposition method, ${\chi }_{i,k|k-1}^{*}$ is the Cubature point generated from states and ${\mathit{z}}_{i,k|k-1}^{*}$ represents the Cubature point generated from measurements.

3. Attitude and Orbit Determination Method with Angle-Only Measurements

In this section, the autonomous attitude and orbit determination problem, with different conditions, are modeled based on the theory in Section 2.

3.1. Case I: Both Orbits of $S_1$ and $S_2$ Are Known, and the Attitude of $S_1$ Is Unchanged

In this case, the navigation problem (19) is simplified as follows: the orbits of both $S_1$ and $S_2$ are known, and the attitude of $S_1$ is unknown but unchanged (i.e., $\mathbf{\omega }={\left[000\right]}^{\mathrm{T}}$). Thus, the state vector can be written as:

(33)$\mathit{X}=\mathit{q}={[q_0,q_1,q_2,q_3]}^{\mathrm{T}}$

The state model is given by:

(34)$\dot{\mathit{X}}=\mathit{F}\left(\mathit{X}\right)=\dot{\mathit{q}}={\mathit{g}}_1\left(\mathit{q}\right)=\frac{1}{2}\mathsf{\Omega }\left(\mathbf{\omega }\right)\mathit{q}=\frac{1}{2}\left[\begin{array}{cc}0& -{\mathbf{\omega }}^{\mathrm{T}}\\ \mathbf{\omega }& -{\mathbf{\omega }}^{\times }\end{array}\right]\mathit{q}={[0,0,0,0]}^{\mathrm{T}}$

The measurement is shown in Equation (17) and then the partial derivative of the intersatellite angle to the quaternion of spacecraft $S_1$${\mathbf{H}}_k$ is of the form:

(35)$\begin{array}{cc}\hfill & \hfill {\mathbf{H}}_k={\left.\frac{\partial \mathit{y}}{\partial \mathit{q}}\right|}_{t=t_k}={\left.\frac{\partial {\mathbf{R}}_E^B}{\partial \mathit{q}}\right|}_{t=t_k}\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥}\hfill \\ \hfill & \hfill =\left[{\left.\frac{\partial {\mathbf{R}}_E^B}{\partial q_0}\right|}_{t=t_k}{\left.\frac{\partial {\mathbf{R}}_E^B}{\partial q_1}\right|}_{t=t_k}{\left.\frac{\partial {\mathbf{R}}_E^B}{\partial q_2}\right|}_{t=t_k}{\left.\frac{\partial {\mathbf{R}}_E^B}{\partial q_3}\right|}_{t=t_k}\right]\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥}\hfill \end{array}$

For the state model in Equation (34), the Jacobi matrix is given as:

(36)$\frac{\partial \mathit{F}\left(\mathit{X}\right)}{\partial \mathit{X}}=\frac{\partial \mathit{F}\left(\mathit{q}\right)}{\partial \mathit{q}}={\mathbf{0}}_{4\times 4}$

Therefore, the STM has the following format:

(37)$\mathsf{\Phi }(t_k,t_0)={\mathbf{I}}_{4\times 4}$

3.2. Case II: The Orbit of $S_1$ Is Known, and the Attitude of $S_1$ Is Unchanged

In this case, the orbit of $S_1$ is known, and the attitude of $S_1$ is unknown but unchanged; hence, the state to be estimated is given by:

(38)$\dot{\mathit{X}}=\mathit{F}\left(\mathit{X}\right)={[{\dot{\mathit{X}}}_2^{\mathrm{T}},{\dot{\mathit{q}}}^{\mathrm{T}}]}^{\mathrm{T}}={[\mathit{f}{\left({\mathit{X}}_2\right)}^{\mathrm{T}},{\mathit{g}}_1{\left(\mathit{q}\right)}^{\mathrm{T}}]}^{\mathrm{T}}={[\mathit{f}{\left({\mathit{X}}_2\right)}^{\mathrm{T}},[0,0,0,0]]}^{\mathrm{T}}$

The corresponding state model is as follows:

(39)$\dot{\mathit{X}}=\mathit{F}\left(\mathit{X}\right)={[{\dot{\mathit{X}}}_2^{\mathrm{T}},{\dot{\mathit{q}}}^{\mathrm{T}}]}^{\mathrm{T}}={[\mathit{f}{\left({\mathit{X}}_2\right)}^{\mathrm{T}},{\mathit{g}}_1{\left(\mathit{q}\right)}^{\mathrm{T}}]}^{\mathrm{T}}={[\mathit{f}{\left({\mathit{X}}_2\right)}^{\mathrm{T}},[0,0,0,0]]}^{\mathrm{T}}$

The partial differential matrix ${\mathbf{H}}_k$ of measurement (17) is written as:

(40)$\begin{array}{cc}\hfill & \hfill {\mathbf{H}}_k={\left.\frac{\partial \mathit{y}}{\partial \mathit{X}}\right|}_{t=t_k}=\left[{\left.\frac{\partial \mathit{y}}{\partial {\mathit{X}}_2}\right|}_{t=t_k}{\left.\frac{\partial \mathit{y}}{\partial \mathit{q}}\right|}_{t=t_k}\right]=\left[{\left.\frac{\partial \mathit{y}}{\partial {\mathit{r}}_2}\right|}_{t=t_k}{\left.\frac{\partial \mathit{y}}{\partial {\mathit{v}}_2}\right|}_{t=t_k}{\left.\frac{\partial \mathit{y}}{\partial \mathit{q}}\right|}_{t=t_k}\right]\hfill \\ \hfill & \hfill =\left[{\mathbf{R}}_E^B{\left.\frac{\partial }{\partial {\mathit{r}}_2}\left(\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥}\right)\right|}_{t=t_k}{\mathbf{R}}_E^B{\left.\frac{\partial }{\partial {\mathit{v}}_2}\left(\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥}\right)\right|}_{t=t_k}{\left.\frac{\partial {\mathbf{R}}_E^B}{\partial \mathit{q}}\right|}_{t=t_k}\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥}\right]\hfill \end{array}$

The propagation of the STM for state model (39) is then given by:

(41)$\left\{\begin{array}{c}\hfill \mathsf{\Phi }(t_0,t_0)={\mathbf{I}}_{10\times 10}\hfill \\ \hfill \dot{\mathsf{\Phi }}(t_k,t_0)={\left.\frac{\partial \mathit{F}\left(\mathit{X}\right)}{\partial \mathit{X}}\right|}_{t=t_k}\mathsf{\Phi }(t_k,t_0)\hfill \end{array}\right.$

Hence, the STM $\mathsf{\Phi }(t_k,t_0)$ has the following format:

(42)$\mathsf{\Phi }(t_k,t_0)=\left[\begin{array}{ccc}{\mathsf{\Phi }}_{rr}(t_k,t_0)& {\mathsf{\Phi }}_{rv}(t_k,t_0)& {\mathbf{0}}_{3\times 4}\\ {\mathsf{\Phi }}_{vr}(t_k,t_0)& {\mathsf{\Phi }}_{vv}(t_k,t_0)& {\mathbf{0}}_{3\times 4}\\ {\mathbf{0}}_{4\times 3}& {\mathbf{0}}_{4\times 3}& {\mathbf{I}}_{4\times 4}\end{array}\right]$

3.3. Case III: Both the Orbits of $S_1$ and $S_2$ Are Unknown, and the Attitude of $S_1$ Is Unchanged

In this case, the orbits of $S_1$ and $S_2$ as well as the attitude of $S_1$ are estimated. The only information for the navigation system is that $\mathbf{\omega }={\left[000\right]}^{\mathrm{T}}$. Therefore, the state to be estimated is:

(43)$\mathit{X}={[{\mathit{X}}_1^{\mathrm{T}},{\mathit{X}}_2^{\mathrm{T}},{\mathit{q}}^{\mathrm{T}}]}^{\mathrm{T}}={[x_1,y_1,z_1,{\dot{x}}_1,{\dot{y}}_1,{\dot{z}}_1,x_2,y_2,z_2,{\dot{x}}_2,{\dot{y}}_2,{\dot{z}}_2,q_0,q_1,q_2,q_3]}^{\mathrm{T}}$

The state model is then given by:

(44)$\dot{\mathit{X}}=\mathit{F}\left(\mathit{X}\right)={[{\dot{\mathit{X}}}_1^{\mathrm{T}},{\dot{\mathit{X}}}_2^{\mathrm{T}},{\dot{\mathit{q}}}^{\mathrm{T}}]}^{\mathrm{T}}={[\mathit{f}{\left({\mathit{X}}_1\right)}^{\mathrm{T}},\mathit{f}{\left({\mathit{X}}_1\right)}^{\mathrm{T}},{\mathit{g}}_1{\left(\mathit{q}\right)}^{\mathrm{T}}]}^{\mathrm{T}}={[\mathit{f}{\left({\mathit{X}}_1\right)}^{\mathrm{T}},\mathit{f}{\left({\mathit{X}}_1\right)}^{\mathrm{T}},[0,0,0,0]]}^{\mathrm{T}}$

Similar to Equation (40), the partial differential matrix of measurement is given as:

(45)$\begin{array}{cc}\hfill & \hfill {\mathbf{H}}_k={\left.\frac{\partial \mathit{y}}{\partial \mathit{X}}\right|}_{t=t_k}=\left[{\left.\frac{\partial \mathit{y}}{\partial {\mathit{X}}_1}\right|}_{t=t_k}{\left.\frac{\partial \mathit{y}}{\partial {\mathit{X}}_2}\right|}_{t=t_k}{\left.\frac{\partial \mathit{y}}{\partial \mathit{q}}\right|}_{t=t_k}\right]=\left[{\left.\frac{\partial \mathit{y}}{\partial {\mathit{r}}_1}\right|}_{t=t_k}{\left.\frac{\partial \mathit{y}}{\partial {\mathit{v}}_1}\right|}_{t=t_k}{\left.\frac{\partial \mathit{y}}{\partial {\mathit{r}}_2}\right|}_{t=t_k}{\left.\frac{\partial \mathit{y}}{\partial {\mathit{v}}_2}\right|}_{t=t_k}{\left.\frac{\partial \mathit{y}}{\partial \mathit{q}}\right|}_{t=t_k}\right]\hfill \\ \hfill & \hfill ={\left.\left[{\mathbf{R}}_E^B\frac{\partial }{\partial {\mathit{r}}_1}\left(\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥}\right),{\mathbf{R}}_E^B\frac{\partial }{\partial {\mathit{v}}_1}\left(\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥}\right),{\mathbf{R}}_E^B\frac{\partial }{\partial {\mathit{r}}_2}\left(\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥}\right),{\mathbf{R}}_E^B\frac{\partial }{\partial {\mathit{v}}_2}\left(\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥}\right),\frac{\partial {\mathbf{R}}_E^B}{\partial \mathit{q}}\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥}\right]\right|}_{t=t_k}\hfill \end{array}$

where:

(46)$\frac{\partial }{\partial {\mathit{r}}_1}\left(\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥}\right)=-\frac{1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥}\left[{\mathbf{I}}_{3\times 3}-\frac{({\mathit{r}}_2-{\mathit{r}}_1){({\mathit{r}}_2-{\mathit{r}}_1)}^{\mathrm{T}}}{{∥{\mathit{r}}_2-{\mathit{r}}_1∥}^2}\right]$

(47)$\frac{\partial }{\partial {\mathit{v}}_1}\left(\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥}\right)={\mathbf{0}}_{3\times 3}$

The STM $\mathsf{\Phi }(t_k,t_0)$ in this case has the following format:

(48)$\mathsf{\Phi }(t_k,t_0)=\left[\begin{array}{ccccc}{\mathsf{\Phi }}_{r_1{r}_1}(t_k,t_0)& {\mathsf{\Phi }}_{r_1{v}_1}(t_k,t_0)& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 4}\\ {\mathsf{\Phi }}_{v_1{r}_1}(t_k,t_0)& {\mathsf{\Phi }}_{v_1{v}_1}(t_k,t_0)& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 4}\\ {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathsf{\Phi }}_{r_2{r}_2}(t_k,t_0)& {\mathsf{\Phi }}_{r_2{v}_2}(t_k,t_0)& {\mathbf{0}}_{3\times 4}\\ {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathsf{\Phi }}_{v_2{r}_2}(t_k,t_0)& {\mathsf{\Phi }}_{v_2{v}_2}(t_k,t_0)& {\mathbf{0}}_{3\times 4}\\ {\mathbf{0}}_{4\times 3}& {\mathbf{0}}_{4\times 3}& {\mathbf{0}}_{4\times 3}& {\mathbf{0}}_{4\times 3}& {\mathbf{I}}_{4\times 4}\end{array}\right]$

3.4. Case IV: Both the Orbits of $S_1$ and $S_2$ Are Unknown, and the Angular Velocity of $S_1$ Is Known

In this case, the attitude of $S_1$ is unknown and changed. For spacecraft $S_1$, the angular velocity $\mathbf{\omega }={\left[{\omega }_x{\omega }_y{\omega }_z\right]}^{\mathrm{T}}$ is foreknown (assume that ${\omega }_x$, ${\omega }_y$ and ${\omega }_z$ are constant and satisfy ${\omega }_x^2+{\omega }_y^2+{\omega }_z^2\ne 0$). Hence, the state vector to be estimated is the same as that of case III, with the form of Equation (43). The state model is written as:

(49)$\dot{\mathit{X}}=\mathit{F}\left(\mathit{X}\right)={[{\dot{\mathit{X}}}_1^{\mathrm{T}},{\dot{\mathit{X}}}_2^{\mathrm{T}},{\dot{\mathit{q}}}^{\mathrm{T}}]}^{\mathrm{T}}={[\mathit{f}{\left({\mathit{X}}_1\right)}^{\mathrm{T}},\mathit{f}{\left({\mathit{X}}_1\right)}^{\mathrm{T}},{\mathit{g}}_1{\left(\mathit{q}\right)}^{\mathrm{T}}]}^{\mathrm{T}}={[\mathit{f}{\left({\mathit{X}}_1\right)}^{\mathrm{T}},\mathit{f}{\left({\mathit{X}}_1\right)}^{\mathrm{T}},{\mathit{q}}^{\mathrm{T}}{\mathsf{\Omega }}^{\mathrm{T}}\left(\mathbf{\omega }\right)/2]}^{\mathrm{T}}$

The STM is in the form of:

(50)$\mathsf{\Phi }(t_k,t_0)=\left[\begin{array}{ccccc}{\mathsf{\Phi }}_{r_1{r}_1}(t_k,t_0)& {\mathsf{\Phi }}_{r_1{v}_1}(t_k,t_0)& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 4}\\ {\mathsf{\Phi }}_{v_1{r}_1}(t_k,t_0)& {\mathsf{\Phi }}_{v_1{v}_1}(t_k,t_0)& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 4}\\ {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathsf{\Phi }}_{r_2{r}_2}(t_k,t_0)& {\mathsf{\Phi }}_{r_2{v}_2}(t_k,t_0)& {\mathbf{0}}_{3\times 4}\\ {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathsf{\Phi }}_{v_2{r}_2}(t_k,t_0)& {\mathsf{\Phi }}_{v_2{v}_2}(t_k,t_0)& {\mathbf{0}}_{3\times 4}\\ {\mathbf{0}}_{4\times 3}& {\mathbf{0}}_{4\times 3}& {\mathbf{0}}_{4\times 3}& {\mathbf{0}}_{4\times 3}& e^{\frac{\mathsf{\Omega }\left(\omega \right)}{2}}\end{array}\right]$

3.5. Case V: Both the Orbits of $S_1$ and $S_2$ Are Unknown, and the Angular Velocity of $S_1$ Is Unchanged

In this case, the angular velocity $\mathbf{\omega }={\left[{\omega }_x{\omega }_y{\omega }_z\right]}^{\mathrm{T}}$ also exists, unknown but unchanged. Therefore, the state vector is obtained by Equation (13) and the state model is given as Equation (14).

The partial derivative of the intersatellite angle to the state vector $\mathit{X}$ is as follows:

(51)$\begin{array}{cc}\hfill {\mathbf{H}}_k=& \hfill {\left.\frac{\partial \mathit{y}}{\partial \mathit{X}}\right|}_{t=t_k}=\left[{\left.\frac{\partial \mathit{y}}{\partial {\mathit{X}}_1}\right|}_{t=t_k}{\left.\frac{\partial \mathit{y}}{\partial {\mathit{X}}_2}\right|}_{t=t_k}{\left.\frac{\partial \mathit{y}}{\partial \mathit{q}}\right|}_{t=t_k}{\left.\frac{\partial \mathit{y}}{\partial \mathbf{\omega }}\right|}_{t=t_k}\right]\hfill \\ \hfill =& \hfill \left[{\left.\frac{\partial \mathit{y}}{\partial {\mathit{r}}_1}\right|}_{t=t_k}{\left.\frac{\partial \mathit{y}}{\partial {\mathit{v}}_1}\right|}_{t=t_k}{\left.\frac{\partial \mathit{y}}{\partial {\mathit{r}}_2}\right|}_{t=t_k}{\left.\frac{\partial \mathit{y}}{\partial {\mathit{v}}_2}\right|}_{t=t_k}{\left.\frac{\partial \mathit{y}}{\partial \mathit{q}}\right|}_{t=t_k}{\left.\frac{\partial \mathit{y}}{\partial \mathbf{\omega }}\right|}_{t=t_k}\right]\hfill \\ \hfill =& \hfill \left[{\mathbf{R}}_E^B\frac{\partial }{\partial {\mathit{r}}_1}\left(\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥}\right),{\mathbf{R}}_E^B\frac{\partial }{\partial {\mathit{v}}_1}\left(\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥}\right),{\mathbf{R}}_E^B\frac{\partial }{\partial {\mathit{r}}_2}\left(\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥}\right)\right.,\hfill \\ & \hfill {\left.\left.{\mathbf{R}}_E^B\frac{\partial }{\partial {\mathit{v}}_2}\left(\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥}\right),\frac{\partial {\mathbf{R}}_E^B}{\partial \mathit{q}}\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥},\frac{\partial {\mathbf{R}}_E^B}{\partial \mathbf{\omega }}\frac{{\mathit{r}}_2-{\mathit{r}}_1}{∥{\mathit{r}}_2-{\mathit{r}}_1∥}\right]\right|}_{t=t_k}\hfill \end{array}$

4. Numerical Simulation

In this section, a series of numerical results for several types of scenarios is presented, with the following three objectives: (1) to demonstrate whether the system is observable or unobservable for each scenario and (2) to validate the observability indices by comparing their estimations to the quality of the solution of the state estimation problem using the fifth-order CKF.

4.1. Simulation Background

An example with two spacecraft in circle orbits is considered. The nominal orbit elements are listed in Table 1, and the corresponding orbits are shown in Figure 2. The elements h, e, i, $\mathsf{\Omega }$, $\omega$, and n denote the orbit altitude, eccentricity, inclination, longitude of the ascending node, and argument of the periapse and true anomaly, respectively. The semimajor axis a is computed by $a=h+R_e$, where $R_e=6378.137\mathrm{km}$ is the radius of the Earth.

The initial position and velocity errors (if unknown) are set to 10 km and 10${}^{-3}$ km/s, respectively. The initial attitude of spacecraft $S_1$ is set to be ${\mathit{q}}_{13}={[0,0,0]}^{\mathrm{T}}$ and $q_0=\sqrt{1-{\mathit{q}}_{13}^{\mathrm{T}}{\mathit{q}}_{13}}=1$. For all the cases above, the initial estimation of the quaternion is given by:

(52)$\widehat{\mathit{q}}={\left[cos\frac{\mathsf{\Theta }}{2}\frac{\sqrt{3}}{3}sin\frac{\mathsf{\Theta }}{2}\frac{\sqrt{3}}{3}sin\frac{\mathsf{\Theta }}{2}\frac{\sqrt{3}}{3}sin\frac{\mathsf{\Theta }}{2}\right]}^{\mathrm{T}}$

where $\mathsf{\Theta }=5^{\circ }$.

The initial angular velocity of spacecraft $S_1$, if existing (i.e., for case IV and case V), is set to be $\mathbf{\omega }={[0.01,0.01,0.01]}^{\mathrm{T}}{}^{\circ }/s$. The initial estimation of angular velocity is considered as $\widehat{\mathbf{\omega }}={[0.09,0.09,0.09]}^{\mathrm{T}}{}^{\circ }/s$ if $\mathbf{\omega }$ is to be estimated (only for case V).

The initial covariance matrix ${\mathbf{P}}_0$ is given as ${\mathbf{P}}_{{\mathit{r}}_1{\mathit{r}}_1,0}={\mathbf{P}}_{{\mathit{r}}_2{\mathit{r}}_2,0}=100{\mathbf{I}}_{3\times 3}\mathrm{k}{\mathrm{m}}^2$, ${\mathbf{P}}_{{\mathit{v}}_1{\mathit{v}}_1,0}={\mathbf{P}}_{{\mathit{v}}_2{\mathit{v}}_2,0}={\mathbf{I}}_{3\times 3}{m/s}^2$, ${\mathbf{P}}_{q_0,0}={(cos{2.5}^{\circ }-1)}^2{\mathbf{I}}_{1\times 1}$, ${\mathbf{P}}_{q_1,0}={\mathbf{P}}_{q_2,0}={\mathbf{P}}_{q_3,0}={(sin{2.5}^{\circ }/\sqrt{3})}^2{\mathbf{I}}_{1\times 1}$ and ${\mathbf{P}}_{\omega ,0}={0.001}^2{\mathbf{I}}_{3\times 3}({}^{\circ }/s{)}^2$.

Suppose spacecraft $S_1$ could track $S_2$ using optical equipment, where possible, with a 1 s measurement. The angle-measurement error in Equation (17) is assumed to be Gaussian white noise with a standard deviation (STD) of ${0.01}^{\circ }$ (equal to ${0.6}^{\prime }$).

4.2. Results and Discussion

During the observability test, the state transfer matrix is obtained using MATLAB function ode45, the rank of OM $\mathrm{rank}\left(\mathbf{N}\right)$ is calculated using function rank and the CN of the system $\mathrm{cond}\left(\mathbf{N}\right)$ is obtained by function cond. All the following operations are executed on MATLAB R2018b [40]. In addition, the particular situations for the autonomous attitude and orbit determination with two spacecraft are simulated to verify the observability analysis. The nominal orbit elements are given in Table 1. The estimation problems are solved using the traditional fifth-order CKF. For convenience, the quaternion of spacecraft $S_1$ is transferred into the form of Euler angles using the MATLAB function quat2angle. The results are given as follows.

4.2.1. Case I: Both Orbits of $S_1$ and $S_2$ Are Known, and the Attitude of $S_1$ Is Unchanged

Figure 3 displays the observability results of case I, using the two-body dynamics (state model as Equation (3)). For case I, 60 measurements (equivalent to 60 s) are executed during the navigation process. As shown in Figure 3, the upper stacked subplot illustrates the rank of the observability matrix (OM), where the red line represents the unobservable period (i.e., $\mathrm{rank}\left(\mathbf{N}\right)<4$) and the blue line represents the observable period (i.e., $\mathrm{rank}\left(\mathbf{N}\right)\ge 4$). Moreover, the lower stacked subplot demonstrates the condition number (CN) of OM (note that in Figure 3, the y-label represents the reciprocal of CN). In addition, the logarithmic value of the reciprocal of CN is also given in the lower stacked subplot, making the change curve of CN more obvious. Figure 3 shows that the CN deceases with respect to the observation time, which means that the observability of the system continuously improves as the number of measurements increases.

The detailed value of the corresponding index in Figure 3 is selectively illustrated in Table 2. For case I, the system is completely observable after only two measurements (at epochs $t_0$ and $t_k$, respectively).

Figure 4 depicts the estimation errors of the attitude of spacecraft $S_1$. It is shown that the errors converge to zero at approximately 50 s when the initial Euler angle errors are set to be ${2.5}^{\circ }$.

4.2.2. Case II: The Orbit of $S_1$ Is Known, and the Attitude of $S_1$ Is Unchanged

The observability results of case II are illustrated in Figure 5 and Table 3. As indicated in Figure 4, the system of Equation (39) is observable after approximately 94 measurements (from epoch $t_0$ to epoch $t_{93}$). Compared with case I, the system corresponding to case II is much more difficult to observe. Note that as containing the state vector of target spacecraft $S_2$, the dimension of the state to be estimated in case II is higher than that of case I, which implies more effort in measurements.

The estimation results are illustrated in Figure 6 and Figure 7. As shown in Figure 6, the attitude converges much faster than the orbit state, as the previous item converges within only a few minutes, while it takes approximately five hours for the orbit state of spacecraft $S_2$ to converge. Note that measurement model (17) is much more sensitive to the attitude of observer $S_1$ than the orbits of the observer and the target. Therefore, the attitude converges before the orbit state converges. Note that the observation interval has no effect on the observability, so the influence of shadow or light conditions are not considered in this paper. In most cases, the long-term observation in Figure 6 is impossible to realize due to eclipses, but this simplification is reasonable considering the research content of this paper.

4.2.3. Case III: Both the Orbits of $S_1$ and $S_2$ Are Unknown, and the Attitude of $S_1$ Is Unchanged

In this case, three kinds of dynamics (state model as Equation (3) and Equation (4), respectively) are considered, as illustrated in Figure 8 and Figure 9. It can be concluded that for dynamics with nonperturbation and $J_2$ perturbation, the systems are unobservable (as shown in Figure 8).

For case III, the state vector with 16 variables is estimated. However, when considering two-body dynamics, the observable states are 13-dimensional, which means that only 13 variables (or variable combinations) of the 16-dimensional state vector can be observed. In addition, when taking the $J_2$ perturbation into consideration, one more state variable (or variable combination) could be observed, suggesting that the system is still unobserved.

Fortunately, as depicted in Figure 6 and Table 4, the system is observable when $J_2$, $J_3$ and $J_4$ perturbations are considered. In this situation, a total of 16 variables are observed after 490 angle-only measurements. It should be noted that strictly two-body dynamics, or dynamics with a particular perturbation, are, of course, unlikely. This is because aspheric perturbation of the celestial body (i.e., Earth in this paper) contains higher order terms. Moreover, the solar pressure, atmospheric drag and gravitational perturbation of the third body could also influence the orbiting spacecraft. However, the significance of the observability analysis is to state (as might be expected) that when the dynamics are very close to the two-body dynamics (e.g., when an orbit is high, for example, high-orbit GPS satellites), it is difficult to estimate the orbit because the influence of aspheric perturbation is weak.

The simulation results of the subcase with nonperturbation and $J_2$ perturbation are illustrated in Figure 10. It was shown that estimation quality is poor and appears to be diverging (the $v_z$ for $S_1$ and $v_x$ for $S_2$ are clearly divergent with time), meaning that the system is unobservable under the given dynamics and initial conditions, which validates the observability analysis listed in Table 4.

Figure 11 and Figure 12 depict the autonomous attitude and orbit determination results of the two spacecraft considering $J_2$, $J_3$, and $J_4$ perturbations. Although exhibiting an obvious oscillation, the system still succeeds in converging. However, compared with the results in case II, the convergence is much weaker, meaning that this form is not stable enough. It is noted that although we have proven that case II under dynamics (5) is observable, this does not contradict the estimation results obtained here but implies that the system is higher-order locally weakly observable [24].

4.2.4. Case IV: Both the Orbits of $S_1$ and $S_2$ Are Unknown, and the Angular Velocity of $S_1$ Is Known

In case IV, the results are demonstrated in Figure 13 and Figure 14 and Table 5. As shown in Figure 13 and Table 5, when applying a known attitude maneuver to observer spacecraft $S_1$, the system is completely observable even under the simplest two-body dynamics.

As indicated in Table 5, only 344 measurements are needed to observe the 16 variables, which is less than that of case III (situation considering $J_2$, $J_3$, and $J_4$ perturbations in Table 4). Furthermore, at time epoch $t_{489}$, which is difficult to observe in case III, the CN of the system is larger than that of case IV (for case III, as shown in Table 4, $\frac{1}{\mathrm{cond}\left(\mathbf{N}\right)}=2.3821\times {10}^{-13}$, while for case IV, listed in Table 5, $\frac{1}{\mathrm{cond}\left(\mathbf{N}\right)}=3.4291\times {10}^{-12}$), indicating that compared to the perturbation acceleration, a suitable attitude maneuver is more likely to attach obvious improvement to the observability of the two-spacecraft system.

Figure 8 compares the influence of different dynamics. In Figure 14, the blue solid line, the red dashed line and the orange dash-dotted line represent the dynamics (3) and (4), respectively. It is illustrated that, for all three conditions, the system becomes observable around epoch $t_{350}$, which means that the dynamics make no difference to the observability of the system (note that the subtle difference could be recognized as the outcome of numerical calculation during the usage of ode45, rank and cond).

The autonomous attitude and orbit determination results of case IV are shown in Figure 15 and Figure 16, respectively. The results show strong convergence within five hours, implying that the corresponding system is completely observable. In addition, it is observed that the state estimations of Case IV present a significantly better stability than those of Case III. With a known attitude maneuver executed, the navigation system is more stable, and the certainty of the estimates improves compared with the results shown in Figure 12. In conclusion, by comparison, the attitude maneuver makes the system more observable and the estimation more accurate.

4.2.5. Case V: Both the Orbits of $S_1$ and $S_2$ Are Unknown, and the Angular Velocity of $S_1$ Is Unchanged

The test results of case V are illustrated in Figure 17 and Table 6. When not considering any perturbation, the system is observable at epoch $t_{1028}$. Compared with the situation in which the angular velocity of $S_1$ is known, it is slightly more difficult to estimate the system when the angular velocity needs to be determined.

As shown in Figure 18, the conclusions are summarized that the perturbations have almost no influence on the observability of the system, although the angular velocity is unknown.

Note that the system is sensitive to the attitude and the angular velocity of $S_1$; hence, for case V, the numerical observability measures are computed from the initial simulation epoch to the end of the simulation (within a total of 4 h) with a measurement frequency of 1 per 0.5 s. The results are illustrated in Figure 19, Figure 20 and Figure 21. Figure 19 and Figure 21 show that the attitude and the angular velocity converge at approximately 1 h, while as expected, the orbit states converge at approximately 2 h (Figure 20).

Even under the situation in which the angular velocity is to be estimated, the estimation quality of case V is still healthier than that of case III, indicating that the attitude maneuver is superior to the complex perturbation with respect to system observability.

The above discussion of the observability of the attitude and orbit is summarized in Table 7. Note that in Table 7, the symbol ‘✓’ denotes that the corresponding item is known, while the symbol ‘✗’ states that the item is unknown. Moreover, one item is suggested to not exist if it is marked ‘none’. For example, subcase II of case III indicates that the $J_2$ perturbation is considered, both the orbits of $S_1$ and $S_2$ and the attitude of $S_1$ are estimated, and no attitude maneuver is executed.

5. Conclusions

In this paper, the autonomous attitude and orbit determination problem of a two-spacecraft system using angle-only measurements is studied. The observability of the system is analyzed based on the theory of observability matrix. Five cases are analyzed and the observability analysis results are as follows:

• •. When the orbits of both observer and target are known, and the attitude of the target is unknown and unchanged, the navigation system is observable.

• •. When the orbit of observer is known, and the attitude and orbit of the target are unknown, the navigation system is observable.

• •. When the orbits of observer and target and the attitude of the target are unknown, the navigation system is unobservable in the two-body dynamics. The navigation system becomes observable when considering high-order perturbations.

• •. When the orbits of observer and target and the attitude of the target are unknown, and the attitude of the target is changed, the navigation system is observable.

In addition, the observability analysis and the filter results both verify that compared to the perturbation acceleration, a suitable attitude maneuver is more likely to attach obvious improvement to the observability of the two-spacecraft system.

Author Contributions

Conceptualization, J.W., C.L., Y.W. and Q.X.; methodology and software, X.Z. and Y.S. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Diagram of a two-spacecraft system.

Figure 2. Orbits of spacecraft [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 3. Observability simulation of case I (nonperturbation).

Figure 4. Attitude determination errors of case I (nonperturbation) and the corresponding larger plot.

Figure 5. Observability simulation of case II (nonperturbation).

Figure 6. Attitude determination errors of case II (nonperturbation) and the corresponding larger plot.

Figure 7. Orbit determination errors of spacecraft [Forumla omitted. See PDF.] in case II (nonperturbation).

Figure 8. Observability simulation of case III (nonperturbation and with [Forumla omitted. See PDF.] perturbation).

Figure 9. Observability simulation of case III (with [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] perturbations).

Figure 10. Orbit determination errors of case III (nonperturbation and with [Forumla omitted. See PDF.] perturbation). (a) Orbit estimation errors of spacecraft [Forumla omitted. See PDF.]. (b) Orbit estimation errors of spacecraft [Forumla omitted. See PDF.].

Figure 11. Attitude determination errors of case III (with [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] perturbations) and the corresponding larger plot.

Figure 12. Orbit determination errors of case III (with [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] perturbations. (a) Orbit estimation errors of spacecraft [Forumla omitted. See PDF.]. (b) Orbit estimation errors of spacecraft [Forumla omitted. See PDF.].

Figure 13. Observability simulation of case IV (nonperturbation).

Figure 14. Observability comparisons of case IV.

Figure 15. Attitude determination errors of case IV (nonperturbation) and the corresponding larger plot.

Figure 16. Orbit determination errors of case IV (nonperturbation). (a) Orbit estimation errors of spacecraft [Forumla omitted. See PDF.]. (b) Orbit estimation errors of spacecraft [Forumla omitted. See PDF.].

Figure 17. Observability simulation of case V (with nonperturbation).

Figure 18. Observability comparisons of case V.

Figure 19. Attitude determination errors of case V (nonperturbation) and the corresponding larger plot.

Figure 20. Orbit determination errors of case V (nonperturbation). (a) Orbit estimation errors of spacecraft [Forumla omitted. See PDF.]. (b) Orbit estimation errors of spacecraft [Forumla omitted. See PDF.].

Figure 21. The estimation errors of angular velocity of case V (nonperturbation).

Table 1

Nominal orbit elements of the circle orbits.

Spacecraft	h/km	e	i/deg	$\mathsf{\Omega }$/deg	$\mathit{\omega }$/deg	n/deg
$S_1$	500	0	45.05	29.93	132.9	−107.74
$S_2$	1000	0	45	94.8	199.0	−54.13

Table 2

Observability test results of case I (nonperturbation).

Epoch ${\mathit{t}}_{\mathit{k}}$	$\mathbf{Rank}\left(\mathbf{N}\right)$	Observability	$1/\mathbf{cond}\left(\mathbf{N}\right)$	$\mathbf{ln}[1/\mathbf{cond}\left(\mathbf{N}\right)]$
0	3	False	-	-
1	4	True	0.0005	−7.4740
2	4	True	0.0009	−6.9835
3	4	True	0.0012	−6.6692

Table 3

Observability test results of case II (nonperturbation).

Epoch ${\mathit{t}}_{\mathit{k}}$	$\mathbf{Rank}\left(\mathbf{N}\right)$	Observability	$1/\mathbf{cond}\left(\mathbf{N}\right)$	$\mathbf{ln}[1/\mathbf{cond}\left(\mathbf{N}\right)]$
0	3	False	-	-
1	5	False	-	-
2	7	False	-	-
3	9	False	-	-
⋯	⋯	⋯	⋯	⋯
92	9	False	-	-
93	10	True	$5.1763\times {10}^{-14}$	−30.5921
94	10	True	$5.3998\times {10}^{-14}$	−30.5498

Table 4

Observability test results of case III (with $J_2$, $J_3$, and $J_4$ perturbations).

Epoch ${\mathit{t}}_{\mathit{k}}$	$\mathbf{Rank}\left(\mathbf{N}\right)$	Observability	$1/\mathbf{cond}\left(\mathbf{N}\right)$	$\mathbf{ln}[1/\mathbf{cond}\left(\mathbf{N}\right)]$
0	3	False	-	-
1	5	False	-	-
2	7	False	-	-
3	9	False	-	-
⋯	⋯	⋯	⋯	⋯
488	15	False	-	-
489	16	True	$2.3821\times {10}^{-13}$	−29.0656
490	16	True	$2.4174\times {10}^{-13}$	−29.0509

Table 5

Observability test results of case IV (nonperturbation).

Epoch ${\mathit{t}}_{\mathit{k}}$	$\mathbf{Rank}\left(\mathbf{N}\right)$	Observability	$1/\mathbf{cond}\left(\mathbf{N}\right)$	$\mathbf{ln}[1/\mathbf{cond}\left(\mathbf{N}\right)]$
0	3	False	-	-
1	6	False	-	-
2	8	False	-	-
3	10	False	-	-
⋯	⋯	⋯	⋯	⋯
342	15	False	-	-
343	16	True	$1.9378\times {10}^{-13}$	−29.2720
344	16	True	$1.9865\times {10}^{-13}$	−29.2472
⋯	⋯	⋯	⋯	⋯
489	16	True	$3.4291\times {10}^{-12}$	−26.3987
490	16	True	$3.4823\times {10}^{-12}$	−26.3833

Table 6

Observability test results of case V.

Epoch ${\mathit{t}}_{\mathit{k}}$	$\mathbf{Rank}\left(\mathbf{N}\right)$	Observability	$1/\mathbf{cond}\left(\mathbf{N}\right)$	$\mathbf{ln}[1/\mathbf{cond}\left(\mathbf{N}\right)]$
0	3	False	-	-
1	6	False	-	-
2	9	False	-	-
3	11	False	-	-
⋯	⋯	⋯	⋯	⋯
1027	18	False	-	-
1028	19	True	$5.7892\times {10}^{-13}$	−28.1776
1029	19	True	$5.8186\times {10}^{-13}$	−28.1725

Table 7

Observability under different situations.

Case	I	II		III		IV	V
Perturbation	none	none	none	$J_2$	$J_2+J_3+J_4$	none	none
Dynamics	Equation (3)	Equation (3)	Equation (3)	Equation  (4)	Equation (4)	Equation (3)	Equation (3)
Quaternion of $S_1$	✗	✗	✗	✗	✗	✗	✗
Orbit state of $S_1$	✓	✓	✗	✗	✗	✗	✗
Orbit state of $S_2$	✓	✗	✗	✗	✗	✗	✗
Angular velocity of $S_1$	none	none	none	none	none	✓	✗
Observability	Y	Y	N	N	Y	Y	Y

References

1. Qin, T.; Qiao, D.; Macdonald, M. Relative orbit determination using only intersatellite range measurements. J. Guid. Control Dyn.; 2019; 42, pp. 703-710. [DOI: https://dx.doi.org/10.2514/1.G003819]

2. Vishwajeet, K.; Singla, P.; Jah, M. Nonlinear Uncertainty Propagation for Perturbed Two-Body Orbits. J. Guid. Control Dyn.; 2014; 37, pp. 1415-1425. [DOI: https://dx.doi.org/10.2514/1.G000472]

3. Vishwajeet, K.; Singla, P. Adaptive split/merge-based Gaussian mixture model approach for uncertainty propagation. J. Guid. Control Dyn.; 2018; 41, pp. 603-617. [DOI: https://dx.doi.org/10.2514/1.G002801]

4. Qiao, D.; Zhou, X.; Zhao, Z.; Qin, T. Asteroid Approaching Orbit Optimization Considering Optical Navigation Observability. IEEE Trans. Aerosp. Electron. Syst.; 2022; 58, pp. 5165-5179. [DOI: https://dx.doi.org/10.1109/TAES.2022.3167653]

5. Yang, Z.; Luo, Y.; Lappas, V.; Tsourdos, A. Nonlinear analytical uncertainty propagation for relative motion near J2-perturbed elliptic orbits. J. Guid. Control Dyn.; 2018; 41, pp. 888-903. [DOI: https://dx.doi.org/10.2514/1.G003071]

6. Li, X.; Qiao, D.; Li, P. Bounded trajectory design and self-adaptive maintenance control near non-synchronized binary systems comprised of small irregular bodies. Acta Astronaut.; 2018; 152, pp. 768-781. [DOI: https://dx.doi.org/10.1016/j.actaastro.2018.09.028]

7. Yang, Z.; Luo, Y.; Zhang, J. Nonlinear semi-analytical uncertainty propagation of trajectory under impulsive maneuvers. Astrodynamics; 2019; 3, pp. 61-77. [DOI: https://dx.doi.org/10.1007/s42064-018-0036-7]

8. Zhou, X.; Cheng, Y.; Qiao, D.; Huo, Z. An adaptive surrogate model-based fast planning for swarm safe migration along halo orbit. Acta Astronaut.; 2022; 194, pp. 309-322. [DOI: https://dx.doi.org/10.1016/j.actaastro.2022.02.006]

9. Wu, W.; Tang, Y.; Zhang, L.; Qiao, D. Design of communication relay mission for supporting lunar-farside soft landing. Sci. China Inf. Sci.; 2018; 61, [DOI: https://dx.doi.org/10.1007/s11432-017-9202-1]

10. Tang, Y.; Wu, W.; Qiao, D.; Li, X. Effect of orbital shadow at an Earth-Moon Lagrange point on relay communication mission. Sci. China Inf. Sci.; 2017; 60, pp. 1-10. [DOI: https://dx.doi.org/10.1007/s11432-016-9069-9]

11. Qin, T.; Qiao, D.; Macdonald, M. Relative orbit determination for unconnected spacecraft within a constellation. J. Guid. Control Dyn.; 2021; 44, pp. 614-621. [DOI: https://dx.doi.org/10.2514/1.G005424]

12. Julier, S.J.; Uhlmann, J.K. Unscented filtering and nonlinear estimation. Proc. IEEE Inst. Electr. Electron. Eng.; 2004; 92, pp. 401-422. [DOI: https://dx.doi.org/10.1109/JPROC.2003.823141]

13. Park, R.S.; Scheeres, D.J. Nonlinear semi-analytic methods for trajectory estimation. J. Guid. Control Dyn.; 2007; 30, pp. 1668-1676. [DOI: https://dx.doi.org/10.2514/1.29106]

14. Zhou, X.; Wang, S.; Qin, T. Multi-Spacecraft Tracking and Data Association Based on Uncertainty Propagation. Appl. Sci.; 2022; 12, 7660. [DOI: https://dx.doi.org/10.3390/app12157660]

15. Qin, T.; Macdonald, M.; Qiao, D. Fully Decentralized Cooperative Navigation for Spacecraft Constellations. IEEE Trans. Aerosp. Electron. Syst.; 2021; 57, pp. 2383-2394. [DOI: https://dx.doi.org/10.1109/TAES.2021.3060734]

16. Han, H.; Qiao, D.; Chen, H.; Li, X. Rapid planning for aerocapture trajectory via convex optimization. Aerosp Sci Technol; 2019; 84, pp. 763-775. [DOI: https://dx.doi.org/10.1016/j.ast.2018.11.009]

17. Chen, Q.; Qiao, D.; Wen, C. Orbital Element Reachable Set After Gravity Assists of Planets in Elliptical Orbits. J. Guid. Control Dyn.; 2020; 43, pp. 989-997. [DOI: https://dx.doi.org/10.2514/1.G004611]

18. Zhou, X.; Qin, T.; Meng, L. Maneuvering Spacecraft Orbit Determination Using Polynomial Representation. Aerospace; 2022; 9, 257. [DOI: https://dx.doi.org/10.3390/aerospace9050257]

19. Abusali, P.A.; Tapley, B.D.; Schutz, B.E. Autonomous navigation of global positioning system satellites using cross-link measurements. J. Guid. Control Dyn.; 1998; 21, pp. 321-327. [DOI: https://dx.doi.org/10.2514/2.4238]

20. Opromolla, R.; Fasano, G.; Rufino, G.; Grassi, M. A review of cooperative and uncooperative spacecraft pose determination techniques for close-proximity operations. Prog. Aerosp. Sci.; 2017; 93, pp. 53-72. [DOI: https://dx.doi.org/10.1016/j.paerosci.2017.07.001]

21. Zhang, L.; Yang, H.; Lu, H.; Zhang, S.; Cai, H.; Qian, S. Cubature Kalman filtering for relative spacecraft attitude and position estimation. Acta Astronaut.; 2014; 105, [DOI: https://dx.doi.org/10.1016/j.actaastro.2014.09.007]

22. Li, Y.; Zhang, A. Observability analysis and autonomous navigation for two satellites with relative position measurements. Acta Astronaut.; 2019; 163, pp. 77-86. [DOI: https://dx.doi.org/10.1016/j.actaastro.2019.02.030]

23. Duong, N.; Winn, C.B. Orbit determination by range-only data. J. Spacecr. Rockets; 1973; 10, pp. 132-136. [DOI: https://dx.doi.org/10.2514/3.61858]

24. Hu, Y.; Sharf, I.; Chen, L. Three-spacecraft autonomous orbit determination and observability analysis with inertial angles-only measurements. Acta Astronaut.; 2020; 170, pp. 106-121. [DOI: https://dx.doi.org/10.1016/j.actaastro.2020.01.005]

25. Hu, Y.; Sharf, I.; Chen, L. Distributed orbit determination and observability analysis for satellite constellations with angles-only measurements. Automatica; 2021; 129, 109626. [DOI: https://dx.doi.org/10.1016/j.automatica.2021.109626]

26. Psiaki, M.L. Absolute orbit and gravity determination using relative position measurements between two satellites. J. Guid. Control Dyn.; 2011; 34, pp. 1285-1297. [DOI: https://dx.doi.org/10.2514/1.47560]

27. Lee, S.; Mortari, D. Design of constellations for earth observation with intersatellite links. J. Guid. Control Dyn.; 2017; 40, pp. 1263-1271. [DOI: https://dx.doi.org/10.2514/1.G001710]

28. Hill, K.; Born, G.H. Autonomous interplanetary orbit determination using satellite-to-satellite tracking. J. Guid. Control Dyn.; 2007; 30, pp. 679-686. [DOI: https://dx.doi.org/10.2514/1.24574]

29. Hill, K.; Born, G.H. Autonomous orbit determination from lunar halo orbits using crosslink range. J. Spacecr. Rockets; 2008; 45, pp. 548-553. [DOI: https://dx.doi.org/10.2514/1.32316]

30. Perez, A.C.; Geller, D.K.; Lovell, T.A. Non-iterative angles-only initial relative orbit determination with J2 perturbations. Acta Astronaut.; 2018; 151, pp. 146-159. [DOI: https://dx.doi.org/10.1016/j.actaastro.2018.06.033]

31. Ou, Y.; Zhang, H.; Li, B. Absolute orbit determination using line-of-sight vector measurements between formation flying spacecraft. Astrophys. Space Sci.; 2018; 363, 76. [DOI: https://dx.doi.org/10.1007/s10509-018-3293-2]

32. Luo, Y.; Qin, T.; Zhou, X. Observability Analysis and Improvement Approach for Cooperative Optical Orbit Determination. Aerospace; 2022; 9, 166. [DOI: https://dx.doi.org/10.3390/aerospace9030166]

33. Mao, X.; Visser, P.; van den IJssel, J. Absolute and relative orbit determination for the CHAMP/GRACE constellation. Adv. Space Res.; 2019; 63, pp. 3816-3834. [DOI: https://dx.doi.org/10.1016/j.asr.2019.02.030]

34. Geller, D.K.; Klein, I. Angles-only navigation state observability during orbital proximity operations. J. Guid. Control Dyn.; 2014; 37, pp. 1976-1983. [DOI: https://dx.doi.org/10.2514/1.G000133]

35. Qian, Y.; Li, C.; Jing, W.; Hwang, I.; Wei, J. Sun-Earth-Moon autonomous orbit determination for quasi-periodic orbit about the translunar libration point and its observability analysis. Aerosp Sci Technol; 2013; 28, pp. 289-296. [DOI: https://dx.doi.org/10.1016/j.ast.2012.11.009]

36. Sun, X.; Geng, C.; Deng, L.; Chen, P. Geolocation of Formation-Flying Spacecraft Using Relative Position Vector Measurements. J. Guid. Control Dyn.; 2022; 45, pp. 764-773. [DOI: https://dx.doi.org/10.2514/1.G006377]

37. Han, H.; Qiao, D. Optimization for the aeroassisted orbital plane change with the synergetic maneuver using the hp-adaptive pseudospectral method. J. Aerosp. Eng.; 2017; 30, 04017076. [DOI: https://dx.doi.org/10.1061/(ASCE)AS.1943-5525.0000788]

38. Han, H.; Qiao, D.; Chen, H. Optimization of aeroassisted rendezvous and interception trajectories between non-coplanar elliptical orbits. Acta Astronaut.; 2019; 163, pp. 190-200. [DOI: https://dx.doi.org/10.1016/j.actaastro.2018.11.028]

39. Chen, Q.; Qiao, D.; Wen, C. Reachable domain of spacecraft after a gravity-assist flyby. J. Guid. Control Dyn.; 2019; 42, pp. 931-940. [DOI: https://dx.doi.org/10.2514/1.G003576]

40. MATLAB. Version R2018b; The MathWorks Inc.: Natick, MA, USA, 2018.