Stochastic noise modelling of kinematic orbit positions in the Celestial Mechanics Approach
Martin Lasser
Ulrich Meyer
Daniel Arnold
Adrian Jäggi
Gravity field models may be derived from kinematic orbit positions of Low Earth Orbiting satellites equipped with onboard GPS (Global Positioning System) receivers. An accurate description of the stochastic behaviour of the kinematic positions plays a key role to calculate high quality gravity field solutions. In the Celestial Mechanics Approach (CMA) kinematic positions are used as pseudoobservations to estimate orbit parameters and gravity field coefficients simultaneously. So far, a simplified stochastic model based on epochwise covariance information, which may be efficiently derived in the kinematic point positioning process, has been applied.
We extend this model by using the fully populated covariance matrix, covering correlations over 50 min. As white noise is generally assumed for the original GPS carrier phase observations, this purely formal variance propagation cannot describe the full noise characteristics introduced by the original observations. Therefore, we sophisticate our model by deriving empirical covariances from the residuals of an orbit fit of the kinematic positions.
We process GRACE (Gravity Recovery And Climate Experiment) GPS data of April 2007 to derive gravity field solutions up to degree and order 70. Two different orbit parametrisations, a purely dynamic orbit and a reduceddynamic orbit with constrained piecewise constant accelerations, are adopted. The resulting gravity fields are solved on a monthly basis using daily orbital arcs. Extending the stochastic model from utilising epochwise covariance information to an empirical model, leads to a – expressed in terms of formal errors – more realistic gravity field solution.
 Article
(683 KB)  Fulltext XML
 BibTeX
 EndNote
Various Low Earth Orbiting (LEO) satellites are equipped with a GNSS (Global Navigation Satellite System) receiver, which may not only be used for navigational purposes and Precise Orbit Determination (POD) but also for gravity field recovery. Since the CHAllenging Minisatellite Payload mission (CHAMP, Reigber et al., 1998) was launched in 2000, further missions dedicated to the determination of the Earth's gravity field were brought into space, such as the Gravity Recovery And Climate Experiment (GRACE, Tapley et al., 2004), the Gravity field and steadystate Ocean Circulation Explorer (GOCE, Drinkwater et al., 2006), Swarm (FriisChristensen et al., 2006), and recently GRACE FollowOn (Flechtner et al., 2013). All of them were equipped with a Global Positioning System (GPS) receiver to support the determination of the long wavelength part of the Earth's gravity field. Kinematic orbit positions are obtained in a precise point positioning process (Švehla and Rothacher, 2005). They are widely used as pseudoobservations in the context of gravity field determination, first demonstrated by Gerlach et al. (2003). In particular in the context of bridging the almost one year gap between GRACE and GRACE FollowOn, gravity field determination using kinematic orbit positions became an essential topic, e.g., Lück et al. (2018) used Swarm kinematic positions to derive timevariable gravity fields and ocean mass changes.
We focus on kinematic positions derived from GRACE GPS data to reconstruct orbits and estimate gravity fields simultaneously by applying the Celestial Mechanics Approach (CMA, Beutler et al., 2010). Studies based on the CMA showed its performance in gravity field determination using kinematic positions, e.g., Prange (2010) or Prange et al. (2009) using CHAMP data, Jäggi et al. (2011a) for CHAMP, GRACE and GOCE, Jäggi et al. (2015) for GOCE or Jäggi et al. (2016) for Swarm data. An accurate description of the stochastic behaviour of the kinematic positions plays a key role to calculate highquality gravity field solutions, not only for its own sake but also for further processing, e.g., a realistic handling of the stochastic noise would enable a rigorous combination of gravity solutions produced by different analysis centres on the level of normal equations (Meyer et al., 2019).
So far, a simplified stochastic model based on epochwise covariance information, derived in the kinematic point positioning process, has been applied in the CMA. We extend this model by using the fully populated covariance matrix as Jäggi et al. (2011b) did, covering correlations over 50 min, and refine it by the use of empirical covariances based on the residuals of an orbit fit of the kinematic positions. In order to assess the mechanism between stochastic models and the orbit parametrisation, we use purely dynamical orbits (six Keplerian elements and additional accelerometer calibration parameters) on one hand and reduceddynamic orbits with additional constrained piecewise constant accelerations set up within intervals of 15 min on the other hand. We take the opportunity GRACE provides and make use of the Kband ranging system between the two satellites to independently validate the reconstructed orbits with highest precision.
The paper is organised in six sections, where Sect. 2 briefly introduces the data and Sect. 3 the methods of gravity field recovery used in this study. In Sect. 4 stochastic noise modelling for kinematic orbit positions is explained together with the results for the different covariance matrices we have employed. First we discuss the simplified stochastic model based on epochwise covariances (Sect. 4.1.1), next the extension to covariances overarching epochs (Sect. 4.1.2) and finally we adopt empirical covariances (Sect. 4.2). Section 5 gives an outline to the stochastic behaviour of undifferenced ambiguity fixed kinematic positions. Section 6 concludes the results of this study.
In order to test our methods of stochastic modelling, we use GRACE data collected in April 2007. We deliberately chose data of good quality to avoid issues with screening methods for a reliable derivation of empirical covariances. Kinematic positions are obtained from the GPS carrier phase observations using the Bernese GNSS Software (Dach et al., 2015). The procedure of kinematic point positioning follows Švehla and Rothacher (2005) using the GPS orbits and satellite clock corrections generated at the Center for Orbit Determination in Europe (CODE, Dach et al., 2009; Bock et al., 2009). The CODE products (Sušnik et al., 2020) are also generated with the Bernese GNSS software, thus providing full consistency. Biased Kband range data are used as independent observations to validate the distance changes between the two GRACE satellites derived from the reconstructed orbits.
In the CMA gravity field recovery is based on the solution of a generalised orbit determination problem for each GRACE satellite (e.g., Jäggi et al., 2011a). We split the approach into two steps. In the first step the kinematic positions together with covariance information are introduced as weighted pseudoobservations and daily arcs of a priori orbits are fitted according to standard leastsquares adjustment by numerically integrating the equation of motion in the force field defined by the models listed in Table 1. This orbit reconstruction process leads to an a priori orbit for each arc.
^{a} Pavlis et al. (2012). ^{b} Savcenko and Bosch (2011). ^{c} Dobslaw et al. (2017). ^{d} Petit and Luzum (2010). ^{e} Folkner et al. (2009). ^{f} Sušnik et al. (2016) and Sušnik et al. (2020).
In the second step the a priori orbits are introduced to set up normal equations along the orbit by using the same orbit parametrisation, and additionally, spherical harmonic (SH) coefficients representing the gravity field. The kinematic positions again serve as pseudoobservations and are weighted according to the covariance information. All daily normal equations are accumulated into one normal equation system covering one month, which is eventually solved for in order to obtain a solution for the correction of the orbit parameters and the SH coefficients. No regularisation is applied to compute the gravity field parameters.
To investigate the robustness and the dependencies of the approach, two types of orbit parametrisations were adopted for each type of covariance in orbit reconstruction and gravity field recovery:

A dynamic orbit represented by only six initial osculating (Keplarian) elements and nine additional parameters characterising the accelerometer, and

a reduceddynamic orbit represented by the same parameters as the dynamic orbit and additional constrained Piecewise Constant Accelerations (PCA) set up over intervals of 15 min in radial, alongtrack and crosstrack direction.
The nine additional accelerometer parameters consist of a bias in radial and crosstrack direction, a polynomial of degree 3 in alongtrack and a scaling factor for each axis. The optimal constraining and sampling of the PCAs for the reduceddynamic orbit are determined empirically. The constraints are set to $\mathrm{1}\times {\mathrm{10}}^{\mathrm{8}}$ m s^{−2}, the sampling is set to 15 min.
As the CMA is only one of several approaches to recover gravity field information from kinematic positions, we refer to Baur et al. (2014) for an overview on different strategies of gravity field determination, including a brief description of stochastic modelling.
The stochastic noise modelling of kinematic orbit positions is a crucial part for highprecision gravity field recovery. The way the noise in the pseudoobservations is treated influences the quality of the result (e.g., in terms of formal errors) significantly, and it might also affect the solution itself. The correct characterisation of noise in the data retains the full signal content and separates signal from noise through a modelling of the latter – provided that the signal component is also adequately modelled.
In case of kinematic positions, which are introduced as pseudoobservations in the CMA, one has to consider the process of calculating these pseudoobservations in the stochastic modelling. The original observations are the GPS dualfrequency carrier phases. The kinematic positions of the LEO satellite are determined at the observation epochs of the GPS carrier phase measurements by a precise point positioning (PPP) approach (Zumberge et al., 1997) where three dimensional coordinates and the receiver clock corrections are estimated. GPS satellite orbits and clocks are introduced as error free, thus, correlations due to the GPS satellite orbits, clocks and underlying network errors are neglected as it is commonly done in most studies. However, the carrier phase ambiguities, which represent the integer number of full cycles between the GPS satellites and the receiver, and additional satellite and receiverspecific biases, correlate the kinematic positions in time, “implying that degraded estimates of carrier phase ambiguities may affect the estimation of the kinematic positions for more than one hour”. (Jäggi et al., 2011b). Consequently, the positions are no longer uncorrelated and the correlation depends on the constellation of the GPS satellites. Initially, carrier phase ambiguities are estimated as float numbers together with the kinematic coordinates and the receiver clock corrections in the leastsquares adjustment. A fixing between satellite differences to the correct integer number, requires in the undifferenced GPS processing the knowledge of satellitespecific bias terms (Schaer et al., 2020). In Sect. 5 an outlook is given on the effect of integerfixed positions on the stochastic behaviour of the kinematic pseudoobservations.
The kinematic positions are introduced as pseudoobservations in a standard leastsquares adjustment to estimate orbit parameters and gravity field parameters in a joint adjustment. The leastsquares adjustment yields
where x is the solution of the linearised system of equations, composed of orbit and gravity field parameters. l denotes the Observed minus Computed (O−C) vector, where the pseudoobservations form the observed vector. A being the first designmatrix and P the weight matrix, which is supposed to describe the stochastic behaviour of the pseudoobservations.
4.1 Formal variance propagation
The process of calculating the pseudoobservations from the original carrier phase measurements allows for a covariance estimation. We derive the covariance matrix ${\mathbf{C}}_{ll}^{\mathrm{formal}}$ from the inverse normal equation matrix, which is the result of a formal variance propagation from the phase observations to the unknown kinematic positions according to
where
Q_{kk} denotes the cofactor matrix of the original phase observations and ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$ is the a priori variance of unit weight. The relation between cofactor and covariance matrix can be established by $\mathbf{C}={\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}\mathbf{Q}$. White noise is generally assumed for C_{kk}. The propagation matrix R stems from the kinematic point positioning with A being the first designmatrix and P the weight matrix derived from ${\mathbf{Q}}_{kk}^{\mathrm{1}}$.
The correlations in ${\mathbf{C}}_{ll}^{\mathrm{formal}}$ feature a twiceperrevolution signal, which is related to the satellite crossings of the pole and the weaker observation geometry there since GPS satellites revolve in inclined orbits at 55^{∘} with respect to the Earth's equator. The better tracking geometry in the equatorial regions leads to a smaller number of interruptions, thus fewer ambiguities, and consequently, to a better determination of each ambiguity and more correlated kinematic positions (Jäggi et al., 2011b). Generally, the correlations are positive and decreasing for a growing distance in time, however jumps may occur due to the set up of new ambiguities, which indicates changes in the observed constellation or data problems, e.g., cycle slips.
The GPS carrier phase ambiguities are the only parameters in A in Eqs. (2) and (3) connecting different epochs when using GPS phase data. Consequently, deficiencies in the modelling of the GPS phase observations may be propagated through the ambiguities over several epochs, which is also reflected by ${\mathbf{C}}_{ll}^{\mathrm{formal}}$. Since ${\mathbf{C}}_{ll}^{\mathrm{formal}}$ only depends on the observation scenario but not on the actual observations (since white noise is assumed for them), any degradation of positions due to GPS data quality issues (including issues with GPS orbits and clocks) is not reflected in this type of covariance information. An outlier in the carrier phase observations, e.g., affects the individual kinematic position and through the ambiguities neighbouring positions as well, but not the formal variance propagation.
The weight matrix P for the leastsquares adjustment of the gravity field and orbit parameters (Eq. 1) is determined from the inverse covariance matrix and the a priori variance of unit weight ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$ as
The a priori variance of unit weight ${\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}$ is chosen based on the expectation about the uncertainty of the original carrier phase observations and is usually set between 1 and 2 mm for GPS L1 and L2 data. However, one is free to chose an arbitrary number, limits may only be given by the numerical precision of computers.
Introducing the pseudoobservations together with the full covariance matrix ${\mathbf{C}}_{ll}^{\mathrm{formal}}$ in the leastsquares adjustment for the gravity field recovery would be equivalent to starting directly with carrier phase observations (shown for the application of carrier phase observations and kinematic positions as pseudoobservations in Jäggi et al., 2011b). Facing the large number of kinematic pseudoobservations, e.g., 3×8640 for one day with a 10 s sampling of GPS data, this would require a huge computational effort and would demand high storage requirements. Two ways how to deal with these issues are presented in the next subsections.
4.1.1 Epochwise covariance information
The epochwise covariance information is a subset of the covariance matrix from the formal variance propagation, which may be easily derived in the kinematic point positioning process using preelimination and backsubstitution techniques. Focusing only on the kinematic positions, the epochwise covariances, assembled in Eq. (5), consist of six distinct elements, which comprise the variances of the three coordinates (${c}_{xx},{c}_{yy},{c}_{zz}>\mathrm{0}$) and the offdiagonal elements ${c}_{xy},{c}_{xz},{c}_{yz}$, which contain information about the correlation between the x,y,z axes in the Earthfixed coordinate system,
Any correlations between subsequent epochs are neglected. The epochwise covariances are only little memory consuming and the 3×3matrix is invertible to derive the (epochwise) weight matrices according to Eq. (4).
The epochwise covariance information mainly features the twiceperrevolution behaviour for polar orbiting satellites (orbital period for GRACE was T=96 min) which reflects the observation scenario. Figure 1 displays all elements of the epochwise covariance matrix for GRACEA of day 091, 2007 in the time and frequency domain. The elements are rotated to a Local Orbit Frame (LOF) based on the unit vectors derived from the reduceddynamic orbit, the axes follow radial (r), alongtrack (a) and crosstrack (c) direction. The twiceperrevolution signal is clearly visible in the radial component at $T/\mathrm{2}=\mathrm{48}$ min where T corresponds to the revolution period of the GRACE satellites.
This information can only be derived from the formal variance propagation and cannot be neglected in the orbit reconstruction and gravity field recovery process without significantly degrading the recovered gravity field solutions using the classical CMA, see Prange (2010).
Impact of epochwise covariance information on orbit reconstruction and gravity field recovery
The residuals of the reconstructed orbit (Fig. 2) with respect to the kinematic positions are characterised by the underlying orbit parametrisation. The reduceddynamic orbit residuals in Fig. 2 (top) show a scatter around zero with few systematics where longperiodic effects related to a deficient force modelling are absorbed by the PCA set up every 15 min. Thus, the reduceddynamic orbit follows the observations to a large extent. In contrast, the residuals of the dynamic orbit shown in Fig. 2 (middle) feature strong systematic signals caused by the deficient force model, which are mainly visible as onceperrevolution signatures.
The Kband validation of the fitted dynamic orbit is smoother than of the fitted reduceddynamic orbit (see Table 2). This reflects that the PCA allow the orbit to follow the excursions in the observations to a certain extent, whereas, the dynamic orbit follows – based on the leastsquares adjustment – the a priori force field. The small highfrequency noise in the Kband validation, visible, e.g., in the peaks of the red curve, can be attributed to the high degrees of the gravity field, which are difficult to estimate using only GPS observations. The Kband validation of the reduceddynamic orbit does not reflect the same level of accuracy as using directly GPS carrier phase observations for orbit determination, which is accounted for by the insufficient stochastic model, see Jäggi et al. (2011b) for a detailed investigation. This means the fitted reduceddynamic orbit is less precise when being computed from kinematic positions with epochwise covariance information than a corresponding orbit being computed directly from GPS carrier phase observations, where no degradation due to neglected correlations is occurring.
The gravity field solutions based on the epochwise covariance information when using the reduceddynamic and the dynamic orbit parametrisations are shown in Fig. 3. The solution based on the dynamic orbit parametrisation is degraded because of the deficient force model. The gravity field solution based on the reduceddynamic parametrisation (red) is of good quality and represents the classical parametrisation used so far in the context of the CMA. In all further degree amplitude figures it will serve as a reference and will always be indicated by red lines. Any further reference such as “classical approach” refers to the use of epochwise covariance information.
In the dynamic and the reduceddynamic case the formal errors (dashed lines) of the solution do not reflect the accuracy assessed by the differences to the (superior) GOCO05s (MayerGürr et al., 2015) gravity field model (solid lines). Hence, the stochastic description of the pseudoobservations needs to be expanded.
4.1.2 Covariance information overarching epochs
Apart from the epochwise covariance information from the formal variance propagation, the covariances connecting a certain number of epochs may be taken into account as well. As discussed at the beginning of Sect. 4.1, their behaviour reflects the presence of the carrier phase ambiguities as solely these parameters connect different epochs in a kinematic point positioning. They comprise the observation geometry, thus the influence of the constellation on the correlation between epochs.
The computation of the covariances is in principle straightforward, however, depending on the number of correlated epochs that are of interest, memory consumption may become an issue (e.g., for one day of 10 s kinematic positions and correlations up to two hours >53 Million elements need to be stored, which consumes in a 64 bit notation about 500 MB memory). The variance/covariance matrix is fully populated and nondegenerate if all epochs are considered. However, when taking only correlations for a certain time interval into account, the matrix is banddiagonal and not necessarily invertible anymore. A simple approach of obtaining an invertible matrix is slicing fully populated blocks of the matrix and inverting these blocks individually (the block diagonal matrix is invertible, cf. Jäggi et al., 2011b). The block length defines the number of epochs that are considered. Correlations between individual blocks are neglected, however, when applying a sufficiently long correlation time (until the correlation becomes negligible) the impact on the solution is insignificant. We use the covariance matrix up to 50 min, which is about twice the time the observed constellation takes to vary completely (the maximum time interval an individual ambiguity is valid is around 25 to 30 min, hence, only ambiguities overlapping each other propagate covariance information for longer time spans with decreasing influence the longer the distance in time is). Figure 4 shows all individual elements of the covariance function for the first 50 min of the day 091, 2007 in Earthfixed frame (panel a) and for a comparison with the subsequent graphs of covariance functions also rotated to a LOF based on a fit using the reduceddynamic orbit parametrisation (panel b). The largest variance (for c_{zz} and c_{rr}) slightly exceeds 50 mm^{2}, but for reasons of clarity it is cut off in Fig. 4. However, this number coincides with the fact that the radial component is about three times less precise than the other two axes (visible in Fig. 1).
The jumps in the covariance function occur due to the setup of new ambiguities (changes in the observed GPS constellation).
Impact of covariances over 50 min on orbit reconstruction and gravity field recovery
For both, the dynamic and the reduceddynamic orbit, covariances covering correlations up to 50 min were applied to determine the weight matrix in Eq. (4). The residuals of the reconstructed orbits with respect to the kinematic pseudoobservations in Fig. 5 exhibit a slightly different behaviour compared to the classical approach shown in Fig. 2. The dynamic orbit is still subject to the deficient a priori force model, however longperiodic variations induced by the observation geometry and ambiguity setup of the original phase observations are no longer allowed to be absorbed by orbit parameters. That accounts for both parametrisations and is visible in particular in the reduceddynamic orbit residuals as superimposed oscillation (see Fig. 5, top).
The Kband range residuals for the reduceddynamic orbit are significantly lower than in the classical approach of weighting (cf. Tables 3 and 2) since longperiodic variations of the pseudoobservations are no longer (erroneously) fitted by the parameters of the orbit model but (correctly) interpreted as a consequence of the ambiguityinduced correlations in time (Jäggi et al., 2011b). In contrast to the epochwise covariance weighting, the Kband range validation is now at the same level as using directly GPS carrier phase observations for orbit determination. Experiments substantiating this behaviour are presented in Jäggi et al. (2011b).
Weighting the kinematic positions including 50 min of correlations leads to slightly more realistic formal errors in the gravity field recovery process (Fig. 6). The very low degrees, however, still reveal deficiencies in the stochastic modelling. The dynamic solution profits most from taking the 50 min of correlations into account (compare Fig. 6 to Fig. 3), however, it cannot fully compete with the reduceddynamic solution in the higher degrees. The very low degrees (best visible for degree 2) are better determined with the dynamic parametrisation as the PCA of the reduceddynamic parametrisation absorb some of the signal.
4.2 Empirical covariances
The residuals of an orbit fit with respect to the kinematic positions reflect all model (functional and stochastic) and data deficiencies. Consequently, deriving covariances from the residuals leads to a stochastic description of the entire physical system. The residuals e are obtained as the difference between the pseudoobservations l (i.e., the kinematic positions) and the estimated observations, which are computed from the estimated model parameters x with
As these residuals depend on the a priori force field, we first estimate a gravity field solution based on epochwise covariances and then reintroduce this solution as new a priori gravity field, to become independent from the a priori gravity field when computing the residuals for the empirical covariances. This procedure could be applied iteratively by taking the solution computed using empirical covariances as new a priori gravity field, however, in our case of starting with epochwise covariances as an initial stochastic model, it turned out that the gravity field solution does not alter after one iteration (not shown). We compute the residuals in a backsubstitution process, thus allowing for efficient solving of the normal equations by prior preelimination of the orbit parameters.
The covariance function for a certain time interval Δt_{k}, with k denoting the sampling interval of the pseudoobservations, is defined by the autocorrelation between the respective estimated residuals e according to
This type of autocorrelation is only legitimate under the assumption that the residuals in the considered interval fulfil the requirements of a (weak) stationary process (Etten, 2005). We assume the noise to be stationary, thus, the more parameters are set up to absorb deficiencies in the modelling, the better this assumption is fulfilled, which can be seen in the orbit residuals of the dynamic and reduceddynamic parametrisation. Moreover, the longer the correlation lengths are that are taken into account, the more coloured is the noise in the residuals (to be clearly seen in Figs. 2 and 5 for the reduceddynamic cases).
The correlations between the axes are obtained through a crosscorrelation between the residuals of the respective axes
Assembling a variance/covariance matrix from the covariance function leads to a block Toeplitz structure, which may be populated by
and
In case of correlating less than all epochs, the matrix is band diagonal. Note that only for Δt_{0} the relation for the offdiagonal elements of c_{xy}=c_{yx}, c_{xz}=c_{zx}, c_{yz}=c_{zy} holds, in case of $\mathrm{\Delta}{t}_{k},k\ge \mathrm{1}$ the covariance block cov(Δt_{k}) is composed of nine independent elements.
We do not want to model short term variations within an arc, therefore, we use the residuals of a whole month to derive a mean covariance function. The use of empirical covariances relies on the quality and the shape of the residuals, in particular outliers will affect the empirical covariances heavily. This leads to an interdependency between the orbit parametrisation and the pseudoobservations. The dependency on the a priori gravity field is mitigated as much as possible by first estimating an independent gravity field solution and obtaining the covariance function from postfit residuals. Consequently, the time until the correlations become negligible is highly dependent on these two factors. Even under the assumption that there are no outliers (a prerequisite) in the data, the parametrisation of the underlying orbit affects the magnitude and time the correlations take to vanish significantly, see Fig. 7 for the two parametrisations adopted in this study.
We introduce the epochwise covariance information from the formal variance propagation into the a priori gravity field and orbit recovery, on which the empirical covariance function is built, because the epochwise covariances transport information about the tracking scenario (inferior constellation over polar regions), which cannot be easily recovered by orbit dynamics only. However, introducing white noise as most simple assumption to weight the kinematic positions in the a priori gravity field and orbit recovery, leads to the same gravity field solution, but requires at least one additional iteration.
Memory consumption is not an issue for the usage of empirical covariances since they are valid over one month and are fully described by one function from which the variance/covariance matrix can be compiled according to Eq. (10).
We expect the nature of the residuals in a LOF to be more stationary, as periodic behaviours due to the satellites' orbit will be reflected much cleaner and easier to interpret, thus, we derive the empirical covariances in a LOF. For the reduceddynamic orbit the largest correlations over epochs occur in the radial direction, which can be attributed to the simultaneous estimation of kinematic positions and receiver clock corrections, while the other directions play a minor role. This is different for the dynamic parametrisation, where the correlations in all three axes are about the same magnitude. Noteworthy is the large and long lasting correlation in crosstrack and between the radial and the alongtrack axis. Thus, the three axis are not independent from each other, even over several epochs.
As we aim to reconstruct the weighting for the observations from estimated residuals e, the relation between the covariance matrix of the residuals C_{ee} and the covariance matrix of the observations is needed. This relation may be established by variance propagation starting from inserting Eq. (1) into Eq. (6)
Setting
the variance propagation for the covariance matrix of the estimated residuals follows
which gives by expansion of the products
A simple rearrangement yields
when deploying the residuals e as estimation for the covariance matrix of the observations C_{ll}. For the processing of kinematic positions with a precision of about 1 cm, the correctional term ${\mathbf{C}}_{\widehat{l}\widehat{l}}={\mathit{\sigma}}_{\mathrm{0}}^{\mathrm{2}}\left[\mathbf{A}({\mathbf{A}}^{\mathrm{T}}\mathbf{PA}{)}^{\mathrm{1}}{\mathbf{A}}^{\mathrm{T}}\right]$ will be negligible with respect to C_{ee} and we may set
even though the estimation of the covariance matrix of the observations is then not unbiased anymore. Figure 8 shows the correctional term ${\mathbf{C}}_{\widehat{l}\widehat{l}}$ for the first minutes of day 091, 2007 for the reduceddynamic (panel a) and dynamic (panel b) orbit parametrisation. The correctional term is given in a LOF and can be directly compared to the empirical covariance functions in Fig. 7. The first minutes of the arc feature the largest correctional term, which is located at the arc's boundaries, a consequence of the leastsquares adjustment. It gets smaller towards to the middle of the arc.
The magnitude of the corrections is much lower than that of the covariance function computed from the residuals (compare graphs in Figs. 7 and 8). The differences in magnitude may be explained by the different nature of C_{ee} and ${\mathbf{C}}_{\widehat{l}\widehat{l}}$. The correctional term is composed of a purely mathematical variance propagation, thus only reflecting formal errors, which uniquely depend on the feasibility of estimating the orbit parameters for the given observation scenario. In case of the dynamic orbit, the 15 parameters which are set up, are well determined with a low RMS and high redundancy, despite the fact that the overall fit exhibits deficiencies (see residuals in Fig. 2). In case of adding PCA to form the reduceddynamic parametrisation, the fit is significantly improved but the redundancy is slightly worse and so is the precision of each parameter, however, in both cases the variance propagation of the parameters' precision to estimated observations, reflects this high degree of redundancy.
The covariance function of the residuals on the other hand reflects all model and data deficiencies, and consequently, its magnitude is significantly larger. Neglecting the correctional term provides the benefit of obtaining a nondegenerate covariance matrix by the auto and crosscorrelations. This simplification might only be allowed for kinematic positions, but not for significantly more precise observables such as ultraprecise KBand observations.
The weight matrix is obtained through Eq. (4). Likewise to the formally propagated covariances, the band diagonal structure is broken up for the sake of simple computation to a blockwise treatment in the leastsquares problem (since a band diagonal version of a positive definite matrix is not necessarily positive definite), where one block is invertible and the blocks are independent from each other. The block length depends on the time interval of the autocorrelation of the residuals. The longer, the more computational effort needs to be carried out for the inversion. If the correlation is as such that only one block remains (e.g., one day for a daily arc), C_{ll} is fully populated and nondegenerate. The blockwise treatment disregards correlations between subsequent blocks including the large correlations between the old block's end and the new block's beginning. Correlations between subsequent blocks do not vanish, however, correlations within a block tend to zero.
Impact of empirical covariances on orbit reconstruction and gravity field recovery
For both orbit parametrisations a weight matrix derived from empirical covariances covering correlations up to 50 min was introduced in the leastsquares adjustment. The empirical covariances for the dynamic orbit were determined utilising the dynamic orbit parametrisation, the empirical covariances for the reduceddynamic orbit from a reduceddynamic a priori orbit fit. Even though the correlations in the dynamic case do not vanish after 50 min, this time interval is sufficiently long to describe the noise behaviour. The orbit residuals of the fit with empirical covariances and the Kband validation (Fig. 9) show almost the same behaviour as in the classical approach (Fig. 2), an indication that the empirical covariances do not introduce systematic behaviour into the physical system. The Kband RMS over one month is 1.8 cm for the reduceddynamic and 0.8 cm for the dynamic parametrisation.
Using the empirical covariances to weight the kinematic positions in the gravity field recovery process, one can find formal errors much closer to the degree variance differences (Fig. 10). The low to medium degree terms, which are prone to systematic errors (e.g., antenna phase centre variations, cf. Jäggi et al., 2009b), still feature discrepancies between the theoretically expected errors and the differences to the superior solution.
Empirical covariances perform better than formally propagated ones because knowledge about model and data deficiencies is both transferred into the leastsquares adjustment. Degrees 10 to 26 are slightly degraded compared to the classical solution (Fig. 10, red curve), which needs further inspection. The influence of the PCA is still visible in the low degrees, however, the dynamic solution is not degraded in the higher degrees anymore (compare with Fig. 3) because the empirical covariances are capable of dealing with noise that is absorbed by the PCA.
In this section we give a short outlook on the use of undifferenced ambiguity fixed kinematic positions for gravity field determination. The process of undifferenced ambiguity fixing can be found in Schaer et al. (2020).
Introducing such kinematic positions potentially helps to improve the gravity field determination process as indicated by the performance of the Kband validation of a classical reduceddynamic orbit (see Fig. 11). The longperiodic noise gets drastically reduced and the fit of the kinematic positions agrees much better with the Kband range observations (RMS over one month: 0.24 cm compared to 1.8 cm). It is lower as if one employs the formally propagated covariance matrix over a sufficient amount of epochs and close to a doubledifference ambiguity fixed baseline solution where submillimetre KBand validations can be achieved (Jäggi et al., 2009a).
In the processing of carrier phase observations as many ambiguities as possible are resolved to their integer values before the leastsquares adjustment for the kinematic positions takes place, and only the unresolved ambiguities remain in the system. Thus, the number of ambiguity parameters is significantly reduced by more than 90 %. This implies that the kinematic positions are almost uncorrelated in time, which is shown in Fig. 12. Note the differences to the covariance from the corresponding ambiguity float solution in Fig. 4. The correlations vanish after less than ten minutes in the ambiguityfixed case.
In theory, these short correlations in time carry the whole information on the stochastic behaviour. However, for a gravity field solution it is not sufficient to introduce only these correlations in the leastsquares adjustment to obtain a reasonably realistic description of the stochastic behaviour; the system behaves almost like employing epochwise covariances, which could be expected from the shape of the covariance function and the Kband validation results. The formal errors of the gravity field solution are at the same level as Fig. 3 shows for float solutions. Thus, the correlations overarching epochs up to 30 min and more, which cannot be reflected in the covariance matrix of the integerfixed positions anymore, transport essential information about the noise characteristics. In other words, the ambiguity fixing process does not obliterate the correlations (it diminishes them slightly) but it narrows the possibilities of mapping them.
Using undifferenced ambiguity fixed positions together with empirical covariances yields results at the same level as float ambiguities and empirical covariances (not shown). Concerning the signal content in the kinematic positions for gravity field recovery, currently the use of ambiguity fixed positions does not outperform results from a float solution.
Kinematic positions are widely used pseudoobservations in gravity field determination and their stochastic behaviour plays a crucial role when assessing the quality of gravity field solutions. Kinematic positions are correlated in time due to the ambiguities in the original GPS phase observations and not considering these correlations degrades a subsequent gravity field solution. Taking these correlations into account can be achieved either by the covariance information of the kinematic positions over longer time spans or by empirical covariances derived from residuals of an orbit fit with respect to the kinematic positions. However, these empirical covariances demand a stationary process and are, therefore, highly dependent on the quality and shape of the residuals.
Formally propagated covariances overarching epochs are in case of a dynamic orbit parametrisation able to absorb parts of the errors caused by a deficient force field, piecewise constant accelerations have similar capabilities for reduceddynamic parametrisations. Empirical covariances support the estimation of physically meaningful formal errors and absorb errors caused by a deficient force field in case of dynamic and reduceddynamic orbit descriptions.
Undifferenced ambiguity fixed kinematic positions feature different stochastic characteristics than ambiguity float kinematic positions – they are almost uncorrelated – however, the possibility of describing the deficiencies that are still in the processing chain is lost through the ambiguity fixing process. Currently, an improvement of a gravity field solutions by using ambiguityfixed kinematic positions as pseudoobservations cannot be substantiated.
Further steps to investigate the described methods of stochastic modelling include the diversification to different satellite missions, such as Swarm or GOCE, and an appropriate treatment of deteriorated data (affected by outliers etc.), for empirical covariances.
Data for GRACE was accessed via GFZ's ISDC data centre ftp://isdcftp.gfzpotsdam.de/gracefo/ (GFZ/JPL, 2020).
ML performed the computational work of gravity field determination, DA contributed the kinematic positions. Concept, framework and interpretation of the study was accomplished by all authors.
The authors declare that they have no conflict of interest.
This article is part of the special issue “European Geosciences Union General Assembly 2019, EGU Geodesy Division Sessions G1.1, G2.4, G2.6, G3.1, G4.4, and G5.2”. It is a result of the EGU General Assembly 2019, Vienna, Austria, 7–12 April 2019.
We thank reviewers for insightful comments which helped to greatly improve the manuscript.
This research has been supported by the Swiss National Science Foundation (SNSF) (grant no. 200021_175942 Assessment of Noise Models for GRACE and GRACEFO).
This paper was edited by Annette Eicker and reviewed by Akbar Shabanloui and one anonymous referee.
Baur, O., Bock, H., Höck, E., Jäggi, A., Krauss, S., MayerGürr, T., Reubelt, T., Siemes, C., and Zehentner, N.: Comparison of GOCEGPS gravity fields derived by different approaches, J. Geodesy, 88, 959–973, https://doi.org/10.1007/s0019001407366, 2014. a
Beutler, G., Jäggi, A., Mervart, L., and Meyer, U.: The celestial mechanics approach: theoretical foundations, J. Geodesy, 84, 605–624, https://doi.org/10.1007/s0019001004017, 2010. a
Bock, H., Dach, R., Jäggi, A., and Beutler, G.: Highrate GPS clock corrections from CODE: support of 1 Hz applications, J. Geodesy, 83, 1083, https://doi.org/10.1007/s0019000903261, 2009. a
Dach, R., Brockmann, E., Schaer, S., Beutler, G., Meindl, M., Prange, L., Bock, H., Jäggi, A., and Ostini, L.: GNSS processing at CODE: status report, Jo. Geodesy, 83, 353–365, https://doi.org/10.1007/s0019000802812, 2009. a
Dach, R., Lutz, S., Walser, P., and Fridez, F.: Bernese GPS software version 5.2, Documentation, University of Bern, Bern Open Publishing, https://doi.org/10.7892/boris.72297, 2015. a
Dobslaw, H., BergmannWolf, I., Dill, R., Poropat, L., Thomas, M., Dahle, C., Esselborn, S., König, R., and Flechtner, F.: A new highresolution model of nontidal atmosphere and ocean mass variability for dealiasing of satellite gravity observations: AOD1B RL06, Geophys. J. Int., 211, 263–269, https://doi.org/10.1093/gji/ggx302, 2017. a
Drinkwater, M., Haagmans, R., Muzi, D., Popescu, A., Floberghagen, R., Kern, M., and Fehringer, M.: The GOCE gravity mission: ESA's first core explorer, 1–7, ESA SP627, Frascati, Italy, 2006. a
Etten, W. V.: Introduction to Random Signals and Noise, John Wiley & Sons, Hoboken, NJ, USA, 2005. a
Flechtner, F., Morton, P., Watkins, M., and Webb, F.: Status of the GRACE followon mission, Gravity, Geoid and Height Systems, IAG Symposia, 117–121, https://doi.org/10.1007/9783319108377_15, 2013. a
Folkner, W. M., Williams, J. G., and Boggs, D. H.: The Planetary and Lunar Ephemeris DE 421, Interplanetary Network Progress Report, 41–178, Jet Propulsion Laboratory, Pasadena, California, 2009. a
FriisChristensen, E., Lühr, H., Knudsen, D., and Haagmans, R.: Swarm – An Earth Observation Mission investigating Geospace, Adv. Space Res., 41, 210–216, https://doi.org/10.1016/j.asr.2006.10.008, 2006. a
Gerlach, C., Földvary, L., Švehla, D., Gruber, T., Wermuth, M., Sneeuw, N., Frommknecht, B., Oberndorfer, H., Peters, T., Rothacher, M., Rummel, R., and Steigenberger, P.: A CHAMPonly gravity field model from kinematic orbits using the energy integral, Geophys. Res. Lett., 30, https://doi.org/10.1029/2003GL018025, 2003. a
GFZ/JPL: GRACE, available at: ftp://isdcftp.gfzpotsdam.de/gracefo/, last access: 2 September 2020. a
Jäggi, A., Beutler, G., Prange, L., Dach, R., and Mervart, L.: Assessment of GPSonly observables for Gravity Field Recovery from GRACE, Observing our Changing Earth, 133, 113–123, https://doi.org/10.1007/9783540854265_14, 2009a. a
Jäggi, A., Dach, R., Montenbruck, O., Hugentobler, U., Bock, H., and Beutler, G.: Phase center modeling for LEO GPS receiver antennas and its impact on precise orbit determination, J. Geodesy, 83, 1145–1162, https://doi.org/10.1007/s0019000903332, 2009b. a
Jäggi, A., Bock, H., Prange, L., Meyer, U., and Beutler, G.: GPSonly gravity field recovery with GOCE, CHAMP, and GRACE, Adv. Space Res., 47, 1020–1028, https://doi.org/10.1016/j.asr.2010.11.008, 2011a. a, b
Jäggi, A., Prange, L., and Hugentobler, U.: Impact of covariance information of kinematic positions on orbit reconstruction and gravity field recovery, Adv. Space Res., 47, 1472–1479, https://doi.org/10.1016/j.asr.2010.12.009, 2011b. a, b, c, d, e, f, g, h
Jäggi, A., Bock, H., Meyer, U., Beutler, G., and van den IJssel, J.: GOCE: assessment of GPSonly gravity field determination, J. Geodesy, 89, 33–48, https://doi.org/10.1007/s001900140759z, 2015. a
Jäggi, A., Dahle, C., Arnold, D., Bock, H., Meyer, U., Beutler, G., and van den IJssel, J.: Swarm kinematic orbits and gravity fields from 18 months of GPS data, Adv. Space Res., 57, 218–233, https://doi.org/10.1016/j.asr.2015.10.035, 2016. a
Lück, C., Kusche, J., Rietbroek, R., and Löcher, A.: Timevariable gravity fields and ocean mass change from 37 months of kinematic Swarm orbits, Solid Earth, 9, 323–339, https://doi.org/10.5194/se93232018, 2018. a
MayerGürr, T., Kvas, A., Klinger, B., Rieser, D., Zehentner, N., Pail, R., Gruber, T., Fecher, T., Rexer, M., Schuh, W.D., Kusche, J., Brockmann, J.M., Loth, I., Müller, S., Eicker, A., Schall, J., Baur, O., Höck, E., Krauss, S., and Maier, A.: The new combined satellite only model GOCO05s, Geophys. Res. Abstr., EGU201512364, EGU General Assembly 2015, Vienna, Austria, https://doi.org/10.13140/RG.2.1.4688.6807, 2015. a
Meyer, U., Jean, Y., Kvas, A., Dahle, C., Lemoine, J. M., and Jäggi, A.: Combination of GRACE monthly gravity fields on the normal equation level, J. Geodesy, 93, 1645–1658, https://doi.org/10.1007/s00190019012746, 2019. a
Pavlis, N. K., Holmes, S. A., Kenyon, S. C., and Factor, J. K.: The development and evaluation of the Earth Gravitational Model 2008 (EGM2008), J. Geophys. Res.Sol. Ea., 117, B04406, https://doi.org/10.1029/2011JB008916, 2012. a
Petit, G. and Luzum, B.: IERS Conventions (2010), IERS Technical Note No. 36, Verlag des Bundesamts für Kartographie und Geodäsie, Frankfurt am Main, 2010. a
Prange, L.: Global Gravity Field Determination Using the GPS Measurements Made Onboard the Low Earth Orbiting Satellite CHAMP, PhD thesis, Schweizerische Geodätische Kommission, Institut für Geodäsie und Photogrammetrie, Eidg. Technische Hochschule Zürich, Zurich, Switzerland, 2010. a, b
Prange, L., Jäggi, A., Beutler, G., Dach, R., and Mervart, L.: Gravity Field Determination at the AIUB – the Celestial Mechanics Approach, Observing our Changing Earth, 133, 353–362, https://doi.org/10.1007/9783540854265_42, 2009. a
Reigber, C., Lühr, H., and Schwintzer, P.: Status of the CHAMP Mission, Towards an Integrated Global Geodetic Observing System (IGGOS), IAG Symposia, 63–65, https://doi.org/10.1007/9783642597459, 1998. a
Savcenko, R. and Bosch, W.: EOT11a – a new tide model from MultiMission Altimetry, OSTST Meeting, 19–21 October 2011, San Diego, 2011. a
Schaer, S., Villiger, A., Arnold, D., Prange, L., and Jäggi, A.: The CODE ambiguityfixed clock and phase bias analysis and their properties and performance, J. Geodesy, in review, 2020. a, b
Sušnik, A., Dach, R., Villiger, A., Maier, A., Arnold, D., Schaer, S., and Jäggi, A.: CODE reprocessing product series, CODE_REPRO_2015, [Dataset], BORIS, https://doi.org/10.7892/boris.80011, 2016. a
Sušnik, A., Grahsl, A., Arnold, D., Villiger, A., Dach, R., Beutler, G., and Jäggi, A.: Validation of the EGSIEMREPRO GNSS Orbits and Satellite Clock Corrections, Remote Sensing, 12, 2322, https://doi.org/10.3390/rs12142322, 2020. a, b
Švehla, D. and Rothacher, M.: Kinematic precise orbit determination for gravity field determination, A window on the future of geodesy, IAG Symposia, 128, 181–188, https://doi.org/10.1007/3540274324_32, 2005. a, b
Tapley, B. D., Bettadpur, S., Watkins, M., and Reigber, C.: The gravity recovery and climate experiment: Mission overview and early results, Geophys. Res. Lett., 31, https://doi.org/10.1029/2004GL019920, 2004. a
Zumberge, J. F., Heflin, M. B., Jefferson, D. C., Watkins, M. M., and Webb, F. H.: Precise point positioning for the efficient and robust analysis of GPS data from large networks, J. Geophys. Res.Sol. Ea., 102, 5005–5017, 1997. a