Angle-Only Cooperative Orbit Determination Considering Attitude Uncertainty

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.


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

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

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: The state equation of the spacecraft orbit under the two-body dynamics can be expressed in general form as:Ẋ In the case of ideal two-body dynamics (i.e., particle dynamics model), Equation (2) is given by:r In the case of considering J 2 , J 3 , and J 4 perturbations, Equation (3) is rewritten as: where µ e is the Earth's gravitational parameter, and a J 2 , a J 3 and 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: where ϕ is the latitude of the spacecraft's ground trace; The attitude of spacecraft S 1 is represented through the quaternion, defined by: where: wheren is a unit vector corresponding to the principal axis of rotation and α is the angle of rotation. The quaternion kinematics are derived through the spacecraft's angular velocity as follows: where ω = [ω x ω y ω z ] T denotes the angular velocity, J represents the moment of inertia of the spacecraft, M is the sum of external moments, and matrix Ω(ω) is defined as: The four quaternion elements satisfy the following normalization constraint: Combining the elements for the above orbits and attitude, the state vector to be estimated is: The state model is given by:

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 α and β, given as: where ∆x = x 2 − x 1 , ∆y = y 2 − y 1 and ∆z = z 2 − z 1 ; ε α and ε β , 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: where r 1 = [x 1 , y 1 , z 1 ] T and r 2 = [x 2 , y 2 , z 2 ] T ; ε denotes the randomly distributed noise for vector measurement y . 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 y is given by the following: where ε is the corresponding measurement noise vector characterized by a normal distribution with zero mean and covariance R ∈ R 3×3 and R B E is the rotation matrix from ECI Hence, the navigation problem is given by:

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 {y 0 , y 1 , · · · , y k−1 } from time epoch t 0 to t k−1 , the OM is represented as: where Φ(t k , t 0 ) is the state transformation matrix (STM) from t 0 to t k and The differential equation of the STM is as follows: and is initialized by: where n is the dimension of the state vector to be estimated and I n×n is an n-dimensional unit matrix. Note that the differential term ∂F(X)/∂X is also known as the Jacobi matrix of state model (14). An OM with a full rank (i.e., rank(N) = 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 cond(N) = N · N −1 . The smaller the CN, the better the observability [24].

Review of Fifth-Order CKF
In this section, the fifth-order CKF algorithm is briefly summarized. First, consider a discrete nonlinear system as: where x k ∈ R n is the state vector at time epoch k and z k ∈ R m is the measurement. w k−1 ∈ R n and v k ∈ R m denote the independent system and measurement noise, respectively, and are both considered independent and white Gaussian distributions with covariances Q k−1 and 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 is a nonlinear function. The integral with respect to the general Gaussian distribution N(x;x, P) can be further approximated by: where P = SS T , N P is the total number of points, and γ i and W i are the quadrature points and weights, respectively, corresponding to the Gaussian distribution N(x;x, P). Specifically, the integral with respect to N(x; 0, I) can be approximated by the following quadrature rule: According to the cubature rule and Equation (26), in fifth-order CKF, the cubature points ζ i are given by [21]: . . , n ; j = 1, 2, . . . , n(n − 1)/2 (28) where n is the state dimension of the system to be estimated, A k is the Cholesky decomposition of covariance matrix P k at epoch t k and P k = A k A T k , and e i is the i-th column of the n-th identity matrix I n . The point sets s + j and s − j are given as follows.
The corresponding weight w i of each Cubature point ζ i is given by Then, the fifth-order CKF algorithm is summarized as follows: i. Time updating: ii. Measurement updating where Cholesky(·) represents the Cholesky decomposition method, χ * i,k|k−1 is the Cubature point generated from states and z * i,k|k−1 represents the Cubature point generated from measurements.

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.

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., ω = [0 0 0] T ). Thus, the state vector can be written as: The state model is given by: The measurement is shown in Equation (17) and then the partial derivative of the intersatellite angle to the quaternion of spacecraft S 1 H k is of the form: For the state model in Equation (34), the Jacobi matrix is given as: Therefore, the STM has the following format: 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: The corresponding state model is as follows: The partial differential matrix H k of measurement (17) is written as: The propagation of the STM for state model (39) is then given by: Hence, the STM Φ(t k , t 0 ) has the following format: 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 ω = [0 0 0] T . Therefore, the state to be estimated is: The state model is then given by: Similar to Equation (40), the partial differential matrix of measurement is given as: where: ∂ ∂r 1 The STM Φ(t k , t 0 ) in this case has the following format: 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 ω = [ω x ω y ω z ] T is foreknown (assume that ω x , ω y and ω z are constant and satisfy ω 2 x + ω 2 y + ω 2 z = 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: The STM is in the form of: 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 ω = [ω x ω y ω z ] 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 X is as follows:

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.

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, Ω, ω, 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 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 q 13 = [0, 0, 0] T and q 0 = 1 − q T 13 q 13 = 1. For all the cases above, the initial estimation of the quaternion is given by: The initial angular velocity of spacecraft S 1 , if existing (i.e., for case IV and case V), is set to be ω = [0.01, 0.01, 0.01] T • /s . The initial estimation of angular velocity is considered asω = [0.09, 0.09, 0.09] T • /s if ω is to be estimated (only for case V).

Results and Discussion
During the observability test, the state transfer matrix is obtained using MATLAB function ode45, the rank of OM rank(N) is calculated using function rank and the CN of the system cond(N) 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.

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., rank(N) < 4) and the blue line represents the observable period (i.e., rank(N) ≥ 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).   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 Figures 6 and 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.  In this case, three kinds of dynamics (state model as Equation (3) and Equation (4), respectively) are considered, as illustrated in Figures 8 and 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.  Figures 11 and 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].  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,   1 cond(N) = 2.3821 × 10 −13 , while for case IV, listed in Table 5,   1 cond(N) = 3.4291 × 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 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 Figures 15 and 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. 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 Figures 19-21. Figures 19 and 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.

Conclusions
In this paper, the autonomous attitude and orbit determination problem of a twospacecraft 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.