Issue 
A&A
Volume 499, Number 1, May III 2009



Page(s)  321  330  
Section  Celestial mechanics and astrometry  
DOI  https://doi.org/10.1051/00046361/200811509  
Published online  25 March 2009 
Estimating the accuracy of satellite ephemerides using the bootstrap method
J. Desmars ^{1,2}  S. Arlot^{3}  J.E. Arlot^{1}  V. Lainey^{1}  A. Vienne^{1,4}
1  Institut de Mécanique Céleste et de Calcul des Éphémérides, Observatoire de Paris, UMR 8028 CNRS, 77 avenue DenfertRochereau, 75014 Paris, France
2  Univ. Pierre & Marie Curie, 4 place Jussieu, 75252 Paris, France
3  CNRS, Willow ProjectTeam, Laboratoire d'Informatique de l'École Normale Supérieure (CNRS/ENS/INRIA UMR 8548), 45 rue d'Ulm, 75230 Paris, France
4  LAL, Univ. de Lille, 59000 Lille, France
Received 12 December 2008 / Accepted 18 February 2009
Abstract
Context. The accuracy of predicted orbital positions depends on the quality of the theorical model and of the observations used to fit the model. During the period of observations, this accuracy can be estimated through comparison with observations. Outside this period, the estimation remains difficult. Many methods have been developed for asteroid ephemerides in order to evaluate this accuracy.
Aims. This paper introduces a new method to estimate the accuracy of predicted positions at any time, in particular outside the observation period.
Methods. This new method is based upon a bootstrap resampling and allows this estimation with minimal assumptions.
Results. The method was applied to two of the main Saturnian satellites, Mimas and Titan, and compared with other methods used previously for asteroids. The bootstrap resampling is a robust and practical method for estimating the accuracy of predicted positions.
Key words: planets and satellites: individual: Mimas  planets and satellites: individual: Titan  ephemerides  methods: statistical
1 Introduction
To compute the motion of solar system objects, we need a dynamic model including all significant dynamical interactions and nongravitational effects for small bodies. In order to quantify the orbital parameters, we need a set of observations. Fitting the model to the observations allows us to estimate the values of the initial conditions and parameters. We then are able to compute the position and velocity of the studied bodies at any time (either inside or outside the observation period).
The predicted positions include errors which have several causes. First, the quality of the theoretical model gives the internal error or even the precision of the theory. Second, the observations used for the fit of the parameters are the cause of the external error. That depends on the accuracy and the distribution of the observations. The observational errors are the main cause of the global error.
During the observation period, accuracy can be estimated by comparing observed and computed positions. Outside the period, this estimation is somewhat difficult. Many methods have been developed for asteroids in order to recover the ones lost after a few observations. Hence, astronomers have to estimate the accuracy of predicted positions of asteroids. Usually, these methods use few observations whereas for natural satellites, many observations are available. Consequently, methods used for asteroids have to be adapted to these objects.
The purpose of this paper is to show that the bootstrap method (Efron 1979; Efron & Tibshirani 1993) can be used successfully to estimate the accuracy of predicted positions inside and outside the observational period. This method is applied on two Saturnian satellites (Mimas and Titan). After comparison with two methods used for asteroids, we show that the bootstrap appears a robust and pratical method for estimating the accuracy of predicted satellite positions.
2 The dynamical model used: TASS1.7
To test the method of boostrap resampling, we used the orbital model TASS1.7 (Vienne & Duriez 1995). This is a theorical model of the motions of the eight major Saturnian satellites^{}. The main difficulty of this dynamical system comes from the various mean motion resonances: 2:4 in inclinations (MimasTethys), 1:2 in eccentricities (EnceladusDione) and 3:4 in eccentricities (TitanHyperion).
TASS theory has been developed using a much more complete dynamical model than Dourneau (1987) or Harper & Taylor (1993). The physical model takes into account Saturn's oblateness (J_{2}, J_{4} and J_{6}), the mutual interactions and the solar perturbation. It is constructed in a dynamically consistent way in which the satellites are considered together; its only parameters are the initial conditions, the masses of the satellites and the oblateness coefficients of Saturn.
First, the Lagrange equations of the osculating elements were developed in a complete and analytical way. A separation between the short period terms which are easily integrated analyticaly and the critical terms (secular, resonant and solar terms) was performed. The internal precision is controlled down to ten kilometers over one century.
There is an advantage in using TASS1.7 instead of a numerical integration. The numerical integration (see Appendix) may be more accurate than TASS for some sets of observations but the computation time is much longer, which is the main advantage in performing the statistical methods presented in this paper. Nevertheless, some tests have been done with the numerical integration for comparison and are presented in the Appendix.
The choice of the Saturnian system in the present study comes from the varied behavior of the satellites. The results presented here concern Mimas and Titan. The dynamics of Mimas is more complex (low order resonance) than that of Titan. Titan is also easier to observe than Mimas (far from Saturn and the rings). The dynamics of the MimasTethys system is regular over at least thousands of years but the orbital solution is very difficult to fit as the motion is very sensitive to the initial conditions. This fact comes from the large amplitude of the libration related to the inclination resonance type of the system. The resonant argument induces a libration of 70 degrees in the mean longitude of Mimas. A small change in the initial conditions induces a rapid separation between two orbits. Furthermore, the partial derivatives in TASS are fixed so they are not computed again between two adjustments. This high sensitivity of the MimasTethys system can explain the behaviours seen in the results presented in the next sections.
3 The observations used for the fit
Dynamical models have to be fitted to observations to provide accurate ephemerides. Fitting to observations involves determining the optimal parameters (initial conditions) by minimizing the difference between observed and computed positions (OC). The least squares method is usually used (see Sect. 4). For this work, we did not choose to determine the physical parameters (masses, J_{2}, J_{4}, J_{6}) as they are sufficiently well known from spacecraft observations (Voyager, Cassini).
The COSS08 catalogue (Desmars et al. 2009) is used in this paper. This catalogue is an extensive set of astrometric observations of the major Saturnian satellites covering the period 1874 to 2007. Figure 1 represents the distribution of observation nights per opposition of Saturn, for all the major satellites.
Figure 1: Histogram of observation nights per opposition. 

Open with DEXTER 
The distribution of these observations is particulary inhomogeneous. Two large gaps appear between 1930 and 1938, and between 1947 and 1961. Before 1947, observations were generally visual observations (micrometer) whereas since 1961, observations have generally been photographic plates and CCD frames. Consequently, the catalogue can be separated into two sets:
 old observations from 1874 to 1947, mostly visual and a priori with low accuracy;
 recent observations from 1961 to 2007, mostly photographic and a priori with better accuracy.
The main source of the ephemeris errors comes from observational errors that have many causes:
 the observer, reading the measurement of the position;
 the instrument used for the observation;
 the uncertainty of the star catalogue used for the reduction;
 corrections which have or have not been taken into account to determine the position of this object (refraction, aberration, ...);
 the difference between center of mass and photocenter due to the inhomogeneity of the surface of the object (phase, albedo, ...);
 the uncertainty of observation time, especially for old observations.
4 Fitting the observations
Fitting the model to observations consists of determining accurate initial parameters of the model
by minimizing the difference between observed and computed positions. The least squares method (Eichhorn 1993) allows this estimation. Generally, initial parameters
are the initial positionvelocity vectors (or osculating elements) and the masses of the studied objects. Denoting
the flow of the dynamical system projected on the position subspace and
the computed positionvelocity vector, we can write:
Fitting to observations amounts to determining , that is the variation of the initial parameter values, for which observed positions are assumed to be:
Thus, we can write to a first order approximation:
Using matrix notation, , and , the previous relation becomes:
As the observations are correlated and have various accuracies, we have to consider the covariance matrix of the observations . In least squares theory, this matrix is supposed to be known. The least squares method (LSM) allows us to estimate which minimizes with . The LSM solution is given by:
where the normal matrix N and covariance matrix are defined as: and .
The main difficulty is to choose the weighting matrix
.
A natural choice would be the covariance matrix of the observations, if it is known. As in most similar works, we chose in this paper to take:
where is the number of observations and is an estimator of the variance of the kth observation. It amounts to dividing each line of the matrix B by with representing the weight of the kth observation, which seems reasonable even if the observations are correlated.
The choice of the weights is detailed in Vienne & Duriez (1995). The observations are sorted according to the author, the instrument used and the observed satellite. Thus, each set of observations has a specific accuracy.
Finally, LSM remains a good estimator of the parameters but the value of the leastsquares criterion at its minimum underestimates the uncertainty of the leastsquares estimator.
5 Related works for asteroids
The problem of the extrapolation of the errors has been partly studied for asteroids or small satellites. For example, for lost asteroids, usual methods of orbit determination do not provide accurate ephemerides for their recovery. A solution is to determine not only the position but the region of the celestial sphere where the asteroid can be found. Several methods exist to determine this region (domain of possible motions). The classical way is to determine the whole family of the most probable orbits constructed using the LSM solution and the covariance matrix. The initial domain of possible motions can be defined by the LSM solution and by the covariance matrix , following a multivariate Gaussian distribution: .
Milani (1999) deals with the problem of the propagation of the normal and covariance matrix for the recovery of lost asteroids. He showed that the linearization hypothesis and so the LSM failed for asteroids observed over a short arc. In that context, he developed algorithms to approximate the recovery region.
Muinonen & Bowell (1993) suggest a statistical approach to the problem of orbit determination. They used Bayesian methods to determine asteroid orbital elements and developed Monte Carlo techniques for orbit determination. This statistical approach allows the estimation of a posteriori probability densities of orbital elements thanks to a priori information. In Muinonen & Bowell (1993), this information is a uniform distibution of the orbital elements. Likewise, Virtanen et al. (2001) use the present distribution of orbital elements of known asteroids as a priori information to constrain the a posteriori distribution of orbits. This method was successfully tested for lost asteroids.
Bordovitsyna et al. (2001) propose algorithms to determine the evolution of the domains of possible motions. These algorithms are based on the realization of a set of possible orbits thanks to the LSM solution and the covariance matrix. Avdyushev & Banshchikova (2007) use this method to deal with the region of possible motions for new Jovian satellites and show that the orbits of some satellites cannot yet be determined with acceptable accuracy.
For asteroids, the former statistical methods seem to be useful because asteroids are generally not much observed. In that case, the classical determination of orbits fitted to observations by LSM cannot always be satisfactory because not enough observations are available. For the main satellites of giant planets, the problem is quite different. Satellites are much observed, and over a large time period (see Fig. 1). Their motions are consequently well known. This better knowledge requires us to create more and more complex dynamical models which are used to produce ephemerides much further in the future. The least square method is used to determine the initial parameters of the model and to provide accurate ephemerides. Nevertheless, the accuracy outside the period of observations is still hard to estimate but the large number of observations allows to use resampling methods (see Sect. 6.3).
6 Methods for quantifying the extrapolated accuracy
We present three methods that we will use to determine the accuracy of predicted positions. Denoting the observed positions^{} at time respectively, the model provides the orbit using LSM fit to observations. The principle of the three methods is to determine the region of possible motions of the satellites using the set of observations .
The first two methods, which come from the study of asteroids, are a Monte Carlo method using the covariance matrix and a Monte Carlo technique applied to the observations. They have been adapted to study the satellites. The last method is the bootstrap.
6.1 Monte Carlo using the covariance matrix (MCCM)
The determination of the region of possible motions using the covariance matrix is probably the most classical method. It consists of simulating orbits using the covariance matrix and LSM solution (Bordovitsyna et al. 2001; Avdyushev & Banshchikova 2007).
The region of possible solutions can be constructed with K solutions:
for k=1,...,K, where is a pdimensional vector of normally distributed random numbers (where p is the number of parameters of the model), A the triangular matrix for which . A can be obtained by the Cholesky decomposition since is symmetric and positivedefinite.
In practice, the model is fitted to observations using LSM. Then we determined the covariance matrix of the parameters and the matrix A by the Cholesky decomposition. Sets of new parameters were computed with Eq. (7), all inside an hyperellipsoid.
The MCCM assumes that estimated parameters are Gaussian random variables with mean (the true parameters) and covariance matrix given by LSM.
6.2 Monte Carlo method applied to the observations (MCO)
The second method comes from a technique developed by Virtanen et al. (2001) and consists of generating orbits by adding random noise to a set of observations. The method of Virtanen et al. (2001) used for asteroids is summarized as follows:
 a method of orbit determination using two observations is used;
 two observations of an asteroid are randomly chosen in the set;
 a random error is introduced on each observation;
 a new orbit is determined with the two new observations;
 if the orbit gives acceptable positions for all observations dates, the orbit is kept. If not, the process starts again.
Contrary to the asteroid problem, the number of observations of satellites is greater and covers a longer period. Hence, the least squares method provides quite accurate satellite orbits. To determine the region of possible motion of satellites, we have adapted this method for a set of N observations :
 we choose the mean and the standard deviation of the random error;
 we create a new set of observations by adding to each observation , related to the Gaussian distribution : ;
 the model is fitted to the new set of observations and a new orbit is determined.
6.3 Bootstrap resampling (BR) and block bootstrap resampling (BBR)
The bootstrap method was first introduced by Efron (1979) in the context of variance estimation. The bootstrap has since been extended successfully to many other problems, such as estimating the distribution of the error of an estimator (see for instance Efron & Tibshirani 1993). Instead of adding some noise to each of the observed positions as in MCO, the idea of the bootstrap is to mimic the whole sampling process in order to create a new set of observations . This operation is also called ``resampling''. Each is obtained by sampling with replacement among the observations . In particular, some of the observations appear several times in the bootstrap sample, which amounts to giving them a weight corresponding to their number of occurences. Then, the model is fitted to the bootstrap sample through LSM and a new orbit is determined. As for MCO, this process can be repeated as many times as desired.
The bootstrap method applied to the estimation of the extrapolation error can be described as follows:
 generate a random set of independent integers with a uniform distribution in the range [1,N];
 build a new set of observations , the bootstrap sample;
 rit the model to the bootstrap sample which determines an orbit;
 repeat this process as many times as desired.
Note that observation errors are usually not independent. Hence, we have to modify the usual bootstrap for our problem. We use a technique similar to block resampling (Politis 2003) which was introduced in the framework of time series analysis. The data are first grouped into independent blocks and then the bootstrap method is applied to these blocks. The block bootstrap resampling can be described as follows:
 group the observations into B independent blocks with and k=1,...,B;
 generate a random set of independent integers with a uniform distribution in the range [1,B];
 build a new set of observations with and ;
 fit the model to the block bootstrap sample which determines an orbit;
 repeat this process as many times as desired.
7 Validation of the methods
We have tested the different methods for two particular Saturnian satellites: Mimas and Titan. Mimas' period is about 0.942 day whereas Titan's period is 15.945 days. Thus, the case of satellites with fast motion and the case with slow motions are studied.
7.1 Simulated observations
To validate the methods presented in the previous section, we first created simulated observations. The interest of creating simulated observations is to compute the real region of possible motion. The aim is to compare the simulated region of possible motion and the ones derived from the different methods.
7.1.1 Simulation of observations
To simulate observations close to reality, we introduce a dependence between the observations. Simulated observations are created in three steps:
 a set of N=3650 dates each 4 days (from 1960 to 2000) is chosen: t_{1},...,t_{N};
 considering positions given by the model as real positions, the positions (right ascension and declination ) for each observation date t_{i} are computed. This corresponds the ``initial orbit'' plotted in Fig. 2;
 for each month of the period, we compute a random variable related to a Gaussian distribution where and ;
 for each coordinate of each observation i, random noise is added, providing new coordinates: right ascension and declination where are normally distributed and independent. Observations performed during the same month have the same accuracy but observations made during different months may not. A dependence between the observations of the same month is introduced since for example, and are correlated for i and j in the same block.
Figure 2: Determination of the simulated region of possible motion with simulated observations. 

Open with DEXTER 
7.1.2 Simulating the region of possible motion
We create K=100^{} samples of simulated observations. A typical result is represented in Fig. 3 obtained for satellites Mimas and Titan with a hundred samples of 3650 simulated observations from 1960 to 2000. For the hundred orbits created, we plotted the difference in separation between positions of the kth orbit and the initial one:
where and .
Figure 3: Difference in separation between 100 orbits obtained by fitting to 100 different samples of simulated observations from 1960 to 2000 for Mimas and Titan. 

Open with DEXTER 
Figure 3 represents the real region of possible motions, in separation distance on the celestial sphere, for Mimas and Titan after fitting to observations from 1960 to 2000. During the period of observations (from 1960 to 2000) the difference is not large (less than 0.1''). Outside the period, the difference grows and can reach 9'' for Mimas after 200 years but remains less than 0.4'' for Titan.
The difference in the results between Mimas and Titan can be explained by the fast motion of Mimas (period of 0.942 day). Consequently, the uncertainty on its positions leads to an uncertainty on its velocity and so the divergence is greater. On the other hand, Titan is a slow motion satellite (period of 15.945 days) and the divergence is less.
7.1.3 The extrapolated standard deviation
Figure 3 represents the region of possible motions. To summarize the information of all the orbits, we introduce the extrapolated standard deviation which is a measure of the size of the region of possible motions over time.
For a time t and for each orbit k, we computed the distance s_{k}(t) which is the difference between the position given by the orbit k and the position given by the initial orbit.
(s_{k})_{k=1,...,K} are independent copies of a random variable S, and the standard deviation
of S(t) measures the uncertainty on the position at time t. This uncertainty is estimated by:
We call the extrapolated standard deviation associated with the separation. It represents at a time t the mean deviation coming from K orbits. is a measure of the size of the region of possible motions, which is a good indicator of the uncertainty of the position since both the estimated orbit (reference orbit) and the true orbit (initial orbit) belong to the region of possible motions with high probability (see Fig. 2). Thereafter, will be used for comparison of the different methods.
7.2 Comparison of methods
To compare the different methods, we choose one of the sets of simulated observations as the reference set of observations^{}. The reference orbit is the orbit fitted to the reference set of observations. Then we apply one of the three methods for determining the region of possible motions, containing the initial orbit assumed as the real orbit (Fig. 4). We compare the region of possible motions provided by the simulated observations (assumed to be real) and the one provided by one of the three methods, using the parameter (Sect. 7.1.3).
Figure 4: Determination of the region of possible motion with the reference set of observations. 

Open with DEXTER 
Figure 5 represents the comparison of the estimation of the standard deviation ( ) between MCCM (see Sect. 6.1) and simulation. The initial parameters (elliptical elements) of TASS are computed at Julian Epoch J1980 i.e. in the middle of the period of observations. For simulated observations, the covariance matrix of the observations is with (the mean of defined in Sect. 7.1.1).
Figure 5: Comparison of between simulations (in green crosses) and MCCM (in red pluses) for Mimas and Titan in the period 19602000. 

Open with DEXTER 
The method using the covariance matrix gives results very different from the simulation for Mimas. An explanation may be that MCCM relies on the assumption that the error on the initial parameters is Gaussian. As the motion of Mimas is very sensitive to initial conditions and the partial derivatives in TASS were fixed, the variancecovariance of the initial parameters is probably not well estimated by the LSM. Hence this method cannot give a good estimation of the positions for Mimas.
However, for Titan, the MCCM provides a good estimate of the extrapolated standard deviation ( ) because the two results seem correlated. Indeed, the correlation coefficient between the simulated value of the standard deviation ( ) and its estimation by MCCM is .
In pratice,
is computed for a set of P dates
.
For each method, we have a set
with
.
The correlation coefficient between
obtained with simulations
and
obtained with one of the methods
is defined as:
The second method MCO (see Sect. 6.2) has been used. We add independent Gaussian errors on the observations with and (the mean of defined in Sect. 7.1.1). The extrapolated standard deviation given by this method and the simulation are compared in Fig. 6.
Figure 6: Comparison of between simulations (in green crosses) and MCO (in red pluses) for Mimas and Titan in the period 19602000. 

Open with DEXTER 
The results obtained with MCO seem very close to those obtained with simulations. We note that the results also seem correlated. The correlation coefficient between simulation and MCO is for Mimas and for Titan. The two correlation coefficients close to 1 mean that the difference between the two methods is only a multiplicative factor, depending on the satellite. The MCO is based on more realistic hypotheses than the ones of MCCM. MCO assumes that observation errors are independent and have the same Gaussian distribution which is true in this particular case of simulated observations, except that they are weakly dependent. Thus, it is not surprising that MCO gives good results with such observations. Errors of real observations are not fully Gaussian and their standard deviations depend on many parameters (see Sect. 3) and obviously are not constant.
The bootstrap (BR) then was applied. The results are similar and the two curves are still correlated (Fig. 7). The correlation coefficients are for Mimas and for Titan. The results appear slightly more accurate than MCO, but without using the knowledge that observation errors are Gaussian.
Figure 7: Comparison of between simulations (in green crosses) and BR (in red pluses) for Mimas and Titan in the period 19602000. 

Open with DEXTER 
To deal with dependent data, a solution is to apply bootstrap block resampling (BBR, see Sect. 6.3). The observations are grouped into independent blocks. For simulated observations, the natural blocks are months of observations. The results appear in Fig. 8. The results between simulation and BBR are also correlated; the correlation coefficients are for Mimas and for Titan. The best method for estimating the extrapolated standard deviation seems to be BBR. However, BR and MCO give very similar results.
Figure 8: Comparison of between simulations (in green crosses) and BBR (in red pluses) for Mimas and Titan in the period 19602000. 

Open with DEXTER 
A first result after the comparison of the methods on simulated observations is that the method of the covariance matrix did not allow us to obtain good estimates because the partial derivatives in TASS are fixed. The three other methods (MCO, BR and BBR) allow us to obtain a good estimation of the region of possible motions.
Nevertheless, to deal with real observations, the bootstrap appears as the best method. Two points lead us to adopt it. The first point is that BR and BBR are ``nonparametric'' methods because no hypothesis on the distribution of errors is made. On the contrary, the two first methods (MCCM and MCO) are ``parametric'', so that they are not accurate whenever the hypotheses made (e.g., Gaussian errors with constant noise level for MCO) are not satisfied. Since the distribution of real observation errors is unknown (and certainly not Gaussian with constant noise level), it seems necessary to use nonparametric methods which are robust because they rely on fewer assumptions.
The second point concerns the implementation of the method. With simulated observations, the implementation is quite easy for the last two methods (MCO and BR). But with real observations, two problems could appear when determining the random error value for MCO:
 observations are given in different formats (absolute, differential coordinates or position angle and separation). Sometimes, especially for observations in the late 19th century and in the early 20th century, only one coordinate of the observation was available. Consequently introducing a random error on the observations can become difficult because the random errors have to be homogeneous. In particular, for the position angle given in degrees, the random errors added have to be homogeneous with the random error added to other coordinates (generally given in arcsec);
 the estimation of the value of the standard deviation of the random error. This value depends on the residuals themselves but also on the way of computing the residuals. If the observation is in intersatellite coordinates (observed satellite compared with reference satellite), the residual of the observation depends on the residuals of the positions of the two satellites. So the standard deviation of the random error will be different if we deal with intersatellite positions (the value will depend on the residuals of the two satellites) or if we deal with absolute coordinates (the value will a priori depend on the residual of the single satellite).
Figure 9: Extrapolated standard deviation of positions for Mimas and Titan after fitting to old observations (18741947) with BBR (in red pluses) and BR (in green crosses). 

Open with DEXTER 
We emphasize that the value of the standard deviation depends on the model. In the Appendix, we tested the bootstrap method with a numerical integration. The test of the methods shows that BR and BBR allow a good estimation of the extrapolated standard deviation. However, results are quite different, particulary in magnitude. This shows that the extrapolated error estimated in this paper depends on the model used. Nevertheless, we have shown that the bootstrap gives a good estimation of the standard deviation, whatever the model is.
8 Estimation of extrapolated errors
The bootstrap allows us to estimate the extrapolated standard deviation of the positions after fitting to real observations, assuming their independence. As we explained in Sect. 3, we applied the bootstrap to two sets of real observations (old and recent ones). This separation will allow us to estimate the extrapolated error with two different periods of obervations: 1874 to 1947 with a priori low accurate observations and 19612007 with a priori good accurate observations.
Figure 10: Extrapolated standard deviation of positions for Mimas and Titan after fitting to recent observations (19612007)with BBR (in red pluses) and BR (in green crosses). 

Open with DEXTER 
However, the distribution of the observations is not homogeneous. In fact, for some years, like 1995, many observations are available. For example, a satellite was observed over a hundred times on the same night. It is obvious that all these observations are not independent. The main hypothesis of BR is precisely the independence of the observations. The similar technique adopted is to group observations into independent groups. The choice of a block of independent observations is not as natural as it seems to be. In fact, we have to consider the cause of dependence between observations. We can reasonably think that observation errors mainly depend on the night of observation because the instrument used, the observer reading the measurent and the observation conditions are probably similar during the night. Consequently, we choose to group observations by night.
The estimation of the standard deviation of the positions was realized with a bootstrap without grouping observations into blocks (BR) and with the block bootstrap, grouping the observations by night (BBR).
This estimation, after fitting to old observations (18741947), is plotted in Fig. 9 for Mimas and Titan. The results given by the two methods are different in value. This difference reveals that there is probably a dependance between observations of the 18741947 period. As the majority of them are micrometric, they probably mainly depend on the observer and grouping observations by nights is not probably natural.
Furthermore, during the observation period, the standard deviation is quite small (about 0.05'' for Mimas and 0.02'' for Titan). Outside this period, the extrapolated standard deviation quickly diverges, particulary for Mimas, and is quite similar to simulations.
Figure 10 represents the estimation of the standard deviation after fitting to recent observations (19612007) with a simple bootstrap and with block resampling. The difference between the two methods is minimal. It appears that the observations between 1961 and 2007 are probably less dependent.
Nevertheless, the divergence for Mimas is more important after fitting to recent observations than after fitting to old observations. This is unexpected since recent observations are a priori better than old ones. This result can be explained because the old observation period stretches from 1874 to 1947 (73 years) whereas the recent observation period stretches from 1961 to 2007 (46 years). The 18741947 observations cover the period of the main term of the mean longitude of Mimas (70.56 years), which is not the case of recent observations. Thus, the mean longitude of Mimas is better estimated with old observations. Consequently, for a good accuracy outside the period of observations, a short period with accurate observations is not systematically better than a long period of average observations.
9 Conclusion
The bootstrap is a quite interesting method to estimate the accuracy of satellite positions over time. The advantages is the robustness, because the only restrictive hypothesis is the independence of observation errors without assumptions on the distribution of these errors, contrary to MCCM and MCO. This hypothesis can be avoided by grouping observations into independent blocks (to be defined according to each dataset). The implementation of BR is also quite easy and practical because no initial information is necessary.
The main constraint is probably the number of observations. To obtain an accurate estimate, the number of bootstrap resamples has to be high. In our particular framework, the restrictive point is the computation time of the fit. Consequently, only a hundred samples can reasonably be created. On the other hand, for asteroids, the number of observations has to be sufficient to allow enough bootstrap resamples to be drawn. In theory, with N observations, N^{N} bootstrap samples can be created.
The bootstrap can easily be adapted to most asteroids since many observations are available. This point will be the subject of a subsequent paper.
Acknowledgements
The second author was partially financed by Univ ParisSud (Laboratoire de Mathematiques, CNRS  UMR 8628) while the others were financed by IMCCE, CNRS  UMR 8028.
Appendix: Results of bootstrap with numerical integration
Numerical integration of the motion of satellites has been done to compare bootstrap results with those computed with TASS. The numerical software is the one used in Lainey et al. (2004a) but adapted to Saturnian satellites.
Equations of motions including perturbations (like J_{2}, J_{4} and J_{6}) are numerically integrated. The variational equations are simultaneously integrated with the equations of motion.
The fit to observations is similar to Lainey et al. (2004b). The positions are compared to observations and the new parameter values can be determined using a least square method (LSM). As an iterative process, the equations of motion and variational equations are integrated again with these new parameters. In practice, the process converges after three or four iterations. This numerical integration (called NUMINT) has been fitted using TASS theory. The numerical accuracy is about a hundred meters over 100 years.
Figure 11: Comparison of between simulations (in green crosses) and BR (in red pluses) for Mimas and Titan in the period 19602000 for NUMINT. 

Open with DEXTER 
Figure 12: Comparison of between simulations (in green crosses) and BBR (in red pluses) for Mimas and Titan in the period 19602000 for NUMINT. 

Open with DEXTER 
As in Sect. 7.1, 3650 simulated observations were created from 1960 to 2000. Only thirty sets of simulated observations were created because of computation time and so thirty new orbits after fit to observation set. It allows us to estimate the extrapolated standard deviation of the satellite positions associated with NUMINT. One of the simulated observation sets was chosen as the reference set of observations. We then applied the bootstrap to the reference set. Figure 11 represents the comparison of the extrapolated standard deviation between the simulation and BR for Mimas and Titan. Figure 12 represents the same comparison between simulation and BBR for Mimas and Titan.
BR and BBR give an estimation of the extrapolated standard deviation close to simulations with correlation coefficients (for Mimas) and (for Titan) for BR, and (for Mimas) and (for Titan) for BBR. However, compared to the TASS model, the value of this standard deviation is different. The accuracy of the predicted positions clearly depends on the model. The predicted positions are more accurate with the numerical integration, as suspected in Sect. 2. Nevertheless, the bootstrap remains a good method to estimate the accuracy of predicted positions, whatever the model used.
References
 Avdyushev, V. A., & Banschikova, M. A. 2007, SoSyR, 41, 413 [NASA ADS] (In the text)
 Bordovitsyna, T., Avdyushev, V., & Chernitsov, A. 2001, CMDA, 80, 227 [NASA ADS] [CrossRef] (In the text)
 Desmars, J., Vienne, A., & Arlot, J.E. 2009, A&A, 493, 1183 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Dourneau, G. 1987, Ph.D., Université de Bordeaux (In the text)
 Duriez, L., & Vienne, A. 1997, A&A, 324, 366 [NASA ADS]
 Eichhorn, H. 1993, Celestial Mechanics and Dynamical Astronomy, 56, 337 [NASA ADS] [CrossRef] (In the text)
 Efron, B. 1979, Bootstrap methods: another look at the jacknife Ann. Statist., 7(1), 1 (In the text)
 Efron, B., & Tibshirani, R. J. 1993, An Introduction to the Bootstrap, in Monographs on Statistics and Applied Probability (In the text)
 Harper, D., & Taylor, D. B. 1993 A&A, 268, 326 (In the text)
 Lainey, V., Duriez, L., & Vienne, A. 2004a, A&A, 420, 1171 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Lainey, V., Arlot, J. E., & Vienne, A. 2004b, A&A, 427, 371 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Milani, A. 1999, Icarus, 137, 269 [NASA ADS] [CrossRef] (In the text)
 Muinonen, K., & Bowell, E. 1993, Icarus, 104, 255 [NASA ADS] [CrossRef] (In the text)
 Peng, Q., Vienne, A., & Shen, K. X. 2002, A&A, 383, 296 [NASA ADS] [CrossRef] [EDP Sciences]
 Politis, D. N. 2003, Statist. Sci., 18, 219 [CrossRef] (In the text)
 Vienne, A., & Duriez, L. 1995, A&A, 297, 588 [NASA ADS] (In the text)
 Vienne, A., Thuillot, W., Veiga, C. H., Arlot, J. E., & Vieira Martins, R. 2001, A&A, 380, 727 [NASA ADS] [CrossRef] [EDP Sciences]
 Virtanen, J., Muinonen, K., & Bowell, E. 2001, Icarus, 154, 412 [NASA ADS] [CrossRef] (In the text)
Footnotes
 ... satellites^{}
 Mimas, Enceladus, Tethys, Dione, Rhea, Titan, Hyperion and Iapetus.
 ... positions^{}
 Here, means coordinate of positions and can be right ascension, declination or differential coordinates, etc.
 ...=100^{}
 We are limited by the computation time of the fitting procedure.
 ... observations^{}
 Similar results were obtained with other simulated observation sets.
All Figures
Figure 1: Histogram of observation nights per opposition. 

Open with DEXTER  
In the text 
Figure 2: Determination of the simulated region of possible motion with simulated observations. 

Open with DEXTER  
In the text 
Figure 3: Difference in separation between 100 orbits obtained by fitting to 100 different samples of simulated observations from 1960 to 2000 for Mimas and Titan. 

Open with DEXTER  
In the text 
Figure 4: Determination of the region of possible motion with the reference set of observations. 

Open with DEXTER  
In the text 
Figure 5: Comparison of between simulations (in green crosses) and MCCM (in red pluses) for Mimas and Titan in the period 19602000. 

Open with DEXTER  
In the text 
Figure 6: Comparison of between simulations (in green crosses) and MCO (in red pluses) for Mimas and Titan in the period 19602000. 

Open with DEXTER  
In the text 
Figure 7: Comparison of between simulations (in green crosses) and BR (in red pluses) for Mimas and Titan in the period 19602000. 

Open with DEXTER  
In the text 
Figure 8: Comparison of between simulations (in green crosses) and BBR (in red pluses) for Mimas and Titan in the period 19602000. 

Open with DEXTER  
In the text 
Figure 9: Extrapolated standard deviation of positions for Mimas and Titan after fitting to old observations (18741947) with BBR (in red pluses) and BR (in green crosses). 

Open with DEXTER  
In the text 
Figure 10: Extrapolated standard deviation of positions for Mimas and Titan after fitting to recent observations (19612007)with BBR (in red pluses) and BR (in green crosses). 

Open with DEXTER  
In the text 
Figure 11: Comparison of between simulations (in green crosses) and BR (in red pluses) for Mimas and Titan in the period 19602000 for NUMINT. 

Open with DEXTER  
In the text 
Figure 12: Comparison of between simulations (in green crosses) and BBR (in red pluses) for Mimas and Titan in the period 19602000 for NUMINT. 

Open with DEXTER  
In the text 
Copyright ESO 2009
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.