{'fraction_non_alphanumeric': 0.032818309510038834, 'fraction_numerical': 0.012310997273403288, 'mean_word_length': 4.628871732862059, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analytical Process Noise Covariance Modeling for Absolute and Relative Orbits

This paper develops new analytical process noise covariance models for both absolute and relative spacecraft states. Process noise is always present when propagating a spacecraft state due to dynamics modeling deficiencies. Accurately modeling this noise is essential for sequential orbit determination and improves satellite conjunction analysis. A common approach called state noise compensation models process noise as zero-mean Gaussian white noise accelerations. The resulting process noise covariance can be evaluated numerically, which is computationally intensive, or through a widely used analytical model that is restricted to an absolute Cartesian state and small propagation intervals. Moreover, mathematically rigorous, analytical process noise covariance models for relative spacecraft states are not currently available. To address these limitations of the state of the art, new analytical process noise covariance models are developed for state noise compensation for both Cartesian and orbital element absolute state representations by leveraging spacecraft relative dynamics models. Two frameworks are then presented for modeling the process noise covariance of relative spacecraft states by assuming either small or large interspacecraft separation. The presented techniques are validated through numerical simulations.


Introduction
In Kalman filtering, process noise represents deviations of the true dynamics from the modeled dynamics. Process noise is always present in orbit determination due to imperfect knowledge of the numerous complex forces acting on a spacecraft and limited computational resources for modeling those forces. Computational capabilities are especially constrained for onboard orbit determination, often necessitating reduced order dynamics models. Accurately capturing dynamics modeling deficiencies through modeled process noise is essential for optimal sequential estimation. Inaccurate process noise modeling increases estimation error and can result in filter inconsistency and divergence. 1, 2 Realistic process noise modeling improves the fidelity of the spacecraft state error covariance, which increases the accuracy of satellite conjunction analysis. In particular, analytical process noise models are advantageous because they are computationally efficient and provide insight into system behavior. This paper addresses gaps in the state of the art through the development of new analytical process noise covariance models for absolute and relative spacecraft states parameterized using both Cartesian coordinates and orbital elements. Notably, analytical process noise covariance models for absolute spacecraft states are obtained by exploiting relative spacecraft dynamics models to describe the state error dynamics.
A variety of approaches have been established to account for unmodeled accelerations in orbit determination. Here the term unmodeled accelerations refers to the difference between the true and modeled spacecraft accelerations, which is due to forces that are omitted in the modeled dynamics as well as errors in the modeled forces. As demonstrated by the TOPEX/Poseidon and Jason-1 missions, unmodeled accelerations can be approximated using trigonometric functions with fixed frequencies and estimated coefficients. 3,4 However, this technique is best suited to a posteriori high-precision orbit determination. 5 In batch estimation, it is also possible to estimate piecewise-constant empirical accelerations. 6,7 In Kalman filtering, a common ad hoc method is to manually tune a diagonal process noise covariance. Such an approach does not scale the process noise covariance according to the length of the propagation interval, so the covariance must be retuned for different propagation interval lengths. 8 Furthermore, this technique does not capture the cross-covariance between elements of the spacecraft state, and does not provide any dynamical consistency between the magnitudes of the process noise covariance of the different spacecraft state components. A second approach in Kalman filtering is to model the effect on the spacecraft state error covariance of specific forces. For example, the effects of auto-correlated central body gravity modeling errors of both commission and omission can be modeled. 9,10 However, these models are rarely used onboard due to computational cost. 8 The process noise due to maneuver execution errors can also be explicitly modeled. 8,11 Olson et al. 12 describe how consider covariance analysis techniques can be used to precompute a process noise covariance profile along a reference trajectory to capture the effects of known model parameter uncertainties on the spacecraft state error covariance. Although this approach is highly accurate provided that the assumed parameter uncertainties are representative and the actual spacecraft trajectory is close to the predefined reference trajectory, the numerical precomputation process is computationally expensive. Another approach called dynamic model compensation directly estimates the unmodeled accelerations by augmenting the estimated state vector with empirical accelerations. 13,14,15,16,17 Augmenting the state vector increases filter computation time, and the performance of dynamic model compensation depends on the tuning of parameters in the dynamical model of the empirical accelerations such as the empirical acceleration time correlation constants. 18,15 Advances in process noise modeling have primarily focused on an absolute Cartesian spacecraft state. However, the most appropriate state representation depends on the application. The spacecraft state can be parameterized as absolute or relative Cartesian coordinates as well as absolute or relative orbital elements. An absolute state describes the motion of a single spacecraft. In contrast, a relative state is a function of the absolute states of two spacecraft and describes the motion of one spacecraft relative to the other. Cartesian coordinates are most frequently used because of their more direct relation to the modeled measurements and because they do not have any singularities. Although orbital element states suffer from singularities, they have several advantages over Cartesian coordinates. Since most orbital elements vary more slowly in time than position and velocity, orbital elements can often be numerically integrated more efficiently. 19 Orbital elements also enable analytical perturbation models 20,21 and provide greater geometric insight than a Cartesian state. See Hintz 22 and Sullivan et al. 23 for comprehensive reviews of absolute and relative state parameterizations respectively.
One of the most widely used process noise models in sequential orbit determination is state noise compensation (SNC). 14,8,24,25 Although the unmodeled accelerations are generally correlated in time to some degree, SNC treats the unmodeled accelerations as zero-mean Gaussian white process noise. The resulting process noise covariance can be evaluated through numerical integration, which is computationally intensive. Alternatively, an analytical SNC model is commonly used that assumes kinematic motion with zero nominal spacecraft acceleration. However, this model is restricted to an absolute Cartesian state and is only valid for small propagation intervals. 8,24,25 The authors are not aware of analytical SNC models for absolute orbital element state representations or mathematically rigorous, analytical process noise covariance models for relative spacecraft states.
To overcome these limitations, this paper develops new analytical SNC process noise covariance models for both absolute and relative spacecraft states. The following section reviews relevant background information on SNC and spacecraft relative motion dynamical modeling. Exploiting relative dynamics models, the subsequent section develops new analytical process noise covariance models for absolute Cartesian and orbital element states by assuming either a circular orbit or a small propagation interval. These models are then extended to long propagation intervals for perturbed, eccentric orbits. The next two sections present frameworks for modeling the process noise covariance of relative spacecraft states by assuming either small or large interspacecraft separation.
The large interspacecraft separation framework is very general, and can leverage any absolute process noise covariance model, not just SNC models. The presented analytical models are validated through numerical simulations in the penultimate section. Finally, conclusions are provided based on the numerical results.

State Noise Compensation for an Absolute Spacecraft State
Let x α (t) ∈ R 6 denote the absolute state of a spacecraft parameterized in either Cartesian coordinates or a set of orbital elements. The specific state representation is indicated by the subscript α. The mean state estimate at time step i after processing all the measurements through time step j isx α,i|j . Assuming an unbiased estimator such that E[x α,i|j ] = x α (t i ), the error covariance or formal covariance of the mean state estimate at time step i after processing all the measurements through time step j is P An essential task in sequential orbit determination as well as satellite conjunction analysis is to propagate the error covariance to some future time t k , which involves determining P α,k|k−1 givenx α,k−1|k−1 and P α,k−1|k−1 . Note that t k can be any time of interest greater than t k−1 .
In general, the dynamical model of the spacecraft state is given by some nonlinear functionẋ α (t) = f (x α (t), u(t), t) where u(t) is the control input. The approach taken in a linearized framework such as an extended Kalman filter is to linearize the dynamical model over the propagation interval [t k−1 , t k ] about the mean state estimate taking into account all measurements through time t k−1 . This linearization results in a linear time-varying systeṁ where X α (t) = x α (t) −x α (t) is the state estimate error, andx α (t) is the mean state estimate at time t taking into account all measurements through time t k−1 such thatx α (t k ) =x α,k|k−1 . Here A α is the plant matrix, B α is the control input matrix, and Γ α is the process noise mapping matrix. These matrices are defined by the partial derivatives evaluated at the mean state estimatex α (t).
The continuous-time process noise describes stochastic deviations from the nominal spacecraft dynamics. If the dynamics were truly linear, the only sources of process noise would be numerical error and unmodeled accelerations. Throughout this paper physically represents unmodeled accelerations since they generally create orders of magnitude more process noise than numerical error. Here ∈ R 3 is modeled as a zero-mean white Gaussian process with autocovariance where δ(·) is the Dirac delta function and Q ∈ R 3×3 is the process noise power spectral density, which describes the strength of the unmodeled accelerations. This approach to process noise modeling is called SNC. 14,8 The stochastic vector of unmodeled accelerations can be inertial or radial-transverse-normal (RTN) coordinates, which will be denoted by the superscripts I and R respectively.
is the rotation matrix from the RTN frame to the inertial frame at time t. The power spectral density of I and R are Q I and Q R respectively. In the RTN frame, also referred to as the Hill frame, the radial and normal axes are aligned with the chief spacecraft position and angular momentum vectors respectively. The transverse axis completes the right-handed triad and is positive in the direction of the chief velocity. Typically, the spacecraft dynamics are highly nonlinear, and Q should be large enough to also accommodate errors due to dynamical nonlinearities in the propagation of the mean state estimate and associated error covariance.
∂xα(t k−1 ) evaluated at the mean state estimate is the state transition matrix. The discrete-time process noise is w α,k ∼ N (0, Q α,k ), and the process noise covariance is Due to the structure of Eq. (5), Q α,k is guaranteed symmetric and positive semi-definite provided that Q is symmetric and positive semi-definite. The propagated or time-updated error covariance is where The process noise covariance integral in Eq. (5) can be evaluated numerically provided that the state transition matrix and process noise mapping matrix are integrable functions. However, analytical solutions are desirable because they significantly reduce computation time and provide insight into system behavior. One common analytical approximation of Eq. (5) considers an absolute, inertial Cartesian state where r and v are the inertial position and velocity vectors respectively. Assuming kinematic motion where the nominal spacecraft acceleration is zero, where ∆t k = t k − t k−1 is the length of the propagation interval [t k−1 , t k ]. The identity matrix and matrix of zeros with three rows and columns are I 3×3 and 0 3×3 respectively. Modeling in the inertial frame, the process noise mapping matrix is Substituting Eqs. (9-10) and Q = Q I into Eq. (5) and then evaluating the integral yields 8, 25 Assuming the orientation of the RTN frame relative to the inertial frame is constant over a small propagation interval [t k−1 , t k ], the relation can be substituted into Eq. (11). 8 Another approach to simplify Eq. (5) is to model as constant over the interval [t k−1 , t k ]. 14 However, the covariance of must then be retuned for different propagation interval lengths as observed by Carpenter et. al. 8 Although the analytical process noise covariance model in Eq. (11) is widely used, it is only valid for small propagation intervals because it assumes kinematic motion with zero nominal spacecraft acceleration. 8,25 Furthermore, the authors are not aware of any analytical SNC models for orbital element states or mathematically rigorous, analytical process noise covariance models for relative spacecraft states. These gaps in the state of the art are addressed in this paper.

Hill-Clohessy-Wiltshire Equations
Many spacecraft relative dynamics models have been developed to describe the motion of a deputy spacecraft relative to a chief spacecraft. 23,21,20 Any relative spacecraft dynamics model can be used to model the dynamics of the state error by considering the estimated and true spacecraft states to be the chief and deputy spacecraft respectively as illustrated in Figure 1. One particularly simple relative dynamics model is the Hill-Clohessy-Wiltshire (HCW) equations, which assume two-body motion, a circular chief orbit, and small interspacecraft separation as compared to the chief orbit radius. The solution of the HCW equations was used by Geller 26 to approximate a scalar metric of the uncertainty of the position of one spacecraft relative to another spacecraft due to navigation and maneuver execution errors as well as umodeled accelerations. In the following sections, the solution of the HCW equations is used for the first time to construct an analytical model of the full process noise covariance of both absolute and relative Cartesian states. This section reviews the HCW equations, 27 which can be parameterized in either rectilinear or curvilinear coordinates, and the corresponding discrete-time solution. 28,29 Relative Spacecraft Coordinates Diagram Throughout this paper, δx α (t) ∈ R 6 denotes a relative spacecraft state where the specific relative state representation is denoted by the subscript α. Let be the state of a deputy spacecraft relative to a chief expressed in the chief RTN frame. The inertial Cartesian coordinates, defined in Eq. (8), of the chief and deputy are denoted by the subscripts c and d respectively. The chief orbit may be the orbit of either a physical or virtual spacecraft. Notice that for the relative velocity vector δv, the time derivative has been taken in the chief RTN frame.
Here δr = [δr r δr t δr n ] T and δv = [δv r δv t δv n ] T . The angular velocity vector of the chief RTN frame with respect to the inertial frame expressed in the chief RTN frame is denoted w. The cross product matrix of w is w × , which is a skew-symmetric matrix defined such that For two-body motion, w r = w t = 0 and w n =θ c whereθ c is the time rate of change of the angle between the chief position vector r c and an arbitrary fixed vector in the chief orbital plane.
Consider the spherical coordinates where r c = ||r c ||. Here ρ = ||r d || − r c , θ is the angle between the projection of the deputy position vector onto the chief orbital plane and the chief position vector, and φ is the angle between the deputy position vector and its projection onto the chief orbital plane. The angles θ and φ are measured positive in the directions of the chief velocity and angular momentum vectors respectively as shown in Figure 1. The nonlinear equations of relative motion parameterized using ρ, θ, and φ are derived by Willis et al. 30 considering only two-body motion.
Linearizing these equations about zero interspacecraft separation, assuming a circular chief orbit, and including differential perturbing accelerations yields 30,28 ρ = 3n 2 ρ + 2nr cθ + δd r (17) The chief mean motion is denoted n, and δd = d d − d c = [δd r δd t δd n ] T is the difference between the perturbing accelerations of the deputy and chief expressed in the chief RTN frame.
One particular curvilinear relative state is The discrete-time solution of Eqs. (17)(18)(19) yields the state transition matrix that propagates δx ψ from time 0 to t, which is the well-known solution of the HCW equations. The HCW equations can also be parameterized in the rectilinear coordinates δx R defined in Eq. (12), and the corresponding state transition matrix is identical to that shown in Eq. (21). Parameterizing the state in curvilinear coordinates instead of rectilinear coordinates has been shown to reduce errors incurred by the linearization when using the state transition matrix to propagate the relative state. 28, 29

Absolute Process Noise Covariance
In general, it is difficult to obtain an exact analytical solution of the integral in Eq. (5) that is reasonably concise even when only considering two-body motion. For an absolute Cartesian state, the complexity of the integral in Eq. (5) is principally due to the state transition matrix. The predominant approach in the literature is to simplify the Cartesian state transition matrix by assuming kinematic motion with zero nominal acceleration, which results in the analytical model in Eq. (11) that is valid for small propagation intervals. An alternate approach is developed in the following subsection by assuming a circular orbit and leveraging the solution of the HCW equations. The resulting analytical process noise covariance model is accurate for much longer propagation intervals for near-circular orbits than the widely used model in Eq. (11). In contrast to a Cartesian state, the complexity of the integral in Eq. (5) for an orbital element state is primarily due to the process noise mapping matrix, which is given by the Gauss Variational Equations. New analytical process noise covariance models for absolute orbital element states are derived in the following subsections by simplifying the process noise mapping matrix through the assumption of either a circular orbit or a short propagation interval. These absolute Cartesian and orbital element process noise covariance models are then extended to long propagation intervals for perturbed, eccentric orbits by exploiting spacecraft relative dynamics models.
Typically, it is advantageous to model spacecraft acceleration uncertainty in the RTN frame since the amount of acceleration uncertainty in each axis of the RTN frame generally varies less than in each axis of an inertial frame. For example, consider a case where atmospheric drag is the largest contributor to spacecraft acceleration uncertainty. Throughout the orbit, the acceleration uncertainty is consistently large in the along-track direction and relatively small in the radial and normal directions. In contrast, the level of acceleration uncertainty in each axis of an inertial frame will change significantly throughout the orbit. Consequently, the unmodeled acceleration power spectral density, Q, will be modeled in the RTN frame throughout the rest of this paper. For simplicity, it is assumed that the power spectral density matrix is diagonal such that Here Q r , Q t , and Q n describe the strength of the unmodeled accelerations in the radial, transverse, and normal directions respectively. A diagonal Q R enables the expression for Q α,k in Eq. (5) to be written as Here is the contribution to Q α,k due to unmodeled accelerations in axis i of the RTN frame where i ∈ {r, t, n}. The vectors Γ α,r (t k , τ ), Γ α,t (t k , τ ), and Γ α,n (t k , τ ) are the first, second, and third columns respectively of The particular values of Q r , Q t , and Q n are based on a priori knowledge of the dynamical environment. For example, these values can be selected in order to match elements of Q α,k to corresponding estimates 25 or in order for the filter error covariance to match an empirical approximation of the true error covariance. 8 Alternatively, the power spectral density can be estimated online using the recently developed adaptive state noise compensation algorithm. 31

Cartesian State and Circular Orbit
This section leverages the solution of the HCW equations to develop a process noise covariance model for the absolute, inertial Cartesian state of a spacecraft x I for near-circular orbits considering only two-body motion. The state error of an inertial Cartesian state X I (t) = x I (t) −x I (t) can be thought of as the state of a deputy relative to a chief spacecraft. The deputy is the true spacecraft state x I , and the chief is the estimated statex I , which is a virtual orbit (See Figure 1). This relative motion is described in a linearized framework through Eqs. (1) and (4). The relative motion can be equivalently parameterized using the curvilinear state δx ψ defined in Eq. (20). The parameterization of X I in curvilinear coordinates will be denoted δX ψ . Assuming an unbiased estimator, the expectation of the estimate error is zero (i.e., E[δX ψ (t)] = 0). The dynamics of δX ψ can be linearized about δX ψ = 0, resulting in a linear time-varying system similar to Eq. (1). Following the same approach used to derive Eqs. (4-7), the covariance of δX ψ is Here Φ ψ,k is the state transition matrix of δX ψ . Assuming two-body motion and a circular orbit, Φ ψ,k is equivalent to Eq. (21). In this case, n in Eq. (21) refers to the mean motion computed from the mean state estimate. The unmodeled accelerations are modeled in the RTN frame of the mean state estimate trajectory. The process noise mapping matrix is then deduced from Eqs. (17)(18)(19) to be The matrix Q ψ,k is obtained by substituting Eqs. (21)(22) and Eq. (28) into Eq. (5). The solution to Eq. (5) is computed by evaluating Eq. (24) and substituting the result into Eq. (23). Analytically evaluating Eq. (24) for each axis yields Each matrix in Eqs. (29)(30)(31) is symmetric, and only the unique elements in the lower triangular portion of each matrix are specified for brevity. The auxiliary variables in Eqs. (29)(30)(31) are Recalling the assumption of an unbiased filter such that E[X can be written in terms of its curvilinear parameterization through a first-order Taylor series expansion about X I (t k ) = 0 as where Here δx I = x d I − x c I is the difference between the absolute, inertial Cartesian states of a deputy and chief. Recall that in this case the chief and deputy are the estimated and true states respectively (See Figure 1). Eq. (34) can be expanded using the chain rule as The first matrix of partial derivatives in Eq. (35) is deduced from Eq. (12) to be Since the relative Cartesian and curvilinear coordinates are equivalent to first order at zero separation for a circular orbit, ∂ δx Using the linearization in Eq. (33), the time updated error covariance of the inertial Cartesian state can be written as Substituting Eq. (27) (7), it is clear that the process noise covariance of x I at time step k is it is equivalent to compute the process noise covariance integral parameterized in the rectlinear coordinates δx R , and then obtain Q I,k through a linear mapping similar to Eq. (41). Thus a higher order mapping from Q ψ,k to Q I,k is required to glean the benefits of spherical coordinates over rectilinear coordinates.

Orbital Element State and Circular Orbit
This section describes a framework for modeling the process noise covariance of an absolute orbital element state by assuming a circular orbit. The set of equinoctial elements is considered as an example, but the developed approach can be applied to any orbital element state representation. The equinoctial elements are defined as in terms of the classical Keplerian orbital elements a, e, i, Ω, w, and M . The state x E is singular for i = 180 • and e = 1. Considering only two-body motion, the state transition matrix of The time derivatives of the equinoctial elements are given by the Gauss Variational Equations, which can be written in matrix form as where d ∈ R 3 is the perturbing accelerations expressed in the RTN frame.
Treating the process noise as unmodeled accelerations expressed in the RTN frame, G(x E ) ∈ R 6×3 is the process noise mapping matrix. This matrix is 32, 33 whereĀ,B, . . . ,M are each a scalar function of x E , and their specific definitions are provided in Appendix A. The matrix Γ E contains trigonometric functions of the true longitude l = ν +Ω+w where ν is the true anomaly. For a near-circular orbit, the approximations f = g = 0 and can be made wheren is the average spacecraft angular rate over the propagation interval [t k−1 , t k ]. These approximations become exact as the orbit eccentricity approaches zero. Hereθ s is the time rate of change of the angle between the spacecraft position vector and an arbitrary fixed vector in the orbital plane. The variable ∆θ s,k = θ s (t k ) − θ s (t k−1 ) is the angle traversed by the spacecraft over the interval [t k−1 , t k ]. In the case of zero eccentricity,n k is equal to the mean motion. For other orbital element state representations, any quickly varying orbital elements such as the true anomaly or argument of latitude can be approximated similar to Eq. (46).
The process noise covariance of the state x E at time step k, Q E,k , is computed by evaluating Eq. (24) and substituting into Eq. (23). Analytically evaluating Eq. (24) for each axis using Eq. (46) and assuming f = g = 0 yields M n = 0 0 0 1 2 Here p = a(1−e 2 ) is the orbit semi-parameter, L = √ µp is the magnitude of the specific angular momentum vector, and l k = l(t k ). Since the orbit is assumed circular, the semi-parameter can be replaced with the semi-major axis.

Orbital Element State and Small Propagation Interval
A concise analytical approximation of Eq. (5) that is valid for eccentric orbits can be obtained by assuming that the propagation interval [t k−1 , t k ] is small such that any quickly varying orbital elements can be approximated as constant in Γ α over that interval. Then Γ α is constant when considering only two-body motion. In the case of the equinoctial elements, it is assumed that the true longitude is constant over the interval [t k−1 , t k ]. For a constant Γ E and the state transition matrix in Eq. (43), the contributions to Q E,k in each axis shown in Eq. (24) evaluate to The first, second, and third columns of Γ E are denoted Γ E,r , Γ E,t , and Γ E,n respectively. Eq. (54) is evaluated using the true longitude halfway through the propagation interval [t k−1 , t k ] at time (t k−1 + t k )/2. The true longitude at time (t k−1 + t k )/2 can be approximated through two-body propagation of the true longitude from time t k−1 .
Both of the developed equinoctial element process noise covariance models provide important insights. For example, notice from Eqs. (44-45) that the time evolution of the equinoctial elements h and k are only affected by perturbing accelerations in the normal direction. As a result, only unmodeled accelerations in the normal direction contribute to the process noise covariance of h and k as is shown in the model described by Eqs. (48-53) and the model in Eq. (54). In both of these models, the process noise covariance of h is small when the true longitude is near l = (2n−1)π 2 throughout the interval [t k−1 , t k ] for any integer n. This occurs because perturbing accelerations have no effect on h at l = (2n−1)π 2 sinceĪ in Eq. (45) is zero. Similarly, the process noise covariance of k is small when the true longitude is near l = nπ throughout the interval [t k−1 , t k ].

Long Propagation Intervals for Perturbed, Eccentric Orbits
In the case of a long propagation interval for an orbit that is eccentric or highly perturbed, the analytical models in the previous sections may not provide sufficient accuracy. This section describes a methodology for such cases. The Using Eq. (7) to sequentially propagate the error covariance over each subinterval, it is evident that Eq. (5) can be equivalently written as is the process noise covariance computed over the subinterval [t i−1 ,t i ].
Formulating the process noise covariance according to Eq. (55) is advantageous because an analytical approximation that is valid over small propagation intervals can be used for each Q α (t i ,t i−1 ) provided that the length of each subinterval is sufficiently small. Then an analytical state transition matrix Φ α (t k ,t i ) with the accuracy required for the specified application can be used to map the effects of each Q α (t i ,t i−1 ) on Q α,k . In this way, each term of the sum in Eq. (55) may neglect certain effects such as perturbations or orbit curvature over the subinterval [t i−1 ,t i ], but these effects can be accounted for throughout the remainder of the propagation interval [t i , t k ] through the selected Φ α (t k ,t i ). Furthermore, Eq. (55) provides the flexibility to use different values of Q for different subintervals. Eq. (55) can also be used to extend many process noise covariance models from literature to perturbed, eccentric orbits such as the dynamic model compensation process noise covariance model developed by Cruickshank,13 which assumes a small propagation interval. 8 Since Φ α (t k ,t i ) describes the evolution of the state estimate error (i.e., the true state relative to the mean state estimate), it can be formulated using any relative spacecraft state transition matrix (See Figure 1). Specifically, where the partial derivatives are evaluated at zero interspacecraft separation, δx α = x d α − x c α is the difference between the states of a deputy and chief, and ∂δx β (t k−1 ) evaluated at zero separation is the state transition matrix of some relative state δx β . A comprehensive review of state transition matrices for relative states available in literature is provided by Sullivan et al. 23 including models for orbits of any eccentricity and models that incorporate a variety or perturbations. Additional state transition matrices for relative orbital elements are provided by Koenig et al. 21 as well as Guffanti and D'Amico 20 that include the dominant effects of earth oblateness, atmospheric drag, and solar radiation pressure. In order to compute each Φ α (t k ,t i ) and Q α (t i ,t i−1 ) in Eq. (55), the spacecraft state at eacht i is required. If the reference trajectory or mean state estimate is numerically integrated, thet i can be chosen as times where the state is already available from the reference trajectory numerical integration.
As an example application of Eq. (55), consider a spacecraft in a nearcircular orbit and a propagation interval that is multiple orbits. To compute each Q α (t i ,t i−1 ), Eqs. (29)(30)(31)(32) and (41) can be used for a Cartesian state or a model similar to Eqs. (48-53) can be used for an orbital element state. The subintervals [t i−1 ,t i ] could be about one orbit or less in length since, as will be shown later, the effects of perturbations on the process noise covariance are relatively small over one orbit. Then an analytical Φ α (t k ,t i ) can be employed in Eq. (55) that accounts for perturbations to a sufficient degree to satisfy mission requirements.
As a second example, consider a spacecraft in an eccentric orbit and a propagation interval that is about one orbit. For an equinoctial element state, each Q α (t i ,t i−1 ) can be modeled through Eq. (54). Since the effects of perturbations are typically small over one orbit, the two-body state transition matrix in Eq. (43) can be used for each Φ α (t k ,t i ). For an inertial Cartesian state, each Q α (t i ,t i−1 ) can be modeled through Eq. (11). The state transition matrix can be computed through Eq. (57) using any of the solutions of the Tschauner-Hempel equations such as the Yamanaka-Ankerson 34 state transition matrix. This matrix is formulated for a state δx R = [δ r T δ v T ] T that is a normalized form of the relative state δx R defined in Eq. (12) where the time derivative is taken with respect to true anomaly. The transformations between these two states are 30 The required partial derivatives in Eq. (57) are evaluated at zero separation where and ∂δx I ∂δx R is defined in Eq. (36). For near-circular orbits, true anomaly and argument of periapsis are poorly defined and can undergo abrupt changes due to perturbations. Thus, it is advantageous to use the Yamanaka-Ankerson state transition matrix formulated using the argument of latitude u = w + f developed by Willis and D'Amico. 35 The substitutions e cos f = e x cos u + e y sin u and e sin f = e x sin u − e y cos u can also be made in Eqs. (58) and (61-62) where the eccentricity vector components e x = e cos w and e y = e sin w are well defined for near-circular orbits. For both Cartesian coordinates and equinoctial elements, the accuracy of the resulting Q α,k increases as the length of the subintervals decreases since Eq. (11) and Eq. (54) assume a small propagation interval. Thus, the number of subintervals N is chosen to balance accuracy and computational cost.

Relative Process Noise Covariance For Small Separations
If the interspacecraft separation is small compared to the orbit radius, then the nonlinear dynamical model of the relative state δx α can be linearized about zero interspacecraft separation over the propagation interval [t k−1 , t k ], resulting in a linear time-varying system The plant matrix, control input matrix, and process noise mapping matrix of the relative state are evaluated at zero interspacecraft separation. The subscript δ indicates a matrix corresponds to a relative state, and the particular relative state is specified by the subscript α. Because the system is linearized about zero separation, the chief and deputy RTN frames are identical for the nominal or reference trajectory. For consistency, the RTN frame will be referred to as that of the chief throughout this section. The differential control inputs δu = u d − u c are the difference between the control inputs of the deputy and chief expressed in the chief RTN frame. Similarly, the differential unmodeled accelerations δ = d − c are the difference between the deputy and chief unmodeled accelerations expressed in the chief RTN frame with autocovariance is the power spectral density of δ . The power spectral densities of the chief and deputy unmodeled accelerations are Q c and Q d respectively. The cross power spectral densities are Q dc and Q cd . Following a similar derivation to that of Eqs. (5 -7), the propagated error covariance of δx α for the linearized system in Eq. (64) is where the process noise covariance of δx α at time step k is The state transition matrix Φ δα,k = Φ δα (t k , t k−1 ) propagates the relative state δx α from time t k−1 to time t k . Due to the structure of Eq. (69), Q δα,k is guaranteed symmetric and positive semi-definite when Q δ is symmetric and positive semi-definite. Throughout this paper, Q δ is expressed in the chief RTN frame, denoted Q R δ . For simplicity, it is assumed that Q R δ is diagonal where

Relative Cartesian State
This section models the process noise covariance of the relative state δx R defined in Eq. (12). First the process noise covariance is constructed for the relative state δx I = x d I − x c I , which is the difference between the absolute, inertial Cartesian states of the deputy and chief. Linearizing the system about zero interspacecraft separation, the state transition matrix and process noise mapping matrix match those of the absolute inertial Cartesian state estimate error X I (t) = x I (t) −x I (t). As a result, the expression for the process noise covariance of δx I , denoted Q δI,k , is equivalent to Q I,k except that Q δI,k considers relative unmodeled accelerations. Thus, Q δI,k can be approximated using Eq. (11) for a small propagation interval, Eqs. (23) and (29)(30)(31)(32) for a near-circular orbit, or Eq. (55) by making the substitutions Q I,k = Q δI,k and Q R = Q R δ where Q R δ is defined in Eq. (70). Then, similar to Eqs. (37)(38)(39)(40)(41), the process noise covariance of δx R is is given in Eq. (63).

Relative Orbital Element State
A framework for modeling the process noise covariance of a relative orbital element state is demonstrated using the set of relative orbital elements 21 as an example. However, this approach can be applied to any set of relative orbital elements. The relative state in Eq. (72) is a function of the equinoctial orbital elements of the chief and deputy, which are defined in Eq. (42). Consider a similar set of relative orbital elements where x c E and x d E are the equinoctial elements of the chief and deputy respectively. Linearizing the system about zero interspacecraft separation and considering only two-body motion, the state transition matrix of δx E matches that shown in Eq. (43). Since the system is linearized about zero interspacecraft separation, the semi-major axis and mean motion in the state transition matrix can refer to the chief or deputy. For consistency, the chief parameters will be used throughout this section whenever there is an arbitrary choice between the parameters of the chief and deputy.
The time derivative of δx E can be expanded using the chain rule as The partial derivatives of δx E with respect to Substituting Eqs. (44) and (75) into Eq. (74) and assuming zero interspacecraft separation yields where δd = d d − d c is the difference between the deputy and chief perturbing accelerations expressed in the chief RTN frame. The matrix c→d R ∈ R 3×3 rotates vectors from the chief RTN frame to the deputy RTN frame and is the identity matrix for zero interspacecraft separation. Thus the process noise mapping matrix is Γ δE = G(x c E ), which matches that of the absolute equinoctial orbital elements shown in Eq. (45).
Since the state transition matrix and process noise mapping matrix of δx E match that of the absolute equinoctial elements, the expressions for Q δE ,k and Q E,k are equivalent except that the expression for Q δE ,k considers differential unmodeled accelerations. Thus, Q δE ,k can be analytically approximated through Eq. (23) and Eqs. (48-53) for a near-circular chief, Eqs. (23) and (54) for a small propagation interval, or Eq. (55) by making the substitutions Q E,k = Q δE ,k and Q R = Q R δ . Then following the approach in Eqs. (37)(38)(39)(40)(41), the process noise covariance of δx E is

Relative Process Noise Covariance for Large Separations
If the interspacecraft separation is large, the linearization about zero interspacecraft separation utilized in the previous section may introduce significant errors. Instead of assuming small separations, this section constructs the process noise covariance of a relative spacecraft state by assuming that the interspacecraft separation is large enough that the unmodeled accelerations of the chief are uncorrelated with those of the deputy. This framework is flexible in that the relative state process noise covariance can be formed using any absolute state process noise covariance model, not just SNC models.
In general, where f (·) is some function. The absolute states of the deputy and chief are x d β and x c β respectively where β denotes the specific absolute state representation. Eq. (80) can be approximated through a first-order Taylor series expansion as which is valid for small deviations of the estimated chief and deputy absolute states from their corresponding true states. Here the partial derivatives are evaluated at the mean state estimate. Using the linearization in Eq. (81), the propagated error covariance of δx α can be written as where is process noise covariance of δx α at time step k. The discrete-time process noise of the absolute chief and deputy states are w c β,k and w d β,k respectively with associated covariances Q c β,k and Q d β,k . In the particular case of SNC, where the superscripts c and d indicate matrices corresponding to the chief and deputy states respectively. Recall that Q cd is the cross power spectral density of the chief and deputy unmodeled accelerations where E[ c (t) d (τ ) T ] = Q cd δ(t−τ ). Note that c and d are expressed in the RTN frames of the chief and deputy respectively. If the interspacecraft separation is large, it can reasonably be assumed that the chief unmodeled accelerations are uncorrelated with those of the deputy such that E[w c β,k w d T β,k ] = 0. Using this assumption, Eq. (85) simplifies to Eq. (87) can be applied regardless of how the process noise of the chief and deputy states are modeled provided that E[w c β,k w d T β,k ] is small. Thus, Q d β,k and Q c β,k can be computed using any process noise covariance models for absolute spacecraft states, not just SNC models. Furthermore, Eq. (87) guarantees Q δα,k is positive semi-definite when Q d β,k and Q c β,k are positive semi-definite. As an example, consider the relative Cartesian state δx R defined in Eq. (12), which is a function of the absolute, inertial Cartesian states of the chief and deputy denoted x c I and x d I respectively. The process noise covariance of δx R can be modeled through Eq. (87) where the process noise covariances of the chief and deputy inertial Cartesian states, Q c I,k and Q d I,k , are separately computed using any absolute process noise covariance model such as Eq. (11) for a small propagation interval, Eqs. (23) and (29-32) for a near-circular orbit, or Eq. (55). The required partial derivatives in Eq.
As another example, the process noise covariance of the relative orbital element state δx E defined in Eq. (72) can be modeled by Eq. (87) where Q c E,k and Q d E,k are each computed through any absolute process noise covariance model such as Eqs. (23) and (48-53) for a near-circular orbit, Eqs. (23) and (54) for a small propagation interval, or Eq. (55). In this case, the required partial derivatives in Eq. (87) are , which is equivalent to Eq. (79).

Numerical Validation
This section validates the developed analytical process noise covariance models by comparing them to numerical solutions of the integrals they approximate for Earth-orbiting spacecraft trajectories. The truth orbit numerical integration scheme and perturbation models are summarized in Table 1. These perturbation models and integration scheme are also used to numerically evaluate the process noise covariance as well as the required state transition matrices using the relations These numerical solutions of the process noise covariance are considered the reference truth. The simulated spacecraft mass is 100 kg, cross-section area is 1 m 2 , drag coefficient is 1, and radiation pressure coefficient is 1.2. Results are also provided for a Keplerian reference truth where the perturbations in Table  1 are omitted. The presented results for a Keplerian reference truth hold for any orbit semi-major axis because the time intervals associated with the results are normalized by the orbit period. The accuracy of the analytical process noise covariance models is quantified using two error metrics. The first error metric is (90) which is the smallest propagation interval length for which the absolute fractional error in any modeled process noise standard deviation is greater than or equal to 0.1. The element in the i th row and column of the process noise covariance as determined numerically and analytically are Q ii andQ ii respectively. The second error metric is the maximum absolute fractional error in any modeled process noise standard deviation for any propagation interval length up to one orbit, which is defined as

Absolute Process Noise Covariance
The analytical process noise covariance models of the absolute Cartesian state x I and equinoctial element state x E are each compared against the corresponding numerical solution of Eq. (5). While the initial osculating eccentricity is varied, the orbit inclination and right ascension of the ascending node are each 45 • . The initial mean longitude is 90 • . The initial periapsis radius is fixed at 6900 km in order to ensure the trajectory is feasible and does not pass through the Earth for any of the considered eccentricities. Thus, the semi-major axis and orbit period vary depending on the eccentricity. The small periapsis radius leads to large perturbing accelerations due to nonspherical gravity and atmospheric drag. The spacecraft trajectory nominally starts at periapsis, although periapsis is not well defined for small eccentricities, in order to maximize the effects of eccentricity.
The absolute power spectral density is as shown in Eq. (22). In one scenario, Q r , Q t , and Q n are all equal to a single value denoted Q * . In a second scenario, Q r = Q n = Q * and Q t = 10 Q * , which is representative of a case where atmospheric drag is the largest source of unmodeled accelerations. Notice in Eq. (5) that each element of the process noise covariance is a linear function of the elements of the power spectral density. Since the nonzero elements of the power spectral density are each a linear function of Q * , each process noise standard deviation is some scalar times ( Q * ) 1 2 for both the numerical and considered analytical solutions. Thus the fractional error of each modeled process noise standard deviation is the same regardless of the simulated value of Q * .
The new model for equinoctial elements in Eq. (54) and the widely used model for an inertial Cartesian state in Eq. (11) both assume a small propagation interval ∆t k . To determine the range of propagation interval lengths for which these models can be accurately applied, ∆t min is plotted in Figure 2 as a function of eccentricity. The fractional error in the process noise standard deviation of the equinoctial element h is as large as 0.14 for ∆t < 0.13 orbits. However, since the initial true longitude is 90 • , the process noise standard deviation of h is very small for ∆t < 0.13 orbits. Thus, the absolute error is very small, and the process noise standard deviation of h for ∆t < 0.13 orbits is neglected in Figure 2. The spacecraft angular velocity is greatest at periapsis and grows with increasing eccentricity, more quickly breaking the assumptions of the models in Eqs. (54) and (11). Thus, ∆t min tends to decrease for larger eccentricities in Figure 2 for both models, and both models generally perform better when the propagation interval is not near periapsis. The equinoctial element model outperforms the Cartesian model. Interestingly, ∆t min for the Cartesian model decreases significantly when the noise is predominantly in the transverse axis. The results in Figure 2 are nearly identical to when the reference truth is computed assuming Keplerian motion.  The new analytical models for absolute equinoctial elements and Cartesian coordinates shown in Eqs. (48-53) and Eq. (41) respectively assume a circular orbit. In order to quantify the maximum orbit eccentricity for which these models can be applied, δ max is plotted in Figure 3 as a function of eccentricity for the case of equal noise in each axis. The results of the transverse dominant noise scenario were similar. The equinoctial model is valid for larger eccentricities than the Cartesian model. For example, δ max is less than 0.1 for eccentricities less than about 0.09 and 0.02 for the equinoctial and Cartesian models respectively for both the Keplerian and perturbed reference truths. Since δ max is the largest fractional error over a full orbit, the fractional errors are smaller on average. Depending on the particular application, δ max ≤ 0.1 may provide sufficient accuracy. The presented model for equinoctial elements in Eq. (54) and the widely used kinematic model for an inertial Cartesian state in Eq. (11) are also included in Figure 3 for reference. Since both of these models assume a small propagation interval and δ max considers intervals up to an entire orbit, it is not surprising that δ max is large for these models. Figure 3 also shows the performance of Eq. (55), which breaks the propagation interval into smaller subintervals, applied to absolute inertial Cartesian coordinates and equinoctial elements using either four or eight subintervals. For Cartesian coordinates, Q α (t i ,t i−1 ) is computed using Eq. (11), and Φ α (t k ,t i ) is modeled using the partial derivatives in Eqs. (61-62) and the Yamanaka-Ankerson state transition matrix formulated using argument of latitude. 35 For equinoctial elements, Eq. (54) is used to compute Q α (t i ,t i−1 ), and Eq. (43) is used for Φ α (t k ,t i ). For simplicity, the subintervals are evenly spaced in time. However, further accuracy could likely be achieved with the same number of subintervals by concentrating more of the subintervals toward the beginning of the propagation interval since the earlier terms of the sum in Eq. (55) tend to be larger due to the propagation of Q α (t i ,t i−1 ) over a larger interval [t i , t k ]. For eccentric orbits, it may also be beneficial to concentrate more of the subintervals near periapsis since that is where the spacecraft angular velocity is greatest, more quickly violating the assumptions of Eqs. (11) and (54). As a worst case scenario, the states at each timet i , which are required to compute each Φ α (t k ,t i ) and Q α (t i ,t i−1 ) in Eq. (55), are computed through a Keplerian propagation of the state from the beginning of the propagation interval. Figure 3 shows that increasing the number of subintervals N increases accuracy. The equinoctial elements model is more accurate than the Cartesian model given the same number of subintervals because Eq. (54) is accurate over longer intervals than Eq. (11) as shown in Figure 2. Even for the considered highly perturbed orbits, the models that use Eq. (55) and neglect perturbations provide sufficient accuracy for many applications for propagation intervals of an orbit or less. The accuracy of these models can be further improved by using Φ α (t k ,t i ) that incorporate perturbations in Eq. (55), which is important for propagation intervals of multiple orbits.

Relative Process Noise Covariance
Consider two spacecraft where the chief has an initial osculating eccentricity of 10 −3 and an inclination and right ascension of the ascending node of 45 • . The initial chief osculating equinoctial elements are [a, f, g, h, k, λ] = [6900 km, 0, 10 −3 , 0.2929, 0.2929, 135 • ]. The initial deputy mean equinoctial elements all match those of the chief with the exception of the mean longitude, which is varied. As a result, the two spacecraft are primarily separated in the along-track direction as described by the initial mean relative mean longitude, δλ 0 , and the secular drift in δλ is minimal. The transformation between mean and osculating orbital elements is done using the first-order J 2 mapping presented by Schaub and Junkins. 37 The unmodeled acceleration power spectral density for each spacecraft is as shown in Eq. (22) where each diagonal element is equal to Q * for both spacecraft. The presented results hold for any Q * . In reality, the degree of correlation between the unmodeled accelerations of two spacecraft depends on the fidelity of the modeled dynamics as well as the orbit geometry. Here, the cross power spectral density of the chief and deputy unmodeled accelerations is modeled as Q cd = (0.9e −δλ0/10 −3 ) Q c . As a result, Q cd decays exponentially as δλ 0 increases. The rate of decay is dictated by the constant 10 −3 , which was chosen for illustrative purposes. The constant 0.9 models the fact that the unmodeled accelerations of the two spacecraft are not perfectly correlated near zero separation due to small differences in attitude and spacecraft physical properties. The reference truth process noise covariance of the relative states δx R and δx E are computed through Eq. (85). For δx R , the partial derivatives ∂δx R (t k ) ∂x c I (t k ) in Eq. (85) are approximated through central finite differences using the exact angular velocity vector of the RTN frame with respect to the inertial frame expressed in the RTN frame w = [r c (r c ·n)/L 0 L/r 2 c ] T , which holds for arbitrary spacecraft motion. 38 The reference truth process noise covariances of the absolute chief and deputy states as well as their cross covariance, which are employed in Eq. (85), are obtained through numerical integration of Eqs. (5) and (86) according to Table 1 both with and without perturbations.
Two approaches were taken to develop analytical process noise covariance models for relative spacecraft states by assuming either small or large interspacecraft separation. For the process noise covariance of δx R , the absolute model in Eq. (41) is applied in Eq. (71) for the small separation framework and Eq. (87) for the large separation framework. The analytical partial derivatives in Appendix B are used in Eq. (87). For the process noise covariance of δx E , the absolute model in Eqs. (48-53) is applied in Eq. (78) for the small separation framework and Eq. (87) for the large separation framework. Since the employed absolute models in Eq. (41) and Eqs. (48-53) assume circular orbits, the modeled relative process noise covariance incurs some error due to eccentricity whether the small or large separation approach is applied. The large separation framework errors increase for small separations due to the assumption that the chief unmodeled accelerations are uncorrelated with those of the deputy. On the other hand, the small separation framework fully accounts for the correlation between the chief and deputy unmodeled accelerations as shown in Eq. (66). However, the small separation framework errors grow as interspacecraft separation increases due to the linearization about zero interspacecraft separation.
In order to quantify these errors, δ max is plotted in Figure 4 for the considered analytical models. These results demonstrate the inherent trade-off in the large and small separation frameworks, which are each valid for a different range of interspacecraft separations. Most notably, the small separation framework can be used for much larger interspacecraft separations for equinoctial elements than for Cartesian coordinates. Specifically, δ max ≤ 0.1 for the Cartesian and equinoctial small separation models for δλ ≤ 2 × 10 −3 and δλ ≤ 1 × 10 −1 respectively for both the Keplerian and perturbed reference truths. These angles can be transformed to approximate interspacecraft distances of 13.8 km and 690 km respectively through multiplication with the orbit semi-major axis. Even for the perturbed reference truth, δ max < 0.013 can be achieved for any δλ 0 for an equinoctial element state by using either the small or large separation framework. On the other hand, the lowest δ max that can be achieved for a Cartesian state using the small and large separation frameworks is as high as 0.1 for the considered δλ 0 .

Conclusions
Accurate process noise modeling is essential for optimal orbit determination in a discrete-time Kalman filtering framework and can improve satellite conjunction analysis. A common approach to process noise modeling called state noise compensation (SNC) treats the process noise as zero-mean Gaussian white noise unmodeled accelerations. The resulting process noise covariance can be evaluated numerically. However, analytical solutions are desirable because they are computationally efficient and provide insight into system behavior. One analytical SNC model is widely used that assumes kinematic motion, but it is restricted to an absolute Cartesian state and small propagation intervals. Furthermore, analytical SNC models for absolute orbital element states are not currently available and there has been little work on process noise covariance modeling for relative spacecraft states. This paper addresses these gaps in the state of the art.
First, a new analytical SNC process noise covariance model is developed for an absolute, inertial Cartesian spacecraft state by leveraging the well-known solution of the Hill-Clohessy-Wiltshire equations. For small orbit eccentricities, this new model can be used for significantly longer propagation intervals than the widely used kinematic model. An approach is presented to obtain analytical SNC process noise covariance models for absolute orbital element states for small eccentricities, and then a second approach is presented that is valid for eccentric orbits over small propagation intervals. The developed absolute Cartesian and orbital element models are then extended to long propagation intervals for perturbed, eccentric orbits by leveraging spacecraft relative dynamics models. Subsequently, two frameworks are established for modeling the process noise covariance of relative spacecraft states by assuming either small or large interspacecraft separations. The large interspacecraft separation framework is particularly flexible in that it can leverage any absolute state process noise covariance model, not just SNC models such as those developed in this paper. All the presented process noise covariance models are analytical, making them much more computationally efficient than numerical solutions, and are guaranteed positive semi-definite. Overall, the orbital element models tend to outperform the Cartesian models. Nevertheless, the newly developed analytical process noise covariance models for absolute and relative spaceraft states provide sufficient accuracy for many applications for both Cartesian and orbital element state representations.
wherer,t, andn are the unit vectors of the RTN frame expressed in the inertial frame. The angular velocity vector of the RTN frame with respect to the inertial frame expressed in the RTN frame is w = [r c (r c ·n)/L 0 L/r 2 c ] T , which holds for arbitrary spacecraft motion. 38 Neglecting out-of-plane accelerations, which are typically very small compared to the in-plane accelerations, the angular velocity vector is w = [0 0 L/r 2 c ] T . This approximation avoids any partial derivatives of the chief inertial accelerationr c . The chief specific angular momentum vector is denoted L = r c ×v c , and its magnitude is L = ||L||. Now the partial derivatives in Eq. (B.1) can be deduced from Eq. (12), which yields Recall that the superscript × indicates a cross-product matrix as illustrated in Eq. (13). The difference between the deputy and chief inertial position and velocity vectors are δr I = r d − r c and δv I = v d − v c . The auxilliary variables K 1 , K 2 ∈ R 3×3 are In order to account for out-of-plane accelerations, the partial derivatives of ([ rc L (r c ·n) 0 0] T ) × R I→R δr I with respect to r c and v c would be subtracted from the right hand sides of Eqs. (B.5) and (B.6) respectively.
in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. This work is also supported by the NASA Small Spacecraft Technology Program cooperative agreement number 80NSSC18M0058 for contributions to the Autonomous Nanosatellite Swarming (ANS) Using Radio-Frequency and Optical Navigation project. Additionally, the authors wish to thank the Achievement Rewards for College Scientists (ARCS) Foundation for their support.