Array-Aided Precise Orbit and Attitude Determination of CubeSats using GNSS

2 COMBINED ARRAY-AIDED PRECISE ORBIT AND ATTITUDE DETERMINATION MODEL

In this section, a combined array-aided model for precise orbit and attitude determination is developed. The sizes of the main matrices defined here are given in Table A-1 in the Appendix. The sizes of the remaining matrices can be computed via linear algebra.

Let us begin by forming a vectorial representation of the inter-satellite single-differenced (SD) code and phase observations for frequencies 1 to f. These observations are obtained from satellites 1 to m, collected by arbitrary onboard antennae r. Considering and for arbitrary frequency f, where s serves as the reference GNSS satellite, we have the following model:

1

Here, E(.) represents the expected value, and ⊗ indicates the Kronecker product (Schott, 2016). e_f = [1,⋯,1]_1×f, and contains the difference of unit vectors to satellites i and s (for i = 1,…, m). This term is computed by using the orbital vectors of GNSS satellites (x^s and xⁱ) at the time of transmission (t – τ) and the antenna coordinates (x_r) at the time of reception (t) of the signals, considering τ as the signal travel time for each observation. b_r contains the unknown orbital components of the CubeSat in the Earth-centered, Earth-fixed (ECEF) coordinate system. is the ionospheric coefficient for frequency f, where λ_f indicates the wavelength of that frequency. I_m−1 is the identity matrix with size m −1. and denote SD ionospheric delays and phase ambiguities, respectively. To obtain a more compact model, we include the SD observations in y_r = [P_r, Φ_r]^T and define and A = [0^T, Λ^T] ⊗ I_m−1, where Λ = diag[λ₁,⋯, λ_f] and . The GNSS observation model is then expressed as follows:

2

Now, let us consider an array of a antennae onboard a CubeSat. For attitude determination, the vectorial form of the linearized double-differenced (DD) GNSS observations between each arbitrary pair of antennae 1 and 2 can be formulated as follows:

3

where is the DD ambiguity and is defined as follows:

4

which can be replaced by denoting the baseline vector as b₁₂ and defining x₁(t) = x₂(t) + b₁₂ as follows:

5

To linearize Equation (5), a Taylor series expansion is applied at b₁₂ = 0 as follows:

6

Finally, by defining and and forming the geometry component of the design matrix as , we obtain the following vectorial forms of the linearized DD GNSS observations:

7

where , , and . The matrix comprises the DD ambiguities, transformed to range via multiplication by the frequency wavelength. Note that the elements of G₁₂ are similar to those of G_r.

It is assumed that the covariance matrices of undifferenced observations are and , which are formed by multiplying the standard deviation σ by the cofactor matrix C computed from elevation-angle-based or SNR-based weighting functions (Allahvirdi-Zadeh & El-Mowafy, 2021). By considering the DD operator as , combining the undifferenced observations of the two antennae as and , and expressing the covariance matrices of these observations as and , we obtain the following for the covariance matrices of the DD observables of the two antennae:

8

The next step is to collect all DD observations and ambiguities of all of the base-lines in and Z = [(N₁₂)^T,⋯, (N_1(a–1))^T]^T. We then define A = [0^T, Λ^T]^T ⊗ I_m–1 and to link the DD data to the ambiguities and geometry of the satellites. Utilizing the multivariate Gauss–Markov equation and arranging the baselines in matrix B = [b₁₂,⋯, b_1(a–1)], we derive the following multivariate formulation of the model for determining the attitude of the CubeSat:

9

There are two critical points to consider in this model:

1. Because of the small size of CubeSats (typically 10×10×10 cm³ per unit), the onboard antennae are mounted very close to each other. Therefore, the unit line-of-sight vectors to the same satellites are considered identical, giving a negligible error. Hence, G₁₂ = G₁₃ =… = G_1(a–1) = G_r.

2. The unknown parameters are DD ambiguities in Z and baselines in B. Because the array geometry is known in the SRF as B₀, an estimation of the baselines will enable us to determine the CubeSat attitude by applying the attitude matrix R as B = RB₀. This attitude matrix is used to transform the SRF to the ECEF frame and can be expressed by using rotation angles or quaternions.

To further strengthen the model and fully exploit the known geometry of the baselines, the two following constraints are applied:

• The DD ambiguities are integers (Z ∈ Z^{f(m–1)×(a–1)}, where a is the total number of antennae in the array).

• This attitude matrix is orthonormal (R^TR = I R ∈ O^3×q, where q is the number of unknown baselines).

The first constraint ensures that highly accurate phase observations are used, which, in turn, enable precise orbit and attitude determination for CubeSats. The second constraint provides increased reliability in the estimation by enabling a higher success rate in determining the ambiguities compared with the unconstrained model.

Taking Equations (2) and (9) with the aforementioned constraints and using the vec operator to transform the m × n matrix into an mn-sized vector, we can define the following combined multivariate constrained orbit and attitude determination model for CubeSats, called the array-aided precise orbit and attitude determination (A-POAD) model:

10

In A-POAD, two distinct models for the orbit and attitude determination of CubeSats are combined; each can be solved independently. However, as explained below, correlations exist between observations in the combined model, which should be taken into account. Otherwise, the optimality of the solution will be reduced.

The observation set in the A-POAD model is reached from the SD observations by using the following equation (Li & Teunissen, 2014):

11

where c₁ = [1,0,⋯,0]^T, is the differencing matrix and Y_r = [y₁,⋯, y_a]. The covariance matrix of the A-POAD observations is calculated by applying the error propagation law to Equation (11) as follows:

12

where indicates the precision contribution of the antenna, Q_a is the covariance matrix of the array, Q_f is the contribution of the frequencies, and the matrices Q_p and Q_φ represent the undifferenced code and phase precision contributions, respectively. The term blkd indicates the block-diagonal matrix. Equation (12) includes non-diagonal values that indicate the correlation between the two sets of observations (y_r and Y) in our combined model.

To obtain a rigorous solution for the A-POAD model, the observations are decorrelated via the invertible decorrelating transformation defined by Teunissen (2012b) as follows:

13

Finally, we obtain the following decorrelated A-POAD model:

14

where the decorrelated observations are derived from and the combined weighted least-squares solution of all antenna positions is . The covariance matrix of the observations in Equation (14) is a block-diagonal matrix, confirming that the observations of the A-POAD model are now decorrelated. For a proof of the covariance matrix calculation, see the work by Teunissen (2012b) and Li & Teunissen (2014).

We can now obtain a robust solution, which is achieved by dividing the A-POAD model into two distinct components: 1) array-aided precise orbit determination (A-POD) and 2) array-aided precise attitude determination (A-PAD). These two components can be solved independently, allowing us to effectively leverage the entire covariance matrix of the observations in the A-POD component. According to Teunissen (2012b), it has been proven that the decorrelated observations have smaller variances than the original observations. The new variance is proportional to the number of antennae in the array as . This theoretical insight also validates the improvement of POD achieved by using the array-aided concept.

The first component of our solution is the A-POD model, defined as follows:

15

We can have two solutions for this A-POD model:

1. A-POD solution without integer ambiguities: In this scenario, the focus is solely on reconstructing the A-POD observations using Y. When the antenna configuration is symmetrical, the A-POD model provides observations of the center of gravity of the antenna array. It is worth noting again that these observations have improved accuracy compared with the primary POD observations, as discussed above. This improved accuracy is particularly valuable for CubeSats equipped with commercial off-the-shelf (COTS) sensors and GNSS receivers, which generally exhibit higher noise levels than the advanced ADCSs employed in larger LEO satellites. The impact of these noise sources on POD outputs was investigated in our earlier work (Allahvirdi-Zadeh, Awange et al., 2022).

2. A-POD solution with integer ambiguities: In this case, we initially utilize observations Y to reconstruct the A-POD observations. Subsequently, we employ Z to address the integer ambiguities within the following A-POD observations:

   16

   indicates the observations with fixed integer ambiguities, highlighting the effectiveness and benefit of our A-POD model. It is important to emphasize that the integer ambiguities are the output of the A-PAD model, which is defined and discussed in detail below.

The second component of our solution is the following multi-constrained A-PAD model:

17

The solution for this model is obtained in two steps. First, we solve the least-squares normal equation to obtain the unconstrained float estimations for ambiguities Ẑ and attitude :

18

Second, we solve the following minimization problem while considering the aforementioned constraints:

19

In this minimization problem, we first estimate the attitude matrix , which is conditioned on the estimation of the ambiguities, as follows:

20

Equation (19) shows that the estimation of the fixed ambiguity solution involves minimizing the objective function F(Z) while conditioning the attitude matrix estimation . This inclusion of a nonlinear constraint in the A-PAD model makes the search space for the integer ambiguities non-ellipsoidal and complex. To address this challenge, constrained searching methods such as search-and-expand and search-and-shrink algorithms are defined to bound the minimization objective function and the designated searching space (Giorgi et al., 2008; Giorgi, 2017). For instance, in the search-and-expand method, a lower bound is defined for the minimization objective function, and a search space for that lower bound is considered. The integer search is performed: if any candidates are found, the minimization objective function is evaluated for the candidates, and the candidate that returns the smallest value is the integer minimizer. If there are no candidates in this search space, the lower bound is expanded in another round. The fixed attitude matrix is estimated by substituting the integer ambiguities in Equation (20) with , which leads to a mixed ambiguity–attitude minimization problem.

Because of the complexity of the search space in the multi-constrained A-PAD model for populated arrays with more than three antennae, a set of linear constraints was implemented, leading to the application of the affine-constrained concept (Teunissen, 2012a), which can be modeled as follows:

21

This affine definition comes from the orthonormal matrix parametrization, which allows us to use the basis matrix (S) of the null space of B and update our A-PAD model with the equivalent implicit form of B = RB₀ as BS = 0. By doing this, we ignore the quadratic constraints and consider only the linear constraints. This replacement ensures a high probability of success for integer ambiguity resolution. The procedure of the affine-constrained A-PAD model is similar to that of the multi-constrained model. However, by defining the projector P_S = S(S^TPS)⁻¹S^TP, we can write the ambiguity objective function for the affine-constrained model as follows:

22

In contrast to the objective function in Equation (19), the above objective function provides an ellipsoidal search space with lower complexity and less computational burden. Thus, the affine-constrained A-PAD model is more efficient for CubeSats with limited power and processing resources.

To apply this approach in real-time mode and make it suitable for onboard processing, A-POD and A-PAD modules were implemented in the LeoPod software (Allahvirdizadeh, 2022). Figure 2 presents a flow diagram summarizing the proposed module for the ADCS based on our models. The observations obtained from the array are utilized to estimate the attitude angles and integer ambiguities within the A-PAD module. This estimation process can be conducted using either a multi-constrained or affine-constrained approach. Subsequently, these observations are decorrelated to improve their accuracy, and then, the fixed integer ambiguities are taken into account to reconstruct the A-POD observations. Finally, the POD procedure can be initiated, utilizing the reconstructed A-POD observations. During this process, the estimated attitude angles are employed as required, facilitating an accurate determination of the CubeSat’s precise orbit. Data quality control is performed by applying statistical testing using the detection, identification, and adaptation method inside the filter (Teunissen, 2006).

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 2

Flow diagram of A-POD and A-PAD modules implemented in the LeoPod software

By integrating the A-PAD and A-POD modules, our approach leverages the benefits of both precise attitude and precise orbit determination, resulting in improved accuracy and performance. To validate the A-POD and A-PAD algorithms, an antenna array is required. However, to our knowledge, no CubeSat is currently equipped with more than one GNSS antenna. Therefore, a simulated antenna array for a CubeSat using actual onboard observations is used in our tests to demonstrate the proposed method.

3 TESTING

3.1 Test Description

In this section, we present a simulation of an antenna array, illustrated in Figure 3. The size of the array and the proximity of the antennae within it are critical design factors, primarily dictated by the constraints of the CubeSat structure. The close proximity of the antennae (less than half of the signal wavelength) introduces mutual coupling and pattern distortion effects, which can impact signal correlation and the effectiveness of conventional signal processing and models designed for independent or uncorrelated signals. Meanwhile, pattern distortion can potentially cause frequency shifts that compromise proper functionality within the designated frequency band. To mitigate such effects, it is crucial to consider an appropriate distance between the antennae during the array design phase. The optimal distance between antenna elements in an array is more than half the wavelength of the designated frequencies, which, for our case in which GPS L1 is the higher frequency used, is roughly 10 cm.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 3

Layout of the simulated antenna array

The left CubeSat depicted in Figure 3 is a three-unit (3U) Spire CubeSat that we selected for the testing. The size of the top panel is 10×10 cm², where placing the array on this panel would not be feasible because of the mutual coupling and pattern distortions it would introduce. Alternatives such as placing antennae on the side panels could be explored to avoid these effects; however, this approach presents challenges regarding the number of satellites tracked by the receiver, which requires further investigation. Additionally, there are template constraints for each satellite based on its mission and the capacity of the bus and payloads, which are beyond the scope of this paper.

To address these considerations, we assumed that the 3U CubeSat is expanded to a 12U CubeSat and simulated an array of four antennae on the top panel of this CubeSat, as depicted in Figure 3. Antenna #1 serves as the main POD antenna, whereas the remaining antennae are simulated based on their coordinates in the SRF, as explained in the next section.

It is worth mentioning that the near-field environment, including other available sensors, may further amplify the aforementioned effects. In addition to constraints on the size of the array and the distances between antennae, this complexity necessitates pre-flight tests and in-flight calibration of both absolute and relative antenna phase patterns (Montenbruck, 2017).

Table 2 provides the coordinates of the simulated antennae in the SRF. The SRF’s orientation for the Spire CubeSats is defined as follows (personal communication with the Spire team):

• The Z-axis points towards the zenith direction.

• The Y-axis points towards the negative orbit normal.

• The X-axis is approximately aligned with the velocity direction but orthogonal to the Y-axis.

View this table:

• View inline
• View popup

TABLE 2

Coordinates of the Center of the Simulated Antennae in the SRF

These axes are joined at the center of mass (CoM) of the CubeSat, as shown in Figure 3.

3.4 Testing Results

In the initial test, we utilize simulated GPS observations obtained from an antenna array to estimate attitude angles. The estimation is performed twice, first with dual-frequency and then with single-frequency observations. The latter is important for CubeSats, as they often face power constraints that may restrict them to single-frequency observations. The primary objective of this test is to examine the effects of this constraint on data collection and its impact on the overall analysis.

Figure 4 illustrates the roll and pitch angles while ambiguities are resolved as float and fixed values via the A-PAD model for 24 h of 1-Hz GPS observations. These solutions are called “float” and “fixed” attitude angles, respectively. The left and right columns of the figure depict the dual-frequency and single-frequency cases, respectively. It is evident from the figure that ambiguity resolution has a notable impact on the estimated roll and pitch angles. Utilizing the antenna array in the A-PAD model not only enhances the precision of the estimated rotation angles but also renders the estimations more precise compared with magnetometers and sun sensors. The Spire CubeSat is equipped with a combination of these sensors in the ADCS (Spire, 2023), and their outputs are shown in the figure to visually confirm the accuracy of the A-PAD model with the antenna array. Additionally, this array configuration is more cost-effective, further emphasizing its advantage over alternative sensor options (see Table 1). Figure 4 consists of two rows, where the top row corresponds to the application of the multi-constrained (MC A-PAD) model, and the bottom row represents the affine-constrained (AC A-PAD) model. In the MC A-PAD case, the array utilizes antennae 1, 2, and 3, whereas the AC A-PAD model employs the entire array. The inclusion of an additional antenna and the consideration of the relevant constraint in the AC A-PAD model result in improved attitude angle estimates, regardless of whether single- or dual-frequency observations are used. This enhancement is the primary advantage of the AC A-PAD attitude model.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 4

Roll and pitch attitude angles estimated from dual-frequency (left column) and single-frequency (right column) observations for MC A-PAD (top row) and AC A-PAD (bottom row)

Given that the antennae are positioned on the top panel perpendicular to the Z-axis of the SRF, fixing ambiguities is not expected to have a substantial impact on the estimated yaw angle unless the CubeSat undergoes rotation around its Z-axis. This specific scenario is illustrated in Figure 5 for one of the testing cases, i.e., MC A-PAD with dual-frequency observations. The outcomes of the other tests exhibit similar patterns, further supporting the notion that fixing ambiguities primarily affects the roll and pitch angles. This result arises because the baselines rotate around the Z-axis simultaneously. Exploring alternative antenna configurations, such as configurations in which one antenna is placed on the side panels or bottom of the CubeSat, may help ascertain the impact of fixed ambiguity on the yaw angle. Such exploration is limited by the structural constraints of CubeSats and the potential for signal blockage from GNSS satellites and will be considered in a future study.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 5

Yaw attitude angle, with magnified results shown for a duration of approximately 20 min on the right

As mentioned earlier, the attitude angles in this CubeSat are a combination of magnetometer and sun sensor outputs. However, the low precision of these kinds of sensors, as indicated in Table 1 and evident from Figures 4 and 5, prevents them from being considered the best reference, in contrast to star trackers. Considering this limitation, we will apply the estimated attitude angles from the A-PAD model in the final POD test. A successful POD output implicitly validates the accuracy of the estimated attitude angles.

The final step involves evaluating the A-POD model. The reconstructed observations from Equation (16) are processed in a kinematic POD, as illustrated in Figure 2. A comparison between these orbits and the reference reduced-dynamic orbit from the Bernese software verifies that our A-POD model significantly enhances the accuracy of all orbital components (3D root mean square [RMS] = 4.1 cm), surpassing the results obtained when utilizing a single antenna (3D RMS = 11.8 cm) with quaternions. Table 5 provides a summary of the RMS values for the discrepancies between the orbit from the A-POD model and the reference orbit. It is important to highlight that this test also serves as an evaluation of the impact of the attitude matrix used in the estimated orbits. This evaluation involves applying the transformation in various parts of the POD process, including the consideration of antenna sensor offsets. By conducting this evaluation, we can gain insights into the performance and accuracy of the attitude matrix in POD. Future investigations will also explore the application of the attitude matrix in a wider range of operational situations, such as pointing operations and formation flying maneuvers.

View this table:

• View inline
• View popup

TABLE 5

RMS Values of the Differences Between Estimated Orbits and the Reference Orbit

The final stage of evaluation involves conducting an observation residual analysis. As depicted in Figure 6, the residuals of all observations in the A-POD model show significantly smoother behavior compared with conventional POD with a single antenna. The RMS value of the residuals for the A-POD model is 6 mm, which is approximately half the value obtained from conventional POD using one antenna. This substantial RMS reduction confirms the superior performance and optimality of the proposed A-POD model in comparison to typical POD processing. Importantly, the outcomes of the A-POD model remain consistent, regardless of whether an MC or AC A-PAD approach is utilized in deriving the observations outlined in Equation (16).

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 6

Observation residuals from conventional POD with one antenna (top) and A-POD with a full array (bottom); UTC: coordinated universal time

It is worth noting that the performance of the developed model will be affected if frequent signal interruptions and initializations occur, as the A-POAD concept relies on the kinematic POD approach, where disruptions in the collection of GNSS data have the potential to trigger a recalculation reset. Naturally, the GNSS observation noise and quality will affect any positioning method.