Relative navigation for autonomous formation flying satellites using the state-dependent Riccati equation filter

Available online at www.sciencedirect.com

ScienceDirect Advances in Space Research xxx (2015) xxx–xxx www.elsevier.com/locate/asr

Relative navigation for autonomous formation flying satellites using the state-dependent Riccati equation filter Han-Earl Park ⇑, Young-Rok Kim Korea Astronomy & Space Science Institute, Daejeon 34055, South Korea Received 30 May 2015; received in revised form 5 October 2015; accepted 6 October 2015

Abstract A relative navigation method for autonomous formation flying using the state-dependent Riccati equation filter (SDREF) is presented. In the SDREF, nonlinear relative dynamics, including J2 perturbation, are parameterized into a state-dependent coefficient (SDC) form without any loss of nonlinearity. The relative navigation algorithm is established based on the carrier-phase differential GPS (CDGPS) and single-frequency GPS data, in which the SDREF is used as a nonlinear estimator. To evaluate the SDREF performance, two different extended Kalman filters (EKFR1 and EKFR2) are introduced. The dynamic models of all the filters are based on relative motion including J2 perturbation. However, the SDREF and the EKFR1 use linear state propagation, whereas EKFR2 employs nonlinear state propagation. The navigation simulation is performed for each filter using live GPS signals simulated by a GPS signal generator, and the result is analyzed in terms of estimation accuracy and computational load. As a result, the SDREF provides a relative navigation solution with 3-D RMS accuracies of 6.0 mm and 0.153 mm/s for position and velocity, respectively, for a separation of 50 km with a computation time of approximately 34 s. The simulation results demonstrate that the SDREF estimates the relative states as rapidly as the EKFR1 and as accurately as the EKFR2, which means that the developed SDREF combines the strong points of EKFR1 and EKFR2 and overcomes their disadvantages. Ó 2015 COSPAR. Published by Elsevier Ltd. All rights reserved.

Keywords: State-dependent Riccati equation filter; Relative navigation; Spacecraft formation flying

1. Introduction For spacecraft formation flying, it is essential to determine precisely the position and velocity of a deputy satellite relative to a chief satellite. In particular, it is important to precisely estimate the velocity because accurate knowledge of the velocity reduces orbit control error (Tillerson et al., 2002), which thereby increases fuel efficiency and the mission lifetime. Additionally, it is desirable to autonomously estimate the relative position and velocity using onboard sensors to conduct a mission in real-time or promptly cope with unpredictable situations (i.e., collision avoidance). ⇑ Corresponding author. Tel.: +82 42 869 5812; fax: +82 42 861 5610.

E-mail address: [email protected] (H.-E. Park).

GPS sensors can be applied advantageously to such an autonomous relative navigation system with inter-satellite communication links (Montenbruck et al., 2002). Differential processing using GPS carrier-phase measurements allows for precise estimation of the relative position and velocity. The estimation performance is closely related to the characteristics of the filter. For spacecraft formation flying, the extended Kalman filter (EKF) is one of the most important optimal tools, particularly in real-time applications. Many studies in this field have been conducted using the EKF (Leung and Montenbruck, 2005; Gill et al., 2007; Ardaens et al., 2009). The procedure of the discrete-time EKF can be essentially divided into a prediction step and an update step. During the prediction step, the state and error covariance

http://dx.doi.org/10.1016/j.asr.2015.10.009 0273-1177/Ó 2015 COSPAR. Published by Elsevier Ltd. All rights reserved.

Please cite this article in press as: Park, H.-E., Kim, Y.-R. Relative navigation for autonomous formation flying satellites using the state-dependent Riccati equation filter. Adv. Space Res. (2015), http://dx.doi.org/10.1016/j.asr.2015.10.009

2

H.-E. Park, Y.-R. Kim / Advances in Space Research xxx (2015) xxx–xxx

are propagated forward in time. For the state propagation, both a linear state propagation using a state transition matrix and a nonlinear state propagation using nonlinear dynamics equations with numerical integration are available (Alfriend et al., 2010). The former has a small computational burden but reduces the estimation accuracy in a highly nonlinear system. On the contrary, the latter requires a relatively large computational effort but can estimate the parameters more precisely. Therefore, we need to determine which approach to use by considering the trade-off between the computational burden and the estimation accuracy according to the objective of the application. It is worthwhile to note that the computational burden is no less important than the estimation accuracy in practical applications, particularly when using an onboard system with limited resources available and requiring real-time implementation, such as autonomous formation flying using small spacecraft. In the current study, the state-dependent Riccati equation filter (SDREF) is proposed as an alternative, which may satisfy both the computational burden and estimation accuracy. The SDREF was first introduced by Mracek et al. (1996) and is based on the state-dependent Riccati equation (SDRE) technique (Cloutier, 1997). In the field of spacecraft formation flying, the SDRE has been continuously used as a nonlinear controller for relative orbit and/ or attitude (Chang et al., 2009; Park et al., 2011; Lee et al., 2012; Massari and Zamaro, 2014). The SDREF is also a nonlinear filter, which estimates states in a similar methodology to the EKF. The only difference is that the EKF is based on linearization, whereas the SDREF is based on a parameterization that transforms the nonlinear system to a linear structure. This allows us to easily deal with a nonlinear system as if it were a linear system. Moreover, there are no linearization effects, and the Jacobian matrix is also not required in the prediction step (Wang, 2011). As mentioned above, to use the SDREF, a nonlinear system must be converted into a linear structure that has a state-dependent coefficient (SDC) form, which is referred to as SDC parameterization or SDC factorization. Note that in a multivariable case, there are an infinite number of ways to convert the nonlinear equation into the SDC form (Cimen, 2012). However, care should be taken so as not to lose the nonlinearity or induce any singularity during the SDC parameterization. In the current study, the nonlinear dynamic model, including J2 orbital perturbation, which is expressed in ECEF (Earth-centered, Earth-fixed) frame, was parameterized into the SDC form without a loss of nonlinearity or any inducement of singularities. In previous formation flying research, SDC forms were typically derived from the nonlinear dynamics expressed in Hill frame (Irvin and Jacques, 2002; Won and Ahn, 2003; Park et al., 2011), which are suitable for control applications for a relative orbit. However, in navigation applications using GPS, because the reference frame (WGS-84) of the GPS is based

on ECEF frame, it is more convenient to use equations based on the ECEF frame rather than those based on the Hill frame. Furthermore, no previous SDC form has included J2 perturbation for relative orbital dynamics except the SDC form of Park et al. (2011), which required the use of the linearized transformation matrix (Gim and Alfriend, 2003) to consider J2 perturbation. The newly derived SDC form not only practically preserves full nonlinearity but is also less complex than the SDC form of Park et al. (2011). The relative navigation simulations are performed for formation spacecraft in low Earth orbit (LEO) with the circular formation introduced by Sabol et al. (2001). To validate the developed navigation algorithm in a realistic environment, the orbital dynamics, including full perturbations, are considered and raw GPS data is obtained from the HILS developed by Park et al. (2013), which includes a GPS signal generator and a spaceborne GPS L1 receiver. In addition, two different extended Kalman filters (EKFR1 and EKFR2) are introduced to evaluate the estimation performance of the SDREF through a comparison analysis amongst the filters. The dynamic models of all the filters used in the current study are based on relative motion including J2 perturbation. However, the SDREF and the EKFR1 use linear state propagation, whereas EKFR2 employs nonlinear state propagation. The numerous navigation simulations are performed based on various simulation variables, such as the separations between two user satellites and filter design parameters, in which the optimal process noise values for the respective filters are determined. Then, the simulation results employing the optimized design parameters are analyzed for each filter in terms of the estimation accuracy and computation time. The current study addresses multiple challenges, which are described below. First, the relative navigation algorithm for autonomous formation flying via the SDREF is developed. Second, the SDC form and the linearized equations of relative motion are formulated under J2 perturbation in the ECEF frame. In particular, the SDC form not only preserves virtually full nonlinearity but also does not have any singularity or ignored terms. Third, to obtain realistic raw GPS data, the HILS consisting of a GPS signal generator and spaceborne GPS L1 receiver is used, which provides a ground-based testing environment to reproduce possible actual spacecraft operations. Fourth, numerous relative navigation simulations are performed and analyzed with a variety of variables, such as various separations between two spacecraft, process noise of a filter, and poor GPS geometry. Lastly, the estimation performance of the SDREF is compared with those of two EKFs to gain insights into the advantages of the SDREF. The remainder of this paper is organized as follows. Chapter 2 introduces the theory of the SDREF and presents the procedure of an SDC parameterization in detail. Chapter 3 introduces the overall simulation system and presents the simulation environment and configurations.

Please cite this article in press as: Park, H.-E., Kim, Y.-R. Relative navigation for autonomous formation flying satellites using the state-dependent Riccati equation filter. Adv. Space Res. (2015), http://dx.doi.org/10.1016/j.asr.2015.10.009

H.-E. Park, Y.-R. Kim / Advances in Space Research xxx (2015) xxx–xxx

Chapter 4 explains the overall navigation algorithm, and Chapter 5 presents the simulation results. Discussion and conclusions are presented in Chapter 6. 2. State-dependent Riccati equation filter The SDRE technique was introduced by Cloutier (1997) and is used to create a linear structure for a nonlinear system for optimal and control problems. Expanding the area, it was introduced as a new technique for nonlinear estimation (Mracek et al., 1996). The new nonlinear filter is referred to as the SDREF or SDRE-based estimator. In the current study, the SDREF is used as a nonlinear estimator for relative navigation for spacecraft formation flying.

xkjk ÞPkjk UT ð^ xkjk Þ þ Qk Pkþ1jk ¼ Uð^

^kþ1jk þ K kþ1 ½ykþ1  Hð^ ^kþ1jkþ1 ¼ x xkþ1jk Þ^ xkþ1jkk  x

x_ ¼ f ðxÞ þ w

ð1Þ

y ¼ hðxÞ þ v

ð2Þ

where x and y denote the system state vector and measurement vector, respectively, w is the Gaussian zero-mean white process noise vector, and v is the Gaussian zeromean white measurement noise vector. If f ð0Þ ¼ 0, hð0Þ ¼ 0, and they are differentiable, the nonlinear system can be converted into a linear structure having an SDC form as follows (Mracek and Cloutier, 1998): x_ ¼ FðxÞx þ w

ð3Þ

y ¼ HðxÞx þ v

ð4Þ

Eqs. (3) and (4) have a linear structure; however, they are exactly or essentially equal to the nonlinear Eqs. (1) and (2). This property enables us to manipulate the nonlinear system like a linear system without any linearization effect. This system can be converted into a discrete-time system with the same method as the classical Kalman filter (KF) (Gelb, 2001). xkþ1 ¼ Uðxk Þxk þ wk

ð5Þ

yk ¼ Hðxk Þxk þ vk

ð6Þ

where the subscript k indicates the kth time step, U is the state transition matrix, and H is the observation matrix. In the SDREF, U is given by ð7Þ

Unlike with the EKF, the Jacobian of f ðxÞ evaluated at xk is not required to calculate U; moreover, Fðxk Þ has already been formulated through the SDC parameterization (Eq. (3)). The SDREF estimates states using the next two steps. In the first step, the states and covariance matrix are predicted as follows: ð8Þ

ð10Þ

xkþ1jk ÞPkþ1jk ½I  K kþ1 Hð^ xkþ1jk ÞT Pkþ1jkþ1 ¼ ½I  K kþ1 Hð^ þ K kþ1 Rkþ1 K Tkþ1

Consider a continuous-time nonlinear system,

^kþ1jk ¼ Uð^ x xkjk Þ^ xkjk

ð9Þ

where subscript kjk indicates the estimate at time k; k þ 1jk indicates the predicted value from time k to time k þ 1; Qk is the covariance matrix for the system process noise, which is defined as Qk ¼ Eðwk wTk Þ. This methodology is similar to the EKF apart from the fact that the state transition matrix is calculated with the state vector in every step. The state can also be predicted with nonlinear equation (1) using numerical integration instead of using Eq. (8) (Alfriend et al., 2010). In the second step, the states and covariance matrix are updated with the measurement as follows:

2.1. SDREF derivation

Uðxk Þ ¼ eFðxk ÞDt

3

ð11Þ

where subscript k þ 1jk þ 1 is the updated estimate with the measurement at time k þ 1; Rkþ1 is the covariance matrix for the measurement noise at time k þ 1, which is defined as Rkþ1 ¼ Eðvkþ1 vTkþ1 Þ. Eq. (11) is the Joseph form of the covariance update, which guarantees a more stable numerical calculation by maintaining a symmetric form (Bucy and Joseph, 1968). The gain matrix K is given by 1

K kþ1 ¼ Pkþ1jk H T ð^ xkþ1jk Þ½Hð^ xkþ1jk ÞPkþ1jk H T ð^ xkþ1jk Þ þ Rkþ1 

ð12Þ As explained above, the SDREF does not approximate a dynamic model and/or measurement model by linearization. Moreover, the SDREF does not require the Jacobian matrix when calculating the state transition matrix used in the prediction of the state vector or error covariance matrix (Wang, 2011). Thus, the SDREF will have a better estimation performance compared with that of the EKF. 2.2. SDC parameterization The equation of motion, including J2 perturbation, is given by €r ¼ 

l r þ aJ 2 r3

ð13Þ T

where r is the position vector (r ¼ ½ x y z  ) in the Earth Centered Inertial (ECI) frame, and x, y, and z are the state variables. l is the gravity parameter of the Earth, r ¼ jrj, and aJ 2 indicates the acceleration vector caused by the J2 perturbation, which is given by (Vallado, 2007)  3 2
2 1  5zr2 x 6
 7 7 3 l 6 5z2 6 aJ 2 ¼  5 f 6 1  r 2 y 7 7 2 r 4
 5 5z2 3  r2 z

ð14Þ

where f  J 2 R2E , J 2 is the coefficient of the J2 perturbation, and RE is the radius of the Earth.

Please cite this article in press as: Park, H.-E., Kim, Y.-R. Relative navigation for autonomous formation flying satellites using the state-dependent Riccati equation filter. Adv. Space Res. (2015), http://dx.doi.org/10.1016/j.asr.2015.10.009

4

H.-E. Park, Y.-R. Kim / Advances in Space Research xxx (2015) xxx–xxx

The relative position vector between the chief satellite and deputy satellite is defined as DrAB  rB  rA

ð15Þ

where the subscripts ‘‘A” and ‘‘B” indicate the chief satellite and deputy satellite, respectively. Differentiating Eq. (15) twice yields D€rAB ¼ €rB  €rA

ð16Þ

Substituting Eqs. (13) and (14) into Eq. (16) yields  3 2
5z2 1  r2B xB B 6  7 7 l l 3 l 6
5z2B 1  yB 7 D€rAB ¼  3 rB  5 f6 2 6 7 þ r 3 rA r B 2 rB 4
rB 5 A  5z2 3  r2B zB B 2
 3 ð17Þ 5z2A 1  r 2 xA 7 6 A
 7 3 l 6 5z2A 7 6 þ 5 f6 1  r 2 y A 7 A 2 rA 6
 7 5 4 5z2A 3  r2 zA A

From Eq. (15), it is clear that xB ¼ xA þ DxAB ;

y B ¼ y A þ Dy AB ;

zB ¼ zA þ DzAB ;

rB ¼ jrA þ DrAB j

ð18Þ

Substituting Eq. (18) into the first equation (the x-axis) of Eq. (17) yields l D€xAB ¼  ðxA þ DxAB Þ 3 jrA þ DrAB j ! 2 3 l 5ðzA þ DzAB Þ  f 1 ðxA þ DxAB Þ 2 2 jrA þ DrAB j5 jrA þ DrAB j     l 3 l 5z2   3 xA  5 f 1  2A xA ð19Þ 2 rA rA rA First,

1 jrA þDrAB jn

can be rewritten as

1 1 n2 n ¼ n ð1  aÞ jrA þ DrAB j rA

ð20Þ

where a is defined as a 

2xA DxAB þ 2y A Dy AB0 þ 2zA DzAB r2A Dx2AB þ Dy 2AB þ Dz2AB r2A

1 1 n ¼ n ð1 þ cðnÞaÞ jrA þ DrAB j rA

ð24Þ

Substituting Eq. (24) into Eq. (19) and rearranging the equation yields   l 3 1 15 1 D€xAB ¼  3 1 þ 2 f  fg DxAB 2 rA 2 r4A rA 15 l þ fxA ð2zA þ DzAB ÞDzAB 2 r7A   l 3 1 15 1  3 cð3Þ þ 2 fcð5Þ  fgcð7Þ ðxA þ DxAB Þa 2 rA 2 r4A rA ð25Þ where, g ¼ z2A þ 2zA DzAB þ Dz2AB . Using the same methodology as the parameterization of the x-axis, the equations of the y-axis and z-axis are also derived as   l 31 15 1 fg Dy AB D€y AB ¼  3 1 þ 2 f  2 rA 2 r4A rA 15 l þ fy ð2zA þ DzAB ÞDzAB 2 r7A A   l 31 15 1  3 cð3Þ þ 2 fcð5Þ  fgcð7Þ ðy A þ Dy AB Þa 2 rA 2 r4A rA 



ð26Þ

l 9 1 15 1 1þ 2 f fg DzAB 2 rA 2 r4A r3A 15 l þ fzA ð2zA þ DzAB ÞDzAB 2 r7A   l 9 1 15 1  3 cð3Þ þ 2 fcð5Þ  fgcð7Þ ðzA þ DzAB Þa 2 rA 2 r4A rA

D€zAB ¼ 

ð27Þ Because Eqs. (25)–(27) are expressed in the ECI frame, the Earth’s rotation should be considered to be transformed into the ECEF frame as follows: 2 3 D€xAB 6 7 D€rAB ¼ 4 D€y AB 5 þ C ECEF ð28Þ D€zAB where C ECEF is the correction for the rotating reference frame, which is given by (Vallado, 2007)

ð21Þ

Using the negative binominal theory, the right-hand side of Eq. (20) can be expressed in series form, n n   þ1 2 1 1 n n2 2 2 a þ ... ð1  aÞ ¼ n 1 þ a þ ð22Þ 2! rnA rA 2 Let us define n n þ1 n cðnÞ  þ 2 2 a þ ... 2 2!

Then, Eq. (20) becomes

ð23Þ

C ECEF ¼ 2xE  D_rAB  xE  ðxE  DrAB Þ

ð29Þ

where xE is the rotation rate of the Earth. As a result, the SDC form, i.e., Eq. (3), is formulated in this study as follows: x_ pv ¼ F SDC ðxpv Þxpv

ð30Þ

where the subscript ‘‘pv” indicates the states for the relative T position and velocity. xpv ¼ ½ DrTAB D_rTAB  ; DrAB ¼ ½ DxAB Dy AB DzAB T ; D_rAB ¼ ½ D_xAB D_y AB D_zAB T ; the components of matrix F SDC are listed in Appendix A. The

Please cite this article in press as: Park, H.-E., Kim, Y.-R. Relative navigation for autonomous formation flying satellites using the state-dependent Riccati equation filter. Adv. Space Res. (2015), http://dx.doi.org/10.1016/j.asr.2015.10.009

H.-E. Park, Y.-R. Kim / Advances in Space Research xxx (2015) xxx–xxx

reference frame of these expressions is the ECEF frame. The SDC form was derived from the full nonlinear equations of relative motion concerning J2 perturbation, i.e., Eq. (17), with virtually no loss of nonlinearity. We can set the SDC form to be practically the same as the full nonlinear equations with the summation of the sufficient terms in the series expansion of Eq. (23). Moreover, the SDC form does not include any singularities and is also simple and easy to use compared with using the previous SDC form (Park et al., 2011) requiring the knowledge of the orbital elements. 3. Simulation setup 3.1. Navigation simulation system All of the navigation simulations are performed on a desktop PC in offline mode, where GPS data have been prepared in advance using a hardware-in-the-loop simulator (HILS). Fig. 1 shows the top-level architecture diagram of the navigation simulation. On the left side of Fig. 1, the HILS, which was developed by Park et al. (2013), can demonstrate orbit determination and closed-loop control as well as open-loop control for spacecraft formation flying. The core components include a GPS signal generator and a GPS receiver. The former is the GSS6560 model manufactured by Spirent Inc., which can generate a GPS L1 signal for up to 12 channels. The latter is the AsteRx1 PRO manufactured by Septentrio NV, which is a spaceborne GPS L1 singlefrequency receiver with 12 tracking channels. The HILS also includes the precise orbit propagator, including full orbit perturbations, to generate true motion data of the formation spacecraft by numerical integration. The details of the adopted force models are listed in Table 1. Note that in view of the objective of this study, orbit control was not performed when preparing the GPS data using the HILS.

5

In addition, the GPS constellation of the GPS signal generator is based on the YUMA almanac, which contains 32 GPS satellites. The details of the GPS configuration are also listed in Table 1. The GPS broadcast ephemeris error, GPS clock error, and ionospheric delay are taken into account when generating the GPS RF signal (Park et al., 2013). According to the GPS configurations in Table 1, the GPS observation condition is worse than what is typical due to five unhealthy GPS satellites. This might induce unfavorable GPS satellite geometry. As mentioned above, all of the navigation simulations are performed on a high-performance desktop PC (IntelÒ i7-3770 CPU at 3.40 GHz), as shown on the right side of Fig. 1. The navigation simulation is performed in three steps (error correction – absolute navigation – relative navigation). This navigation algorithm will be presented later in detail. In the relative navigation, the SDREF is used as a nonlinear estimator. In addition, two different EKFs are also applied in the relative navigation using the same scenarios to compare with the SDREF. 3.2. Formation flying scenarios Twin satellites equipped with a spaceborne GPS L1 receiver orbit in LEO with the circular formation. The orbit of the chief satellite is circular with a semi-major axis of 7000 km and an inclination of 60°. The orbit of the deputy satellite is defined such that its relative orbit is circular. Fig. 2 illustrates the circular formation for a 50-km separation in the Hill frame. The equations of the nominal relative orbit are given by (Sabol et al., 2001) xðtÞ ¼ A0 cosðnt þ h0 Þ yðtÞ ¼ 2A0 sinðnt þ h0 Þ

ð31Þ

zðtÞ ¼ B0 cosðnt þ h0 Þ

Fig. 1. Top-level architecture diagram of the navigation simulation. Orbit control part (dot line) is not used in the current research.

Please cite this article in press as: Park, H.-E., Kim, Y.-R. Relative navigation for autonomous formation flying satellites using the state-dependent Riccati equation filter. Adv. Space Res. (2015), http://dx.doi.org/10.1016/j.asr.2015.10.009

6

H.-E. Park, Y.-R. Kim / Advances in Space Research xxx (2015) xxx–xxx

Table 1 Simulation environment configurations. Force models Earth’s gravity Air drag Solar radiation pressure 3rd body GPS configurations Number of GPS satellites Broadcast ephemeris error Clock error Ionospheric model

Joint gravity model 3 (JGM3) up to 20  20 Harris–Priester density model O Sun and Moon Total 32 (healthy: 27, unhealthy: 5) 1.713 m 5 ns Constant model (TEC: 1.0  1017 electrons/m2)

Fig. 2. Nominal relative orbit for the circular formation with 50-km separation in the Hill frame.

pffiffiffi where, B0 ¼ 3A0 , and x, y, and z are the radial, alongtrack, and cross-track components in the Hill frame, respectively. n indicates the mean motion of the chief satellite, and h0 is constant. For a 50-km separation, A0 is set at 25 km. Note that the distance between the two satellites in Fig. 2 is not constantly maintained in the current study due to various orbital perturbations. The navigation simulations are performed for six scenarios, where each scenario has its own separation. The length of the separation, which is closely related with the nonlinearity of the relative dynamics, was set to 1, 10, 20, 30, 40, and 50 km. Table 2 lists the configuration of the formation flying scenarios in detail. Note that the chief and deputy satellites have the same hardware specifications. 4. Navigation algorithms The overall strategy for the relative navigation is presented, which is based on the carrier-phase differential GPS (CDGPS) technique. The relative navigation can be largely divided into three steps (error correction, absolute navigation, and relative navigation), which is performed

on the desktop PC with a 1-s interval. Fig. 3 presents the flow chart of the overall navigation strategy. All of the symbols and equations in Fig. 3 are presented throughout this chapter. 4.1. Error corrections Error sources of GPS data should first be eliminated to achieve precise relative navigation. The expected error sources consist of GPS satellite clock error, GPS broadcast ephemeris error, ionospheric delay, GPS receiver clock error, and the GPS receiver clock jump effect. As shown in Fig. 3, the receiver clock jump effect and GPS satellite clock error are eliminated in the error correction part. The remaining errors are corrected or estimated in the relative navigation step. The AsteRx1 PRO GPS receiver used in the current paper executes a clock jump to synchronize the receiver clock with the GPS time within a predefined threshold (e.g., 0:5 ms), which induces two main effects on carrierphase measurements. One effect is the common bias for all measurements, which is naturally corrected during the estimation of the receiver clock offset or doubledifference. The other effect is the individual error for each observation because it is related with the range rate between a GPS satellite and receiver (Zhang et al., 2013). The latter effect must be corrected first because it is not eliminated through differential processing. Moreover, the magnitude of the error caused by the effect is too large (approximately 4 cycles) to precisely estimate the relative states (Zhang et al., 2013). The receiver clock jump is corrected using the receiver clock offset and Doppler measurements with the following equations (Kim and Langley, 2002): y P ;corr ¼ y P  y D  kL1  srcv

ð32Þ

y u;corr ¼ y u  y D  srcv

ð33Þ

where, y P , y u , and y D indicate the code measurement, carrier phase, and Doppler measurement, respectively. kL1 is the wavelength of the L1 frequency, and srcv is the receiver clock offset. As shown in Fig. 3, srcv is obtained from the solution of the least-squares (LSQ). Therefore, the ‘‘Receiver clock jump”-‘‘GPS satellite clock”-‘‘LSQ”-‘‘Re ceiver clock jump” procedure should be iterated at least once. The GPS satellite clock error is also corrected using the GPS clock coefficients of the broadcast navigation message. However, the corrected GPS measurements are only used for the absolute solution. In the case of the relative navigation, the error is eliminated through the doubledifference. 4.2. Absolute navigation The absolute navigation is performed for each user spacecraft by two steps, as shown in Fig. 3. First, coarse absolute solutions are obtained at every step using the

Please cite this article in press as: Park, H.-E., Kim, Y.-R. Relative navigation for autonomous formation flying satellites using the state-dependent Riccati equation filter. Adv. Space Res. (2015), http://dx.doi.org/10.1016/j.asr.2015.10.009

H.-E. Park, Y.-R. Kim / Advances in Space Research xxx (2015) xxx–xxx

7

Table 2 Formation flying scenario configurations. Epoch time (UTC)

2014-01-07 00:00:00

Simulation time/step size (second)

6000/1.0

Formation configuration Formation geometry Six separations (km)

Circular formation 1, 10, 20, 30, 40, and 50

Initial position and velocity vectors for 50-km separation (ECI frame) Chief satellite Deputy satellite x (m) 7,000,000.0 7,025,000.0 y (m) 0.0 37,500.0 z (m) 0.0 21,650.6 Vx (m/s) 0.0 0.0 Vy (m/s) 3773.0 3759.6 Vz (m/s) 6535.1 6511.7 Spacecraft specifications Mass (kg) Cross-sectional area (m2) Drag coefficient (C D ) Solar radiation pressure coefficient (C R ) Equipment

100 1 2.2 1.3 Spaceborne GPS L1 receiver

Fig. 3. Flow chart of the overall navigation algorithm.

LSQ technique, where the absolute position (rA , rB ), velocity (_rA , r_ B ), receiver clock offset (srcv;A , srcv;B ), and drift (_srcv;A , s_ rcv;B ) are coarsely estimated. The absolute states are used as a priori states of EKF and are used to calculate the elevation angle of the GPS satellites. The coarse solutions provide an absolute position accuracy of approximately 10 meters. Although relative navigation is not sensitive to the accuracy of the absolute solution (Leung and Montenbruck, 2005), it is desirable to reduce the absolute position error down to the level of a few meters to precisely

and stably estimate relative states. However, it is difficult to achieve such accuracy in this step because ionospheric delay is not corrected. In the next step, the EKF and ionosphere-free group and phase ionospheric calibration (GRAPHIC) measurements (Yunck, 1993) are employed for the final absolute solutions. In addition, the linearized dynamic model based on Eq. (13), including only J2 perturbation, is used in the EKF. Henceforth, we refer to the EKF used in the absolute navigation as EKFA. Because the GRAPHIC measure-

Please cite this article in press as: Park, H.-E., Kim, Y.-R. Relative navigation for autonomous formation flying satellites using the state-dependent Riccati equation filter. Adv. Space Res. (2015), http://dx.doi.org/10.1016/j.asr.2015.10.009

8

H.-E. Park, Y.-R. Kim / Advances in Space Research xxx (2015) xxx–xxx

ments are obtained from the arithmetic mean of the code and carrier phase, a bias caused by carrier-phase ambiguity occurs in the GRAPHIC measurements. y G ¼ ðy P þ y u Þ=2 ¼ q þ cðsrcv  sGPS Þ þ bG

ð34Þ

where y G is the ionosphere-free GRAPHIC measurement, q is the geometric range to the GPS satellite, c is the speed of light, and bG indicates the bias of the GRAPHIC measurement. Consequently, the EKFA estimates the biases (bG ) of the GRAPHIC measurements as well as other parameters, such as the absolute position (rA , rB ), absolute velocity (_rA , r_ B ), and the receiver clock offset (srcv;A , srcv;B ), in real-time. 4.3. Relative navigation The relative navigation is performed using the CDGPS technique to precisely estimate relative states, and the SDREF is employed as the nonlinear estimator. The CDGPS technique uses double-difference carrier-phase measurements, which makes it possible to resolve the double-difference ambiguity vector into an integer separated from biases because the initial phases and clock offsets of the receivers of the two satellites cancel out. This improves the precision of the relative position estimation (Kroes et al., 2005; Psiaki and Mohiuddin, 2007). The SDREF estimates the relative position (DrAB ) and velocity (D_rAB ), ionospheric delay in the zenith direction and the double-difference ambiguity vector (aiono L1 ), (rDN ijAB ). When processing GPS data based on doubledifference measurements using the filtering algorithm, the estimated ambiguities rDN ijAB may be floating point numbers; these float ambiguities should be fixed to an exact integer in the filter algorithm. In the current study, the least-squares ambiguity decorrelation adjustment (LAMBDA) method is used as an integer ambiguity resolution algorithm. The LAMBDA method, introduced by Teunissen (1995) and developed by the Mathematical Geodesy and Positioning of the Delft University of Technology, can provide accurate integer ambiguity combinations, and by updating the state vectors with fixed integer ambiguities, the relative navigation performance can be improved. In relative navigation, the two different EKFs (EKFR1 and EKFR2) are also used instead of the SDREF. The two EKFs are distinguished by their approach to time propagation. The EKFR1 employs the linear state propagation, whereas the EKFR2 employs the nonlinear state propagation. Note that the SDREF uses the linear approach to time propagation. Table 3 lists the characteristics of the filters in detail. All simulation results of the SDREF, EKFR1, and EKFR2 are compared to evaluate the performance of the filters. 4.3.1. Dynamic models As mentioned before, the SDREF and the two EKFs use dynamic models to predict the state and covariance matrix. All of the dynamic models regarding the relative

motion of the three filters used in relative navigation are based on Eq. (13) (i.e., Equations of motion, including J2 perturbation). More specifically, the SDREF uses the SDC form of Eq. (30), and the EKFR1 uses the linearized equations of Eq. (17), which can be found in Appendix B. Note that both Eqs. (17) and (30) are derived from Eq. (13). In addition, the EKFR2 employs the ‘‘pseudo” relative dynamic model that is established by independently propagating the motion of each user spacecraft (Kroes, 2006). Each dynamic model is based on Eq. (13). In other words, the ‘‘pseudo” relative dynamic model is driven by the direct difference between the accelerations described by Eq. (13) for each satellite. It is worthwhile to remind the reader that the dynamic models used in the relative navigation filters are simpler than the dynamic model used to simulate the GPS observations (e.g., J2 only vs. full perturbations, including a 20  20 gravity field model). Nevertheless, the errors caused by these differences will be compensated in the filters through the optimized design parameters, such as process noise. By including not only the relative state (xpv ) mentioned above but also the zenith ionospheric delay (amk  1iono L1 ) and the double-difference ambiguity vector (rDN ijAB ), the dynamic model (Eq. (5)) used in the relative navigation can be formulated as follows: 2 3 2 3 06m0k Udyn;66 061 xpv 6 iono 7 6 7 01m0k 4 aL1 5 ¼ 4 016 Uion;11 5 0m0kþ1 6 0m0kþ1 1 Uamb;ðm0kþ1 m0k Þ rDN ijAB kþ1 ð35Þ 2 3 xpv 6 7  4 aiono 5 þ wk L1 rDN ijAB

k

where Udyn can be calculated by substituting F SDC ðxpv Þ of Eq. (30) into Eq. (7) and defining Dt of Eq. (7) as one second. Uion is defined as I 11 in the current study, and Uamb is determined in accordance with the combination of the used GPS satellites in the double-difference. For example, if there is no change in the combination between the kth and k + 1th steps, Uamb becomes I m0k m0k . mk is the number of used GPS measurements at the kth step, and m0k is defined as. Finally, wk is given by T wk ¼ ½ rr rr rr rv rv rv riono rAMB1 rAMB2 . .. rAMBm0kþ1 

ð36Þ

where, rr , rv , riono , and rAMB are the process noise values for the relative position, velocity, zenith ionospheric delay, and the double-difference ambiguity, respectively. 4.3.2. Measurement model Double-differential code measurements and carrierphase measurements are modeled for precise relative navigation. Unknowns include the relative position and velocity (DrAB , D_rAB ), zenith ionospheric delay (aiono L1 ), and

Please cite this article in press as: Park, H.-E., Kim, Y.-R. Relative navigation for autonomous formation flying satellites using the state-dependent Riccati equation filter. Adv. Space Res. (2015), http://dx.doi.org/10.1016/j.asr.2015.10.009

H.-E. Park, Y.-R. Kim / Advances in Space Research xxx (2015) xxx–xxx

9

Table 3 Three filters for relative navigation. Filter

SDREF

EKFR1

EKFR2

Dynamic model State propagation

SDC form based on Eq. (13) (Appendix A) Linear approach xkþ1 ¼ Uðxk Þxk

Linearized relative dynamic model based on Eq. (13) (Appendix B) Linear approach xkþ1 ¼ Uk xk

Measurement model

Eq. (40)

Eq. (40)

‘‘Pseudo” relative dynamic model based on Eq. (13) Nonlinear R t approach xkþ1 ¼ tkkþ1 f ðx; tÞdt Eq. (40)

Fig. 4. Number of common GPS satellites (left) and the GDOP (right).

Fig. 5. Absolute position errors of the chief satellite in ECEF.

Please cite this article in press as: Park, H.-E., Kim, Y.-R. Relative navigation for autonomous formation flying satellites using the state-dependent Riccati equation filter. Adv. Space Res. (2015), http://dx.doi.org/10.1016/j.asr.2015.10.009

10

H.-E. Park, Y.-R. Kim / Advances in Space Research xxx (2015) xxx–xxx

Fig. 6. 3-D RMS values of relative solutions for process noise of a filter.

Please cite this article in press as: Park, H.-E., Kim, Y.-R. Relative navigation for autonomous formation flying satellites using the state-dependent Riccati equation filter. Adv. Space Res. (2015), http://dx.doi.org/10.1016/j.asr.2015.10.009

H.-E. Park, Y.-R. Kim / Advances in Space Research xxx (2015) xxx–xxx

double-differential ambiguity of the carrier phase (rDN ijAB ). Because the focus of this study is relative state estimation, the absolute position and velocity are only obtained from the estimation result of the absolute navigation. The receiver clock offset was canceled because of the usage of only double-differential measurements. Because the ionospheric delay is not completely canceled for a long baseline (>10 km) because it is related with the signal path between a GPS satellite and receiver, the ionospheric delay should be considered with each elevation angle of the GPS satellite. The measurement models can be described as follows (Psiaki and Mohiuddin, 2007): rDy ijP;AB

¼

rD^eijAB DrAB

þ

ij aiono L1 rDM AB

ð37Þ

"

rDy ijP;AB rDy iju;AB

2

#

¼4 k

rD^eijAB;ðm0 3Þ k

0m0k 3

rD^eijAB;ðm0 3Þ 0m0k 3 k 2 3 DrAB 6 D_r 7 AB 7 6  6 iono 7 þ vk 4 aL1 5 rDN ijAB

11 rDM ijAB;ðm0 1Þ

0m0k m0k

rDM ijAB;ðm0 1Þ

kL1;ðm0k m0k Þ

k

k

3 5

k

ð40Þ

where vk is defined as vk ¼ ½ rPR1

rPR2

...

rPRm0

k

ru1

ru2

. . . rum0 T k ð41Þ

where rPR and ru are the measurement noise values for the code and carrier phase, respectively. 5. Simulation results

ij ij rDy iju;AB ¼ rD^eijAB DrAB  aiono L1 rDM AB þ kL1 rDN AB

ð38Þ

for the code and carrier phase, respectively. Here, rDðÞijAB is the double-differential operator, which is defined as ij i i j j rDðÞAB  ðÞA  ðÞB  ðÞA þ ðÞB , and the superscripts i and j indicate the ith and jth GPS satellite, respectively. ^e is a line of sight vector, f L1 is the frequency of the GPS L1 signal, and M is a mapping function related with the ionospheric delay as follows (Leung and Montenbruck, 2005): iono gps I gps rcv ¼ aL1 M rcv ;

aiono L1 ¼

40:3N e ; f 2L1

2:037 M gps rcv ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi gps sin2 elgps rcv þ 0:076 þ sin elrcv

ð39Þ

where I gps rcv is the ionospheric delay, N e is the total electron content, and elgps rcv is the elevation angle of the GPS satellite. The GPS satellite ephemeris errors and clock errors as well as the GPS receiver clock offset are canceled by the doubledifferential operation. Other sources of measurement errors, such as the GPS receiver’s thermal noise, are considered as measurement noise. Note that in the SDC parameterization, the measurement model is also converted into the SDC form (Eq. (4)). However, as described above, we used the purely linearized measurement model in the current study, which is similar to the general measurement models used in the GPS navigation applications. Unlike the dynamic model, using a nonlinear measurement model hardly influences the estimation accuracy in GPS applications because the measurement system is highly linear. In addition, the other EKFs and the SDREF use the same measurement model. As a result, the measurement model, which corresponds to Eq. (6), used in the relative navigation is formulated as follows:

To evaluate the estimation performance of a filter, relative navigation simulations were numerically performed for the six formation flying scenarios. This chapter first presents the information on the GPS satellite geometry because it is closely related with the estimation performance. Then, in sequential order, the absolute and relative navigation results are presented. The left graph of Fig. 4 shows the number of common GPS satellites simultaneously observed by both user spacecraft during the simulation time, where the number is, in general, maintained for more than six satellites except approximately 2000 s and from 3800 s to 5200 s. During those times, the number of GPS satellites decreases to six or five. Because of this, the three evident peaks of the GDOP occurs approximately 2000, 4000, and 5000 s, as shown in the right graph of Fig. 4. In particular, the GDOP increases to approximately 11 around 4,000 s, which indicates poor GPS satellite geometry. 5.1. Absolute solutions The final absolute positions were estimated using the EKFA with the ionosphere-free GRAPHIC measurements. Fig. 5 shows the absolute position errors of the chief satellite during the simulation time, where the blue line is the Table 4 Design parameter configurations of all the filters. SDREF Process noise rr ðmÞ rv ðm=sÞ riono ðmÞ rAMB ðcycleÞ

EKFR1 1=2

ð1  107 Þ 1=2 ð1  109 Þ 1 1=2 ð1  10 Þ 1=2 ð1  101 Þ

Measurement noise rPR ðmÞ 0.4 ru ðmÞ 0.002

EKFR2 1=2

1=2

ð1  103 Þ 1=2 ð1  105 Þ 1 1=2 ð1  10 Þ 1=2 ð1  101 Þ

ð1  107 Þ 1=2 ð1  109 Þ 1 1=2 ð1  10 Þ 1=2 ð1  101 Þ

0.4 0.002

0.4 0.002

Please cite this article in press as: Park, H.-E., Kim, Y.-R. Relative navigation for autonomous formation flying satellites using the state-dependent Riccati equation filter. Adv. Space Res. (2015), http://dx.doi.org/10.1016/j.asr.2015.10.009

12

H.-E. Park, Y.-R. Kim / Advances in Space Research xxx (2015) xxx–xxx

Fig. 7. Relative navigation accuracy for the 50-km separation.

Please cite this article in press as: Park, H.-E., Kim, Y.-R. Relative navigation for autonomous formation flying satellites using the state-dependent Riccati equation filter. Adv. Space Res. (2015), http://dx.doi.org/10.1016/j.asr.2015.10.009

H.-E. Park, Y.-R. Kim / Advances in Space Research xxx (2015) xxx–xxx

13

Table 5 Relative navigation accuracy comparisons of the three filters for six different separations.

Position (mm)

3-D RMS

Max 3-D error

Velocity (mm/s)

3-D RMS

Max 3-D error

Separation (km)

SDREF

EKFR1

EKFR2

1 10 20 30 40 50 1 10 20 30 40 50

1.4 2.0 2.8 3.9 5.7 6.0 5.2 7.6 11.4 15.7 19.3 19.9

4.1 4.6 5.1 6.1 8.1 9.1 27.8 26.1 26.2 34.9 48.6 80.3

1.4 2.0 2.8 3.9 5.7 6.0 5.2 7.6 11.5 15.7 19.4 19.9

1 10 20 30 40 50 1 10 20 30 40 50

0.089 0.095 0.102 0.116 0.139 0.153 0.248 0.254 0.271 0.332 0.425 0.492

0.40 0.50 1.06 2.17 3.76 5.81 2.40 2.05 3.09 4.28 6.80 9.36

0.089 0.095 0.102 0.116 0.140 0.154 0.248 0.256 0.276 0.338 0.428 0.499

accuracy is not highly precise but is acceptable when considering the size of the simulated broadcast ephemeris errors (refer to Table 1) and when compared with other studies using the GRAPHIC measurement (e.g., Leung and Montenbruck, 2005). The final absolute solution will be utilized as the reference point in the subsequent relative navigation, and its accuracy will be helpful to securely and precisely estimate relative states. 5.2. Filter performance and process noise

Fig. 8. Computational load for the three filters.

coarse solutions of the LSQ, and the red line indicates the final solution of the EKFA. The errors were calculated from the difference between the estimated absolute positions and the true motion data described in Section 3.1. The errors of the coarse solution using LSQ increased occasionally by more than 20 m due primarily to the ionospheric delay. Meanwhile, the final solution using the EKFA was estimated with a 3-D RMS error of 1.39 m. Note that the value was calculated using data after the initial convergence (500 s). The accuracy of the final absolute solution was significantly improved due to the ionospheric correction by the GRAPHIC measurements. This absolute

The process noise in filter design parameters significantly influences the estimation performance, particularly in a dynamic filter. To find the optimized process noise of each filter, the 3-D RMS errors of the relative position and velocity were investigated through numerous relative navigation simulations using various configurations. Fig. 6 shows the 3-D RMS errors as a function of the separation and the process noise value (r2 ) for each filter. Note that 3-D RMS errors of more than 10 mm are not expressed in the graphs for the relative position, and those above 8 mm/s are not expressed in graph (d) of Fig. 6. As a result, 9.1% of all of the simulations using the SDREF or EKFR2 were eliminated by a threshold of 10 mm, whereas 40.9% of all of the simulations using EKFR1 were eliminated by the same threshold. For relative velocity, no result exceeded the threshold of 0.8 mm/s in the simulations using the SDREF or EKFR2; however, 31.8% of all of the simulations using the EKFR1 were eliminated even though a threshold of 8 mm/s, which is 10 times larger than one of the SDREF or EKFR2, were considered.

Please cite this article in press as: Park, H.-E., Kim, Y.-R. Relative navigation for autonomous formation flying satellites using the state-dependent Riccati equation filter. Adv. Space Res. (2015), http://dx.doi.org/10.1016/j.asr.2015.10.009

14

H.-E. Park, Y.-R. Kim / Advances in Space Research xxx (2015) xxx–xxx

As seen in (a) and (e) of Fig. 6, it is difficult to distinguish the difference between the relative position results of SDREF and EKFR2. Both demonstrated the best esti1=2

mation performance at a process noise of ð1  107 Þ m. As seen in (c) of Fig. 6, the EKFR1 demonstrated the best performance with small process noise for only short separations. If we must determine the optimized process noise available regardless of separations, process noises of 1=2

1=2

ð1  103 Þ m or ð1  102 Þ m are suitable. As seen in (b) and (f) of Fig. 6, it is still difficult to recognize the difference between the two graphs for the relative velocity results of SDREF and EKFR2. The best estimation performance occurred with a process noise of 1=2

ð1  109 Þ m/s. In graph (d) of Fig. 6, the 3-D RMS errors are much larger than the others and significantly increase as the separation increases. Consequently, a

EKFR1 was primarily caused by the linearization. The reason for the larger EKFR1 noise is because its set process noise value is larger than that of SDREF or EKFR2 (see Table 4). However, as mentioned above, this is already the optimized value for the EKFR1. Incidentally, despite the inaccurate velocity, the position was estimated relatively accurately, which is related with the process noise in the filter. In other words, the process noise of the EKFR1 was set to larger values in comparison with those of the SDREF or EKFR2, which thus minimized the impact of the inaccurate velocity when estimating the relative position. Table 5 summarizes the 3-D RMS values of the position errors and velocity errors for all of the relative navigation simulations

1=2

process noise of ð1  105 Þ m/s is always available. As a result, the optimized process noise values are listed in Table 4, where the other parameters (riono , rAMB , rPR , ru ) are set to the same value for each parameter because all of the filters use the same ionospheric model, ambiguity resolution, and GPS measurements. As seen in Table 4, the EKFR1’s process noises (rr and rv ) concerning the relative position and velocity are much larger than those of the SDREF or EKFR2. It might appear that the EKFR1 receives unfair treatment because a smaller process noise may induce better precision in general. However, an excessively small process noise that is inadequate for filter performance prevents accurate and precise estimation, as shown in Fig. 6. In other words, larger approximations are made in the EKFR1, which can be compensated for through a larger process noise. 5.3. Relative solutions The relative position and velocity were estimated with the filter design parameters of Table 4 for the respective filters. Fig. 7 shows the relative navigation results of the three filters for the 50-km separation, which includes six graphs that show the position and velocity errors for the radial, along-track, and cross-track axes. As seen in (a), (c), and (e) of Fig. 7, the results of the SDREF and EKFR2 are extremely similar to each other, whereas the position error of the EKFR1 is not only noisier but also occasionally unstable with the evident growth of the error. The unstable time of the EKFR1 is coincident with the time of the high GDOP in Fig. 4, which means that EKFR1 is more sensitive to the effect of a poor GPS satellite geometry than the SDREF or EKFR2. In the case of velocity error, as seen in (b), (d), and (f) of Fig. 7, there is nearly no difference between the errors of the SDREF and EKFR2, whereas the EKFR1 error is not only noisier but also larger, with any pattern caused by the discrepancy between the dynamic model of the filter and the true dynamics. The dynamic discrepancy of the

Fig. 9. Analysis of the relative position accuracy for the separations.

Please cite this article in press as: Park, H.-E., Kim, Y.-R. Relative navigation for autonomous formation flying satellites using the state-dependent Riccati equation filter. Adv. Space Res. (2015), http://dx.doi.org/10.1016/j.asr.2015.10.009

H.-E. Park, Y.-R. Kim / Advances in Space Research xxx (2015) xxx–xxx

As mentioned before, the computation speed of each filter was also analyzed. Fig. 8 shows the computation times for the respective filters. The times were measured over the execution of the entire navigation algorithm (see Fig. 3) on the desktop PC. The SDREF required the shortest computation times for all of the simulations. However, the time differences between the SDREF and the EKFR1 were insignificant. On the contrary, the EKFR2 required much longer computation times by a factor of approximately four compared with that of the SDREF. Specifically, the SDREF and EKFR1 required similar times of 34 and 36 s, respectively, whereas the EKFR2 required a much longer time of approximately 121 s. The long computation time of the EKFR2 was mostly caused by the two orbit propagations. The EKFR2 independently propagates the

15

motion of each user spacecraft through numerical integration. Fig. 9 shows the 3-D RMS and maximum errors of the relative position as a function of the separation for the three filters. Regardless of the separation, the estimation performance of the SDREF and the EKFR2 was essentially the same and better than that of the EKFR1. The overall accuracy was reduced when the separation increased because the effect of common error cancelation decreased. For a 50-km separation, the relative positions of the SDREF and the EKFR2 were estimated with a 3-D RMS error of 6.0 mm, whereas those of the EKFR1 were estimated with a 3-D RMS error of 9.1 mm. In addition, the maximum errors of the SDREF and EKFR2 were 19.9 mm and that of the EKFR1 was 80.3 mm. Fig. 10 shows the 3-D RMS and maximum errors of the relative velocity as a function of the separation for the three filters. For the velocity errors, the difference between the SDREF and EKFR1 appeared more evidently. For the 50-km separation, the relative velocities of the SDREF, EKFR1, and EKFR2 were estimated with a 3-D RMS error of 0.153 mm/s, 5.81 mm/s, and 0.154 mm/s, respectively. For the maximum errors, the SDREF, EKFR1, and EKFR2 were 0.492 mm/s, 9.36 mm/s, and 0.499 mm/s, respectively. As seen in Figs. 9 and 10, the SDREF and EKFR2 achieved the best estimation accuracy out of all the simulations. The SDREF was more accurate than the EKFR1. In particular, for the velocity, the accuracy difference between the two filters increased to approximately 38 times for the 50-km separation. 6. Conclusions

Fig. 10. Analysis of the relative velocity accuracy for the separations.

The SDREF was developed for relative navigation, in which the SDC form was derived from the nonlinear relative dynamics concerning J2 perturbation, which is expressed in ECEF frame, without any loss of nonlinearity. The SDC form not only does not include any singularities but is also simpler and easier to use in GPS relative navigation applications than the previous SDC forms based on Hill frame. The navigation simulations were performed for formation spacecraft in LEO with the circular formation, where GPS data were received using the HILS every second. To evaluate the developed SDREF performance, the same simulations were performed using the two different EKFs and the SDREF; then, the results were compared with each other in terms of the estimation accuracy and the computation load. In the simulations, the SDREF provided precise estimation accuracies of 6.0 mm and 0.153 mm/s for position and velocity, respectively, for the 50-km separation, with a relatively short computation time of approximately 34 s. Compared with the EKFR1, the accuracies improved by 34% and 3797% for position and velocity, respectively. Compared with the EKFR2, the computation time reduced by 72%. In conclusion, the SDREF estimated the relative states as accurately as the EKFR2 while requiring considerably

Please cite this article in press as: Park, H.-E., Kim, Y.-R. Relative navigation for autonomous formation flying satellites using the state-dependent Riccati equation filter. Adv. Space Res. (2015), http://dx.doi.org/10.1016/j.asr.2015.10.009

16

H.-E. Park, Y.-R. Kim / Advances in Space Research xxx (2015) xxx–xxx

less computational effort and as quickly as the EKFR1 but with better estimation accuracy, particularly in the velocity estimation. This means that the developed SDREF combined the strong points of EKFR1 and EKFR2 and overcame their disadvantages. This property of the SDREF will be highly helpful in practical applications, particularly in the case in which there is an onboard system with limited resources available and real-time implementation is required, such as autonomous formation flying using small spacecraft. Acknowledgment This Research has been performed as a collaborative research project of Building Response System for National-wide Issues Based on High performance Supercomputer supported by the Korea Institute of Science and Technology Information (KISTI). Appendix A. SDC form, including J2 perturbation x_ pv ¼ F SDC ðxpv Þxpv   F 11 F 12 F SDC ðxÞ ¼ F 21 F 22 2 3 2 3 0 0 0 1 0 0 6 7 6 7 F 11 ¼ 4 0 0 0 5; F 12 ¼ 4 0 1 0 5; 0 0 0 0 0 1 2 3 F 21 ð1; 1Þ F 21 ð1; 2Þ F 21 ð1; 3Þ 6 7 F 21 ¼ 4 F 21 ð2; 1Þ F 21 ð2; 2Þ F 21 ð2; 3Þ 5; F 21 ð3; 1Þ F 21 ð3; 2Þ F 21 ð3; 3Þ 2 3 0 2XE 0 6 7 0 05 F 22 ¼ 4 2XE 0 0 0   l 31 15 1 F 21 ð1;1Þ ¼  3 1 þ 2 f  fg 2 rA 2 r4A rA   l 31 15 1 þ 5 cð3Þ þ 2 fcð5Þ  fgcð7ÞÞðxA þ DxAB 2 rA 2 r4A rA  ð2xA þ DxAB Þ   l 3 1 15 1 F 21 ð1; 2Þ ¼ 5 cð3Þ þ 2 fcð5Þ  fgcð7Þ 2 rA 2 r4A rA

  l 3 1 15 1 1 þ f  fg 2 r2A 2 r4A r3A   l 3 1 15 1 þ 5 cð3Þ þ 2 fcð5Þ  fgcð7Þ rA 2 rA 2 r4A

F 21 ð2; 2Þ ¼ 

 ðy A þ Dy AB Þð2y A þ Dy AB Þ F 21 ð2; 3Þ ¼

 ðy A þ Dy AB Þð2zA þ DzAB Þ   l 9 1 15 1 F 21 ð3; 1Þ ¼ 5 cð3Þ þ 2 fcð5Þ  fgcð7Þ 2 rA 2 r4A rA  ðzA þ DzAB Þð2xA þ DxAB Þ   l 9 1 15 1 F 21 ð3; 2Þ ¼ 5 cð3Þ þ 2 fcð5Þ  fgcð7Þ rA 2 rA 2 r4A  ðzA þ DzAB Þð2y A þ Dy AB Þ   l 9 1 15 1 F 21 ð3; 3Þ ¼  3 1 þ 2 f  fg 2 rA 2 r4A rA 15 l þ fzA ð2zA þ DzAB Þ 2 r7A   l 3 1 15 1 fgcð7Þ þ 5 cð3Þ þ 2 fcð5Þ  2 rA 2 r4A rA  ðzA þ DzAB Þð2zA þ DzAB Þ Appendix B. Linearized equation of relative motion, including J2 perturbation x_ pv ¼ F EKF R1 xpv  F 11 F EKF R1 ¼ F 21 2 0 0 6 F 11 ¼ 4 0 0 2

15 l fxA ð2zA þ DzAB Þ 2 r7A   l 3 1 15 1 fgcð7Þ þ 5 cð3Þ þ 2 fcð5Þ  rA 2 rA 2 r4A

 ðxA þ DxAB Þð2zA þ DzAB Þ   l 3 1 15 1 F 21 ð2; 1Þ ¼ 5 cð3Þ þ 2 fcð5Þ  fgcð7Þ rA 2 rA 2 r4A  ðy A þ Dy AB Þð2xA þ DxAB Þ

0

0

F 12

2

F 22

0

0

1

0 0

F 21 ð1; 1Þ

F 21 ð1; 2Þ

0

nz ¼ 3

F 21 ð2; 2Þ

F 21 ð3; 1Þ

0 6 ¼ 4 2XE

nxy ¼ 3



F 22 3 2 0 1 0 6 7 0 5; F 12 ¼ 4 0 1

6 F 21 ¼ 4 F 21 ð2; 1Þ

 ðxA þ DxAB Þð2y A þ Dy AB Þ F 21 ð1; 3Þ ¼

15 l fy ð2zA þ DzAB Þ 2 r7A A   l 3 1 15 1 fgcð7Þ þ 5 cð3Þ þ 2 fcð5Þ  rA 2 rA 2 r4A

F 21 ð3; 2Þ 3 2XE 0 7 0 05 0

3

7 0 5;

F 21 ð1; 3Þ

3

7 F 21 ð2; 3Þ 5; F 21 ð3; 3Þ

0

l 15 fl 105 fl 2 þ  z r5A 2 r7A 2 r9A A

l 45 fl 105 fl 2 þ  z r5A 2 r7A 2 r9A A

F 21 ð1; 1Þ ¼ 

l 3 fl 15 fl 2  þ z þ nxy x2A 2 r7A A r3A 2 r5A

Please cite this article in press as: Park, H.-E., Kim, Y.-R. Relative navigation for autonomous formation flying satellites using the state-dependent Riccati equation filter. Adv. Space Res. (2015), http://dx.doi.org/10.1016/j.asr.2015.10.009

H.-E. Park, Y.-R. Kim / Advances in Space Research xxx (2015) xxx–xxx

F 21 ð1; 2Þ ¼ nxy xA y A F 21 ð1; 3Þ ¼

30 1 flxA zA 7 þ nxy xA zA 2 rA

F 21 ð2; 1Þ ¼ nxy y A xA F 21 ð2; 2Þ ¼  F 21 ð2; 3Þ ¼

l 3 fl 15 fl 2  þ z þ nxy y 2A 2 r7A A r3A 2 r5A

30 1 fly A zA 7 þ nxy y A zA 2 rA

F 21 ð3; 1Þ ¼ nz zA xA F 21 ð3; 2Þ ¼ nz zA y A F 21 ð3; 3Þ ¼ 

l 9 fl 15 fl 2 30 1  þ z þ flz2A 7 þ nxy z2A 2 r7A A 2 rA r3A 2 r5A

References Alfriend, K.T., Vadali, S.R., Gurfil, P., How, J.P., Breger, L.S., 2010. Spacecraft Formation Flying: Dynamics, Control and Navigation. 1st ed., Elsevier Astrodynamics Series. Ardaens, J.-S., D’Amico, S., 2009. Spaceborne autonomous relative control system for dual satellite formations. J. Guidance Control Dyn. 32 (6), 1859–1870. Bucy, R., Joseph, P., 1968. Filtering for Stochastic Processes with Applications to Guidance. Interscience Publishers, New York. Chang, I., Park, S.-Y., Choi, K.-H., 2009. Decentralized coordinated attitude control for satellite formation flying via the state-dependent Riccati equation technique. Int. J. Non-Linear Mech. 44 (8), 891–904. Cimen, T., 2012. Survey of state-dependent Riccati equation in nonlinear optimal feedback control synthesis. J. Guidance Control Dyn. 35 (4), 1025–1047. Cloutier, J.R., 1997. State-dependent Riccati equation techniques: an overview. In: The American Control Conference. IEEE Publ., pp. 932– 936. Gelb, A., 2001. Applied Optimal Estimation, 16th ed. The M.I.T. Press. Gill, E., D’Amico, S., Montenbruck, O., 2007. Autonomous formation flying for the PRISMA mission. J. Spacecraft Rockets 44 (3), 671–681. Gim, D.W., Alfriend, K.T., 2003. State transition matrix of relative motion for the perturbed noncircular reference orbit. J. Guidance Control Dyn. 26 (6), 956–971. Irvin, D.J., Jacques, D.R., 2002. A study of linear versus nonlinear control techniques for the reconfiguration of satellite formations. Adv. Astronaut. Sci. 109, 589–608. Kim, D., Langley, R.B., 2002. Instantaneous real-time cycle-slip correction for quality control of GPS carrier-phase measurements. J. Inst. Navig. 49 (4), 205–222. Kroes, R., Montenbruck, O., Bertiger, W., Visser, P., 2005. Precise GRACE baseline determination using GPS. GPS Solutions 9, 21–31.

17

Kroes, R., 2006. Precise relative positioning of formation flying spacecraft using GPS (Ph.D. Dissertation), Delft University of Technology, Delft, Netherlands. Lee, D., Cochran Jr., J., No, T., 2012. Robust position and attitude control for spacecraft formation flying. J. Aerosp. Eng. 25 (3), 436– 447. Leung, S., Montenbruck, O., 2005. Real-time navigation of formationflying spacecraft using global-positioning-system measurements. J. Guidance Control Dyn. 28 (2), 226–235. Massari, M., Zamaro, M., 2014. Application of SDRE technique to orbital and attitude control of spacecraft formation flying. Acta Astronaut. 94 (1), 409–420. Montenbruck, O., Ebinuma, T., Lightsey, E.G., Leung, S., 2002. A realtime kinematic GPS sensor for spacecraft relative navigation. Aerosp. Sci. Technol. 6, 435–449. Mracek, C.P., Cloutier, J.R., D’Souza, C.A., 1996. A new technique for nonlinear estimation. In: The 1996 IEEE International Conference on Control Applications. IEEE Publ., pp. 338–343. Mracek, C.P., Cloutier, J.R., 1998. Control designs for the nonlinear benchmark problem via the state-dependent Riccati equation method. Int. J. Robust Nonlinear Control 8, 401–433. Park, H.-E., Park, S.-Y., Choi, K.-H., 2011. Satellite formation reconfiguration and station-keeping using state-dependent Riccati equation technique. Aerosp. Sci. Technol. 15 (6), 440–452. Park, H.-E., Park, S.-Y., Kim, S.-W., Park, C., 2013. Integrated orbit and attitude hardware-in-the-loop simulations for autonomous satellite formation flying. Adv. Space Res. 52 (12), 2052–2066. Psiaki, M.L., Mohiuddin, S., 2007. Modeling, analysis, and simulation of GPS carrier phase for spacecraft relative navigation. J. Guidance Control Dyn. 30 (6), 1628–1639. Sabol, S., Burns, R., McLaughlin, C.A., 2001. Satellite formation flying design and evolution. J. Spacecraft Rockets 38 (2), 270–278. Teunissen, P.J.G., 1995. The least-squares ambiguity decorrelation adjustment: a method for fast GPS integer ambiguity estimation. J. Geodesy 70, 65–82. Tillerson, M., Inalhan, G., How, J.P., 2002. Co-ordination and control of distributed spacecraft systems using convex optimization techniques. Int. J. Robust Nonlinear Control 12, 207–242. Vallado, D.A., 2007. Fundamentals of Astrodynamics and Applications, third ed., Space Technology Library. Wang, X., 2011. Nonlinear control and estimation with general performance criteria (Ph.D. Dissertation), Marquette University, Wisconsin, US. Won, C.H., Ahn, H.S., 2003. Nonlinear orbital dynamic equations and state-dependent Riccati equation control of formation flying satellites. J. Astronaut. Sci. 51 (4), 433–449. Yunck, T.P., 1993. Coping with the atmosphere and ionosphere in precise satellite and ground positioning. In: Valance-Jones, A. (Ed.), Environmental Effects on Spacecraft Trajectories and Positioning, AGU Monograph. Zhang, X., Guo, B., Guo, F., Du, C., 2013. Influence of clock jump on the velocity and acceleration estimation with a single GPS receiver based on carrier-phase-derived Doppler. GPS Solutions 17, 549–559.

Please cite this article in press as: Park, H.-E., Kim, Y.-R. Relative navigation for autonomous formation flying satellites using the state-dependent Riccati equation filter. Adv. Space Res. (2015), http://dx.doi.org/10.1016/j.asr.2015.10.009