Jet transport-based nonlinear state and parameter estimation for geostationary spacecraft

Jet transport-based nonlinear state and parameter estimation for geostationary spacecraft Jianlin Chena,b,c,∗, Josep J. Masdemontc, Gerard G´omezd, Jianping Yuana,b,∗ aNational Key Laboratory of Aerospace Flight Dynamics, Xi’an, Shaanxi, 710072, China bSchool of Astronautics, Northwestern Polytechnical University, Xi’an, Shaanxi, 710072, China cIEEC &Departament de Matem`atiques, Universitat Polit`ecnica de Catalunya, Diagonal 647, 08028 Barcelona dIEEC &Departament de Matem`atiques i Inform`atica, Universitat de Barcelona, Gran Via de les Corts Catalanes 585, 08007 Barcelona Abstract Based on the jet transport technique, this paper proposes a novel nonlinear Kalman filter for simultaneously estimating the spacecraft state vector and uncertain parameters, either physically related with the spacecraft or with the measurement procedure. Two different coordinate representations, including Cartesian and hybrid geostationary orbital elements, are exploited in the new nonlinear estimators. The performance and sensitivity analyses of the proposed jet transport-based nonlinear estimators are assessed by numerical simulations and compared with the classical extended Kalman filter. Keywords: jet transport, parameter estimation, nonlinear estimator, geostationary spacecraft 1. Introduction Geostationary or geosynchronous Earth Orbit (GEO) regime has the valuable characteristic that keeps the spacecraft at rest relative to the Earthcentered Earth-fixed (ECEF) frame. This particular feature reduces the difficulty for providing long-time continuous direct broadcast, communication, ∗Corresponding author Email addresses: [email protected] (Jianlin Chen), [email protected] (Jianping Yuan) and observation services for a certain fixed region of the Earth’s surface. A precise initial GEO orbit determination (IOD) is of extreme importance for these space-related missions and, in particular, emerging high-accuracy GEO applications are being proposed to monitor the land stability subjected to natural hazards, like volcanic activities and earthquakes (future GEOSAR missions [1]), as well as to implement meso-scale measurements of sea altimetry using TV signals from GEO satellites which, being much stronger than the GPS ones, provide several advantages [2]. Currently, various filtering algorithms have been proposed and extensively exploited in the IOD problems [3–7]. A complete discussion of the pros and cons of these filtering algorithms can be found in [3, 8]. The well-known extended Kalman filter (EKF) behaves well in the linear or weak nonlinear estimation problem [9]. However, it loses accuracy, or even diverges, in some certain cases such as: when the problem of state estimation is highly nonlinear, in systems with low frequency measurements, or when the initial state estimation errors are large. The unscented Kalman filter (UKF) overcomes this drawback by removing the linearization assumption but with a higher computational cost [6]. A better nonlinear extended Kalman filtering algorithm was proposed by Park and Scheeres [7] by means of a semi-analytic uncertainty prediction technique to describe the localized nonlinear motion around the nominal trajectory, and carrying out a high order Taylor expansion of the nonlinear measurements. The bottleneck of these high order Kalman filters (HEKF) is the computational complexity of the required high order derivatives that limits its applications. Making use of the fact that Differential Algebra (DA) techniques determine high order derivatives in an accurate and efficient way, Valli et al. [3] developed a high order EKF algorithm with the use of DA for orbit determination in two simple cases: Sun-Earth halo orbit determination, and Keplerian orbit determination around the Earth. Further usages of this high order EKF in the angles-only initial orbit determination problem [4] and spacecraft relative state estimation problem [5] also validated its effectiveness. However, most of the references only consider the state estimation and leave aside the estimation of parameter uncertainties in inaccurate dynamical models or in measurement models, such as the spacecraft area-to-mass ratio, or the uncertain location of a removable ground tracking station. A simultaneous estimation of these additional parameters is meaningful to further improve the estimation accuracy of trajectories; for instance reference [10] 2 underlined that the solar radiation pressure (SRP) modeling error is the largest dynamical error source in the process of the geosynchronous navigation satellites orbit determination problem, a precise area-to-mass ratio estimation can provide good a priori information for the subsequent accurate trajectory determination. The purpose of this paper is to implement a precise geostationary or geosynchronous orbit determination, taking into account the dominant perturbations in the GEO regime. At the same time, the procedure enables to estimate additional uncertain parameters in the equations of motion or in the measurement equations, which allows to further improve the accuracy of the models and of the orbit determination solution. The methodology is systematically implemented and can be easily extended for any parameter estimation or orbit determination problem in another regime. The paper is organized as follows: Section 2 describes the equations of motion used in this work, including a Cartesian representation and a hybrid GEO element representation. In Section 3 we develop the jet transport-based HEKF algorithms for a spacecraft physical parameters estimation, and for a ground tracking station position estimation, using both kinds of coordinates. Section 4 contains some of the simulations done to validate the efficiency of the estimators. The conclusions are given in Section 5. 2. Dynamical model Several coordinate representations were proposed in the past to describe the motion of GEO satellites [11]. However, the particularity of GEO orbits leads specially to two main representations for the description of its motion without singularities. They are the usual Cartesian representation and the particular GEO element representation [12], that will be both considered in this work. 2.1. Cartesian representation Consider the motion of a spacecraft close to the GEO regime subjected to the central gravitational acceleration of an homogeneous spherical Earth acen =−µ r3x, −µ r3y, −µ r3zTand to the four main perturbing accelerations: anon,aSRP ,as, and am, that respectively denote perturbing Earth’s nonspherical gravitational acceleration, solar radiation pressure, Sun and Moon 3 gravitational accelerations [9]. The equations of motion in Cartesian coordinates are ¨ r=acen +ap, ap=anon +aSRP +as+am,(1) where rindicates the spacecraft position vector in the Earth-centered inertial (ECI) frame, and ris its norm, the notation apindicates the total perturbing acceleration. The adopted formulae for the dominant perturbing accelerations can be found in detail in [9]. The SRP, Sun and Moon gravitational accelerations can be expressed in the ECI reference frame as, aECI SRP =−νP Cr A m r r3  AU2,(2) aECI M=GM rM−r |rM−r|3−rM |rM|3,(3) where ris the relative position vector pointing from the Sun to the spacecraft, Pindicates the SRP force at a distance of 1 AU, νdenotes the Earth’s shadow function and Mrepresents the mass of the Sun or Moon, rMis the position vector of the mass Mrelative to the Earth. Additional parameters, including the radiation pressure coefficient Cr, the spacecraft mass m, and the cross sectional area A, are related to the individual spacecraft properties. Note that all these spacecraft physical parameters can be affected by some uncertainties, and can be estimated by the nonlinear filters introduced in this work. The adopted Earth’s non-spherical gravity is expressed in the ECEF frame as, aECEF non =X n,m ax,nm +µ r3x, X n,m ay,nm +µ r3y, , X n,m az,nm +µ r3z, , 4 where ax,nm =µ 2R2 ⊕{(−CnmVn+1,m+1 −SnmWn+1,m+1)+ (n−m+ 2)! (n−m)! (CnmVn+1,m−1+SnmWn+1,m−1)}, ay,nm =µ 2R2 ⊕{(−CnmWn+1,m+1 +SnmVn+1,m+1)+ (n−m+ 2)! (n−m)! (−CnmWn+1,m−1+SnmVn+1,m−1)}, az,nm =µ R2 ⊕{(n−m+ 1)(−CnmVn+1,m −SnmWn+1,m)}. In these equations the Earth’s equatorial radius R⊕, and the geopotential coefficients Cnm,Snm are provided by the 5 ×5 EGM96S gravity model. In particular, m=n= 0 corresponds to the situation in which only the central gravitational acceleration is considered. The values of Vnm and Wnm are recursively given by V00 =R⊕ r, Vm−1,m = 0 , Vmm = (2m−1) R⊕ r2{xVm−1,m−1−yWm−1,m−1}, Vnm =2n−1 n−mzR⊕ r2Vn−1,m −n+m−1 n−mR2 ⊕ r2Vn−2,m, W00 = 0 , Wm−1,m = 0 , Wmm = (2m−1) R⊕ r2{xWm−1,m−1+yVm−1,m−1}, Wnm =2n−1 n−mzR⊕ r2Wn−1,m −n+m−1 n−mR2 ⊕ r2Wn−2,m. Besides, the time-dependent coordinate transformation from the ECEF frame to the ECI frame can be obtained via the use of average angular speed of the Earth’s rotation ωe= 7.292115 ×10−5rad/s, thus TECI ECEF =  cos (ωet)−sin (ωet) 0 sin (ωet) cos (ωet) 0 0 0 1  . 5 Therefore, the perturbing Earth’s non-spherical gravitational acceleration is expressed in the ECI reference frame as aECI non =TECI ECEF ·aECEF non .(4) 2.2. GEO representation In [13], Tombasco introduced a non-dimensional GEO element set in terms of the classical Keplerian elements {a, e, i, ω, Ω, θ}, which consists of the Earth-fixed sub-spacecraft longitude λ, the longitudinal drift rate δ¯a, two eccentricity vector components (ex, ey), and two equinoctial elements (Q1, Q2). The expressions of the GEO elements, as well as four intermediate variables, can be written as follows, λ,(ω+Ω+θ)−GA (t), δ¯a,a−An An , ex,ecos (ω+ Ω) , ey,esin (ω+ Ω) , Q1,tan i 2sin (Ω) , Q2,tan i 2cos (Ω) , r=An(δ¯a+ 1) 1−e2 x−e2 y 1 + excos s+eysin s, s=λ+GA (t) = ω+Ω+θ, p=An(δ¯a+ 1) 1−e2 x−e2 y, h=√pµ, where An= 42164.2 km indicates a nominal GEO semi-major axis, GA (t) = GA(t0)+ωe(t−t0) stands for the Greenwich sidereal angle at t. The notation sdenotes the spacecraft sidereal angle, prepresents the semi-latus rectum of the orbit, and his the norm of the angular momentum. With the above 6 variables, the equations of motion using GEO element representation are ˙ λ=h r2+r hQ2sin s−Q1cos sah−ωe, δ˙ ¯a=2(δ¯a+ 1)2 hAn(exsin s−eycos s)ar+p raθ, ˙ex=r hnp rsin s·ar+ex+ (1 + p r) cos saθ+eyQ1cos s −Q2sin saho, ˙ey=r hn−p rcos s·ar+ey+ (1 + p r) sin saθ−exQ1cos s −Q2sin saho, ˙ Q1=r 2h(1 + Q2 1+Q2 2) sin s·ah, ˙ Q2=r 2h(1 + Q2 1+Q2 2) cos s·ah. (5) We note that, using Cartesian coordinates, the perturbing accelerations in (1) can be easily computed through (2), (3), and (4). However, when we implement the nonlinear estimation using GEO elements, a series of conversions have to be done for computing ap= (ar, aθ, ah)Tin (5). It is as follows: i) at a certain epoch transform the known GEO elements into Cartesian coordinates (the transformation is given in [14]); ii) compute the dominant perturbing accelerations using the Cartesian coordinates through (2), (3), and (4); iii) project all the perturbations into the Local Vertical Local Horizontal reference frame (LVLH). For the transformation of accelerations from the ECI into to the LVLH frame we make use of the cosines matrix TLV LH ECI =  I·i I ·j I ·k J·i J ·j J ·k K·i K ·j K ·k ,(6) where I,J,Kand i,j,kare the unitary vector sets defining the ECI and LVLH reference frames respectively. Using the spacecraft state (r,v) in the ECI frame, we have, i=r krk,k=r×v kr×vk,j=k×i. Thus, the four dominant perturbations we consider can be rotated into the LVLH frame by means of aLV LH p=TLV LH ECI ·aECI p.(7) 7 Finally, a jet transport-based numerical integrator is developed using κth order Taylor expansions (JT-HEKF-κalgorithm) in subsection 3.2. It is used to propagate both dynamic models and produce the spacecraft state and parameter predictions at a new epoch; then, the JT-HEKF-κfiltering algorithms do the state and parameter updates by incorporating the new measurements. The proposed JT-HEKF-κfiltering algorithms make an intensive use of this procedure to estimate the state vector and parameters in the GEO element space. 3. Jet transport technique The DA technique, as a way for doing efficient symbolic computations, was first proposed by Berz and Makino [15] in 1999 for the study of particle beam accelerators, and implemented by the same authors in the COSY Infinity package. With similar ideas, P´erez-Palau [16] studied the application of JT technique in the structure detection of astrodynamics problem. Essentially, either DA method or JT method is an automatic differentiation technique that enables to derive high order Taylor expansions (or other kind of expansions) of general nonlinear functions and, in particular, of the flow associated to an ODE in an accurate and efficient way. For convenience of the reader, some vital informations on jet transport technique are reviewed in this section. More fundamental details about Jet Transport formulation, as well as its computer implementation can be found in [16, 17]. 3.1. Jet transport nonlinear expansion Consider an ODE system ˙x=f(x,p, t), such as the ones given by (1) or (5), with a n-dimensional state vector x, and a l-dimensional parameter vector p, as well as the initial conditions x(t0) = x0, and p(t0) = p0. The JT method performs the propagation of a neighborhood N0around x0and p0from t0up to a final time td. A sketch of the procedure is the following: first parameterize the initial neighborhood N0by a pair of polynomial vectors [x0] = ¯x0+δx0and [p0] = ¯p0+δp0, where ¯x0and ¯p0are the initial nominal values, while δx0and δp0are their initial uncertainties. Then, the initial state and parameter uncertainty neighborhood N0is integrated in the accurate nonlinear vectorfield, that includes the Earth potential and perturbations, from the initial epoch t0up to a final epoch td, using a JT-based eighth order variable step Runge 8 Kutta method, or by means of any other numerical integrator. The result is aκ-th order state polynomial vector, [x1]κ=¯xd+Pκ xd(δx0, δp0), which approximates the state uncertainty neighborhood Φ(td;t0,x0+δx0,p0+δp0), where Φdenotes the flow associated to the differential equation, and a similar representation, [pd] = ¯pd+Pκ pd(δp0), that approximates the parameter uncertainty neighborhood at time td. In particular, we note that the truncated state Taylor series at time tdcan be written as [xd]κ=¯xd+Pκ xd(δx0, δp0) =P 0≤γ1+···+γn+l≤κ aγ1···γn+lδxγ1 0,1···δxγn 0,n δpγn+1 0,1···δpγn+l 0,l ,(8) where δx0= (δx0,1, . . . , δx0,n)T,δp0= (δp0,1, . . . , δp0,l)T, while aγ1...γn+l=1 γ! ∂γΦxd ∂δxγ1 0,1···∂δxγn 0,n ∂δpγn+1 0,1···∂δpγn+l 0,l , are the coefficients of the Taylor expansion of the solution flow, and γ= γ1+···+γn+l. Clearly, the accuracy of (8) is mainly affected by the expansion order κ, the size of δx0and δp0, and the total propagation time ∆t=td− t0. In principle, for a given uncertainty neighborhood, the approximation accuracy can reach a given desired tolerance by tuning the expansion orders; a detailed analysis about the procedure is given in [14]. Note that if κ= 1, the JT method degrades to the linear approximation of the flow around the nominal orbit, and the usual state transition matrix is given by the coefficients aγ1...γn+lassociated to the first order of the expansion (γ1+···+ γn+l= 1). 3.2. JT-based Runge Kutta integrator In order to perform a JT-based numerical integration, all the arithmetic operations of the usual numerical propagator must be replaced with the corresponding polynomial algebra in the JT scheme. After each time step, the state xand the right hand side of the ODE system ˙ x=f(x, t) is expressed as a kth-order expansion around a certain point of the nominal solution. At this point, it is worth to remark a particularity of the JT-based numerical integrators with adaptive step size control, such as the Runge-Kutta-Fehlberg procedures. The Runge-Kutta methods implemented in the JT shceme (JTRK) are the only ones we consider in our study. These methods can be divided in two 9 zd −,k+1 =X 0≤γ1+···+γn+l≤κ cd γ1...γn+lEδXγ1 k,1···δXγn+l k,n+l, Px,i1i2 −,k+1 =P 1≤γ1+···+γn+l≤κP 1≤˜γ1+···+˜γn+l≤κ ai1 γ1...γn+lai2 ˜γ1...˜γn+l· Enδxγ1+˜γ1 k,1···δxγn+˜γn k,n δpγn+1+˜γn+1 k,1···δpγn+l+˜γn+l k,l o −δmx,i1 k+1δmx,i2 k+1 +Ewi1 kwi2 k, i1, i2∈[1, n], Pxp,ij −,k+1 =P 1≤γ1+···+γn+l≤κP 1≤˜γn+1+···+˜γn+l≤κ ai γ1...γn+lbj ˜γn+1...˜γn+l· Enδxγ1 k,1···δxγn k,nδpγn+1+˜γn+1 k,1···δpγn+l+˜γn+l k,l o −δmx,i k+1δmp,j k+1, i ∈[1, n], j ∈[1, l], Pp,j1j2 −,k+1 =P 1≤γn+1+···+γn+l≤κP 1≤˜γn+1+···+˜γn+l≤κ bj1 γn+1...γn+lbj2 ˜γn+1...˜γn+l· Enδpγn+1+˜γn+1 k,1···δpγn+l+˜γn+l k,l o−δmp,j1 k+1δmp,j2 k+1 +Evj1 kvj2 k, j1, j2∈[1, l], where δmx,i k+1 =¯xi k+1 −xi −,k+1, δmp,i k+1 =¯pi k+1 −pi −,k+1, and the covariance matrix of the augmented state is P−,k+1 =Px −,k+1 Pxp −,k+1 Ppx −,k+1 Pp −,k+1 . •Update equations: incorporting the new measurement zk+1 at time tk+1 into the estimation algorithm, Pd1d2 zz,k+1 =P 1≤γ1+···+γn+l≤κP 1≤˜γ1+···+˜γn+l≤κ cd1 γ1...γn+lcd2 γ1...γn+l· Enδxγ1+˜γ1 k,1···δxγn+˜γn k,n δpγn+1+˜γn+1 k,1···δpγn+l+˜γn+l k,l o −δnd1 k+1δnd2 k+1 +Eud1 k+1ud2 k+1, d1, d2∈[1, m], (28) Pid xz,k+1 =P 1≤γ1+···+γn+l≤κP 1≤˜γ1+···+˜γn+l≤κ ai γ1...γn+lcd ˜γ1...˜γn+l· Enδxγ1+˜γ1 k,1···δxγn+˜γn k,n δpγn+1+˜γn+1 k,1···δpγn+l+˜γn+l k,l o −δmx,i k+1δnd k+1, i ∈[1, n], d ∈[1, m], (29) 16 Pjd pz,k+1 =P 1≤γn+1+···+γn+l≤κP 1≤˜γ1+···+˜γn+l≤κ bj γn+1...γn+lcd ˜γ1...˜γn+l· Enδx˜γ1 k,1···δx˜γn k,nδpγn+1+˜γn+1 k,1···δpγn+l+˜γn+l k,l o −δmp,j k+1δnd k+1, j ∈[1, l], d ∈[1, m], (30) where δnd k+1 =¯zd k+1 −zd −,k+1. Using (29) and (30), we can obtain an augmented matrix PXz,k+1 = [PT xz,k+1 PT pz,k+1]T. Finally, we can update the estimated state and parameter vectors, as well as the covariance matrix, incorporating the measurement information at time tk+1. So, at time tk+1 Kk+1 =PXz,k+1 (Pzz,k+1)−1, X+,k+1 =X−,k+1 +Kk+1 (zk+1 −z−,k+1), P+,k+1 =P−,k+1 −Kk+1Pzz,k+1KT k+1. It is clear that, due to the nonlinearity of the system, at time tk+1 the distributions of the state and parameter vectors in the JT-HEKF-κfilters (excluding JT-HEKF-1) are no longer Gaussian, even that their distributions at time tkwere Gaussian. Although reference [7] showed that the Gaussian assumption is accurate enough in the orbit determination application using the proposed JT-HEKF-κfilters, in a forthcomming work we will remove this assumption. Under the Gaussian assumption, Isserlis’s formula [18] provides an analytical way to compute the associated expectation values in the JTHEKF-κalgorithms. As it has already been said, the purpose of this paper is to show how some parameters can be accurately estimated jointly with the spacecraft state. The nature of the parameters can be very different, from physical ones related with the spacecraft to the uncertain ones associated with the measurement procedure. Without loss of generality, and in order to show the behavior of the developed methodology, here we only present two illustrative cases: one physically affected by the time-varying area-to-mass ratio, and another one related with an uncertain tracking station position. 3.3.1. Case study A: spacecraft state and physical parameter estimation Although many spacecraft physical parameters are pre-designed and precisely manufactured for concrete missions, still some of them may be determined or adjusted during the mission, and need to be estimated in real time. Examples of such ones could be the spacecraft mass and the illuminated 17 cross sectional area. The uncertainty on the mass is mainly originated from unaccurate operations of the propulsion system, while the uncertainty on the illuminated cross sectional area is affected by the shape and attitude of the spacecraft, as well as the relative position between the spacecraft and the Sun. In general, the modeling error of the SRP acceleration is the largest dynamical error source in the GEO regime [10]. Morever, due to maneuvers, the area-to-mass ratio η=A/m changes with time and needs refitting. In this case study we assume that we are inside a period of time without maneuvers, but we have a time varying ηbecause the non-spherical GEO spacecraft is rotating with respect to the Sun in a diurnal basis, and so, its cross sectional area varies accordingly. The nonlinear variation law considered in the simulation is assumed to be ˙η=η0sin(ωet) where ωestands for the average angular speed of the Earth’s rotation. 3.3.2. Case study B: spacecraft state and tracking station position estimation For usual data collection techniques, underlying systematic biases and unmodeled observation errors significantly degrade the measurement quality. These facts limit the orbit determination accuracy, especially for the GEO case, where one arcsecond of angular error from the ground tracking station position, corresponds to approximately 200 meters of the spacecraft position observational error. Therefore, the precise estimation of these systematic biases and unmodeled observation errors, as well as the offset to remove them, is of significant interest and importance. In particular, the accurate estimation of the height of a removable ground tracking station is meaningful, since it is difficult to be accurately determined via the GPS services. In this case study we attempt to precisely estimate the uncertain position of a ground tracking station. In general, the position of the ground tracking station is provided by specifying its east longitude Λ, geodetic latitude φ, and an elevation Habove the WGS84 ellipsoidal surface, and the position vector has to be converted into the ECI frame. Unlike the previous case, (19) is now fixed and simplifies into (25). Therefore, assuming that the height of the ground tracking station is denoted by Hand its corresponding uncertainty by δH, the parameter equation can be expressed as Hk+1 =Hk+δHk+vk.(31) 18 3.3.3. Sensitivity analysis Complementing the previous discussions for the case studies A and B, we focus on the sensitivity analysis of the JT-HEKF-κfilters with respect to initial estimation errors, measurement acquisition frequencies, and the observation geometry between the spacecraft and the ground tracking station (i.e., ∆Λr=λ−Λ). In general, the classical EKF filtering algorithms lose accuracy in the low measurement acquisition frequency case. Therefore, three different measurement acquisition frequencies are considered in order to produce some significant conclusions. On the other hand, we analyze the sensitivity of JT-HEKF-κalgorithms with respect to the initial simulation conditions subjected to a large estimation error distribution (i.e., large values of σr,σv, and σp) via Monte Carlo (MC) simulations. It is worth mentioning that we omit the simulation with small initial estimation error distribution case, since both the usual EKF and the JT-HEKF proposed in this paper work very well for the simultaneous estimation of the spacecraft state and additional parameters in the GEO regime. Finally, a sensitivity analysis with respect to the observational geometry (25 uniformly distributed sampling points are taken from the interval ∆Λr∈[−60◦,−60◦]) is also implemented via MC simulations. In order to assess the sensitivity tests of the JT-HEKF-κfiltering algorithms, some statistical indices are defined as follows, κ¯εNp i=PNp j=1 εj i Np , τκ2 κ1=κ2¯εNp i κ1¯εNp i , ¯ε κσNp i=sPNp j=1(εj i−κ¯εNp i)2 Np−1, ζκ2 κ1= ¯ε κ2σNp i ¯ε κ1σNp i . where εj i, i =r, v, p indicates the position error, velocity error and parameter error of the jth Monte Carlo simulation calculated at the steady stage of the filtering process, κ¯εNp iand ¯ε κσNp idenote the mean and the standard deviation of εj iaccounting for the expansion order κ,Npis the total number of Monte Carlo simulations. Note that τκ2 κ1reveal the error ratio between JT-HEKF-κ2 and JT-HEKF-κ1filtering algorithms. if 0 < τκ2 κ1<1, JT-HEKF-κ2algorithm features the better accuracy than JT-HEKF-κ1algorithm, otherwise is worse. The ratio ζκ2 κ1shows the dispersion of the estimation error around 19 the mean. If ζκ2 κ1>1, JT-HEKF-κ2algorithm possesses a higher dispersion than JT-HEKF-κ1algorithm, otherwise is smaller. 4. Results Currently, optical telescopes are typically exploited to track the GEO objects. As a passive data collection technique, optical telescopes enable to efficiently track numerous objects surreptitiously. The recent development of more powerful wide-field-of-view optical sensors possesses sub-arcsecond observation accuracy for the near-geosynchronous objects [12]. In particular, the state of the art accuracy of the most sophisticated optical sensors enables to reach 20 milli-arcseconds [19]. In the following simulations, we assume that the measurement errors on the topocentric right ascension (RA, α) and declination (δ) are Gaussian white noises with the same standard deviation 66.6 milli-arcseconds (corresponding to σm= 3.232 ×10−7rad and 14 m in the GEO regime). The position of the ground tracking station is known by specifying its east longitude Λ, geodetic latitude φ, and an elevation Habove the WGS84 ellipsoidal surface. The local sidereal time of the ground tracking station is so= Λ + GA(t). After the coordinate transformation, the position of the ground tracking station is expressed in the ECI frame as ro= (xo, yo, zo) = (Rccos φcos so, Rccos φsin so, Rssin φ) , where roindicates the inertial position vector of the tracking station, while Rcand Rsare two intermediate values depending on the equatorial radius, Re, flattening, f, and height of the tracking station, H. The spacecraft position vector relative to the ground tracking station is expressed as %=r−ro= (x−xo, y −yo, z −zo). Therefore, the measurement vector z= (α, δ) can be computed by the corresponding measurement equations,        α= arctan y−yo x−xo+u1, δ= arcsin z−zo ||%||2+u2. (32) In practice, these measurements are provided by optical telescopes or radars; here we numerically generate such inputs according to the following procedure: we assume a true nominal trajectory of the higher fidelity dynamical model (including 10 ×10 EGM96S gravity model, SRP, albedo, Sun and Moon gravitational perturbations) that is propagated forward to time tk, 20 Table 1: True initial values including spacecraft state, parameters, and the ground tracking station position Parameter Nominal value Initial state Nominal value Epoch 15 November 2015, 0:0:0 UTC x024487.8 km η=A/m 74/3300 m2/kg y034324.4 km Cr1.3 z00 km Λ 42.0516528◦˙x0-2.50298 km/s φ0.7293472◦˙y01.78568 km/s H1620 m ˙z00 km/s where the observational measurement is done; then, numerical measurement at this epoch is generated by adding some Guassian white noises into the spacecraft topocentric right ascension RA (α) and declination (δ), which is calculated in terms of the true spacecraft state from the above state propagation and the true position of the tracking station (32). Table 1 shows the true initial spacecraft state vector and parameters, as well as the true tracking station position considered in the simulations. In order to compare the performances of the Cartesian and GEO representations, the JT-HEKF-κfiltering algorithms are implemented in both models with the same simulation conditions. The performance of a JTHEKF-κfilter includes the estimation accuracy, computational burden and its robustness. In particular, the estimation error is defined as the difference between the estimated state and parameter vectors (i.e., x+,k,p+,k), computed by the JT-HEKF-κfiltering algorithms, and the true state and parameter vectors (i.e., xk,pk), obtained by the numerical integration of the state and parameter equations, this is, εx=x+,k −xk, and εp=p+,k −pk. Furthermore, we define the position error εr=pε2 x+ε2 y+ε2 zand velocity error εv=qε2 ˙x+ε2 ˙y+ε2 ˙z. For reasons of space and clarity, only the sensitivity analysis results for case study A are given, the same conclusions are also valid for case study B. All the codes are written in C++ compiled with gcc version 5.5.0 and implemented with a laptop processor Intel(R) Core(TM) i5-7300HQ under Linux. 4.1. Spacecraft physical parameter estimation In this study case we consider all the spacecraft physical parameter equations, the state and measurement equations are nonlinear. The performances 21 of the JT-HEKF-κfiltering algorithms, implemented in both coordinate representations, are assessed by the simultaneous estimation of the state xand the time-varying area-to-mass ratio η. We assume that the initial guess of the area-to-mass ratio is 0.1 m2/kg off from the true value given in Table 1, and its initial standard deviation is ση= 0.1 m2/kg. The initial estimation errors are of 100 km in the position vector components, and of 0.5 m/s2in the velocity vector components. The adopted initial covariance matrix is P+,0=Px +,00 0σ2 η,where Px +,0=1010 I3×30 00.25 I3×3. The total simulation time is 4 days, provided that 7 measurements can be done per night (10 hours per day), separated by regular time intervals. Let the angle measurement noise be a zero mean white noise with standard deviation σm= 3.232×10−7rad (corresponding to 14 m in the GEO regime). The measurement noise matrix Ris R=σ2 m0 0σ2 m.(33) Figures 1 and 2 exhibit the results of the simultaneous estimation of the area-to-mass ratio and the spacecraft state, respectively implemented in the Cartesian and GEO representations. Both figures show that the JT-HEKF-2 filter outputs a better accuracy than the JT-HEKF-1 filter (i.e., the classic EKF filter). For both JT-HEKF-2 and JT-HEKF-3 filters, the position and velocity estimation errors are respectively less than 10 m and 5 ×10−4m/s. The estimation error of the area-to-mass ratio converges to 2 ×10−4m2/kg, which is around 0.1% of the nominal value. In contrast, the JT-HEKF-1 filter produces the worse estimation errors of the spacecraft position, velocity, and area-to-mass ratio: respectively, 100 m, 3 ×10−3m/s and 3 ×10−3m2/kg (14%). In essence, the JT-HEKF-1 algorithm fails in the estimation of the area-to-mass ratio, as is shown in the bottom line sub-figures of Figs. 1 and 2. Both sub-figures illustrate that the area-to-mass ratio estimation error is almost 14% of the initial error, so the filter does not work properly. From this simulation, it follows that we can at least obtain one order of magnitude accuracy gain replacing the linear EKF (i.e., JT-HEKF-1 filter) by the JTHEKF-2 filter. Furthermore, due to normality hypothesis, the expectations of the third order terms vanish, thus no accuracy gain is obtained replacing the JT-HEKF-2 algorithm by the JT-HEKF-3 one. 22 10-4 10-2 100 102 104 0 12 24 36 48 60 72 84 96 εr (km) JT-HEKF-1 JT-HEKF-2 JT-HEKF-3 10-6 10-4 10-2 100 102 0 12 24 36 48 60 72 84 96 εv (m/s) 10-6 10-4 10-2 100 0 12 24 36 48 60 72 84 96 εA/m (m2/kg) time (h) Figure 1: Case A: State and area-to-mass ratio estimation errors with a measurement frequency 7 times/night. Implemented in the GEO representation. 10-4 10-2 100 102 104 0 12 24 36 48 60 72 84 96 εr (km) JT-HEKF-1 JT-HEKF-2 JT-HEKF-3 10-6 10-4 10-2 100 102 0 12 24 36 48 60 72 84 96 εv (m/s) 10-6 10-4 10-2 100 0 12 24 36 48 60 72 84 96 εA/m (m2/kg) time (h) Figure 2: Case A: State and area-to-mass ratio estimation errors with a measurement frequency 7 times/night. Implemented in the Cartesian representation. 23 Comparing Fig. 1 with Fig. 2, it follows that the accuracies of both representations are of the same order of magnitude but, in general, the accuracy of GEO representation is a little better than that of the Cartesian representation if the algorithms are convergent (i.e., JT-HEKF-2 and JT-HEKF-3). This phenomenon originates from the more accurate prediction of the GEO representation in the time-update step. Although the Cartesian representation behaves better than GEO representation in the JT-HEKF-1 algorithm, it makes no sense since the 1st order algorithm is essentially failing. Another important issue to be pointed out here is related with the computational efficiency. Table 2 shows the computational cost of the simulations corresponding to Figs. 1 and 2 (without including the computational cost for generating the measurements). It underlines that the GEO representation is much faster than the Cartesian one, since the GEO elements vary much slower than the classical Cartesian ones in the propagation of GEO orbits. Then, a larger integration step size is adopted in the prediction process, resulting in less computational time and in a higher efficiency. Therefore, one can conclude that GEO representation possesses a better performance (including a superior estimation accuracy and less computational burden) than the Cartesian representation. This conclusion is also validated in the case study B. But, for brevity, in what follows we only show the results of the GEO representation for the case study B, while the ones corresponding to the Cartesian representation are omitted. Table 2: Computational time (s) Model JT-HEKF-κ κ= 1 κ= 2 κ= 3 GEO 24 44 114 Cartesian 144 149 197 4.1.1. Sensitivity analysis with respect to the observational geometry Clearly, in order to estimate the spacecraft state and area-to-mass ratio, the GEO spacecraft must be in the visible region from the ground tracking station. The sensitivity analysis of the JT-HEKF-κalgorithms with respect to the observational geometry is implemented, in particular with respect to the longitude difference ∆Λrbetween the spacecraft and the tracking station. For this purpose, we take 25 uniform sampling points in the interval 24 10-4 10-2 100 102 104 0 12 24 36 48 60 72 84 96 εr (km) JT-HEKF-1 JT-HEKF-2 10-6 10-4 10-2 100 102 0 12 24 36 48 60 72 84 96 εv (m/s) 10-6 10-4 10-2 100 0 12 24 36 48 60 72 84 96 εA/m (m2/kg) time (h) Figure 3: Case A: Sensitivity analysis with respect to the observational geometry with a measurement frequency of 7 times/night. Implemented in the GEO representation. Table 3: Case A: Sensitivity analysis to the observational geometry Indices i=r i =v i =η Unit m 10−4m/s 10−4m2/kg JT-HEKF-1 1¯ε25 i36.7 27.6 7.96 ¯ε 1σ25 i24.1 421.9 438.5 JT-HEKF-2 2¯ε25 i7.88 6.45 1.50 ¯ε 2σ25 i2.08 76.8 81.7 Ratio τ2 1=2¯ε25 i/1¯ε25 i0.215 0.234 0.188 ζ2 1=¯ε 2σ25 i/¯ε 1σ25 i0.086 0.182 0.186 25 transport-based high order filtering algorithms could be suitable to carry out future applications related with high precise geostationary or geosynchronous orbit determination. Appendix To incorporate more measurement information, the range from the spacecraft to the tracking station is assumed to be available. The standard deviation of the measurement error for the range is supposed to be 1 meter. A Monte Carlo simulation is implemented to assess its performance. To compare the results with Figure 4, the JT-HEKF-κalgorithms are fed with the same measurements (excluding the range information), and initialized with the same initial estimates of the state and area-to-mass ratio. Figure 8 shows that the supplement of the range measurement information decreases the performance of JT-HEKF-1 filter, while improves the performances of the JT-HEKF-2 and JT-HEKF-3 filters. Tables 4 and 5 provide more detailed information. Their comparison further reveals that more measurement information makes the original divergent algorithm, this is JT-HEKF-1, diverge now more severely, but it improves the estimation accuracy of the original convergent methods (i.e. the JT-HEKF-2 and JT-HEKF-3 filters). Table 5: Case A: Sensitivity analysis with respect to initial estimation errors incorporating the range measurement information Indices i=r i =v i =η Unit m 10−3m/s 10−4m2/kg JT-HEKF-1 1¯ε25 i204.5 288.0 83.2 ¯ε 1σ25 i77.5 1953.0 2147.0 JT-HEKF-2 2¯ε25 i1.22 0.80 0.18 ¯ε 2σ25 i0.07 11.64 12.27 Ratio τ2 1=2¯ε25 i/1¯ε25 i0.006 0.003 0.002 ζ2 1=¯ε 2σ25 i/¯ε 1σ25 i0.001 0.006 0.006 32 10-4 10-2 100 102 104 0 12 24 36 48 60 72 84 96 εr (km) 10-6 10-4 10-2 100 102 0 12 24 36 48 60 72 84 96 εv (m/s) 10-6 10-4 10-2 100 102 0 12 24 36 48 60 72 84 96 εA/m (m2/kg) time (h) Figure 8: Case A: Sensitivity analysis with respect to initial estimation errors with a measurement frequency of 7 times/night. Implemented in the GEO representation and incorporating the range measurement information. Acknowledgments J.C. thanks the support of Doctorate Foundation of Northwestern Polytechnical University and the Chinese Scholarship Council. He also thanks National natural science foundation of China for the grant 11572248, MINECOFEDER for the grant MTM2015-65715-P and Key Laboratory of Aerospace Flight Dynamics for the grant 2015afdl015. J.J.M. thanks MINECO-FEDER for the grant MTM2015-65715-P and the Catalan government for the grant 2017SGR-1049. G.G. thanks the Catalan grant government 2017SGR-1374, and MINECO-FEDER for the grant MTM2016-80117-P. References [1] R. M. Fuster, M. F. Us´on, A. B. Ibars, Interferometric orbit determination for geostationary satellites, Science China Information Sciences 60 (6) (2017) 060302 (2017). doi:10.1007/s11432-016-9052-y. [2] S. Ribo, J. C. Arco, E. Cardellach, S. Oliveras, A. Rius, C. Buck, Preliminary results of digital satellite tv opportunity signals scattered on 33 the sea-surface, in: Workshop on Reflectometry using GNSS and Other Signals of Opportunity (GNSS+R), Purdue University, West Lafayette, USA, October 10-11, 2012 (October 10-11, 2012). [3] M. Valli, R. Armellin, P. di Lizia, M. R. Lavagna, Nonlinear filtering methods for spacecraft navigation based on differential algebra, Acta Astronautica 94 (1) (2014) 363 – 374 (2014). doi:10.1016/j.actaastro.2013.03.009. [4] R. Armellin, P. Di Lizia, R. Zanetti, Dealing with uncertainties in anglesonly initial orbit determination, Celestial Mechanics and Dynamical Astronomy 125 (4) (2016) 435–450 (2016). doi:10.1007/s10569-016-9694-z. [5] F. Cavenago, P. Di Lizia, M. Massari, A. Wittig, On-board spacecraft relative pose estimation with high-order extended kalman filter, Acta Astronautica 158 (2019) 55–67 (2019). [6] E. A. Wan, R. Van Der Merwe, The unscented kalman filter for nonlinear estimation, in: Adaptive Systems for Signal Processing, Communications, and Control Symposium 2000. AS-SPCC. The IEEE 2000, Ieee, 2000, pp. 153–158 (2000). doi:10.1109/ASSPCC.2000.882463. [7] R. S. Park, D. J. Scheeres, Nonlinear semi-analytic methods for trajectory estimation, Journal of guidance, control, and dynamics 30 (6) (2007) 1668–1676 (2007). doi:10.2514/1.29106. [8] J. L. Crassidis, J. L. Junkins, Optimal estimation of dynamic systems, CRC press, Boca Raton, FL, 2004, Ch. 5, pp. 243–320 (2004). [9] O. Montenbruck, E. Gill, Satellite orbits: models, methods and applications, Springer Science & Business Media, 2012, Ch. 3, pp. 53–81 (2012). [10] C. Tang, X. Hu, S. Zhou, R. Guo, et al., Improvement of orbit determination accuracy for beidou navigation satellite system with two-way satellite time frequency transfer, Advances in Space Research 58 (7) (2016) 1390–1400 (2016). doi:10.1016/j.asr.2016.06.007. [11] E. M. Soop, Handbook of geostationary orbits, Kluwer Academic, Dordrecht, The Netherlands, 1994, Ch. 2, pp. 13–26 (1994). 34 [12] J. Tombasco, Orbit estimation of geosynchronous objects via groundbased and space-based optical tracking, Ph.D. thesis, University of Colorado at Boulder (2011). [13] J. Tombasco, P. Axelrad, M. Jah, Specialized coordinate representation for dynamic modeling and orbit estimation of geosynchronous orbits, Journal of guidance, control, and dynamics 33 (6) (2010) 1824–1836 (2010). doi:10.2514/1.48903. [14] J. Chen, J. J. Masdemont, G. G´omez, J. Yuan, High accuracy state and parameter estimation of geo-stationary satellites using jet transport based nonlinear filtering algorithm, in: 69th International Astronautical Congress, Bremen, Germany, 2018 (2018). [15] M. Berz, Modern map methods in particle beam physics, Vol. 108, Academic Press, 1999, Ch. 2, pp. 81–111 (1999). [16] D. P´erez-Palau, Dynamical transport mechanisms in celestial mechanics and astrodynamics problems, Ph.D. thesis, Universitat de Barcelona (2016). [17] D. P´erez-Palau, J. J. Masdemont, G. G´omez, Tools to detect structures in dynamical systems using jet transport, Celestial Mechanics and Dynamical Astronomy 123 (3) (2015) 1–24 (2015). [18] L. Isserlis, On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables, Biometrika 12 (1/2) (1918) 134–139 (1918). doi:10.2307/2331932. [19] A. T. Tokunaga, New generation ground-based optical/infrared telescopes, in: Encyclopedia of the Solar System (Third Edition), Elsevier, 2014, pp. 1089–1105 (2014). doi:10.1016/B978-0-12-415845-0.00051-7. 35