Aleš Bezděk
email: bezdek@asu.cas.cz

Calibration of accelerometers aboard GRACE satellites by comparison with POD-based nongravitational accelerations

Keywords: Space accelerometers, Nongravitational forces, Generalized least squares, Autoregressive processes

The proposed calibration method uses the precise kinematic positions derived from the data of the GPS receivers aboard the twin GRACE satellites (POD, Precise Orbit Determination). The total satellite accelerations are obtained numerically as a second derivative of the kinematic positions, from these the modelled forces of gravitational origin are subtracted. The resulting nongravitational accelerations then serve as a calibration standard for the uncalibrated accelerometer data. The calibration parameters for the GRACE accelerometers have already

been published using other methods. The aim of our study was to obtain not only the calibrated accelerometer measurements, but also a statistically correct estimate of their uncertainty.

The main problem in the application of a numerical derivative to observational data is the amplification of noise, especially at high frequencies. Besides, the filter of the numerical derivative introduces the correlation structure in the noise, which complicates the uncertainty estimates using the ordinary least squares. We succeeded in solving both of these problems by using the generalized least squares (GLS) method.

A c c e p t e d M a n u s c r i p t

Introduction

The wealth of quality data from the two GRACE satellites (launched in 2002), and also from its predecessor satellite CHAMP (launched in 2000), has substantially contributed to the improved modelling of the global Earth's gravity field, its static part as well as its temporal variations [START_REF] Reigber | CHAMP Mission 5 Years in Orbit[END_REF][START_REF] Schmidt | Static and Time-Variable Gravity from GRACE Mission Data[END_REF]. As the orbital altitude of these satellites is very low (below 550 km), they are equipped with space accelerometers, whose purpose is to measure the nongravitational accelerations. When processing the measurements from the CHAMP and GRACE missions to produce the gravity field models, the measurements from the onboard accelerometers have to be calibrated. The gravity on the ground is so much larger than the nongravitational accelerations measured in space that the electronic properties of a space accelerometer do not allow it to be calibrated before the launch. Many scientific teams using the CHAMP and GRACE data for the gravity field modelling therefore calibrated the accelerometer measurements Page 3 of 64 A c c e p t e d M a n u s c r i p t [START_REF] Frommknecht | Integrated Sensor Analysis GRACE[END_REF]Reigber et al., 2003Reigber et al., , 2005a)). Ideally, the calibrated accelerometer measurements should be accompanied with correct uncertainty estimates, but this is usually impossible, because the accelerometer calibration parameters constitute only a tiny part of the fitted parameters. Moreover, to stabilize the solution of large regression equations in the gravity field studies, one must usually use some regularization scheme, but then the regularized solution is biased and the bias could be much larger than the computed confidence intervals [START_REF] Aster | Parameter Estimation and Inverse Problems[END_REF]. Over the years it was found that the accelerometer calibration parameters can vary a lot depending on the analysis methods and the context of data usage (Bettadpur, 2004a).

Besides the gravity field modelling, the accelerometer measurements may also be used for the analysis of sources of the nongravitational forces themselves, especially for studying problems related to thermospheric density and winds [START_REF] Doornbos | Air density models derived from multisatellite drag observations[END_REF][START_REF] Flury | Precise accelerometry onboard the GRACE gravity field satellite mission[END_REF]. Specific for the attitude stabilized satellites CHAMP and GRACE are the firings of the attitude control thrusters, which show up in the linear accelerometer measurements, mainly because of thruster misalignments [START_REF] Frommknecht | Integrated Sensor Analysis GRACE[END_REF]. As real forces, the thruster firings are properly registered by the onboard accelerometers so that the full nongravitational signal can later be eliminated in the gravity field determination; this is one of the reasons why space accelerometers are useful in geodetic missions. However, from the point of view of aeronomy studies, the magnitude of the thruster events is often at the same order as that of the external nongravitational accelerations themselves, especially in the cross-track and radial directions. On the other hand, thermospheric density is derived from air drag, which is the dominant nongravitational acceleration in the along-track direction. The use of the properly calibrated ac-A c c e p t e d M a n u s c r i p t celerometer data for aeronomy studies was the main motive for writing this paper.

From the several calibration methods currently in use, we chose the satellite acceleration approach. The basic idea of the acceleration approach is to derive satellite accelerations by double numerical differentiation of the satellite positions along the precise orbit. Newton's second law of motion then links the resulting acceleration vectors to the forces acting upon the satellite. The successful implementation of this technique with results comparable to the classical, dynamic approach was enabled by the fact that kinematic orbits can nowadays be determined at a few cm accuracy. The satellite acceleration approach has been used by several scientific teams for the modelling of the geopotential (e.g., [START_REF] Ditmar | DEOS CHAMP-01C 70': a model of the Earth's gravity field computed from accelerations of the CHAMP satellite[END_REF][START_REF] Reubelt | Harmonic analysis of the Earth's gravitational field from kinematic CHAMP orbits based on numerically derived satellite accelerations[END_REF][START_REF] Švehla | From kinematic orbit determination to derivation of satellite velocity and gravity field[END_REF]. Numerically, the accelerometer calibration is much simpler compared to the gravity field determination, where one needs an inversion of normal matrix with tens of thousands unknowns and millions of measurements, a difficult computational problem, which requires special techniques to be applied. In this study, the calibration standard, the vector of the POD-based nongravitational acceleration, is projected into the accelerometer reference frame, where each component is directly compared with the uncalibrated accelerometer data and the linear least-squares calibration model may be used. In the ideal case, the residuals should be approximately independent and normally distributed to enable statistical inference concerning the regression results. On the other hand, the acceleration approach has the problem with the amplified noise.

The double numerical differentiation increases the noise in the positions proportionally to the squared frequency, and, therefore, the high-frequency noise will be amplified very significantly.

The prime motivation of this paper is that for a proper use of the accelerometer Page 5 of 64 A c c e p t e d M a n u s c r i p t measurements, and more generally of any observational data, one needs not only the measurement result, a point estimate of the true value, but also an estimate of the uncertainty of the result, a realistic error bar, which is a quantitative statement about where the true value 'really' is, with a given probability (cf., [START_REF] Taylor | Guidelines for Evaluating and Expressing the Uncertainty of NIST Measurement Results[END_REF]. Without error bars it is not possible to assess the quality of observations in question, to compare two competing theories using the observational data, to properly combine measurements from different sources, or to correctly combine the measurements even from the same experiment, if they have nonconstant variance.

Gravitational and nongravitational accelerations

GRACE project and SuperSTAR accelerometer

The Gravity Recovery And Climate Experiment (GRACE) is a joint US/German satellite mission [START_REF] Tapley | The gravity recovery and climate experiment: Mission overview and early results[END_REF] designed to very accurately map variations in the Earth's gravity field. The two almost identical GRACE satellites were launched in March 2002 into a near polar orbit at about 500-km altitude, separated by approximately 200 km. Each spacecraft carries a science payload consisting of microwave ranging system, GPS receiver, star cameras and accelerometer. Based on data from this mission, the most recent global Earth gravitational field models were published [START_REF] Förste | EIGEN-GL05C -A new global combined high-resolution GRACEbased gravity field model of the GFZ-GRGS cooperation[END_REF][START_REF] Pavlis | An Earth Gravitational Model to Degree 2160[END_REF][START_REF] Tapley | The GGM03 Mean Earth Gravity Model from GRACE[END_REF].

The SuperSTAR accelerometer on board of the GRACE satellites is a three axis capacitive accelerometer with two sensitive and one less sensitive axes. The sensitive axes point in the flight and radial directions, the less sensitive axis points in the cross-track direction. The precision of the sensitive axes is specified to be 10 -10 m s -2 , and that of the less sensitive axis 10 -9 m s -2 , within the bandwidth

5

A c c e p t e d M a n u s c r i p t of 2×10 -4 -10 -1 Hz [START_REF] Flury | Precise accelerometry onboard the GRACE gravity field satellite mission[END_REF]. Compared to the CHAMP accelerometer, the GRACE accelerometers have thermally controlled environment with the temperature variations below 0.1 K/orbit (Tapley and Reigber, 2002). ical properties; in this study, for the GRACE satellites we used the macro model and surface properties from [START_REF] Bettadpur | Product Specification Document[END_REF] and the mass from ISDC/GFZ data centre (http://isdc.gfz-potsdam.de/grace/). General formulae for computing the nongravitational accelerations may be found e.g. in [START_REF] Montenbruck | Satellite Orbits -Models, methods and applications[END_REF] or [START_REF] Milani | Non-gravitational perturbations and satellite geodesy[END_REF], in this study we used the model of neutral thermospheric density DTM-2000 [START_REF] Bruinsma | The DTM-2000 empirical thermosphere model with new data assimilation and constraints at lower boundary: accuracy and properties[END_REF] and the zonal and seasonal models of the Earth's albedo and emissivity [START_REF] Knocke | Earth radiation pressure effects on satellites[END_REF]. sudden spikes in the waveform correspond to the cold-gas thruster firings, which are activated on average every 2.3 minutes by the attitude control system in order to satisfy the pointing requirements of the microwave ranging system [START_REF] Flury | Precise accelerometry onboard the GRACE gravity field satellite mission[END_REF].

Nongravitational accelerations

Figure 3 should be positioned here.

Gravitational vs. nongravitational accelerations

The histograms in accelerations approximately equal in magnitude to the upper limit must start at degree 50-60/80-100/70-90 and go up to degree 125-150 to cover 1 nm s -2 lower limit of the nongravitational acceleration level, or up to degree 150-180 to reach 0.1 nm s -2 . In this study, we used the geopotential harmonic expansion up to degree and order 180 (or the maximum allowable value of a given model).

Figure 4 should be positioned here.

Method of calibration -a general look

In this section, we will explain the proposed method of calibration using the simulated positions and accelerations. To the simulated satellite positions we will add white noise of a known variance, to have an approximate representation of the POD positions. The uncalibrated accelerometer data will be represented by Filtering the positions yields the estimated second derivatives, the POD-based total acceleration vectors, from which the modelled gravitational accelerations are subtracted. In this way, the POD-based nongravitational acceleration vector is obtained, which serves as the calibration standard (etalon). The calibration equation then connects the mean curve, given here by the simulated nongravitational accelerations, with the calibration standard as the observation vector containing a random component. From this simple linear regression model, we find the bias and scale factor as the calibration parameters for each accelerometer component.

When filtering the positions, the filter of the second derivative introduces serial correlation into the random component of the POD-based nongravitational accelerations. While the mean values of the fitted calibration parameters are not much affected, the standard fit error and all the confidence intervals are not correct. The generalized least squares method (GLS) is used to find the correct estimates of the uncertainty in the calibrated nongravitational accelerations.

An important aspect of the presented calibration method is that we use the kinematic orbits, i.e. those determined directly from GPS measurements and not influenced by any force models (cf. [START_REF] Ditmar | DEOS CHAMP-01C 70': a model of the Earth's gravity field computed from accelerations of the CHAMP satellite[END_REF]. This is of concern es-

Filter of the second derivative

We obtain the POD-based total accelerations a POD TOTAL by double differentiation of the positions r. For this purpose we used the Savitzky-Golay or polynomial smoothing filters (e.g., [START_REF] Press | Numerical recipes in Fortran[END_REF]. A polynomial of a chosen order is least-squares fitted to the data points within a running window of a chosen length;

the approximate numerical derivative at the central point is obtained by the differentiation of the fitted polynomial.

We looked for the best agreement between the simulated and POD-based nongravitational accelerations, when no noise in positions is introduced. We started with the first approximation to the numerical second derivative, the simple three--point formula, but we found that such low order derivatives produce too high a bias (10 -6 m s -2 with the time step of 1 sec, 10 -4 m s -2 with 10 sec) between the simulated and POD-based nongravitational accelerations. We then systematically tested many combinations of the polynomial orders and window lengths to find a suitable pair with low values of both parameters that would produce a satisfactory agreement between the simulated and POD-based nongravitational accelerations.

A c c e p t e d M a n u s c r i p t

Finally, we have chosen the combination of the polynomial order 6 with the window length 9; other combinations, e.g. 8/13, 9/11, 9/21 yielded similar results.

The tested combinations comprised also the case with no smoothing, where the window length equals the polynomial order plus one, but again, the bias was too high for our purposes. For later reference, we will symbolically write the filtering of positions as the convolution of the second-derivative filter F and the radiusvector r, a POD TOTAL = F * r.

(1)

POD-based nongravitational accelerations

The calibration standard, the POD-based nongravitational acceleration vector a POD NG , is obtained from the POD-based vector of total accelerations a POD TOTAL by subtracting the modelled accelerations of gravitational origin a SIM GRAV ,

a POD NG = a POD TOTAL -a SIM GRAV , (2)
where the vector a SIM GRAV is the sum of the acceleration vectors caused by the Earth static gravitational field, direct lunisolar perturbations, solid Earth and ocean tides, and relativistic effects (Sec. 2.3). The relatively high degree and order of the geopotential model, which is necessary for the generation of gravitational accelerations of low enough magnitude comparable to that of the calibrated accelerometer measurements, was discussed in Section 2.4.

While the numerical differentiation of the positions is most easily done in the (inertial) celestial reference frame, the POD-based nongravitational accelerations obtained in Eq. ( 2) must be projected into the science reference frame, in which all GRACE Level-1B data products are specified [START_REF] Case | GRACE Level 1B Data Product User Handbook[END_REF]. derivative filter, the 1-cm noise in positions is amplified to high-frequency noise in accelerations with oscillations on the order of 10 -4 m s -2 . The "true" signal a SIM NG of amplitudes 10-500 nm s -2 is buried in noise.

Figure 5 should be positioned here.

Calibration equation

The calibration equation is given by the linear model

a POD NG = B + S a UNCAL ACC + ǫ, (3)
where B is bias, S scale factor, a UNCAL ACC uncalibrated accelerometer data, ǫ statistical error. On the assumption that the accelerometer measures independently in its three axes, we have one independent calibration equation (3) for each accelerometer axis.

In this section, the uncalibrated data a UNCAL ACC are represented by the simulated

13

A c c e p t e d M a n u s c r i p t nongravitational accelerations a SIM NG , which were scaled by S =1.1 and shifted by B=1.2×10 -6 m s -2 .

Problem of autocorrelated noise

The probability model, y = b 0 + b 1 x + ǫ, for which the ordinary least squares (OLS) method of estimation is best suited, relates the error-free predictor variable

x and the random variable y (see Appendix A). In this respect, the calibration equation ( 3) matches well the OLS model: the noise in the simulated nongravitational accelerations a SIM NG ≡x is several orders of magnitude lower than that of the response variable a POD NG ≡y (Fig. B.5). Also the noise in the accelerometer readouts should be, according to the specifications (Sec. 2.1), much lower than that of a POD NG .

The OLS provide correct uncertainty estimates, if the errors ǫ are independent and normally distributed. If the random errors are positively correlated, the uncertainty in the fitted parameters is usually underestimated, thus giving a false sense of accuracy (e.g., [START_REF] Chatterjee | Regression Analysis by Example[END_REF][START_REF] Rawlings | Applied Regression Analysis: A Research Tool[END_REF].

When a digital filter is applied to a data sequence containing a random component, the random errors within the filter window are linearly combined to the new output value; hence the newly formed random vector has components, which are correlated. This happens to the POD-based nongravitational accelerations a POD NG obtained from the positions by applying the second-derivative filter (1) and after subtracting the modelled accelerations of gravitational origin in Eq. ( 2); the noise in positions, which in this section is supposed to be white (Sec. and S , we will use the generalized least squares method; see Appendix B for a short review.

Use of GLS to remove autocorrelation

In fact, the non-diagonal covariance matrix of the random component in a POD NG was created by the action of the second-derivative filter F from the covariance matrix of the white noise Var(Z i )=σ 2 1. Namely,

Var(ǫ) = FVar(Z)F ′ = σ 2 FF ′ , ( 4
)
where F is a square matrix, generated from the coefficients of the filter F and whose multiplication is equivalent to the action of the filter (e.g., [START_REF] Gray | Toeplitz and Circulant Matrices: A review[END_REF].

But the situation, where we know the covariance matrix of the random errors in a linear model, is exactly what the GLS method is suited for. In our case, finding the GLS transformation matrix is straightforward, W=F -1 . After applying W to the calibration equation (3), and solving the transformed equation (Eq. B.3) through the OLS, the residuals become again uncorrelated and the original σ 2 should be recovered. As regards the implementation of the filtering, we throw away the first and last few acceleration points during the filter warm-up phase, and we find the transformation matrix W through the Cholesky decomposition of the covariance matrix FF ′ (cf. Eq. B.2).

Decorrelation of the observations

The results of the GLS transformation of the POD-based nongravitational ac- Figure 6 should be positioned here.

Very high correlation between the calibration parameters

In [START_REF] Chatterjee | Regression Analysis by Example[END_REF][START_REF] Rawlings | Applied Regression Analysis: A Research Tool[END_REF]Weisberg, 2005). i.e. a mean ratio of the oscillatory amplitudes (provided that x and y are in phase, which is true here). But this is not the case of the calibration equation (3); here the parameter b 1 multiplies the predictor x, which is a sum of the constant x and oscillations (xx), but the predictor connected with b 0 is also a constant, hence the collinearity. What makes the correlation between b0 and b1 so high is the very large value of the offset x in the accelerometer readouts compared to the amplitude of the oscillations (xx). For large sample sizes and x2 ≫ σ2

x , where σ2

x is the sample variance of x, we may approximate the expression for the coefficient of correlation (Eq. A.6) by ) 2 m 2 s -4 and that of the oscillatory component σ2

ρ( b0 , b1 ) = -x σ2 x + x2 ≃ - x | x| 1 - 1 2 σ2 x x2 . ( 5
x ≃(5.10 -8 ) 2 /2 m 2 s -4 give ρ( b0 , b1 ) ≃ 0.9995.

The extremely high correlation between the parameters b0 and b1 may be avoided by changing the calibration model (3). From Eq. ( 5), the correlation between the parameters in the simple linear regression is zero, if the predictor x has zero mean. In the notation of Appendix A, a modified calibration model might be

y -x = b ⋆ 0 + b ⋆ 1 (x -x) + ǫ, (6)
b 0 = b ⋆ 0 + x(1 -b ⋆ 1 ), b 1 = b ⋆ 1 . (7)
We might believe that the 'statistically better', completely uncorrelated parameters b ⋆ 0 , b ⋆ 1 and their uncertainties would somehow help b 0 , b 1 to have less correlation -but this does not happen; starting from (7) and using the rules for variances of the linear functions of random variables (e.g., [START_REF] Rawlings | Applied Regression Analysis: A Research Tool[END_REF], we arrive at exactly the same formulae (A.4), (A.6) as before.

In this study, for the regression calculations themselves we used the modified model ( 6

Calibration of the accelerometer data over several revolutions

In this section we will apply the calibration method to the real GRACE data covering several orbital revolutions in order to analyze the calibration results in more detail. As the POD positions, we used the high-quality 10-second kinematic orbits of the GRACE satellites, kindly provided by D. Švehla (TU Munich). The orbits were computed using the zero-difference ionosphere-free phase measurements, the 10-sec orbits are based on the interpolated 30-sec POD satellite clocks [START_REF] Švehla | Kinematic positioning of LEO and GPS satellites and IGS stations on the ground[END_REF].

The simulated gravitational accelerations, needed for obtaining the POD-based nongravitational accelerations (Sec. 3.3), the coordinate transformations and the simulated nongravitational accelerations were calculated by our own orbital propagator NUMINTSAT [START_REF] Bezděk | Simulation of free fall and resonances in the GOCE mission[END_REF]. When working with the real-world data, it has become clear that in contrast to simulations the use of the most up-todate physical models is crucial for obtaining meaningful calibration results. We used: coordinate transformations between ICRF and ITRF systems [START_REF] Mccarthy | IERS Conventions[END_REF], the model of static gravitational field EIGEN-5C to order and A c c e p t e d M a n u s c r i p t degree 180 [START_REF] Förste | EIGEN-GL05C -A new global combined high-resolution GRACEbased gravity field model of the GFZ-GRGS cooperation[END_REF], lunar and solar ephemerides JPL DE405, the model of solid Earth tides (anelastic Earth;[START_REF] Mccarthy | IERS Conventions[END_REF], the model of ocean tides CSR 4.0 (Bettadpur, 2004b).

We obtain the POD-based nongravitational accelerations a POD NG in Eq. ( 2 . This is not surprising, the kinematic orbits are reported to be correlated [START_REF] Švehla | From kinematic orbit determination to derivation of satellite velocity and gravity field[END_REF]. On the other hand, the standard error of the OLS fit σ of a few centimetres as an estimate of the noise in the real kinematic POD positions is a plausible value.

Page 21 of 64

Removing the autocorrelation with an AR model

In this subsection, we will use a general approach for drawing statistical inferences from time series [START_REF] Brockwell | Introduction to Time Series and Forecasting[END_REF][START_REF] Chatfield | The analysis of time series -An introduction[END_REF]. In most practical problems involving time series we see only one realization, but we imagine it to be one of the many sequences that might have occurred. It is necessary to setup a hypothetical probability model to represent the data; after an appropriate family of models has been chosen, it is then possible to estimate parameters, check for goodness of fit to the data, and possibly to use the fitted model.

We suppose that the correlated OLS residuals (middle panel of well with the sample ACF for lags less than 100; in our experience, the order 7 is sufficient to match the correlation structure of the OLS residuals.

We suppose that the OLS residuals may be viewed as a realization of the fitted AR(7) process, in other words, as an output to filtering a white noise input by the corresponding AR filter. Therefore, the covariance matrix of the correlated residuals in Fig. B.8 is now given as that of the fitted AR process. This new covariance matrix replaces the matrix Var(Z) in Eq. ( 4) and the GLS method is applied in the same way as in Section 3.6. We will use the subscript 2 to distinguish the new GLS transformation. The GLS 2 transformation matrix W 2 is obtained numerically by the Cholesky decomposition of the new covariance matrix (Eq. B.2). After transforming the calibration equation (3) using W 2 , and using the OLS 2 estimation to find the calibration parameters, we finally obtain an approximately uncorrelated series of residuals, in the middle panel of Figure B.9.

Indeed, the ACF and PACF (bottom left panel) are negligible except at zero lag.

Figure 9 should be positioned here.

Calibrated accelerometer measurements

On solving the calibration equation ( 3) by the GLS 2 method described in the previous section, we obtained the calibrated accelerometer measurements a CAL ACC and their estimated uncertainty band σ(a CAL ACC ) given by the confidence interval (B.5).

The fact that the GLS 2 residuals appear to be approximately uncorrelated and normal (bottom panels of eter should be located with a high level of confidence within the ±3 σ(a CAL ACC ) band around a CAL ACC . This is in accordance with the usual definition of the 99.7-percent confidence interval, within which we expect the 'true' value of the estimated parameter to be located with the coverage probability of 99.7 %, when the normal distribution is sampled ('three-sigma rule'). For the statement of uncertainties in this study, we used the coverage factor (CF) of 1 ('one-sigma' uncertainty, coverage probability 68.3 %) or that of 3 (coverage probability 99.7 %).

The calibrated accelerometer measurements a CAL ACC together with the 3 σ(a CAL ACC )

uncertainty band for two orbital revolutions are in Figure 10 should be positioned here.

Similar calibration results have been obtained also for the radial component;

the mean uncertainty σ(a CAL ACC ) is around three times larger, but the normality of the GLS 2 residuals is questionable. In the cross-track direction, we have not succeeded to find a suitable AR process to decorrelate the GLS 1 calibration residuals.

So, in the cross-track and radial directions, we found the calibration parameters b0 and b1 , but we are not able to calculate a reliable estimate of the uncertainty of a CAL ACC . From the point of view of the atmospheric density modelling, this is not a problem, by far the strongest signal from the atmospheric drag is in the along-track component and besides, the cross-track and radial components of the accelerometer readouts contain the disturbing signal from the attitude thrusters.

Evolution of calibration parameters over 1.5 years

The presented calibration method has been applied to the accelerometer data of both GRACE satellites within a period of 1. are better suited for a comparison of different gravitational models and calibration algorithms than a few days studies, where chance may play a role.

Long-term values of the obtained uncertainties

In the regression analysis, the squared standard error of the fit σ2 (A. 5) is an estimate of the constant variance of the observations, provided the assumptions of the OLS are met. As a factor, σ then enters the uncertainty estimates (A.4, A.7, A.8).

Although the correlated noise in the POD positions prevents the usual 3-σ interpretation of the OLS 1 residuals (in the middle panel of

Long-term evolution of scale factors and biases

In the long term, the scale factor b1 of the accelerometer data is approximately constant (upper panel of Figure B.12), with the mean value near 1 for both GRACE satellites, with the 3-σ uncertainty of a few percent. Using the fitted value of b1 and Eq. ( 7), the biases b0 are obtained, which we can subsequently fit with a straight line regression model to obtain simple long-term expressions (lower panel of Fig. B.12), similarly to Bettadpur (2004a).

In Figure B.13, there are the results of the same procedure applied to the modelled nongravitational accelerations. While the long-term statistical results of the scale factor b1 are comparable for both accelerometer-based and simulated accelerations, the biases are different: on average, the simulated nongravitational accelerations are very close to the calibration standard a POD NG , the fitted mean value of b0 is less than 0.01 nm s -2 ; but the variation in the straight-line model of the bias is 3-8 times greater in the simulated accelerations than in the accelerometerbased accelerations. In other words, the long-term accelerometer bias is more stable with respect to the calibration standard than the bias of the simulated nongravitational accelerations. This may be attributed to the fluctuating errors in the nongravitational acceleration models, which depend on the orbital conditions.

Figure 12 should be positioned here.

Figure 13 should be positioned here.

A c c e p t e d M a n u s c r i p t

Comparison with the calibration parameters from an independent study

In a technical note, Bettadpur (2004a) states the constant scale factors and simple linear (or quadratic) models of the changes in bias for each accelerometer axis of the GRACE A/B satellites. These estimates were obtained in the GRACE data processing for the precise orbit and gravity field determination, and their limit of applicability is from the launch until 1 November 2003.

Considering the very high correlation between the fitted calibration parameters

Uncertainties for different gravity field models

In Table B.1, there are the long-term means of the estimated RMS of noise in the POD positions σOLS 1 and of the uncertainty in the calibrated accelerometer measurements σ(a CAL ACC ) obtained using selected models of the static gravity field. We calculated the accelerations for degree/order 180 or less, according to the definition of the model (indicated by superscripts).

In the first group, there are the most recent models based also on the data from the GRACE mission: EIGEN-5C [START_REF] Förste | EIGEN-GL05C -A new global combined high-resolution GRACEbased gravity field model of the GFZ-GRGS cooperation[END_REF], EGM08 [START_REF] Pavlis | An Earth Gravitational Model to Degree 2160[END_REF], GGM03C/S [START_REF] Tapley | The GGM03 Mean Earth Gravity Model from GRACE[END_REF]. These models provided the best results; the accelerometer calibration also does not indicate any statistical difference between the results from the combination and satellite-only gravity field models GGM03C and GGM03S.

The second group in Table B.1 are models computed using the CHAMP data, but not those from GRACE: EIGEN-CHAMP03S (Reigber et al., 2005b), DEOS CHAMP-01C 70 [START_REF] Ditmar | DEOS CHAMP-01C 70': a model of the Earth's gravity field computed from accelerations of the CHAMP satellite[END_REF]. To test the influence of including the higher degree terms of the static geopotential models on the proposed accelerometer calibration, we also used the EIGEN-5C model limited to degree/order 70 (Sec. 2.4). From the statistical point of view, the results of this group of models are equivalent to the GRACE models. While there is no visible change in the results pertaining to the along-track component, those of the cross-track and radial components display a slight systematic decrease in the precision for the models with the maximum degree/order 70, which might be attributed to the lower magnitude of the nongravitational accelerations in these directions (Fig. B.3). So for a precise accelerometer calibration it is better to include the higher degree/order gravity terms.

The results based on the pre-CHAMP gravity models EGM96 (Lemoine et al.,

29

A c c e p t e d M a n u s c r i p t 1998) and GRIM5C [START_REF] Gruber | GRIM5-C1: Combination solution of the global gravity field to degree and order 120[END_REF] are worse by a factor of about 4 in the cross-track and radial components. Thus the proposed accelerometer calibration provides an indirect evidence that the gravity missions CHAMP and GRACE have appreciably contributed to improve the higher degree/order terms of the current global static gravity field models.

Table 1 should be positioned here.

The main purpose of including this section was to show that the calibration method does not depend upon a particular gravity model used, in other words that it is plausible to suppose that the errors in the accelerations derived from the gravity field model are negligible compared to those of the accelerations derived from kinematic positions. This is clearly demonstrated by the long-term results in Table 1, where the four most recent gravity models, derived by different groups using different processing schemes, give statistically equivalent results in all three accelerometer components. Besides, if nowadays the best available gravitational model EGM08 goes up to degree/order 2159, and the new EIGEN or GGM models go up to degree/order 360, then we may expect that they are consistent in predicting the geopotential functionals with a relatively low limit of degree/order less than 150 and that they should generate rather close vectors of the gravitational acceleration.

Discussion

As mentioned in Section 1, many scientific teams have calculated the calibration coefficients of the GRACE accelerometers for periods of differing length,

Page 31 of 64

A c c e p t e d M a n u s c r i p t from days to years. The question of the accuracy of the calibrated accelerometer measurements, however, seems not to be discussed very much, as either the primary research objective in other studies is the gravity field modelling, or the complexity of the calibration process prevents the uncertainty estimates from quantifying, e.g. due to regularization.

Van den Ijssel and Visser (2007) estimated the nongravitational accelerations for the CHAMP and GRACE A satellites as piecewise constant empirical accelerations via the reduced-dynamic POD approach. To obtain a solution, regularization was necessary. Only the longer wavelengths were recovered, at best in the along-track direction, with a bias in the cross-track direction. The authors concluded that no meaningful solution could be obtained in the radial direction.

Van Helleputte et al. ( 2009) used the reduced-dynamic POD technique to determine the calibration parameters of the CHAMP and GRACE A/B satellites over a 5-year period. The method needs strong constraints to be set on the a priori bias values in the cross-track and radial direction.

For the derivation of the satellite accelerations from kinematic positions, [START_REF] Reubelt | Harmonic analysis of the Earth's gravitational field from kinematic CHAMP orbits based on numerically derived satellite accelerations[END_REF] used the second derivative of the Gregory-Newton interpolation scheme;

the explicitly stated coefficients of the 9-point filter are the same as those from the second derivative of a 9-point polynomial filter of order 8 (i.e. with no smoothing, cf. Sec. 3.2). The choice of this filter was driven by the aim of the study, which was the determination of the gravity field parameters from 2 years of the CHAMP kinematic orbits without a regularization to guarantee an unbiased solution.

There are several scientific groups, which used the fitted ARMA models when solving the inverse problem of the gravity field determination, but with different aims and details of implementation compared to our method. In the context of A c c e p t e d M a n u s c r i p t

processing the future GOCE gradiometer data, [START_REF] Schuh | The Processing of Band-Limited Measurements; Filtering Techniques in the Least Squares Context and in the Presence of Data Gaps[END_REF] used the discrete linear filters and the GLS method for handling the correlated measurements in the frequency domain. The target was to obtain decorrelated observational equations and to distribute the computational effort to a cluster of computers. A need to treat the huge least-squares problems in the gravity field determination motivated [START_REF] Klees | How to handle colored observation noise in large least-squares problems[END_REF] and [START_REF] Ditmar | Frequency-dependent data weighting in global gravity field modeling from satellite data contaminated by non-stationary noise[END_REF] to study how the coloured noise represented by the ARMA processes might be used as a fast method to solve a Toeplitz system of linear equations. [START_REF] Ditmar | Frequency-dependent data weighting in global gravity field modeling from satellite data contaminated by non-stationary noise[END_REF] points out that the assumption about the stationarity of the noise in the kinematic POD positions may not be realistic in many cases, due to a quickly changing constellation of visible GPS satellites for a LEO satellite, and therefore, the orbit accuracy may vary considerably in time. This might be the reason for the increase in the estimated RMS of the POD positions with longer length of the fitting windows (Sec. 5.1).

The fact that the RMS of noise in the cross-track and radial components of the POD positions is several times worse, when comparing the accelerometer calibration statistics based on the pre-CHAMP gravity field models with those using the recent models including the CHAMP and GRACE data (Sec. 5.4), is in accordance with a similar improvement in the accuracy of the radial orbit component of the altimeter satellites [START_REF] Klokočník | Evaluation of the accuracy of the EIGEN-1S and -2 CHAMP-derived gravity field models by satellite crossover altimetry[END_REF][START_REF] Klokočník | Improvement in the radial accuracy of altimeter-satellite orbits due to the geopotential[END_REF].

Conclusions

In this study it was demonstrated that the proposed method of calibration of the linear accelerometer measurements is capable of finding the point esti- The calibration procedure makes use of the generalized least squares method, which might be useful in other linear regression problems, where one has to deal with the correlated residuals. In the case of the accelerometer calibration, the situation is particularly convenient for the application of the GLS method, as we know exactly the regression mean function, equal to the uncalibrated accelerometer measurements, and we need to shift it to the "right place" determined by the calibration standard.

From the point of view of aeronomy and atmosphere research, the most important is the along-track component of the accelerometer data, where the signal from the atmospheric drag is dominant; moreover, the cross-track and radial components of the accelerometer data contain the relatively strong disturbing signal due to the action of the attitude control thrusters.

Throughout the study, we have also used the modelled nongravitational accelerations, whose waveform matches well that of the accelerometer readouts but is generally smoother, and in the cross-track and radial components it does not contain the spikes caused by the attitude thrusters. After the calibration of the along-track component, the accelerometer data and the modelled nongravitational accelerations have approximately the same mean uncertainty; this is due to the fact that the GLS calibration effectively integrates the acceleration signal, so in the calibration only the longer period waves are actually used. This is closely connected with the fact that the calibration standard is calculated from the orbital positions.

Page 37 of 64 where P * ≡X * (X * ′ X * ) -1 X * ′ is the 'hat matrix' of the transformed model. 6 2002.8 2003 2003.2 2003.4 2003.6 2003.8 2004 2004.2 4 4.5 ) GLS bias (GB) Bettadpur(2004) 2002.6 2002.8 2003 2003.2 2003.4 2003.6 2003.8 2004 2004.2

Figure B. 1

 1 Figure B.1 shows the simulated nongravitational accelerations acting on the GRACE A satellite during one orbital revolution. The projection of the acceleration vectors refers to the satellite local reference frame; the three components are the along-track (A-T; projected to the velocity direction), the cross-track (C-T; direction of angular momentum) and the radial one (RAD; completes the righthanded system). The figure is typical for satellites in low Earth orbits (LEO, altitudes 100-2000 km, mainly 150-800 km): the dominant nongravitational accelerations change with the satellite local reference frame directions; the close-Earth motion makes the satellite to pass through the Earth's shadow, which is visualized by the characteristic jumps. In the along-track component, the main nongravitational driver is the atmospheric drag (DRAG), pointing always in the direction opposite to the satellite's motion. Even in the along-track component, there may appear jumps in the smooth waveform of the drag acceleration caused by the direct solar radiation pressure (DSRP), whose action is dominant in the sunlit part of the cross-track and radial components. In the shadow of Earth, the terrestrial infrared radiation (IR) dominates the radial component. Sometimes, when the satellite passes directly below the Sun, also the signal from the reflected solar radiation (ALB) may be recognizable in the graphs of the nongravitational accelerations. In each panel, there is also the sum of the individual simulated nongravitational accelerations, a SIM NG . The magnitude of the nongravitational forces in each local reference frame direction depends on the satellite shape and its phys-

Figure 1

 1 Figure 1 should be positioned here.

Figure 2

 2 Figure 2 should be positioned here.

Figure B. 2

 2 Figure B.2 displays the Level-1B accelerometer data of GRACE A during the same period as in Figure B.1. There is an apparent similarity between the waveforms of the sum of the simulated nongravitational accelerations (Fig. B.1) and the uncalibrated accelerometer readouts (Fig. B.2). This is typical for all GRACE Level-1B accelerometer data and provides evidence that the smoother simulated nongravitational accelerations and the accelerometer measurements are consistent with each other. On the other hand, if we compare the units on vertical axes of graphs in Figures B.1 and B.2, it is clear that the accelerometer data are not calibrated; for example, it follows from the geometry of the GRACE A motion during the revolution in question that in the radial component the nongravitational acceleration must pass through zero. In the cross-track and radial components, the

  Figure B.3 show the magnitude of the individual accelerations in the satellite local reference frame components. We simulated the orbital evolution of the GRACE A satellite during 1.5 years, every 60 minutes we recorded the magnitudes of the accelerations acting on the satellite and then draw a histogram for each acceleration. We do not show the specific numbers for the histogram counts on the vertical axis, which is linear, as these are only formal depending on the sampling period and would add complexity to the graphs.The dominant acceleration is due to the static gravitational field; the acceleration caused by the central term (GRAV µ/r; 8.5 m s -2 ) is projected mainly in the radial direction because of the almost circular orbit of the GRACE satellites. Then follows the acceleration due to the Earth flattening (GRAV J 2 ) and to the remaining terms of the geopotential (GRAV rest). Considering the range of the nongravitational accelerations (DRAG, DSRP, ALB, IR: 1-500 nm s -2 ), it is clear that for a successful accelerometer calibration also the other accelerations of gravitational origin must be taken into account: direct lunisolar perturbations (LUNISOL), solid Earth tides (SETID), ocean tides (OTID), and relativistic cor-Geopotential -acceleration with respect to its degree The graphs in Figure B.4 show the accelerations produced by the spherical harmonic terms of the geopotential model EGM96 summed over the orders for a given degree. The individual curves correspond to the altitude of a satellite in a circular orbit around the Earth. The histograms in Figure B.3 set the upper limit of the nongravitational accelerations acting on the GRACE satellites to be 500/30/70 nm s -2 for the alongtrack/cross-track/radial components, while the altitude of the satellites decreased from 510 km to 450 km. From Figure B.4 we infer that the geopotential-induced

  accelerations, shifted and scaled by given values. We will look for a linear filter that would realize the second derivative of positions, taking into account the character of the waveforms in question (Figs. B.1 and B.2).

  pecially for modelling the accelerations due to the geopotential, where different geopotential models might give different POD-based nongravitational signals. It is an assumption of the presented method that the noise in the modelled gravitational accelerations is negligible compared to that of the POD-based total acceleration (more on this point in Sec. 5.4). Simulated positions are computed by the numerical integration of the satellite motion using the simulated gravitational a SIM GRAV and nongravitational a SIM NG accelerations (SIM stands for 'simulated' or 'modelled'). The time step of positions and other quantities used in this study is 10 seconds. To these approximately error-free positions, which are given in the celestial reference frame, we added a normally distributed white noise Z, with a variance of σ 2 =1 cm in each position component. The resulting sequence of random vectors r represents the kinematic positions from the POD.

  The axes of the science reference frame are close to those of the satellite local reference Page 13 of 64 A c c e p t e d M a n u s c r i p t frame (Sec. 2.2) to within a few degrees, except for the sign. In this section, we use the exact satellite local reference frame, in Section 4, where the attitude information of the GRACE satellites is used, we perform a simple sign change to have all our calculations and figures in an approximate satellite local reference frame. The motivation for using the satellite local reference frame lies in its clear physical meaning, e.g. the air drag vector always points in the negative alongtrack direction, the terrestrial infrared radiation in the positive radial direction. In Figure B.5 the components of the POD-based nongravitational acceleration vector a POD NG in the satellite local reference frame are shown. Using the second-

  3.1), after filtering becomes a correlated random component of a POD NG . The OLS applied to the calibration equation (3) now enables one to calculate acceptable estimates of B and S , as the point estimates of the regression parameters are usually not much affected by the autocorrelated errors, but it is not possible to correctly estimate the uncertainty Page 15 of 64 A c c e p t e d M a n u s c r i p t of the calibrated accelerations. For a correct estimation of the uncertainties in B

  celerations a POD NG are in Figure B.6; only the solution in the along-track component is shown. As the GLS transformation matrix W=F -1 is actually the inverse A c c e p t e d M a n u s c r i p t to the second-derivative filter, which produces accelerations from positions, the "nongravitational positions" are obtained as a sort of double integral of a POD NG . Effectively, we got back into the positions, but now with the gravitational signal removed. In the upper panel of Figure B.6, the nongravitational positions are shown (y OLS 1 ) as the observations for the OLS estimates, and the fitted function (ŷ OLS 1 ), which is the simulated nongravitational acceleration a SIM NG transformed to positions by W. Several statistics shown in the lower panels confirm the fact that the OLS residuals in the middle panel are uncorrelated normal: autocorrelation function (ACF), partial autocorrelation function (PACF; more about it in Sec. 4.2), normal probability plot and Jarque-Bera test (e.g., Brockwell and Davis, 2002). Through the OLS applied to the transformed linear model (Eq. B.3), apart from the estimates of the calibration parameters b0 and b1 , the original error variance of the nongravitational positions (Sec. 3.1) is estimated by the OLS residual mean square σ2 (labelled as σ iid,est in Fig. B.6).

  Figure B.6 the reader may have noticed that the coefficient of correlation between the fitted calibration parameters b0 and b1 is very close to one, typically, when calibrating the simulated or real accelerometer data, we get ρ( b0 , b1 )≃0.999. . . Of course, such a high correlation is not good for the stability of the fitted parameters. The cause of this situation lies in the collinearity of the predictor variables, Page 17 of 64 A c c e p t e d M a n u s c r i p t one of the standard problems encountered in multiple regression (e.g., Chatterjee

  For simplicity, let us use for the calibration equation (3) the notation of the OLS from Appendix A and calibrate the accelerometer measurements against the simulated nongravitational accelerations, so in this subsection x≡a UNCAL ACC and y≡a SIM NG . We may approximately take both x (Fig. B.2) and y (Fig. B.1) as signals made up by two components, by a constant signal plus an oscillatory component (sum of sinusoids). This is not very far from the truth, as the patterns of one revolution in Figures B.1 and B.2 repeat themselves relatively regularly during a period of weeks or so. From the point of view of Fourier analysis, the constant component x and the oscillatory component (xx) are orthogonal to each other, the same applies to ȳ and (yȳ), so comparing the constants x, ȳ would produce an estimate of an 'intuitive' bias, i.e. a distance between the mean values x and ȳ, and fitting the oscillations (xx) and (yȳ) would estimate the 'scale factor',

  ) Taking the along-track component of a UNCAL ACC in Figure B.2 as a quantitative ex-A c c e p t e d M a n u s c r i p t ample, the power of the constant component x2 ≃(10 -6

  together with the definitions y ⋆ =yx and x ⋆ =xx. The modified model has perfectly uncorrelated parameters b ⋆ 0 and b ⋆ 1 , moreover, one can easily show that b ⋆ 0 =ȳ-x is the 'intuitive' bias mentioned above. The scale factors b 1 , b ⋆ 1 of both models have the same fitted value, b⋆ 1 = b1 , and, perhaps surprisingly, also the same standard error, σ( b1 )= σ( b⋆ 1 ). Only the modified intercept b ⋆ 0 has a substantially smaller standard error, from (A.4), σ( b⋆ 0 )= σ/ √ n. Indeed, the calculated values of the modified intercept b⋆ 0 are much less noisy compared to those of b0 . But on rearranging the terms in (6), y = x + b ⋆ 0b ⋆ 1 x + b ⋆ 1 x + ǫ, one can express the 'old' calibration parameters b 0 and b 1 by means of the modified ones,

  ). During the inversion of the normal equations,MATLAB (2007) indi-condition number, which was caused by a difference of several orders between the magnitudes of the two predictors; a simple solution was to multiply the intercept b ⋆ 0 by 10 -7 . In fact, both these computational modifications are analogous to standardizing the predictor variables in multiple regression or using the MATLAB option 'center and scale X data'. For the sake of comparison of our calibration parameters with those computed by other groups, and because, after all, the calibration models (3) and (6) are equivalent, the final results are given in terms of the original parameters b 0 and b 1 .

  ) using the second-derivative filter (1) and the modelled accelerations of gravitational origin. Figure B.7 shows a typical result for the three accelerometer axes, the amplified noise from the POD positions being roughly of the same order of magnitude as that for the simulated case in Figure B.5. The components shown in Figure B.7 are not exactly 'along-track', 'cross-track' and 'radial', as the accelerometer readouts are now given in the science reference frame (Sec. 3.3).

Figure 7

 7 Figure 7 should be positioned here.

Figure 8

 8 Figure 8 should be positioned here.

  Fig. B.8) are a realization of a stationary process and we want to represent its correlation structure by fitting an appropriate autoregressive moving-average (ARMA) model. This class of linear time series models has the property that any autocovariance function that asymptotically tends to zero can be approximated arbitrarily well by the autocovariance function of some ARMA process. The fact that the sample autocorrelation function (ACF) is negligible for some finite lag q suggests that a moving-average model MA(q) might provide a good representation of the data. Analogously, the partial autocorrelation function (PACF; in cyan, bottom left panel of Fig. B.8) of a causal autoregressive process AR(p) is zero for lags greater than p. Both the ACF and PACF of the OLS residuals are in the bottom left panel of Figure B.8. The sample PACF clearly falling off, we chose the pure AR(7) process to be fitted to the residuals using the Yule-Walker estimation. The ACF of the fitted AR process of order 7 (in green, bottom left panel of Fig. B.8) agrees A c c e p t e d M a n u s c r i p t

  Fig. B.9) for the along-track component permits us to use statistical inference and to assert that the 'true' signal measured by the accelerom-

  Figure B.10. The uncertainty band is wider when the fitted value is farther from the mean, similarly to the usual OLS model (A.7). The sample mean 3 σ(a CAL ACC ) , which we can use to characterize the obtained uncertainty band in the along-track component, is around 25 nm s -2 . In the same way, we can use the calibration equation to fit the simulated nongravitational accelerations and obtain a SIM,CAL NG . As is apparent from Figure B.10, the uncertainty bands of both σ(a CAL ACC ) and σ(a SIM,CAL NG ) are of similar size. But the calibration equation (3) was used in a usual OLS sense, however, after the GLS 2 transformation W 2 was applied. In the bottom panel of Figure B.10 there are the calibrated accelerometer and simulated nongravitational accelerations with their means subtracted and then projected to the W 2 space. It is evident that the W 2 transformation matrix is an integrator, which, inversely to the second-derivative filter (1), effectively filters out the high frequencies from both a CAL ACC and a SIM,CAL NG .Indeed, the estimated frequency response of the filter W 2 shows that only sinusoids of periods longer than 30 minutes are retained. Although the accelerometer waveform give more details in the 'acceleration domain' than the modelled nongravitational accelerations, the calibration in the GLS-induced nongravitational

  scheme is based on the assumption that the calibration parameters vary slowly in time. As the accelerometer data as well the POD positions contain relatively frequent portions of outliers (cf.[START_REF] Flury | Precise accelerometry onboard the GRACE gravity field satellite mission[END_REF], we used a running window covering several satellite revolutions, within which we calibrated the accelerometer readouts. From these calibration results we selected the non-overlapping seg-

  Fig. B.8), in physics and engineering this 'RMS value' σ is widely used to characterize the power of the residual signal. The upper panel of Figure B.11 shows the standard error of the fit σOLS 1 for the 1.5-year period. The label OLS 1 refers to the case, where the GLS transformation is based only on the inverse second derivative filter W=F -1 , and thus the accelerations a UNCAL ACC and a SIM NG are 'integrated' to give the 'nongravitational positions' (Sec. 3.6). This is interesting, because on supposing that the modelled gravitational accelerations have negligible errors, σOLS 1 then estimates the RMS value of the POD positions when compared with the independently measured accelerometer data. The figure shows that the empirical distributions of σOLS 1 for both a CAL ACC and a SIM,CAL NG are very close, with no statistically significant difference, their mean values being equal to around 3 cm with an approximate uncertainty of 1-2 cm. The aim of this paper is to obtain the calibrated accelerometer data together with a realistic error bar. As mentioned in Section 4.3, this can be achieved in A c c e p t e d M a n u s c r i p t the along-track component only. The uncertainty estimates of the calibrated accelerometer and simulated nongravitational accelerations σ(a CAL ACC ) and σ(a SIM,CAL NG ) in the lower panel of Figure B.11 are again statistically equivalent, the mean uncertainty being 8.5±3.0 nm s -2 . This is due to the severe smoothing, when the accelerations a CAL ACC and a SIM,CAL NG are calibrated against the POD positions, as explained in Section 4.3.

Figure 11

 11 Figure 11 should be positioned here.

(

  Sec. 3.8), we can set the scale factors b1 equal to the values specified in Bettadpur (2004a) and expect that the biases will 'adapt' their values accordingly. Indeed, in the three accelerometer axes of GRACE A, Figure B.14 shows a similar time evolution of our biases and those from the report. In this case, the fixed values of the scale factors were 0.961 (along-track), 0.98 (cross-track), 0.94 (radial). We obtained similar results for GRACE B, Figure B.15, for the fixed scale factors 0.947 (along-track), 0.97 (cross-track), 0.92 (radial).

Figure 14

 14 Figure 14 should be positioned here.

Figure

  Figure 15 should be positioned here.

  mates of the calibration parameters in all three accelerometer components for Page 33 of 64 A c c e p t e d M a n u s c r i p t both GRACE A/B satellites. A statistically correct estimate of the accuracy of the calibrated accelerometer measurements have been obtained for the along-track component of the accelerometer data.

  original linear model into a new one, y * = X * b + ǫ * , (B.3)such that the covariance matrix of the transformed errors ǫ * is again a scaled identity matrix. Indeed,Var(ǫ * ) = WVar(ǫ)W ′ = σ 2 T -1 T T ′ T -1′ = σ 2 1. (B.4)In the transformed variables, y * =Wy, X * =WX, the usual OLS are then used to find the regression parameters b of the original problem with correct estimates of their uncertainties. The thus obtained GLS estimator b is also known as the Aitken estimator. By using the inverse transformation matrix T =W -1 , we may obtain the confidence and prediction intervals of the original fitted function ŷ=T ŷ * from (A.7) and (A.8). Namely, the estimated confidence interval of ŷ is expressed in matrix notation as the square root of the diagonal of the covariance matrix Var(ŷ) = σ2 T P * T ′ , (B.5)

Figure 16

 16 Figure 16 should be positioned here.

Figure 17

 17 Figure 17 should be positioned here.

Fig

  Fig. B.17).

Figure B. 1 :

 1 Figure captions

Figure B. 2 :

 2 Figure B.2: Uncalibrated accelerometer data a UNCAL ACC (the same arc as in Fig. B.1).

Figure B. 3 :

 3 Figure B.3: Histograms of gravitational and nongravitational accelerations in the satellite local reference frame components (GRACE A, 08/2002-03/2004).

Figure B. 4 :

 4 Figure B.4: Acceleration due to the spherical harmonic terms of the gravitational model EGM96 grouped according to the degree.

Figure B. 5 :

 5 Figure B.5: The POD-based nongravitational accelerations a POD NG in the satellite local reference frame (derived from the simulated POD positions). Also shown are the simulated nongravitational accelerations a SIM NG .

Figure B. 8 :

 8 Figure B.8: The ordinary least squares applied to "nongravitational positions" (panels as in Fig. B.6). Real data used (GRACE A, 25 Nov 2003, along-track).

Figure B. 9 :

 9 Figure B.9: The ordinary least squares applied to the transformed residuals from Fig. B.8, the transformation matrix is based on the fitted AR(7) process.

Figure

  Figure B.10: Calibrated accelerometer readouts and simulated nongravitational accelerations (upper panel), after centring and the transformation given by W 2 (bottom panel) (GRACE A, 25 Nov 2003, along-track).

Figure B. 11 :

 11 Figure B.11: Time evolution of the standard fit error for the nongravitational positions (upper panels) and the uncertainty of the calibrated accelerations (lower panels) compared for the accelerometer-based and simulated nongravitational accelerations (GRACE A, along-track, window of 2 revs., approx. 2000 values).

Figure

  Figure B.12: Long-term fit of the calibration parameters for the accelerometer measurements (GRACE A, along-track, window of 2 revolutions).

Figure B. 13 :

 13 Figure B.13: Long-term fit of the calibration parameters for the simulated nongravitational accelerations (GRACE A, along-track, window of 2 revolutions).

Figure B. 14 :

 14 Figure B.14: Comparison of the computed bias for GRACE A with that derived independently by Bettadpur (2004a).

Figure B. 15 :

 15 Figure B.15: Comparison of the computed bias for GRACE B with that derived independently by Bettadpur (2004a).

  Figure B.16: Example of a linear model with the errors generated by a stationary AR(7) process: the direct ordinary least squares solution. Upper panel: ŷCI define the confidence interval around ŷ, ŷPI the prediction interval; middle panel: standardized residuals and the fit results; lower panel: autocorrelation function of residuals.

Figure B. 17 :

 17 Figure B.17: Data as in Fig. B.16: the generalized least squares solution.

  Figure14

  Figure15

A c c e p t e d M a n u s c r i p t

We applied the calibration method to the accelerometer data covering a 1. 5year period in 2002-2004. Taking into account the previous experience, we supposed that the calibration parameters, i.e. the scale factors and biases for each accelerometer axis, evolve slowly in time. We used the running window of 2-4 orbital revolutions, within which we calibrated the accelerometer data and finally selected the non-overlapping segments with the best statistical results. The time evolution of the calibration parameters agrees well with that published in an independent report.

On the assumption that the errors in the modelled accelerations of gravitational origin are very small, the GLS calibration method defines a transformation of the accelerometer data, which may be used to estimate the RMS of noise in the kinematic positions. Based on this comparison between the POD kinematic positions with the independently measured accelerometer data set, we found plausible mean values of (3-4; 4-7; 6-12) cm in the (along-track; cross-track; radial) directions.

We compared the long-term calibration results for several models of the Earth static gravity field. The recent models EIGEN-5, EGM08 and GGM03, which are based also on the data from the CHAMP and GRACE missions, gave statistically equivalent results, the mean uncertainty in the along-track component of the calibrated accelerometer data being 6.5-8.5 nm s -2 (one sigma). The same long-term results were also obtained using the EIGEN-CHAMP03 model, which does not contain the GRACE data. The estimated statistical errors produced using the pre-CHAMP gravity models were several times worse in the cross-track and radial components.

Page 35 of 64 A c c e p t e d M a n u s c r i p t

Acknowledgements

The author would like to thank Dražen Švehla from TU Munich for providing the GRACE kinematic orbits. The ISDC online data centre of the GFZ is acknowledged for the GRACE data products used in this study (http://isdc.gfz-potsdam.de/ grace/), the ICGEM centre of the IAG for the coefficients of the global gravity field models (http://icgem.gfz-potsdam.de/ICGEM/). Thanks are also due to Jaroslav Klokočník and Radek Peřestý for useful discussions, and to the reviewers for their relevant comments leading to many improvements in the manuscript.

Appendix A. Ordinary least squares (OLS)

In the ordinary least squares we suppose that the vector of observations y is given as the sum of a deterministic mean function E(y) to which a vector of random errors ǫ with constant variance is added. The probability model of the simple linear regression is

where y is the vector of n observations, b 0 intercept, b 1 slope, x predictor, ǫ statistical error. The OLS estimates b0 and b1 are given by b1 = SXY/SXX, b0 = ȳ -b1 x, (A.2)

where SXY= (x ix)(y iȳ), SXX= (x ix) 2 . Using b0 and b1 we form the fitted function ŷ as the estimate of the mean function E(y)

Under the assumption that the errors ǫ i are independent and normal with constant variance σ 2 , the OLS estimates b0 , b1 are also normally distributed with the stan-

where the standard error of the OLS fit is

In general, the estimated parameters are correlated with the coefficient of correlation (Weisberg, 2005)

We can calculate the uncertainty band around the fitted function ŷ, which is called the confidence interval,

and the prediction interval, the uncertainty of a single (possibly future) observa-

Appendix B. Generalized least squares (GLS)

Defining X as the matrix of predictors and b as the vector of parameters, let