Xiucong Sun

Pei Chen
email: chenpei@buaa.edu.cn

Chao Han

Christophe Macabiau

G Taes-201500387-R1

Autonomous Orbit Determination via Kalman Filtering of Gravity Gradients

Keywords: Autonomous orbit determination, Gravity gradient, Extended Kalman filter, Augmented state filter, GOCE I

Spaceborne gravity gradients are proposed in this paper to provide autonomous orbit determination capabilities for near Earth satellites. The gravity gradients contain useful position information which can be extracted by matching the observations with a precise gravity model. The extended Kalman filter is investigated as the principal estimator. The stochastic model of orbital motion, the measurement equation and the model configuration are discussed for the filter design. An augmented state filter is also developed to deal with unknown significant measurement biases. Simulations are conducted to analyze the effects of initial errors, data-sampling periods, orbital heights, attitude and gradiometer noise levels, and measurement biases. Results show that the filter performs well with additive white noise observation errors. Degraded observability for the along-track position is found for the augmented state filter. Real flight data from the GOCE satellite are used to test the algorithm. Radial and cross-track position errors of less than 100 m have been achieved.

In early researches the gradiometer was considered as a useful aid to the inertial navigation system (INS) [START_REF] Britting | The impact of gradiometer techniques on the performance of inertial navigation systems[END_REF][START_REF] Grubin | Accuracy improvement in a gravity gradiometer-aided cruise inertial navigator subjected to deflections of the vertical[END_REF][START_REF] Metzger | Inertial navigation performance improvement using gravity gradient matching techniques[END_REF]. Two integration methods, gravity disturbance compensation and gravity gradient map matching, were studied. The former made use of gravity gradients combined with the velocity information to estimate gravity disturbance forces, whereas the latter utilized a stored map for observation data matching and then to estimate the trajectory directly. In 1990, Affleck and Jircitano [START_REF] Affleck | Passive gravity gradiometer navigation system[END_REF] presented a passive navigation system, in which the measured gravity gradients were compared with map values and the differences were processed by an optimal filter to correct INS errors. A parametric study was conducted for an airborne INS. Position accuracies of better than 100 m were demonstrated. Gleason [START_REF] Gleason | Passive airborne navigation and terrain avoidance using gravity gradiometry[END_REF] continued the work and developed a fast Fourier transform algorithm to efficiently generate constant altitude grids of reference gravity gradients. More recently, Richeson [START_REF] Richeson | Gravity gradiometer aided inertial navigation within non-GNSS environments[END_REF] presented a comprehensive discussion on the gravity gradient map matching technique for inertial navigation aiding. Monte Carlo simulations of a hypersonic cruise showed that sub-meter position errors could be possible with a future grade gradiometer.

Over the past few decades, gravity models of terrestrial planets, especially the Earth, have been improved dramatically. Derived from the combination of satellite geodetic data with high-resolution gravitational information collected from surface gravimetry, the new developed Earth Gravitational Model 2008 (EGM2008) is complete to degree 2190 and order 2159 [START_REF] Pavlis | The development and evaluation of the earth gravitational model 2008 (EGM2008)[END_REF]. Meanwhile, significant advances in spaceborne gravity gradiometry have occurred due to geophysical activities [START_REF] Difrancesco | Gravity gradiometer systems-advances and challenges[END_REF]. The electrostatic gravity gradiometer (EGG) implemented on ESA's gravity field and steady-state ocean circulation explorer (GOCE) achieved a noise density level of 0.01 E/√Hz within measurement bandwidth (MBW). A further better accuracy of 10 -3 E/√Hz is projected in the future using cold atom interferometers [START_REF] Carraz | A spaceborne gravity gradiometer concept based on cold atom interferometers for measuring Earth's gravity field[END_REF].

The progresses in global gravity field modeling and spaceborne gravity gradiometry provide opportunities for the application of gravity gradiometry to spacecraft navigation. Chen [START_REF] Chen | Gravity gradient tensor eigendecomposition for spacecraft positioning[END_REF] introduced an idea of using full-tensor gravity gradients combined with precise inertial attitudes to determine position. The system's observability was explained from geometry and the effects of possible sources of errors were analyzed. An eigendecompostion based positioning algorithm using the J2 gravity model was developed.

Simulation results showed that a mean position error of 421 m could be possible for a spacecraft with a gradiometer having a noise level of 0.1 E at 300 km altitude. Sun [START_REF] Sun | Test of GOCE EGG data for spacecraft positioning[END_REF] investigated the use of a least squares searching method which employed a high-degree gravity model for map matching. Real flight data from the GOCE satellite were used to test and verify the method and a mean positioning error of 620 m was achieved.

In this study, autonomous orbit determination using gravity gradients within an extended Kalman filter (EKF) is presented. The main advantage of the EKF-based orbit determination is that the measurement noise can be reduced through incorporation of the orbital motion. Moreover, spacecraft velocity can be output simultaneously. Within the EKF, a 120th degree and order EGM2008 model is used to calculate gravity gradients and high-precision inertial attitudes are used to remove the contribution of instrumental orientation. An augmented state EKF is also developed to deal with significant measurement biases. The algorithms are applied to both simulated data and real GOCE data and effects of several important factors are analyzed.

The remainder of this paper is organized as follows. Section II introduces the mathematical model of gravity gradients. Section III presents the EKF-based orbit determination algorithm, including orbital dynamic modeling, the measurement equation, the estimation formulas, and the model configuration. Section IV presents the metrics used for evaluating filter performance.

Simulation results are given in Section V. The effects of initial errors, data-sampling periods, orbital heights, noise levels, and measurement biases are discussed. Section VI describes the GOCE satellite whose flight data have been used to test the algorithm, and the orbit determination results are presented. The major finding of this study is concluded in Section VII.

II. GRAVITY GRADIENT MODELING

The gravity gradients (or gravitational gradients) are the second-order spatial derivatives of the gravitational potential and form a 3 × 3 matrix called the gravity gradient tensor (GGT)

  2 2 2 2 2 2 2 2 2 2 2 2 U U U x y x z x xx x y x z U U U yx yy yz y x y z y U U U zx z y z z z x z y z V V V U V V V V V V                                                    V (1)
The unit of GGT is Eötvös, denoted by the symbol E. One Eötvös equals 10 -9 1/s 2 in SI units. There are only five independent terms in matrix V. The continuity of the gravitational field guarantees that V is symmetric and the Laplace's equation constrains its trace to be zero [START_REF] Hofmann-Wellenhof | Physical Geodesy[END_REF].

The expression of GGT depends on the choice of reference system. Expressions of GGT in two different reference systems have the following relationship

  T b b b a a a b a a a b   V C V C C V C (2)
where the symbols a and b denote the reference frames, the superscript T denotes the transpose operation, and b a C is the transformation matrix from frame a to b.

The common reference frame used in most of the Earth gravity models (EGM96, EGM2008, GGM03, etc.) is the Earth-centered Earth-fixed (ECEF) frame, one realization of which is the International Terrestrial Reference Frame (ITRF) developed by the International Earth Rotation and Reference Systems Service (IERS). The ECEF frame cannot be physically materialized onboard a satellite. The gravity gradients can only be measured in the Gradiometer Reference Frame (GRF), which is defined by the three orthogonal baselines of a gradiometer. The transformation matrix from ECEF to GRF will be specified in Section III for the measurement equation formulation.

A series of spherical harmonics is used in geodesy to model the gravitational potential

TAES-201500387-R1         0 0 , , sin cos sin n n E nm nm nm n m R GM U r P C m S m r r                       (3)
where r, ϕ, and λ are the geocentric distance, latitude, and longitude of the position (expressed in ECEF), GM is the geocentric gravitational constant, E R is the reference equatorial radius of the Earth, n and m are the degree and order of the normalized coefficients nm C and , nm S and nm P is the normalized Legendre function. The values of GM as well as the coefficients, nm C and , nm S are determined from observations such as satellite tracking and surface gravimetry. The gravity gradients can be obtained by double differentiation of Eq. ( 3). The calculation formulas are given in Appendix.

The spherical harmonic model is usually truncated at a maximum degree for the practical computation. The precision of gravity models with respect to GGT truncated at different degrees is shown in Fig. 1, in which the 300th degree and order EGM2008 model is used as a benchmark. The reference height is 300 km. It can be seen that the mean GGT error decreases with increase of the truncating degree. The effects of tidal variations (solid Earth tide, ocean tide, and pole tide) as well as the third body attraction (Sun and Moon) on GGT modeling are insignificant and are not considered.

III. ORBIT DETERMINATION ALGORITHM

Orbital motion of spacecraft is subject to various random perturbations. The orbit determination problem is to estimate the states (position and velocity) of such a stochastic dynamic system from noisy observations, which in this study refers to the gravity gradient measurements. The EKF is a suboptimal filter which gives real-time estimates of the orbital states and the error covariances based on sequential observations. This section describes the dynamic model and the measurement equation and presents details of the estimation algorithm.

A. Stochastic Model of Orbital Motion

The orbital motion is modeled as a nonlinear stochastic dynamic system and is described by the following stochastic differential equation  

, ,

t d t dt                     v 0 r f r v w v ( 4
)
where the state vector comprises the position and velocity (r and v) in the ECI frame, f is a 3-dimensional vector function representing the accelerations due to the deterministic forces, and t w represents the remaining unmodeled perturbation accelerations and is assumed to be a zero-mean stationary white Gaussian noise.

The continuous time dynamic model needs to be discretized before being used for the state prediction in the EKF, since the measurements are taken at discrete instants. The discretized dynamic equation has the following form

  1 1 1 , , k k k k k t t      x φ x w (5)
where k x is the orbital state at time , k t φ is a 6-dimensional vector function relating states at two adjacent instants, and x is called the state transition matrix, denoted as   1 , .

k k t t  Φ   1 , k k t t  Φ
can be obtained by numerically integrating the following differential equation from 1

k t  to k t       3 3 3 3 1 1 1 1 3 3 3 3 , , k k k k d t t t t t t dt                      f r 0 I Φ Φ 0 ( 6
)
with the identity matrix as the initial value  

1 1 3 3 , k k t t     Φ I (7)
The covariance matrix of 1 k  w is called the process noise matrix, denoted as 1 .

k  Q 1 k  Q can be obtained by numerically integrating the following differential equation from 1 k t  to k t     3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 t t T w w d dt                                         f f r r 0 I 0 I 0 0 Q Q Q 0 σ σ 0 0 ( 8
)
with zero initial values. i.e., Eq. ( 6) and Eq. ( 8), are solved in parallel with the integration of the orbit trajectory.

TAES-201500387-R1 6

B. Measurement Equation

A full-tensor gravity gradiometer is considered in this study to provide gravity gradients for orbit determination. The GRF frame is assumed to be always aligned with the satellite RSW (radial, transverse, and normal) reference frame. The X-axis is in the flight direction and is perpendicular to the radius vector (along-track), the Y-axis is normal to the orbit plane (cross-track), and the Z-axis is radially downwards to the Earth (radial). An onboard star tracker is assumed to simultaneously provide high-precision inertial attitudes, which are used to compute the transformation matrix from ECI to GRF. In addition, the IERS Conventions [START_REF] Petit | IERS Conventions 2010[END_REF] defines precise models for coordinate transformation from ECI to ECEF. Thus the transformation matrix from ECEF to GRF can be determined directly without a prior knowledge of position and velocity

g g i e i e  C C C (9)
where the symbols i, e and g denote the frames ECI, ECEF, and GRF, respectively. Referring to Eq. ( 2), the GGT in GRF can be obtained by

 

T g e g e e g

g i e i i e e i g

  V C V C C C V C C ( 10
)
where e V is the GGT in ECEF and can be calculated using the formulas in Appendix A.

The gravity gradients measured by the gradiometer are output as a column vector

      , , , , , , , , , 1 2 1 2
V V V V V V V V V                            V (11)
Thus the measurement equation is expressed as

    , k g k k k k k t      z V ξ h x ξ ( 12
)
where k z is the 6-dimensional measurement vector at time , k t h is a 6-dimensional vector function, and k ξ is the observation error. The covariance matrix of k ξ is called the measurement noise matrix, denoted as .

k k k k k            r v z H H H x (13)
Since the expression of gravity gradients does not explicitly contain velocity, the partials of  

c c c c c c c c c c c c c c c c c c c c c c c                         , , , , , , , , , , , , , , , , , , e xx
V V V x y z V V V x y z V V V x y z e V V V x y z V V V x y z V V V x y z                                                               

T

where g e Π is the transformation matrix from e  V to , error depends on the accuracies of the attitude measurements and the IERS models. From recent geodetic reports, the relative accuracy of the IERS Earth orientation model is better than 10 -8 [START_REF] Petit | IERS Conventions 2010[END_REF]. According to the analysis in Chen [START_REF] Chen | Gravity gradient tensor eigendecomposition for spacecraft positioning[END_REF], the corresponding GGT transformation error is smaller than 0.1 mE (1 mE = 10 -3 E) and can be neglected. By contrast, the effect of attitude errors is much more significant. An attitude error of 1 arcsec can cause a GGT transformation error of about 0.01 E at 300 km height. The vector-valued GGT transformation error can be written as

g  V   , , 1,2,
1 2 3 t k k                  ξ M ( 15
)
where 1 2

, ,     and 3   are Euler angle errors and k M is a 6 × 3 mapping matrix. The derivation of Eq. ( 15) is found in Appendix B. Suppose that the Euler angle errors are white Gaussian noises and the noise standard deviation is denoted as .

 

The covariance of t k ξ can be given by

2 t k T k k    ξ R M M (16)
The gradiometer noise is a function of the accelerometer noise. Take the electrostatic gravity gradiometer for example. Three pairs of accelerometers are mounted at the ends of three orthogonal baselines [START_REF] Cesare | Performance requirements and budgets for the gradiometric mission[END_REF]. Specific forces measured by each pair of accelerometers are differenced to provide gravity gradients along that direction. Suppose that the six accelerometers are identical and the noise standard deviation is denoted as 0 .

 Since the accelerometer noises are mutually independent, the covariance of the

gradiometer noise n k ξ is a diagonal matrix 2 0 2
2 0 0 0 0 0 0 2 0 0 0 0 0 0 2 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1

n k l                       ξ R ( 17
)
where l is the length of the gradiometer baseline. The measurement noise matrix is finally obtained by

t n k k k   ξ ξ R R R (18)

C. Extended Kalman Filter

The extended Kalman filter is a linear minimum mean square error (LMMSE) estimator which is based on linearization of nonlinear systems and consists of a recursive calculation of the approximate conditional mean and covariance of the state. The linearization of the dynamic and measurement equation is implemented by the Taylor series expansion in which the second-and higher-order terms are neglected, yielding a first-order EKF.

The recursive estimation algorithm of the EKF is given as follows

  1 1 ˆ, , k k k k t t    x φ x (19)     1 1 1 1 , , T k k k k k k k t t t t       P Φ P Φ Q (20) T k k k k k   S H P H R (21) 1 k k k k   K P H S (22)   ˆ, k k k k k k t        x x K z h x (23)     ˆT T k k k k k k k k k     P I K H P I K H K R K (24)
where ˆk x and ˆk P are the estimated conditional mean and covariance at k t , k x is the predicted state to k t from 1 ,

k t  k P and k S
are the state and measurement prediction covariances, and k K is the filter gain.

The basic EKF algorithm presented above only deals with the case in which the observation errors are additive white noises. If unknown biases are contained in the observations, an augmented state extended Kalman filter (ASEKF) needs to be developed in order to estimate the biases. The augmented state vector consists of the orbital states as well as the biases

k a k k        x x B ( 25
)
where a k x is the augmented state vector, k B is a 6-dimensional bias vector. The dynamic evolution of the bias is modeled as a random walk

1 , 1 k k k     B B B w (26)
where

, 1 k  B w
is the bias process noise and its covariance B Q is determined based on the characteristics of the bias. Even for estimating constants, a small fictitious process noise has to be used to prevent divergence of the filter.

The measurement equation for the ASEKF is rewritten as

    , , k g k k k k k k k t       z V B ξ h x B ξ ( 27
)
And the augmented measurement partial matrix becomes

    , 6 3 6 6 6 12 6 3 a k k k a k             r z H H 0 I x ( 28
)
where the identity matrix 6 6

 I represents the partials of k z with respect to .

k

B

The recursive estimation algorithm for the ASEKF can be obtained by replacing k x in Eqs. ( 19)-( 24) with the augmented vector .

a k

x In addition, to balance the matrix ,

a k
H the algorithm uses mE as the unit of gravity gradients.

D. Model Configuration

The model configuration of the EKF includes the definition of deterministic forces for orbital motion and the choice of gravity models for GGT computation. The deterministic forces determine the accuracy of the orbital dynamic model. Its configuration should take into account the computational efficiency as well as the accuracy of the measurement system. For the basic EKF, only the central force and the J 2 perturbation are considered. The analytical expression of f is

TAES-201500387-R1         2 2 2 2 2 2 3 2 2 2 3 1 5 s i n 1 2 3 , , 1 5sin 1 2 3 1 5 s i n 3 2 E E E R J r R GM t J r r R J r                                                  f r v r (29)
where J 2 is the second zonal harmonic coefficient.

For the ASEKF processing model, the gravitational forces up to degree 20 and order 20 are used. The accelerations due to higher degree potential coefficients, lunar and solar gravitational attractions, and the non-gravitational forces (atmospheric drag and solar radiation pressure) are included into the process noise .

t w The standard deviation of t w is required to appropriately reflect the accuracy of the dynamic model. Numerical simulations have been conducted to evaluate the model accuracy. It is found that standard deviations of 0.01 m/s 2 and 5 × 10 -4 m/s 2 can be adopted in the basic EKF and the ASEKF, respectively. The accuracy requirement for the state transition matrix evaluation is not stringent. Only the central and J 2 terms are considered for both the two filters.

The choice of gravity models for GGT computation depends on the noise level of gradiometers. A 120th degree and order EGM2008 model is used in this study. As shown in Fig. 1, the mean GGT error of this model is on the order of 0.001 E, which is smaller than the noises of most current gradiometers (ranging from 0.01 E to 10 E). Thus, the GGT model errors can be neglected in the filter design. The accuracy requirement for the measurement partial matrix evaluation is also not stringent and only the central and J 2 terms are considered.

IV. FILTER PERFORMANCE EVALUATION METRICS

Several criteria are used to evaluate the performance of the orbit determination filter. The first consideration is the filter accuracy.

The estimated states can be compared with true values or more accurate values to obtain time-domain error curves, which imply the actual accuracy of the filter. The conditional covariance matrix measures the goodness of the estimate in a probabilistic sense and represents the predicted accuracy of the filter. In this study, the orbital states are estimated in the ECI frame, whereas the position and velocity errors and their predicted variances are expressed in the RSW frame. The total 3-dimensional (3D) position and velocities errors are also used to compare the orbit determination accuracy between different simulation cases.

The filter consistency is another important metric for performance [START_REF] Bar-Shalom | Estimation with Applications to Tracking and Navigation: Theory Algorithms and Software[END_REF]. The consistency consists of two conditions

  k k E   x x 0 (30)    k k k k k E        x x x x P (31)
Eq. ( 30) is the unbiasedness requirement. Eq. ( 31) is the covariance matching requirement, which means that the actual errors should match the filter predicted variances. The consistency is evaluated by the chi-square test. Define the normalized estimation error squared (NEES) as

      1 ˆT k k k k k k      x x P x x (32) Under hypothesis H 0 that the EKF is consistent,   k  is chi-square distributed with x n degrees of freedom   2 ~x n k   (33)
where x n is the dimension of the state vector x. H0 is accepted if   k  does not exceed the upper limit of a probability region

    2 1 x n k      ( 34
)
where α is the significance level and is set to 5% in this study.

Observability can be inferred from the predicted variances of filter. The observability condition guarantees a "steady flow" of information about the state components and leads to steady-state errors. The steady-state variances are not unique and are affected by many factors, such as data sampling rates and measurement noise levels. The goodness of filter observability can also be inferred from the convergence speed of the filter, i.e., the time of convergence to the steady-state errors.

V. SIMULATION RESULTS AND DISCUSSIONS

Numerical simulations have been conducted to test the performance of the EKF for GGT-based autonomous orbit determination.

The basic EKF cases cover a 6-hour data arc starting from October 1, 2014, 12:00:00.0 (UTC), whereas the ASEKF case covers a 40-hour data arc in order to observe the convergence of the biases. The true orbit ephemerides are generated using a high-precision numerical orbit simulator. The Adams-Bashforth-Moulton method for numerical integration, the EGM2008 model truncated at the 120th degree and order for non-spherical gravitational perturbation, the NRLMSISE-00 model for the atmospheric density, and the analytical formulas for the lunar and solar ephemerides are used in the simulator. The true gravity gradients are generated using a 300th degree and order EGM2008 model. Attitude errors and gradiometer noises at different levels are added to simulate the noisy measurements. In addition, constant biases are also simulated for the ASEKF test.

In the baseline case a nearly circular orbit having a height of 300 km is assumed. The initial osculating orbital elements are listed in Table I. The GGT measurements are simulated using a data-sampling period of 30 s. The standard deviations of the Euler angle The position error of the EKF baseline case can also be compared with the results obtained from the eigendecomposition method given in Chen [START_REF] Chen | Gravity gradient tensor eigendecomposition for spacecraft positioning[END_REF]. Under the similar simulation conditions, the mean position error of the eigendecomposition method is 421 m, which is on the same order of magnitude as the steady-state variance of the 3D position error. However, the EKF performs better than the GGT inversion method which used the J 2 gravity model to extract position. The improvement of the position accuracy is attributed to the incorporation of orbital motion and the use of a higher degree reference gravity model. Fig. 3. Radial, along-track, and cross-track velocity errors and the corresponding 3σ boundaries for the baseline case.

Fig. 4. Normalized state estimation error squared (NEES) and the 95% probability region for the baseline case.

A. Effects of Initial Errors and Data-Sampling Periods

Effects of the initial position and velocity errors are tested by comparing the baseline case with another two cases in which larger initial errors of [100 km, 100 km, 100 km, 100 m/s, 100 m/s, 100 m/s] and [1000 km, 1000 km, 1000 km, 1 km/s, 1 km/s, 1 km/s] are assumed, respectively. The other simulation conditions are set to be the same with the baseline case. The evolution of the total 3D position errors is depicted in Fig. 5. Increasing the initial errors slows down the convergence speed of the filter. From Fig. 5, the convergence time is 40 min for the 100 km and 100 m/s case, whereas for the 1000 km and 1 km/s case the convergence time is 80 min. In these two cases the filter consistency at the transient stage is not good. That is because the large initial errors increase the nonlinearity and the corrections at the first few epochs are too small. Nevertheless, the steady-state estimation errors are not affected by the initial errors, revealing that the EKF is asymptotically consistent.

Cases having different data-sampling periods are also simulated. The comparison results are summarized in Table II. It is found that increasing the sampling periods degrades the orbit determination accuracy. From Table II, the total position error of the 300-s sampling interval case is almost three times that of the 10-s interval case. The reason is that the observation information collected over the same time span gets less when using a larger sampling period. In addition, the weight of the orbital motion information in the EKF is also reduced, leading to white noise-like estimation errors. An even larger sampling interval of 1000 s has also been tested and the EKF filter still converges. It is shown that the navigation ability under the condition of sparse gravity gradient measurements can be guaranteed. However, an extremely large sampling period will cause divergence of the filter, since the state prediction error is positively correlated with the sampling period and might exceeds the linear correction range of the EKF.

B. Effects of Orbital Heights

The spherical shape of the Earth's gravitational field indicates that the accuracy of the GGT-based orbit determination is highly dependent of the height and the eccentricity of the orbit but has little relationship with the other orbital elements, such as inclination, true anomaly, etc. Besides the 300 km height case, simulations of circular orbits at heights of 600 km, 1000 km, 2000 km, and 5000 km have also been conducted. The steady-state variances of estimation errors are summarized in Table III. It is seen that the orbit determination accuracy decreases with increasing the orbital height. For the 5000 km case, the steady-state 1σ variance of the total 3D position error is 458 m, which is about 3 times that of the 300 km case (148 m). This phenomenon is due to the fact that the GGT signals and their sensitivity factors with respect to position variation are inversely related to the altitude. Thus the GGT-based orbit determination is more suitable for low-Earth orbiting (LEO) satellites. Another phenomenon is that the oscillation of steady-state variances becomes weaker at higher altitudes. The reason is that the effects of attitude errors become weaker at higher altitudes. As seen from Eq. (55), the GGT transformation errors are positively related to gravity gradients, which decreases with height. The filter consistency is not affected by the orbital height. All the NEES values of the four cases in Table III are below the upper limit of 12.6. The effects of dynamic variation of orbital heights are also investigated by assuming an elliptical orbit (eccentricity = 0.26) having a perigee height of 300 km and an apogee height of 5000 km. The evolution of the total 3D position errors and the corresponding 3σ boundaries are plotted in Fig. 6. The 6-hour simulation covers 2.5 orbital periods. The predicted variances are consistent with the actual errors. The varying position errors are mainly due to the varying height of the elliptical orbit. The maximum position variance (1σ) at the apogee height and the minimum position variance (1σ) at the perigee height are 448 m and Fig. 6. Evolution of the total 3D position errors and the 3σ boundaries for the elliptical orbit.

C. Effects of Noise Levels

The effects of attitude errors are first investigated. A small gradiometer noise of by the attitude errors. The reason is that the radial position is mainly related to GGT eigenvalues has little to do with the transformation error [START_REF] Chen | Gravity gradient tensor eigendecomposition for spacecraft positioning[END_REF]. It is noticed that the position accuracy is comparable to the GPS pseudorange performance when the attitude error is reduced to 0.1 arcsec.

To investigate the effects of gradiometer noise on orbit determination accuracy, different levels of gradiometer noise (

0 2 l  ),
1 mE, 0.01 E, 0.1 E, and 1 E, have been simulated. The 1 E and 0.1 E noise levels represent the precision of most current generation gradiometers, such as the ARKeX's Exploration Gravity Gradiometer and the Gedex's High-Definition Airborne Gravity

Gradiometer [START_REF] Difrancesco | Gravity gradiometer systems-advances and challenges[END_REF]. The 0.01 E noise level represents the precision of latest gradiometers, such as the GOCE's EGG. The 0.001 E noise level represents the precision of future grade gradiometers [START_REF] Carraz | A spaceborne gravity gradiometer concept based on cold atom interferometers for measuring Earth's gravity field[END_REF]. The attitude error is assumed to be 1 arcsec and the other simulation conditions stay the same with the baseline case. The steady-state variances of estimation errors are presented in Table V.

D. Effects of Measurement Biases

Simulated constant biases have been added to true GGT signals to examine the performance of the ASEKF. The biases assumed on the six gravity gradient components are 300 E (bxx), -2500 E (byy), 1500 E (bzz), 420 E (bxy), 900 E (bxz), and -120 E (byz), respectively. These significant biases are on the same order of magnitude as the true GGT signals. The initial errors of the bias estimation are set to 10 E. The standard deviation of the fictitious process noise is set to 1 mE. As mentioned earlier, the total simulation time for the ASEKF has been increased to 40 hours. The other simulation conditions are set to be the same with the baseline case.

The position and velocity estimation errors are depicted in Fig. 7 7 and Fig. 9, the along-track orbit position error is negatively related to the b xz error, indicating that the ASEKF cannot easily distinguish the variations of these two states. This is due to the fact that the along-track position error varies slowly during the orbit evolution and its dynamic behavior resembles a bias. To illustrate this, assume a Kepler circular orbital motion with an initial error in the true anomaly. The errors of the along-track position as well as the radial velocity both stay constant in the orbit prediction, leading to nearly constant errors of the gravity gradient component Vg,xz. In addition, the biasedness exists in the bxz estimation, leading to the biasedness in the estimation of along-track position and radial velocity.

VI. APPLICATION TO REAL GOCE DATA

Real flight data from the GOCE satellite have been used to test the orbit determination algorithm. The satellite was launched into a sun-synchronous orbit on March 17, 2009, with a mean apogee height of 278.65 km and an inclination of 96.7 deg [START_REF] Floberghagen | Mission design, operation and exploitation of the gravity field and steady-state ocean circulation explorer mission[END_REF]. The gradiometer onboard the satellite was originally designed to determine the medium-and short-wavelength parts of the gravity field.

To ensure a stable and quiet environment, a drag-free control mode was employed. The non-gravitational forces (drag and solar radiation pressure) were continuously compensated by an electric propulsion engine. In addition, three advanced star trackers and two dual-frequency GPS receivers were carried onboard to provide high-precision attitude and orbit information capabilities.

The gravity gradients were sampled at a rate of 1 Hz. The gradiometer achieved the best performance within the MBW from 5 mHz to 0.1 Hz. Data analysis has shown that the noise density level inside the MBW is approximately 20 mE/√Hz for Vg,xx, Vg,yy, V g,zz and V g,xz , 350 mE/√Hz for V g,xy , and 500 mE/√Hz for V g,yz . The larger noise levels of V g,xy and V g,yz are due to the less sensitive axes of the accelerometers. Below the MBW the noises increase with 1/f and are superimposed by cyclic distortions [START_REF] Rummel | GOCE gravitational gradiometry[END_REF]. Therefore, the biases of GOCE gradiometer are not constant. They are slowly drifting and have small-amplitude oscillations. To cope with the time-varying biases, a large process noise is assumed. The standard deviation is set to 20 mE.

A 40-hour data arc starting from September 8, 2013, 00:00:00.0 (GPS time) are used for the test. The data arc covers about 27 orbital periods and includes 4800 measurement epochs with a data-sampling period of 30 s. The GOCE Level 1b product EGG_NOM_1b which contains raw GGT measurements (EGG_GGT) and gradiometer inertial attitudes (EGG_IAQ) are used as inputs. The outliers in the raw GGTs have been detected and removed via polynomial fitting. The standard deviation of the attitude errors is set to 2 arcsec [START_REF] Stummer | Upgrade of the GOCE Level 1b gradiometer processor[END_REF]. To assess the accuracy of the orbit determination results, the GPS-derived high-precision orbits with an accuracy of 2 cm are used as references. The reference values of the GGT biases have also been obtained through comparing the actual measurements with GGTs calculated from the 300th degree and order EGM2008 model.

VII. CONCLUSIONS

A geophysical navigation scheme using the Earth's gravitational field is proposed in this paper for satellite autonomous orbit determination. The gravity gradiometer which measures full-tensor gravity gradients is surveyed as a primary sensor.

High-precision attitude sensors are assumed available to remove the orientation contributions in the observations. The novel method does not rely on ground stations or any other satellites and is suitable for fully autonomous navigation.

The Extended Kalman filter is implemented for orbit determination. Several simulations are conducted to test the algorithm. The initial errors only affect the convergence speed but not the steady-state errors. Increasing the sampling periods degrades the filter accuracy. In addition, the orbit determination accuracy is affected by the orbital heights as well as the attitude errors and gradiometer noises. With a future grade gradiometer, accuracy comparable to the GPS pseudorange will be possible. The bias on Additional calibrations could be helpful to improve the orbit accuracy. Nevertheless, this study demonstrates the feasibility of the EKF algorithm for GGT-based orbit determination, which could be used as a backup for near Earth autonomous spacecraft operation, or as a primary navigation system for future planetary exploration.

APPENDIX A. CALCULATION OF GGT IN ECEF

The computational form of the gravitation potential is given as follows

  0 0 n nm nm nm nm n m E GM U CV S W R       (35) with     1 sin cos n E nm nm R V P m r               1 sin sin n E nm nm R W P m r          
The recursive calculation of nm V and nm W is given by

    , 1 , 1 1 , 1 2 2 0, 1 2 1 2 E E m m m m m m m m xR yR V V W r r m                 (36)     , 1 , 1 1 , 1 2 2 0, 1 2 1 2 E E m m m m m m m m xR yR W W V r r m                 (37)               2 , 1 , 2 , 2 2 1 2 1 1 1 2 1 2 3 E E n m n m n m n n n m n m n zR R V V V n m n m n m n m n r r                       (38)               2 , 1 , 2 , 2 2 1 2 1 1 1 2 1 2 3 E E n m n m n m n n n m n m n zR R W W W n m n m n m n m n r r                       (39)
where x, y, and z here refer to the three position components in the ECEF frame. The initial values are

    0, 0 , 0, 0 0 E R V W r   (40)
The diagonal terms mm V and mm W can be first calculated using Eq. (36) and Eq. (37). By fixing m, the terms nm V and nm W (n > m) can be subsequently obtained using Eq. (38) and Eq. (39).

The gravitational force is obtained by differentiation of Eq. ( 35) and is given by 0 0

n nm nm e n m n m n m E e e V W GM C S R                 a r r (41)
where e a and e r refer to the gravitational force vector and position vector in the ECEF coordinate frame. The partial derivatives of nm V and nm W are

   1,1 1, 1 1, 1 0, 1 1 ( 2)( 1)(2 1) , 0 2(2 3) 2 2 1( 2 1 ) 1 ( 2)( 1 )(2 1 ) , 0 2 ( 2 3 ) ( 2 ) ( 2 3
)

n E nm n m n m E m n n n V m R n V x n m n m n n m n m n V V m R n n                                                (42)    1,1 1, 1 1, 1 0, 1 1 ( 2)( 1)(2 1) , 0 2(2 3) 2 2 1( 2 1 ) 1 ( 2)( 1 )(2 1 ) , 0 2 ( 2 3 ) ( 2 ) ( 2 3
)

n E nm n m n m E m n n n W m R n V y n m n m n n m n m n W W m R n n                                                (43)   1, 1 ( 1)(2 1) 1 ( 2 3
)

nm n m E n m n m n V V z R n            (44)    1, 1 1, 1 0, 1 0, 0 2 2 1( 2 1 ) 1 ( 2)( 1 )(2 1 ) , 0 2 ( 2 3 ) ( 2 ) ( 2 3 ) nm n m n m E m m W n m n m n n m n m n W W m x R n n                                          (45)    1, 1 1, 1 0, 1 0, 0 2 2 1( 2 1 ) 1 ( 2)( 1)(2 1) , 0 2 ( 2 3 ) ( 2 ) ( 2 3 ) nm n m n m E m m W n m n m n n m n m n V V m y R n n                                         (46)   1, 1 ( 1)
(2 1) 1 ( 2 3
)

nm n m E n m n m n W W z R n            (47)
The gravity gradient tensor can be further obtained by differentiation of Eq. ( 41) and is given by

2 2 2 2 0 0 n nm nm nm nm n m E e e V W GM C S R                 V r r ( 48
)
where V represents GGT and should be distinguished from nm V . The spatial derivatives of V are

3 3 0 0 n ij nm nm nm nm n m E V V W GM C S x R ijx ijx                       (49) 3 3 0 0 n ij nm nm nm nm n m E V V W GM C S y R ijy ijy                       (50) 3 3 0 0 n ij nm nm nm nm n m E V V W GM C S z R i j z i j z                       (51)
where i and j refer to the indexes x, y, and z. cos sin sin sin cos cos cos sin sin sin 0 cos cos cos sin 0 sin sin sin cos cos sin cos sin cos sin 0

g i                                               C (53)
g i                                                C L 2 3 1 2 3 1 2 3 2 2 1 2 1 2 2 2 3 1 2 3 1 2 3
cos sin sin cos sin cos sin sin sin sin sin cos cos cos cos sin cos cos cos sin cos

g i                                         C L 3 1 2 3 1 3 1 2 3 1 2 3 3 3 3 1 2 3 1 3 1 2 3 1 2 3
sin cos sin cos sin sin sin sin cos cos cos cos 0 0 0 cos cos sin sin sin cos sin sin sin cos cos sin He is currently a professor at the school of astronautics, Beihang University. He is also an editorial board member of Chinese Journal of Aeronautics. The focus of his research activities lies in the area of spacecraft orbit and attitude dynamics, spacecraft guidance, navigation, and control.

g i                                                C L

Fig. 1 .

 1 Fig. 1. The precision of EGM2008 gravity models truncated at different degrees. The benchmark is a 300th degree and order version. The reference height is 300 km.

  discrete process noise. The random sequence   k w is independent and zero-mean Gaussian. The function φ has no explicit expression and is only numerically obtained using an ordinary differential equation (ODE) solver.The first-order partial derivative of k x with respect to 1 k 

w

  In the orbit determination process, the matrix differential equations,

kR

  The measurement equation considering biases is further given in Eq. (27).The first-order partial derivative of  

  e T comprises the spatial derivatives of the gravity gradients. The observation error k ξ includes the GGT transformation error t k ξ and the gradiometer noise n k ξ . The GGT transformation



  is viewed as the equivalent gradiometer noise level.

  error and the gradiometer noise are 10 arcsec and 0.1 E, respectively. Initial errors of 10 km and 10 m/s are added to each component of the ECI position and velocity vectors. The basic EKF algorithm has been implemented for orbit determination. The position and velocity errors in the RSW coordinates are plotted in Fig.2and Fig.3, respectively. The solid lines are the actual errors and the dotted lines are the corresponding plus and minus 3σ boundaries. As seen from the plots, the filter is consistent: the estimation error has a zero mean and the predicted variances match the actual errors. The filter consistency can also be indicated TAES-201500387-R1 from the chi-square test results. The NEES values at each epoch are shown in Fig.4. The upper limit of the probability region is 12.6 for the basic EKF since 6.x n  It is seen that all of the points are inside the region.Orbit observability can be seen from the predicted variances of the estimation errors. The position errors decrease from 10 km to hundreds of meters level at the first epoch and rapidly (t < 30 min) converge to steady states (1σ variance) of 31.5-32.5 m (radial), 72.0-95.5 m (along-track), 104-119 m (cross-track), and 142-148 m (3D). The Root Mean Square (RMS) values of the steady-state actual position errors are 29.3 m (radial), 74.8 m (along-track), 89.2 m (cross-track), and 120 m (3D). The position error reduction at the first epoch is due to the position fixing ability of full-tensor GGTs. The oscillation of the steady-state variances is resulted from the combined effects of the satellite rotation motion and the attitude errors. The oscillation frequency is twice that of the rotation motion. In addition, the steady-state position error along the radial direction is much smaller than those along the other two directions. The reason is that the sensitivity factor for the vertical position is much smaller [16]. The velocity errors converge more rapidly (t < 20 min) to steady states (1σ variance) of 0.172-0.178 m/s (radial), 0.206-0.220 m/s (along-track), 0.241-0.246 m/s (cross-track), and 0.364-0.373 m/s (3D). The RMS values of the steady-state actual velocity errors are 0.111 m/s (radial), 0.099 m/s (along-track), 0.120 m/s (cross-track), and 0.192 m/s (3D). The baseline case demonstrates the feasibility of the EKF for autonomous orbit determination from unbiased gravity gradient measurements. The position accuracy is competitive with that of traditional ground-based station tracking considering the fast convergence of the filter.

Fig. 2 .

 2 Fig.2. Radial, along-track, and cross-track position errors and the corresponding 3σ boundaries for the baseline case.

Fig. 5 .

 5 Fig. 5. Evolution of the total 3D position errors for the EKF with different initial estimation errors.

  The accuracies of position and velocity decrease with increase of gradiometer noise. The RMS values of the actual 3D position errors for these four cases are 20.8 m, 22.1 m, 76.5 m, and 475 m, respectively, and the RMS values of the actual 3D velocity errors are 0.127 m/s, 0.131 m/s, 0.179 m/s, and 0.614 m/s.

  and Fig. 8, respectively. The solid lines are the actual errors and the dotted lines are the corresponding plus and minus 3σ boundaries. Chi-square test has been conducted to evaluate the filter consistency. The NEES values at each epoch are shown in Fig. 9. The upper limit of the probability region is 21. All the points are inside the region. Out of the 4800 points, only 15 are found outside the region. The radial and cross-track position errors rapidly (t < 3 hours) converge to steady states (1σ variance) of 20.0-20.5 m and 30.0-32.0 m, respectively. The along-track position error converges quite slowly and the 1σ variance at the last epoch is about 400 m. The RMS values of the steady-state actual position errors are 21.5 m (radial), 271 m (along-track), 23.0 m (cross-track), and 273 m (3D). The along-track and cross-track velocity errors rapidly (t < 2 hours) converge to steady states (1σ variance) of 0.0198-0.0202 m/s and 0.0355-0.0375 m/s, respectively. The radial velocity error converges quite slowly and 1σ variance at the last epoch is about 0.47 m/s. The RMS values of the steady-state actual velocity errors are 0.291 m/s (radial), 0.0159 m/s (along-track), 0.0257 m/s (cross-track), and 0.293 m/s (3D). The steady-state variances of the radial and cross-track position components are smaller than those in the baseline case due to the more accurate orbital dynamic model. In addition, the error curves seem much smoother since the orbit prediction constrains the jumps in state estimates.

Fig. 7 .

 7 Fig. 7. Radial, along-track, and cross-track position errors and the corresponding 3σ boundaries for the ASEKF case.

Fig. 8 .Fig. 9 .

 89 Fig.8. Radial, along-track, and cross-track velocity errors and the corresponding 3σ boundaries for the ASEKF case.

  the bias estimation errors and the corresponding 3σ boundaries is plotted in Fig. 10. Among the six biases, the b xx , b yy , b zz , b xy , and b yz components rapidly (t < 3 hours) converge to steady states (1σ variance) of less than 11 mE, 11 mE, 13 mE, 10 mE, and 16 mE, respectively. However, the bxz component converges quite slowly. The 1σ variance at the last epoch is about 240 mE. The RMS values of the steady-state actual bias errors are 7.52 mE (b xx ), 8.74 mE (b yy ), 7.26 mE (b zz ), 6.73 mE (b xy ), 158 mE (bxz), and 11.1 mE (byz). The poor observability of bxz accounts for the large errors in the along-track position estimation. As seen from Fig.

Fig. 10 .

 10 Fig.10. Evolution of the bias estimation errors and the corresponding 3σ boundaries for the ASEKF case.

  The ASEKF has been implemented for orbit determination. The initial position and velocity errors are 10 km and 10 m/s. The initial bias errors are 10 E. The position and velocity errors are shown in Fig.11and Fig. 12, respectively. The chi-square test has been conducted to evaluate the filter consistency. The NEES values are shown in Fig. 13. Only 22 points are found outside the probability region. Similar to the ASEKF simulation case, the radial and cross-track position errors converge rapidly (t < 4 hours) to steady states (1σ variance) of 60-66 m and 55-75 m, respectively. The along-track position error converges slowly and 1σ variance at the last epoch is about 900 m. The RMS values of the steady-state actual position errors are 15.6 m (radial), 559 m (along-track), 52.2 m (cross-track), and 562 m (3D). The along-track and cross-track velocity errors converge rapidly (t < 5 hours) to steady states (1σ variance) of 0.050-0.052 m/s and 0.065-0.080 m/s, respectively, whereas the radial velocity error converges slowly with a final 1σ variance value of about 1.05 m/s. The RMS values of the steady-state actual velocity errors are 0.658 m/s (radial), 0.0166 m/s (along-track), 0.0587 m/s (cross-track), and 0.662 m/s (3D). The GGT biases have also been estimated along with the orbits and the errors are depicted in Fig. 14. The bxx, byy, bzz, bxy, and byz components rapidly (t < 3 hours) converge to steady states (1σ variance) of less than 43 mE, 42 mE, 80 mE, 83 mE, and 105 mE, respectively. The b xz component converges very slowly with a final 1σ variance value of about 600 mE. The RMS values of the steady-state actual bias errors are 12.3 mE (bxx), 30.7 mE (b yy ), 21.4 mE (b zz ), 85.5 mE (b xy ), 351 mE (b xz ), and 101 mE (b yz ). The use of bias process noise allows successful estimation of the varying biases. The reference values as well as the estimated values of byy are plotted in Fig. 15 as an example. It is seen that the estimates follow closely the references values and compensate the drift as well as the small periodic fluctuations. Although not shown in a figure, the post-fit measurement residuals have been calculated. The statistical standard deviations of the six components are 12 mE, 11 mE, 16 mE, 335 mE, 15 mE, and 494 mE, respectively, consistent with the flat noise levels in the MBW.

Fig. 11 .

 11 Fig. 11. Radial, along-track, and cross-track position errors and the corresponding 3σ boundaries for the GOCE case.

Fig. 12 .

 12 Fig.[START_REF] Richeson | Gravity gradiometer aided inertial navigation within non-GNSS environments[END_REF]. Radial, along-track, and cross-track velocity errors and the corresponding 3σ boundaries for the GOCE case.

Fig. 13 .

 13 Fig.[START_REF] Pavlis | The development and evaluation of the earth gravitational model 2008 (EGM2008)[END_REF]. Normalized state estimation error squared (NEES) as well as the 95% probability region for the GOCE case.

Fig. 14 .

 14 Fig. 14. Evolution of the bias estimation errors and the corresponding 3σ boundaries for the GOCE case.

Fig. 15 .

 15 Fig. 15. The estimated and reference values of the byy bias of GOCE.

  track-radial gravity gradient components worsens the system's observability and leads to large estimation errors of the along-track position and the radial velocity. The along-track position error after 40 hours is about 350m for the real GOCE case.

and 3 

 3 are the yaw, pitch, and roll angles, respectively. The relationship between

  the vector-valued GGT transformation error can be written as matrix k M can be directly obtained from Eq. (55).Xiucong Sun was born inLaiwu, China, in 1988. He received his B.S. degree in aerospace engineering from Beihang University, Beijing, in 2010. He is a Ph.D. student at the school of astronautics, Beihang University. Since 2014, he has been with ENAC (Ecole Nationale de l'Aviation Civile) TELECOM laboratory as a visiting researcher, funded by the China Scholarship Council. The focus of his current research mainly lies in gravity gradient based spacecraft navigation. Pei Chen was born in Beijing, China, in 1979. He received his Ph.D. degree in aerospace engineering from Beihang University, Beijing, in 2008. He is currently an associate professor at the school of astronautics, Beihang University. His current research activities comprise spacecraft navigation, GNSS application, and astrodynamics and simulation. . Christophe Macabiau graduated as an electronics engineer in 1992 from the ENAC (Ecole Nationale de l'Aviation Civile) in Toulouse, France. He received his Ph.D. in 1997. Since 1994 he has been working on the application of satellite navigation techniques to civil aviation and has been in charge of the TELECOM Lab of the ENAC since 2000. Chao Han was born in Beijing, China, in 1960. He received his M.S. and Ph.D. degrees in Applied Mechanics from Beihang University in 1985 and 1989, respectively.

Fig. 1 .

 1 Fig.1. The precision of EGM2008 gravity models truncated at different degrees. The benchmark is a 300th degree and order version. The reference height is 300 km.

Fig. 2 .

 2 Fig. 2. Radial, along-track, and cross-track position errors and the corresponding 3σ boundaries for the baseline case.

Fig. 3 .

 3 Fig.3. Radial, along-track, and cross-track velocity errors and the corresponding 3σ boundaries for the baseline case.

Fig. 4 .

 4 Fig. 4. Normalized state estimation error squared (NEES) and the 95% probability region for the baseline case.

Fig. 5 .

 5 Fig. 5. Evolution of the total 3D position errors for the EKF with different initial estimation errors.

Fig. 6 .

 6 Fig.6. Evolution of the total 3D position errors and the 3σ boundaries for the elliptical orbit.

Fig. 7 .

 7 Fig. 7. Radial, along-track, and cross-track position errors and the corresponding 3σ boundaries for the ASEKF case.

Fig. 8 .

 8 Fig.8. Radial, along-track, and cross-track velocity errors and the corresponding 3σ boundaries for the ASEKF case.

Fig. 9 .

 9 Fig.9. Normalized state estimation error squared (NEES) as well as the 95% probability region for the ASEKF case.

Fig. 10 .

 10 Fig.10. Evolution of the bias estimation errors and the corresponding 3σ boundaries for the ASEKF case.

Fig. 11 .

 11 Fig. 11. Radial, along-track, and cross-track position errors and the corresponding 3σ boundaries for the GOCE case.

Fig. 12 .

 12 Fig.[START_REF] Richeson | Gravity gradiometer aided inertial navigation within non-GNSS environments[END_REF]. Radial, along-track, and cross-track velocity errors and the corresponding 3σ boundaries for the GOCE case.

Fig. 13 .

 13 Fig.[START_REF] Pavlis | The development and evaluation of the earth gravitational model 2008 (EGM2008)[END_REF]. Normalized state estimation error squared (NEES) as well as the 95% probability region for the GOCE case.

Fig. 14 .

 14 Fig. 14. Evolution of the bias estimation errors and the corresponding 3σ boundaries for the GOCE case.

Fig. 15 .

 15 Fig. 15. The estimated and reference values of the b yy bias of GOCE.

TABLE I INITIAL

 I ORBITAL ELEMENTS OF THE REFERENCE ORBIT IN THE BASELINE CASE

	Orbit Element	Initial Conditions
	height	300 km
	eccentricity	0
	inclination	60 deg
	right ascension of the ascending node	120 deg
	argument of perigee	0 deg
	true anomaly	80 deg

TABLE II STEADY

 II -STATE VARIANCES FOR THE EKF WITH DIFFERENT DATA-SAMPLING PERIODS

	Sampling		Position Variances (1σ, m)			Velocity Variances (1σ, m/s)
	Period (s)	Radial	Along-track	Cross-track	3D	Radial	Along-track	Cross-track	3D
	10	21.4-21.8	50.0-67.5	71.5-82.5	98.5-103	0.148-0.152 0.185-0.202 0.213-0.220 0.324-0.332
	60	40.0-41.0	90.0-117	130-147	177-184	0.190-0.197 0.218-0.230 0.260-0.265 0.392-0.400
	100	46.8-47.8	106-135	152-170	206-213	0.203-0.210 0.226-0.238 0.275-0.281 0.412-0.422
	300	60.4-61.5	146-179	205-227	276-285	0.233-0.243 0.240-0.251 0.316-0.327 0.464-0.476
	1000	67.2-68.5	188-251	265-295	360-380	0.272-0.295 0.263-0.275 0.390-0.410 0.550-0.572

TABLE IV STEADY

 IV -STATE VARIANCES FOR THE EKF WITH DIFFERENT LEVELS OF ATTITUDE NOISE

	Attitude Noise		Position Variances (1σ, m)			Velocity Variances (1σ, m/s)
	(arcsec)	Radial	Along-track	Cross-track	3D	Radial	Along-track	Cross-track	3D
	0.1	0.640-0.649	1.91-2.86	2.86-3.46	4.00-4.24	0.042-0.043 0.061-0.071 0.071-0.076 0.106-0.110
	1	0.640-0.650	9.25-16.5	17.0-20.8	22.8-24.5	0.043-0.046 0.105-0.130 0.132-0.140 0.180-0.191
	10	0.650-0.700	51.5-80.0	96.0-112	123-129	0.070-0.097 0.178-0.195 0.235-0.240 0.305-0.321
	30	0.650-0.850	105-141	195-220	238-251	0.125-0.163 0.200-0.208 0.301-0.315 0.385-0.410

TABLE V STEADY

 V -STATE VARIANCES FOR THE EKF WITH DIFFERENT LEVELS OF GRADIOMETER NOISE

	Gradiometer		Position Variances (1σ, m)			Velocity Variances (1σ, m/s)
	Noise (E)	Radial	Along-track	Cross-track	3D	Radial	Along-track	Cross-track	3D
	0.001	0.640-0.650	9.25-16.5	17.0-20.8	22.8-24.5	0.043-0.046 0.105-0.130 0.132-0.140 0.180-0.191
	0.01	4.83-4.90	12.8-18.4	18.5-22.0	25.9-27.3	0.086-0.087 0.120-0.135 0.136-0.143 0.205-0.211
	0.1	31.3-31.7	57.4-58.5	59.6-60.1	88.7-89.6	0.168-0.169 0.194-0.195 0.199-0.200 0.325-0.326
	1	206-209	292-296	270-272	448-452	0.451-0.453 0.294-0.295 0.357-0.358 0.646-0.648

  The second-and third-order partial derivatives of nm V and nm W can be derived from

	Eqs. (42)-(47) and are not given here.
						APPENDIX B. DERIVATION OF (15)
	Let	g i C denote the rotation uncertainty of	. C The resulting GGT transformation error in the matrix form will be g i
				g    V	g i C C VC C C C VC C i e i g i e e e i g i e e i  	i g	(52)
	The transformation matrix g i C can be expressed as a function of three Euler angles
		3 cos cos	1	2 sin sin sin 3	1	3 cos sin	1	2 sin sin cos 3	1	2 cos sin	3
			2 cos sin	1				2 cos cos	1	sin	2
		3 sin cos	1	2 sin cos sin 3	1	3 sin sin	1	2 sin cos cos 3	1	2 cos cos	3

  L are the coefficient matrices. Therefore, the GGT transformation error in the matrix form becomes

	where 1 , L	2 , L and 3
		g		1	i e	e i e i g	g i i e	e T 11  e i
			 	2 3	i e i e	e i e i g e i e i g	g i i e g i i e	22 e T e i  e T  33 e i

This research was supported by the National Nature Science Foundation of China through cooperative agreement No. 11002008 and has been funded in part by Founded by MOST through cooperative agreement No. 2014CB845303. The European Space Agency is acknowledged for providing the GOCE data.


DSNA/ Nicolas Durand

MODELING THE CONTROLLER'S CONFLICT DETECTION TASK USING FAST TIME SIMULATION

This paper shows how the CATS (Complete Air Traffic Simulator) can model the controller's conflict detection task using a sliding time window. Because of uncertainties, detection is by far the most time consuming and cognitively challenging aspect of air traffic control. The processes governing detection are distinct from resolution actions. Few conflicts detected lead to a resolution maneuver. Because future aircraft positions are uncertain, controllers detect many more conflicts or "potential conflicts" than actually occur. In this paper, we solve conflicts in a real time context and compare different uncertainty scenarios. We show how uncertainties impact the number of theoretical potential conflicts detected and the increase the number of unnecessary maneuver actions required to keep the traffic safe. When we deal with realistic uncertainty values, the majority of calculation time in CATS is used for calculating maneuvers that will usually not be transmitted to pilots. This observation reveals the true nature of the controller's workload today.

Introduction

Much current research models and estimates the controller's workload. Some studies analyse questionnaires filled by the controllers themselves in different situations. Others use talking time measurements between controllers and pilots. In experiments, controllers have been fitted out with sensors to measure their heartbeat or eye movements [START_REF] Averty | Evaluating a new index of mental workload in real control situation using psychophysiological measures[END_REF], [START_REF] Collet | Workload in air traffic controllers estimated by autonomic nervous system activation and subjects self estimation[END_REF]. Existing models often divide the controller's workload into three specific tasks. The coordination task is clearly identified as the action of taking into account new aircraft, establishing contact with the pilot, and delivering the aircraft to the next sector with respect to existing rules. The monitoring task consists in making sure that the aircraft respects the trajectory defined by the flight plan. Modeling the monitoring task is not easy because an objective of trajectory monitoring is conflict detection. The detection-resolution task is the only task in which the controller can have a concrete effect on aircraft trajectory. Most people include conflict detection as part of the controller's conflict resolution task. However, this paper argues that conflict detection must be conceptualized as a specific task. Because of uncertainties, detection is by far the most time consuming and cognitively challenging aspect of air traffic control. The processes governing detection are distinct from resolution actions. Few conflicts detected lead to a resolution maneuver. Because aircraft future positions are uncertain, controllers detect many more conflicts or "potential conflicts" than actually occur. They generally try to avoid unnecessary maneuvers but they sometimes solve "fake" conflicts for two reasons: 1) usually maneuvers must be initiated far in advance of the conflicting zone. By the time a maneuver is decided, the controller cannot estimate precisely the aircraft's future positions; 2) because controllers are dealing with many aircraft at the same time, they often need to decide maneuvers that will definitely solve potential conflicts in order to give them more time to concentrate on other parts of the traffic. The controllers' priority is security rather than efficiency, they must use learned methods and experience and usually do not show initiative when solving a conflict.

In this paper, we show how the CATS fast time simulator models both the conflict detection and resolution task. By using a sliding window over a day of traffic, we can simulate the controller's behaviour in a real time context and compare different uncertainty scenarios. We show how uncertainties impact the number of theoretical potential conflicts detected and the increase of unnecessary maneuver actions required to keep the traffic safe. We used the French upper airspace to compare different scenar-ios on different days of traffic. The simulations show many situations where maneuvers are anticipated and finally not executed because trajectory prediction is updated and potential conflicts disappear as the conflict zone approaches. When we deal with realistic uncertainty values, the most calculation time in CATS is used for calculating maneuvers that will usually not be transmitted to pilots. This observation reveals the true nature of controller's workload today. Many complexity studies have shown that the number of aircraft in a sector explains most of the controller's workload [START_REF] Gianazza | Selection and evaluation of air traffic complexity metrics[END_REF]. The monitoring task is thus very greedy compared to the resolution task itself because uncertainties create many potential conflicts and increase the monitoring and detection work of controllers. Today's new communication, navigation and positioning systems can help predict aircraft future positions accurately. New tools are tested and installed in countries around the world to assist controllers who are still responsible for conflict detection. In order to be effective and gain the trust of controllers, these tools should never fail to detect a conflict and must be better than or as good as controllers. The challenging question is not how close to trajectory prediction should aircraft stick in order to build an efficient detection resolution tool, but how accurate can trajectory prediction be to build tools that could help controllers in their detection task. The different stakeholders of the system sometimes have the feeling that Air Navigation System Providers want to constrain the trajectories. This paper shows the importance of better sharing uncertainty information.

In the first section we describe the general framework of the CATS simulator. The second part details the ERCOS solver that is used to build conflict free trajectories. The third part focuses on the solver algorithm. Part four details the experimental results performed with different uncertainty scenarios. We show how the controller's detection task can be modeled and how future TP improvements could decrease the controller's workload in the future.

I. The Complete Air Traffic Simulator (CATS)

A. General framework

The CAT Simulator was used in many projected conducted by the French DSNA to study either centralized control models [START_REF] Durand | A statistical analysis of the influence of vertical and ground speed errors on conflict probe[END_REF], [START_REF] Granger | A traffic complexity approach through cluster analysis[END_REF], [START_REF] Granger | Optimal resolution of en route conflicts[END_REF] or Free-Flight concepts [START_REF]FACES: a Free flight Autonomous and Coordinated Embarked Solver[END_REF], [START_REF]Token allocation strategy for free-flight conflict solving[END_REF], [START_REF] Durand | Faces: a free flight autonomous and coordinated embarked solver[END_REF].

In this paper, we used real data traffic simulation. The CAT Simulator takes as input flight plans given by airlines or pilots before or after regulation. The simulator uses the BADA (Base of Aircraft DAta) tabulated model for modeling aircraft performances. For an aircraft type, it gives a vertical speed and a ground speed depending on whether the aircraft is climbing, leveled or descending. The BADA performance summary tables are derived from the total energy model of EUROCONTROL.

Aircraft speeds can be modified by a random value to take into account different factors of uncertainty (aircraft load, winds, etc. . . ). These values can be either computed once at aircraft activation and remain the same for all the flight, or can be modified anytime during the flight. Uncertainty modeling for conflict detection and resolution is introduced in the next section. We will see how using uncertainty on trajectory prediction can model the controller's conflict detection-resolution task.

Aircraft follow classical routes (from way-point to way-point). The flight model is simple. An aircraft first climbs up to its RFL (Requested Flight Level), then remains leveled until its top of descent, then descends to its destination.

Aircraft flight is calculated per time step. The time step is always chosen in order to guarantee that the simumlator cannot miss any conflict. For all our simulations, we use a 15s time step.

B. General architecture of the system

We sketch here the general architecture of the simulator. Each part will be detailed in the following sections. The system architecture is presented in figure 1 and 2. The system relies on three main processes P1, P2, and P3:

• P1 is the traffic simulator. P1 sends current aircraft positions and flight plans to process P2. Process P2 builds trajectory forecasts for T w minutes, does conflict detection by pairs and transforms 1-to-1 conflicts in n-aircraft conflicts. Then, process P3 (the problem solver) solves in parallel each cluster, as aircraft in each cluster are independent from aircraft in the other clusters. The problem solver sends to P2 new orders and P2 builds new trajectory forecasts based on these orders. Then P2 once again runs a conflict detection process to check that modified aircraft trajectories do not interfere with aircraft in another cluster, or with new aircraft. If no interference is found, new flight orders can be sent to P1. If there are interferences, interfering clusters are joined and the solver is used again on that (these) cluster(s). The process is iterated until no interference between clusters remains, or no new aircraft is concerned by modified trajectories. The new orders are sent back to the traffic simulator.

The above process is iterated and all trajectories are optimized each δ minutes (3 or 5 minutes in the experiments). However, during the computation time, aircraft are flying and need to know if they must change their route or not. δ should be large enough to compute a solution, send it to the pilot and leave him enough time to begin the maneuver. Consequently, for each aircraft, at the beginning of the current optimization, trajectories are determined by the previous run of the problem solver and cannot be changed for the next δ minutes.

C. Conflict detection and clustering

1) Trajectory forecast and 1-to-1 conflict detection: As described above, the P2 process does trajectory prediction for T w minutes. This trajectory prediction is done again by a simulation on a slightly modified version of the Air Traffic simulator. But, as stated above, we assume that there is an error about the aircraft's future location because of ground speed prediction uncertainties 1 . Climbing and descending rate uncertainties are larger than ground speed uncertainties. Because the conflict free trajectory must be robust regarding these and many other uncertainties, an aircraft is represented by a point at the initial time. The point becomes a line segment in the direction of uncertainty (the speed direction here, see figure 3). The first point of the line "flies" at the maximum possible speed, and the last point at the minimum possible speed.

When changing direction on a beacon, the heading of the line segment's "fastest point" changes as described on figure 3.

To check the standard separation at time t, we compute the distance between the two line segments modeling the aircraft positions and compare it to the standard separation at each time step of the simulation.

In the vertical plane, we use a cylindrical modeling (figure 3). Each aircraft has a mean altitude, a maximal altitude and a minimal altitude. To check if two aircraft are in conflict, the minimal altitude of the higher aircraft is compared to the maximal altitude of the lower aircraft. 2) Clustering: After pair detection, P2 clusters conflicting aircraft. Each equivalence class for the relation "is in conflict with" becomes a cluster.

For example, if aircraft A and B are in conflict in the T w window, and if B is also in conflict with C in the same time window, then A, B, C is the same cluster and will be solved globally by the conflict solver.

The conflict solver sends back to P2 maneuver orders for solving conflicts. Then P2 computes new trajectories for all aircraft and checks if new interferences appear. For example, if the new trajectory given to aircraft B to solve the conflict with A and C interferes with cluster D, E and with aircraft F , then A, B, C, D, E, F will be sent back to the problem solver as one conflict to solve.

The process will always converge. In the worst case, P3 will have to solve a very large cluster including all aircraft present in the next T w minutes. However, this technique is usually efficient because a very large number of clusters can be solved very quickly in parallel.

Human controller's do not use any timestep, time window T w or update period δ to deal with the traffic. We will see in the experimental results that a lot of calculation is done by the simulator but that only a few maneuvers calculated become effective. This is due to uncertainties that the simulator has to take into account. In real life, controllers also have to deal with uncertainties. Counting the number of maneuvers that never become effective can give

II. Solver Modeling

A. Maneuver modeling

In the horizontal plane, classical maneuvers given to aircraft are heading deviation. In the simulator, 10, 20 or 30 degrees deviations will be allowed. The deviation starts on a virtual beacon created on the route (see figure 4). This beacon is defined by the position of the head of the segment at some time t 0 . It ends on a second virtual beacon, position of the head of the segment at time t 1 . An angle criteria is defined to find on which beacon the modified and initial routes should connect.

A maneuver will be determined by:

• t 0 which defines the first virtual beacon B 0 .

• the deviation angle α.

• t 1 which defines the second virtual beacon B 1 .

In the vertical plane, the aircraft trajectory is divided in 4 periods (figure 5):

• Climbing period. In this period, aircraft can be leveled at a lower than requested flight level to solve a conflict. The aircraft climb is stopped at flight level F L 0 and starts again on a virtual beacon B 1 as stated on figure 6. F L 0 and B 1 are defined by the position of the head of the uncertainty segment at time t 0 and t 1 . • Cruising period. When aircraft have reached their desired flight level, they may be moved

0 FL B 0 FL 1 B 0 B 1 B 1 B 1 l l
Climbing Cruising End of Cruising Descending period period period to the nearest lower level to resolve a conflict. Aircraft starts descending when reaching a virtual beacon B 0 and starts climbing at B 1 (α = 0, B 0 and B 1 are defined by the position of the head of the uncertainty segment at time t 0 and t 1 ). An example of maneuver is represented on figure 7. • End of Cruising period. When aircraft are about 50 nautical miles from beginning their descent to destination, they may be moved to a lower level to resolve a conflict. Aircraft start descending on B 0 and are leveled at F L 1 (α = 0) (see figure 8). B 0 and F L 1 are defined by the position of the head of the uncertainty segment at time t 0 and t 1 . • Descending period. During this period no vertical maneuver is possible. No maneuver will be simultaneously done in the horizontal and vertical plane. This model has the great advantage of reducing the size of the problem.

For a conflict involving n aircraft, the dimension of the search space is 3 n. This will allow us to solve

B. Maneuver decision time

Because of uncertainties, some conflicts that are detected too early would not actually occur in the end. Consequently, deciding to move aircraft in such cases would be useless and could even generate other conflicts that would not occur if no maneuver had been decided. This explains why controllers do not solve conflicts too early. When there is no uncertainty, the earlier the maneuver is started, the lower the delay. However, if speed is not strictly maintained, the earlier the conflict is detected, the lower the probability it will actually happen. Thus, a compromise must be reached between the delay generated and the risk of conflict.

Because of uncertainties, maneuvers should be started as late as possible while still respecting aircraft constraints. First, it prevents the system from deciding unnecessary maneuvers. Second, an aircraft that is already maneuvering cannot make another maneuver before resuming its initial speed.

The solver was modeled this way to keep the maneuvers simple to understand and execute. Starting maneuvers as late as possible increases the number of maneuverable aircraft. Duration T w can be changed, but must be at least equal to 2 × δ. A good evaluation of T w is difficult. With a perfect trajectory prediction, the largest T w should be chosen. However, this is not true as soon as uncertainties are included in the model. A large value of T w induces a large number of 1-to-1 conflicts, as the size of segments (modeling aircraft positions) grows quickly with time. Therefore, the conflict solver can become saturated.

When controllers are not too busy, they try to optimize the traffic and wait before giving a maneuver order that is not necessary. They sometimes prefer to solve a conflict ahead of time before it occurs in order to be able to take care of the rest of the traffic. It is very hard to reproduce this behaviour with an automatic solver, but we will show in this paper how the solver calculates many possible maneuvers that do not become effective but correspond to some control taskload.

C. A sliding forecast time window model

In order to limit the size of the problem and to be reactive to uncertainties, only the next T w minutes of the flights are considered. T w represents the lookahead time also called forecast time window. The situtation is revised every δ minutes with δ << T w . δ is the time step used in the model to make the T w time window slide. This approach ensures that the problem can be updated every δ minutes : current aircraft positions are updated which reduces uncertainties.

In figure 9, at t = 0, the aircraft trajectory cannot be modified before t = δ because any maneuver requires advance notice. Any maneuver that would occur between t = δ and t = 2δ would be kept as a constraint for the next optimization run (in the example, no maneuver is decided) because it will then be to close to the current time. In figure 9, the maneuver described on the first line resulting from an optimization at t = 0 (left turn of 30 degrees) is revised at time t = δ (left turn of 20 degrees) and then kept at time t = 2δ. The end of the maneuver can be recalculated until the optimization starting at t = 3δ.

Pilots should only be given maneuver orders that will not be modified; if no conflict occurs, no order will be given. In the example, the pilot will be notified of the beginning of the speed change at time δ and the end at time 3δ.

The size of the forecast time window is an important parameter. If it is too big, the size of the problem will include a very large number of variables and the resolution might be more difficult. If it is too small, the solutions found might be worse and the total delay induced over the day much higher. In order to solve conflicts with small speed adjustments, T w needs to be large enough. Figure 10 gives an example of two aircraft flying at the same level. In this example T w = 20 minutes and δ = 4 minutes. The speed uncertainty used is 5%. At time t = 4 minutes a potential conflict is detected and the solver calculates a maneuver (heading change) for one of the aircraft. Because the maneuver starts late enough, it is not sent to the aircraft. At time t = 8 minutes, the current positions of the aircraft are updated and the conflict is still detected. A new maneuver is calculated (heading change for the other aircraft), but still not sent to the aircraft because it starts late enough. At time t = 12 minutes, ther is a new update and the solver calculates a new maneuver (shorter than the previous one) but still late enough to remain an option. And at t = 16 minutes (see figure 10), after updating the current position of the aircraft the conflict has disappeared. This behaviour mimics the controllers' conflict detection task. When they are not too busy, controllers try to avoid unnecessary manuevers but in dense traffic situations, they might however give a maneuver that solves immediately a potential conflict to concentrate on the rest of the traffic.

III. Solver Algorithm

Classical Evolutionary Computation (CEC) principles such as described in the literature [START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Machine Learning[END_REF], [START_REF] Michalewicz | Genetic algorithms + Data Structures = Evolution Programs[END_REF] is used in the solver.

A. Fitness function

The cost function used in this part is simply the sum of the delays over the aircraft population.

Solutions respecting the separation constraints cannot be built easily. Consequently, we need to include the separation constraint in the fitness function.

The fitness function chosen is:

F = 2 n -nb man -n i=1 δi Tw
1 + n rc where n is the number of aircraft, nb man the number of maneuvers, δ i the delay resulting from the maneuver of aircraft i and n rc is the number of remaining conflicts.

The fitness function increases when the number of remaining conflicts and the delays decrease. It takes its values in [0, 2 n].

B. Crossover operator

The conflict resolution problem is partially separable as defined in [START_REF] Durand | Genetic crossover operator for partially separable functions[END_REF], [START_REF] Durand | Automatic aircraft conflict resolution using genetic algorithms[END_REF]. In order to increase the probability of producing children with a better fitness than their parents, principles applied in [START_REF] Durand | Genetic crossover operator for partially separable functions[END_REF] were used. For each aircraft i of a population element, a local fitness F i value is defined as follows:

F i = 2 -m i -( δi Tw ) 1 + nrc i
where nrc i is the number of remaining conflicts involving aircraft i and m i = 0 if aircraft i has no maneuver and 1 if it has a maneuver.

Figure 12 presents the crossover operator. First two population elements are randomly chosen. For each parent A and B, fitness A i and B i of aircraft i are compared. If A i < B i , the children will take aircraft i of parent A. If B i < A i , the children will take aircraft i of parent B. If A i = B i children randomly choose aircraft A i or B i or even a combination of A i and B i .

C. Mutation operator

For each candidate to mutation, the delay of an aircraft having one of the worst local fitnesses is modified. If every conflict is solved, an aircraft is randomly chosen and its parameters changed. In practice, a number m is randomly chosen in the interval [1, n 2 ] and we pick up m times an aircraft to find the most constrained aircraft among these m trials. The delay of this aircraft is then either locally optimized or randomly modified with a probability of 50%. We may be tempted to always locally optimize the delay of the worst aircraft, but this would make the algorithm become very deterministic and lead to a premature convergence of the algorithm.

The crossover and mutation operators are more deterministic during the first generations because there are many conflicts to solve. They focus on making feasible solutions. When the solutions without conflicts appear in the population, they become less deterministic.

Sharing: The problem is highly combinatorial and may have many local optima. In order to prevent the algorithm from premature convergence, the sharing process introduced by Yin and Germay [START_REF] Yin | A fast genetic algorithm with sharing scheme using cluster analysis methods in multimodal function optimization[END_REF] is used. The complexity of this sharing process has the great advantage to be in n log(n) (instead of n 2 for classical sharing) if n is the size of the population. The distance used to compare two population elements p and q is:

D = n i=1 |δ p i -δ q i | n

D. Parameters

In the experiments, the following parameters were empirically chosen: the size of the population was set to 100, 20% of the population is crossed, 60% is muted, the selection uses the stochastic remainder without replacement. A sharing process is used. As time to solve a problem is limited, the number of generations is limited to 500.

IV. Experimental Results

CATS was used on a busy day of traffic (July 17 th 2010) in the French airspace. We were only interested in the upper airspace (above FL195). The number of flights in this airspace is 8870 for that specific day. Without any control maneuver, 2305 conflicts are detected (26% of the aircraft population). The mean time of flights is 57 minutes and the mean travelled distance in the french airspace is 394 nautical miles. Different experiments were done using CATS with different hypotheses.

Tables I give the results of simulations for T w = 9 minutes and δ =2, 3 or 4 minutes for four sets of horizontal and vertical uncertainties (ranging from no uncertainty to 10% horizontally and 30% vertically). Column 5 gives the number of remaining conflicts. Column 6 gives the number of maneuvered aircraft and colum 7 the number of maneuvers sent to the pilots. Column 8 gives the number of maneuvers that were calculated by the solver but not sent to the pilots because they were supposed to start late enough to be revised at the next detection resolution process. Column 9 gives the number of potential conflicts detected in each scenario and column 10 gives a normalized number of potential conflicts by deviding the number of potential conflicts by Tw δ . This table shows that the number of potential conflicts detected and the number of maneuvers not sent to pilots grow with uncertainty and also with the update frequency. For example, with 2% and 5% uncertainties, the normalized number of potential conflicts is multiplied by 1.6, with 5% and 10% uncertainties by 2.5 and with 10% and 30% uncertainties by 5.

When T w = 9 minutes and δ = 4 minutes the solver cannot solve every conflict which means that T w needs to be big enough compared to δ to make the solver become efficient. The number of remaining conflicts increases with uncertainty. The fact that the number of maneuvers sent to pilots increases with δ reflects that updating the situation is very important. This explains why the detection resolution process is a permanent task for controllers. The high number of maneuvers that are never sent to pilots can also give an idea of the complexity of the detection task. Indeed, many potential conflicts disappear after updating aircraft positions. However, both an automatic solver or a human controller have to deal with these potential conflicts. The solver calculates maneuvers starting late enough to be updated in the future (unsent maneuvers). Human controllers have to keep the traffic safe and send as little unnecessary maneuvers as possible. These experimental results show that the number of maneuvers sent to pilots is a small portion of the total number of maneuvers calculated (for example with T w = 9 minutes, δ = 2 minutes and a (5%, 10%) uncertainties, only 30% of the maneuvers are sent to pilots). This suggests that the human controllers taskload is probably more dedicated to conflict detection than to conflict resolution.

Table II gives the results of the simulation for T w = 12 minutes. The same comments can be made on this table as on table I. The number of maneuvers sent to pilots decreases with the update frequency and increases with uncertainties. For example, with no uncertainty and δ = 2 minutes, the number of potential conflicts is 10166, the number of maneuvers sent to pilots is 2095 and the number of maneuvers unsent is is 6325. When the uncertainty increases to 5% horizontally and 10% vertically the number of maneuvers sent to pilots doubles (4010), the number of unsent maneuvers is multiplied by 2.5 (16087) and the number of potential conflicts is multiplied by 2.6 (26860). With T w = 9 min and δ = 2 minutes the number of maneuvers necessary to solve every conflict is multiplied by 3 when dealing with 10% and 30% uncertainty and the number of maneuvers calculated but not sent is multiplied by almost 4 and the number of potential conflicts is multiplied by almost 5. This is not surprising because unsent maneuvers deal with conflicts that might occur later and are more influenced by uncertainties. This suggests that uncertainty has a great influence on the controller's taskload and that a good Trajectory Prediction tool could probably help controllers to deal with uncertain conflicts and reduce their detection monitoring task.

Results also show that uncertainties have a big influence on the number of maneuvers needed to solve the conflicts. Controllers are used to saying that they often have to solve conflicts that might never happen because of uncertainties. It appears that an automatic tool would encounter the same issue. These experiments show how an automatic solver can model the taskload of human controller. It first needs to calculate many maneuvers to solve potential conflicts but tries to only make decisions that are absolutely necessary. The number of potential conflicts detected and the number of maneuvers needed to solve them give a better idea of the conflict detection complexity than the final number of actions that the controller really had to take.

V. Conclusion

The CAT Simulator cannot reproduce human controllers' behaviour because it uses the computer's power to detect and build solutions in advance and validate them when necessary. Human controllers rely on their expertise and experience on specific sectors to take the right action at the right moment. It is very difficult to reproduce the human behaviour. Researchers studying chess have managed to build softwares able to beat the best players by using the computer's ability to calculate very fast but not by imitating the players' way of thinking. However the amount of calculation necessary to replace the controller's decision on fast time simulation probably gives a better idea of the controller's taskload than the result of the calculation itself (i.e. the number of conflicts detected or maneuvers sent to pilots).

In this article we tried to show how the CATS solver can reproduce some of the controller's task that does not appear by simply counting the number of actions that were taken. The controller's permanent detection process is modeled by a sliding forecast time window that is updated every δ minutes. The smaller δ is, the more efficient the solver is.

Uncertainties are costly regarding the number of potential conflicts and unsent maneuvers that are detected. The controllers detection taskload is hard to measure because only a small proportion of conflicts end up by a maneuver order. The number of potential conflicts detected and unsent maneuvers can give a good idea of the extra amount of work controllers have to deal with because of uncertainties.

In order to complete the work presented in this paper it would be interesting to correlate the number of potential conflicts or unsent maneuvers with the human taskload. This is a difficult task because it would require defining scenarios on which we could both have measures of the human taskload and CATS simulation results. Exising data in France involve TMA (Terminal Control Area) traffic for wich CATS is not yet adapted. We look forward to work on this in the future.

Figure 2 .

 2 Figure 2. Detailed architecture of the prototype

Figure 3 .

 3 Figure 3. Modeling of speed uncertainties.

Figure 4 .

 4 Figure 4. Horizontal maneuver modeling.

Figure 5 .

 5 Figure 5. Vertical maneuver modeling.

Figure 6 .

 6 Figure 6. Vertical maneuver during the climbing period.

Figure 7 .

 7 Figure 7. Vertical maneuver during the cruising period.

Figure 8 .

 8 Figure 8. Vertical maneuver during the end of cruising period.

igure 9 .

 9 The model and real time optimization.

Figure 10 .Figure 11 .

 1011 Figure 10. Two aircraft potential conflict at time t=4,8,12 minutes.

Figure 12 .

 12 Figure 12. Crossover operator

TABLE I

 I Tw = 9 MINUTES:NUMBER OF REMAINING CONFLICTS,NUMBER OF MANEUVERED AIRCRAFT AND GIVEN AND CALCULATED BUT NOT GIVEN MANEUVERS FOR DIFFERENT δ VALUES AND DIFFERENT UNCERTAINTIES.

	T w δ	horiz	vert remain maneuvered maneuvers maneuvers	pot	norm
		uncert uncert conflict	acft	sent	not sent	confs pot confs
	12 2	0	0	0	1756	2095	6325 10166	1694
	12 3	0	0	0	1914	2338	2942	6304	1576
	12 4	0	0	0	2058	2496	1158	4355	1452
	12 2	2	5	0	2384	2924	10887 17200	2867
	12 3	2	5	0	2721	3457	5153 10699	2675
	12 4	2	5	0	3040	4125	2087	7771	2590
	12 2	5	10	0	3127	4010	16087 26860	4477
	12 3	5	10	0	3600	4840	7902 17444	4361
	12 4	5	10	0	4188	6160	2846 12637	4212
	12 2	10	30	0	5047	7809	24634 51917	8653
	12 3	10	30	0	5559	9164	11349 37738	9434
	12 4	10	30	65	6308	11680	3771 29729	9910

TABLE II Tw

 II = 12 MINUTES:NUMBER OF REMAINING CONFLICTS,NUMBER OF MANEUVERED AIRCRAFT AND GIVEN AND CALCULATED BUT NOT GIVEN MANEUVERS FOR DIFFERENT δ VALUES AND DIFFERENT UNCERTAINTIES.

Ground track uncertainties will not be considered, as they do not increase with time and will be included in the standard separation