Research article 15 Aug 2019
Research article  15 Aug 2019
Multiresolution wavelet analysis applied to GRACE rangerate residuals
 ^{1}Graz University of Technology, Institute of Geodesy, Steyrergasse 30/III, 8010 Graz, Austria
 ^{2}Leibniz University Hanover, Institute of Geodesy, Schneiderberg 50, 30167 Hanover, Germany
 ^{3}Jet Propulsion Laboratory, NASA, Pasadena, CA, USA
 ^{*} Invited contribution by Saniya Behzadpour, recipient of the EGU Geodesy Outstanding Student Poster and PICO Award 2018.
 ^{1}Graz University of Technology, Institute of Geodesy, Steyrergasse 30/III, 8010 Graz, Austria
 ^{2}Leibniz University Hanover, Institute of Geodesy, Schneiderberg 50, 30167 Hanover, Germany
 ^{3}Jet Propulsion Laboratory, NASA, Pasadena, CA, USA
 ^{*} Invited contribution by Saniya Behzadpour, recipient of the EGU Geodesy Outstanding Student Poster and PICO Award 2018.
Correspondence: Saniya Behzadpour (behzadpour@tugraz.at)
Hide author detailsCorrespondence: Saniya Behzadpour (behzadpour@tugraz.at)
For further improvements of gravity field models based on Gravity Recovery and Climate Experiment (GRACE) observations, it is necessary to identify the error sources within the recovery process. Observation residuals obtained during the gravity field recovery contain most of the measurement and modeling errors and thus can be considered a realization of actual errors.
In this work, we investigate the ability of wavelets to help in identifying specific error sources in GRACE rangerate residuals. The multiresolution analysis (MRA) using discrete wavelet transform (DWT) is applied to decompose the residual signal into different scales with corresponding frequency bands. Temporal, spatial, and orbitrelated features of each scale are then extracted for further investigations.
The wavelet analysis has proven to be a practical tool to find the main error contributors. Besides the previously known sources such as Kband ranging (KBR) system noise and systematic attitude variations, this method clearly shows effects which the classic spectral analysis is hardly able or unable to represent. These effects include longterm signatures due to satellite eclipse crossings and dominant ocean tide errors.
For more than 15 years, the Gravity Recovery and Climate Experiment (GRACE) satellite mission measured the time variation of Earth's gravity field with high temporal and spatial resolutions (Tapley et al., 2004). The mission was a trailing formation of two satellites, GRACEA and GRACEB, and provided the observation signals of intersatellite ranging, GPS tracking, the satellite attitudes, and nongravitational accelerations, which are required for the gravity field parameter estimation.
Based on these observations, various timevariable gravity models with monthly resolution were published by different analysis centers (e.g., Bettadpur, 2012; Dahle et al., 2012; MayerGürr et al., 2016). The accuracy level of such models has gradually increased in recent years; however, it has not reached the GRACE baseline accuracy computed through prelaunch simulations (Kim, 2000; Kim and Tapley, 2002). This results in an ongoing effort to understand the error content of GRACE observations, as well as any inaccuracies in the physical and stochastic models used for processing GRACE data.
In recent years, significant research efforts have been made to identify and parametrize the systematic errors such as uncertainties in star camera alignment (Bandikova and Flury, 2014; Harvey, 2016) and accelerometer calibration (Klinger and MayerGürr, 2016). Along with the effects of geophysical aliasing and uncertainties in background models, these errors propagate through the numerical estimation of a large number of parameters. These parameters include gravity parameters in terms of spherical harmonic coefficients as well as orbit and sensor calibration parameters (MayerGürr, 2013).
If the calibration parameters are correctly adjusted and the stochastic model fully describes the observation noise, it is expected that all of the mentioned errors are completely contained within the residuals. In reality, however, these errors might affect the gravity parameters due to imperfections in modeling. Therefore, residual analysis becomes a research topic as it is not only a way to study measurement and physical modeling errors, but also helps to evaluate and improve the gravity field solutions.
The studies in this field have been conducted mainly on the theoretical residuals, which are the difference between the actual GRACE ranging observation and simulated observation computed through force models. Ditmar et al. (2012) applied spectral analysis on theoretical residuals and showed that the major contributors to the noise budget at high frequencies are Kband ranging (KBR) sensor noise and inaccuracies in Earth's assumed static gravity field at higher degrees. It also has been shown that uncertainties in background models and errors in computed dynamic orbits contribute to lowfrequency noise.
The main challenge in the spectral analysis of the residuals is that several noisy signals and disturbances are known to be superimposed at each frequency. Furthermore, the analysis is based on the assumption of the stationary behavior of these signals. However, in reality, most of these signals have nonstationary behavior, meaning that they have dynamic frequency components over time. Classical spectral analysis using Fourier transforms only represents the frequency content of such signals (Fig. 1a) and does not provide any information about the time at which a signal at a specific frequency occurred or the duration for which it lasted (Keller, 2004). Consequently, in this framework, it is not possible to localize each component of the residuals in time for further statistical, spatial, or orbital analysis. The drawback of this framework draws our attention to spatiotemporal approaches, which incorporate data analysis as well as geophysical model validation (e.g., Dransch et al., 2010).
In an attempt to consider time variations in the soughtafter signals, time–frequency methods can be applied to identify and localize the content of the nonstationary signals in the time and frequency domains simultaneously. The simplest method is the shorttime Fourier transform (STFT), which is implemented by sliding a window throughout a signal and applying a Fourier transform to each windowed data segment. The squared magnitudes of the STFT coefficients form a spectrogram, representing the variation of the signal's spectrum over time (Fig. 1b). The shape and length of the window function determines the fixed time and frequency resolution of the STFT. Due to this uniform time resolution for all frequencies, the STFT is limited to capturing time information on rapid changes in a signal as well as spectral information in its lowerfrequency components.
To overcome STFT drawbacks, the wavelet analysis was introduced as a more effective technique for representation, decomposition, and reconstruction of nonstationary signals (Keller, 2004). In contrast to STFT, the wavelet transform provides a better tradeoff between time and frequency resolution by using windows with shorter time spans at higher frequencies and windows with longer time spans at lower frequencies. The multiresolution analysis (MRA), introduced by Mallat (1989) and Meyer (1993), is an efficient implementation of a wavelet transform for real signals. MRA can decompose a signal into multiscale components which can describe all timevariable structures in that signal.
The aim of this paper is to exploit the advantages of the wavelet transform to investigate the major contributors to GRACE rangerate residuals and ideally detect nonstationary noise sources in sensors and background models which cannot be observed with traditional spectral analysis. The results of this study will further improve gravity field modeling based on GRACE data. In addition, they will be beneficial for the preparation of GRACE FollowOn data processing infrastructure. To reach this goal, we decompose the residual signal into three groups of scale and compare the characteristics of each group with known or supposed sources.
In the upcoming Sect. 2, we explain how the residual signal is obtained in the frame of computing the ITSGGrace2016 model and review the performed data processing steps in order to introduce potential error sources. Section 3 discusses the methodology of the multiresolution analysis and the wavelet transform. In Sect. 4, results of the employed method on the residuals are described. Finally, Sect. 5 presents the interpretation of results and a discussion.
In this study, we use GRACE rangerate residuals obtained in the course of computing the ITSGGrace2016 (MayerGürr et al., 2016) gravity field model up to degree and order 60. Therefore, in order to introduce the residual signal, we briefly explain the processing chain of the model (Klinger et al., 2016).
In the ITSGGrace2016 gravity field processing, highprecision kinematic orbits (Zehentner and MayerGürr, 2013) with a sampling of 5 min and Kband intersatellite range rates with a sampling of 5 s serve as observations. Using the approach of variational equations, dynamic orbits are computed for each day (Ellmer and MayerGürr, 2017), and normal equations are set up with an arc length of 3 h. The accumulated normal equations are then solved to estimate gravity parameters in terms of spherical harmonic coefficients, spanning from degree 2 up to degree 60. The background models used during the dynamic orbit integration are listed in Table 1.
MayerGürr et al. (2015)Folkner et al. (2009)Savcenko and Bosch (2012)Dobslaw et al. (2013)van Dam and Ray (2010)Petit and Luzum (2010)In the course of the adjustment process, nongravity parameters are also coestimated for each day. These parameters include the initial orbit states of both satellites, accelerometer scale factor matrices, accelerometer biases modeled by cubic splines with 6 h nodes, and daily gravity field variations up to degree and order 40.
It is worth mentioning that unlike in the standard GRACE monthly solutions, in ITSGGrace2016 the correlations between observations within a data block of 3 h are taken into account. For each observation type, a stochastic model of the observation noise is built under the assumption of stationarity. This model is estimated once per month directly from the observation residuals.
The weights for the different frequency components of the observations are determined through the residual power spectral density (PSD). This PSD is iteratively computed directly from the residuals through variance component estimation (VCE) (Koch, 1999). VCE is also used to estimate the relative weights for the combination of different data types, i.e., kinematic orbits and rangerate observations. This modeling approach seems to appropriately separate the complex colored noise in the observations from the gravity signal; therefore, we expect the residuals to contain most of the imperfections caused by the instruments and background models.
The wavelet transform Wf(u,s) of a signal f(t)∈L^{2}(ℝ),
is the decomposition of that signal over a set of scaled and translated versions of a finite energy and normalized function, the mother wavelet $\stackrel{\mathrm{\u203e}}{\mathit{\psi}}$.
For the wavelet transform Wf(u,s), the translation parameter u determines the location of the wavelet in the time domain, while the scale parameter s is related to the frequency location. These parameters are continuous real values, therefore an infinite number of coefficients are needed to describe a signal in this framework. In a practical implementation, it is convenient to discretize these parameters, as the real signals are band limited. The usual choice is to follow a Jscale dyadic discretization based on powers of 2. This transform is then called a discrete wavelet transform (DWT). For a signal with sampling frequency of F_{S}, the resulting coefficients d(j,n) can be interpreted as detailed subsignals at the scale 2^{j} ($\mathrm{1}\le j\le J$), corresponding to the frequency interval $[{F}_{\mathrm{S}}/{\mathrm{2}}^{j+\mathrm{1}},{F}_{\mathrm{S}}/{\mathrm{2}}^{j}]$ :
The approximation of the signal at the scale J, which corresponds to the frequency interval $[\mathrm{0},{F}_{\mathrm{S}}/{\mathrm{2}}^{J+\mathrm{1}}]$ is also given by
where ϕ(J,n) is the scaling function, associated with the wavelet function ψ(j,n) .
The original signal can be reconstructed by adding all layers of details up to decomposition scale J as well as the approximation subsignal:
Mallat (1989) showed that for a discrete signal f[n], any DWT on the orthonormal basis of L^{2}(ℝ) could be characterized by a particular class of digital filters, the conjugate mirror filters. He introduced a fast discrete wavelet transform by implementing a pair of conjugate mirror filters, corresponding to a specific mother wavelet:
Mathematically, the convolution of the filter response with the discrete signal is expressed as follows:
The scaling function, defined by the filter coefficients h[n], provides approximation coefficients a, which are also referred to as lowpass output. The wavelet function, defined by the filter coefficients g[n], provides the detailed coefficients d, or alternatively the highpass output. This decomposition step is followed by a factor 2 downsampling of the output signals. According to Vetterli and Herley (1992), downsampling cancels the aliasing between the resulting coefficients. This is a necessary condition for recovery of the original signal with an inverse DWT.
A fast inverse DWT reconstructs the initial signal f[n] by upsampling and filtering. The upsampling operation is done by inserting zeroes between every other coefficients in the output signals a[n] and d[n]. The zeropadded coefficients $\widehat{a}$ and $\widehat{d}$ are then filtered by the corresponding inverse filters $\stackrel{\mathrm{\u0303}}{h}\left[n\right]$ and $\stackrel{\mathrm{\u0303}}{g}\left[n\right]$:
As described before, the DWT decomposes the original signal into an approximation subsignal and detailed subsignals. The MRA algorithm suggested by Mallat (1989) and Meyer (1993) calls for this decomposition to be repeated on the approximation subsignal, again yielding detailed subsignals and an approximation subsignal. The selection of the decomposition level depends on the initial size of the original signal, and the desired spectral and temporal resolution. Finally, the original signal can be represented by the approximation coefficients of the last decomposition level and the accumulated detailed coefficients of all decomposition levels. Figure 2 shows a threelevel decomposition MRA algorithm.
We applied MRA using a discrete Daubechies wavelet transform with 20 vanishing moments (Daubechies, 1992) to decompose a monthly time series of residuals into eight different scales. The choice of the Daubechies wavelet is due to its usual application in signal detection and classification. The selection of a high vanishing moment is due to a high smoothness property of the resulting mother wavelet, leading to a better frequency localization in the millihertz (mHz) frequency band. Figure 3 shows scaling and wavelet functions for Daubechies20 together with its corresponding decomposition and reconstruction conjugate mirror filters.
As shown in Fig. 4, we merged detailed coefficients into three major groups, defined approximately through three frequency subbands (Fig. 5):
 a.
short timescale details, containing the details at levels 1 to 3, corresponding to the frequency band above 12.5 mHz;
 b.
medium timescale details, containing the details at levels 4 to 5, corresponding to the frequency range from 3.125 up to 12.5 mHz;
 c.
long timescale details, containing the details at levels 6 to 8, corresponding to the frequency range from 0.391 up to 3.125 mHz.
Each group is then reconstructed into a time series of residuals using Eq. (9). Afterward, the time series are analyzed in three different domains. We have chosen the domains in such a way that they highlight specific characteristics of the error sources contained within the residual time series. They are the following.
 Spectral/temporal domain.

As mentioned in the first section, a spectrogram shows the variation of a signal's energy as a function of time and frequency. Another tool which can be used directly on the wavelet coefficients is the scalogram, in which the amplitude of the coefficients are plotted as a function of the scale and transition parameters. In our analyses, we used spectrograms because the interpretation of a signal in terms of frequency is more accessible than in terms of scale (Fig. 6).
 Spatial domain.

Plotting each time series with respect to the satellite ground track is useful to identify any features of geophysical origin in the data (Fig. 7).
 Orbital domain.

Plotting each time series as a function of satellite position and time reveals features related to the orbit geometry or instrument errors caused by orbital conditions. As the GRACE orbits are nearcircular, the position of each satellite can be specified without loss of accuracy by the argument of latitude, ranging from −180 to 180^{∘}. This domain represents the ascending Equator pass of the satellite at 0^{∘}, the north pole at 90^{∘}, the descending Equator pass at 180 or −180^{∘}, and then the south pole at −90^{∘} (Fig. 8).
These analyses are carried out on the whole ITSGGrace2016 time span (April 2002–June 2017). However, due to low data quality before 2004 and several data gaps and degraded quality of the measurements after 2016, these time periods are excluded from the illustrations. Highlights of this analysis are presented in the next section.
To prove whether or not our applied method using the DWT is applicable to detect the error sources, we initially focused on the investigation of known issues. For instance, it is known that the Kband system noise is dominant in the frequency range above 12.5 mHz. This frequency band corresponds to the short timescale details of the residuals. The power of the noise in this band increases linearly with frequency. This is a result of the way the rangerate observations are derived from the range measurements by differentiation. Investigations by Ko et al. (2012) and Harvey et al. (2017) showed that the excessive highfrequency signatures in this band are highly correlated with low signaltonoise (SNR) values of the Kband frequency observation by GRACEB. Figure 9 compares these SNR values with the wavelet short timescale components, revealing this strong correlation.
According to Bandikova et al. (2012), residuals due to errors in the satellite attitude determination and their effects on the computed antenna offset correction are also expected to be found in the millihertz (mHz) frequency band. Our time–frequency analysis shows a similarity between the residuals in medium timescale and the angular acceleration variations derived from star camera observations (Fig. 10). The spatial pattern of the residuals related to the attitude variations appears as horizontal bands (Fig. 7b), consistent with the results presented by Inácio et al. (2015).
These first investigations already show that our applied method is well suited to identify error sources. However, compared to the spectral analysis, the advantage of the implemented method of DWT is a better separation of superimposed signals in frequencies lower than 12.5 mHz. This enabled the identification of (a) systematic errors caused by eclipse crossings of the satellites and (b) dominant ocean tide model errors, which are explained in the following sections.
4.1 Satellite eclipse crossings
Analysis of the medium timescale details throughout the GRACE time span reveals longterm systematic signatures (Fig. 11a). Although the source of these errors is unknown, our investigation revealed a high correlation with the eclipse transit phases of GRACEA and GRACEB.
Each satellite passes through partial or full eclipse phases when it enters Earth's shadow. Occasionally the Moon also casts a shadow on the satellites. The eclipse factor is defined as the fraction of the Sun's light that reaches the satellite. It has a minimum value of zero if the satellite is in the umbra of the occulting body and a maximum value of one if the satellite is in direct sunlight. For a detailed calculation, the reader is referred to Montenbruck and Gill (2000).
The difference between GRACEB and GRACEA eclipse factors indicates if the mission, i.e., one of the satellites, is in a transit mode. Difference values not equal to zero are interpreted as transit events, in which one of the satellites is passing through a partial eclipse phase. Figure 11b compares the medium timescale details of the residuals with the transit events for the complete GRACE time span. Before 2011, the signatures are most obvious when the difference value is negative, meaning that GRACEA is in the shadow and GRACEB is in sunlight.
The GRACE formation mission started with GRACEA as leading and GRACEB as trailing satellite. After 3 years in orbit, the satellites had to exchange their positions to limit the damage on the Kband horn caused by atomic oxygen. This swap maneuver happened at the end of 2005. Before this time, eclipse crossing signature occurs when the pair entered sunlight. After the orbit swap maneuver in December 2005, when GRACEB became the leading satellite, the signatures are visible when the pair enters the shadow area.
However, after the year 2011, these rules cannot explain the eclipse crossing signatures in the residuals as they appear in both entering and leaving shadow conditions with different intensities. The unstable thermal condition due to the disabled thermal controls might be a possible reason.
We compared the temperature measurements obtained from Level1A HighResolution Temperature data (HRT1A) for November 2008 and October 2011 with these signatures. It becomes obvious that there is a high correlation between the GRACEB Kband antenna horn temperature variation and the disturbances during eclipse crossing events (Fig. 12). We suggest that the increasing temperature on the GRACEB antenna horn may produce disturbances in the KBR measurements. This hypothesis can be investigated in more detail once the complete GRACE Level1A datasets become publicly available. From a gravity field recovery point of view, these eclipse crossing signatures can be interpreted as a temporary unmodeled signal in the rangerate measurements.
4.2 Ocean tide model
Errors in the background force models of temporal gravity field variations can be found in the long timescale details. Due to the spatial nature of these errors and the periodicity of satellite passes over their source regions, different model errors are superimposed at the same n cyclesperrevolution frequencies. Therefore, frequency or time–frequency plots cannot differentiate the dominant source from other influences at this detailed scale.
The two main potential error sources at this scale are (a) inaccuracies in the employed ocean tide model EOT11a (Savcenko and Bosch, 2012) and (b) inaccuracies in the employed nontidal atmosphere and ocean mass variation model, AOD1B RL05 (Dobslaw et al., 2013). To better understand the contributions of the individual models, we swap in alternative models of the same forces and studied the resulting differences. This is best done in a closedloop simulation, where other contributors to noise can be controlled. The simulation is carried out for the time period 2008–2009, when GRACE delivered highquality measurements and comparison of the actual data with the output of the simulation is more relevant. The following steps outline our employed simulation process:

Dynamic orbits are computed based on the background models mentioned in Table 1 with two exceptions. First, the FES2014 ocean tide model (Carrere et al., 2015) is substituted for the EOT11a model. Second, the AOD1B RL05 model and the van Dam–Ray atmospheric tide model (van Dam and Ray, 2010) were substituted with the AOD1B RL06 model. Compared to AOD1B RL05, the AOD1B RL06 model (Dobslaw et al., 2017) has undergone several improvements, amongst them a higher temporal resolution and the separation of nontidal and tidal signals, including atmospheric tides with 12 selected frequencies. Therefore, there was no need to consider a dedicated atmospheric tide model in the simulation employing AOD1B RL06.

Errorfree observations for position, velocity, nongravitational accelerations, and the KBand instrument are synthesized from these ideal orbits.

Realistic models of instrument noise are used to degrade synthesized observations. White Gaussian noise with a standard deviation of 3 cm is added to the simulated satellite positions. Accelerometer observations are degraded by white noise with a standard deviation of 0.3 nm s^{−2} in alongtrack and radial directions and 3 nm s^{−2} in the crosstrack direction. Star camera instrument noise is added as white Gaussian noise with a standard deviation of 0.05 mrad to the orientation quaternions. KBR instrument noise is computed by applying a differential filter to white Gaussian noise with a standard deviation of 0.25 µm s^{−1}, which is then added to the simulated rangerate observations.

The final step is to recover a monthly gravity field using the simulated degraded observations. To this end, the dynamic orbits are reintegrated using the artificially degraded accelerometer observations and the separate models under study, each in a dedicated scenario. The respective obtained residuals are then analyzed and compared.
4.2.1 Simulation scenario 1: propagated errors due to instrument noise
In the first scenario, the same background models as mentioned in the first step of the simulation process are used to compute the reintegrated dynamic orbits. Therefore the results only show the effects of instrument noise. As expected, the propagated noise is 1 order of magnitude smaller than the real residuals in frequency range from 0.391 up to 3.125 mHz (Fig. 13b) and obviously cannot explain the errors in the long timescale details. Analyzing the solution in terms of RMS geoid heights per degree with respect to the reference field GOCO05s, it can also be seen that the monthly solution based on instrument noise alone exhibits differences smaller than those of the GRACE baseline (Fig. 13a).
4.2.2 Simulation scenario 2: propagated errors due to AOD1B RL05
The second scenario studies the propagated errors due to inaccuracies of the nontidal mass variation model. In order to recover a gravity field in this scenario, the AOD1B RL05 model and the van Dam–Ray atmospheric tide model (van Dam and Ray, 2010) were substituted for the AOD1B RL06 model. The simulated residual signal is then decomposed, and its long timescale components are compared to those obtained from real data. As shown in Fig. 13b, although the propagated errors have the same spectral behavior at frequency range from 0.391 to 3.125 mHz, their magnitude and spatial structure (Fig. 14a) cannot explain the real residuals.
4.2.3 Simulation scenario 3: propagated errors due to EOT11a
In the third scenario, we study the contribution of the ocean tide model. To recover a gravity field in this scenario, the EOT11a ocean tide model is substituted for the FES2014 model. After decomposition of the simulated residual signal, its long timescale components are compared to the real data. These errors have comparable magnitude and spatial pattern (Fig. 14b) as those in the real data (Fig. 14c). This leads to the conclusion that the ocean tide model is the dominant error source at the long timescale detailed level.
These results showcase the capability of wavelet analysis in studying the signals due to geophysical processes in GRACE rangerate residuals. The implemented method efficiently finds structures in the signal which are not explicitly apparent in the PSD of the residuals. The wavelet analysis proves to be an efficient tool in decomposing the background model errors and finding the most prominent sources.
The results presented in this paper show the advantages of using a DWT in analyzing the rangerate residuals from the ITSGGrace2016 gravity field model. Several improvements in ITSGGrace2016 resulted in a cumulative noise reduction of 20 %–40 % compared to its predecessor ITSGGrace2014. The proposed analysis framework confirms known and reveals previously unknown systematics in the residuals that allow for a specifically tailored parametrization in the gravity field retrieval.
We showed that the short timescale details of the residuals, equivalent to frequencies above 12.5 mHz, are dominated by KBR system noise. This is in agreement with the results presented by Ko et al. (2012) and Ditmar et al. (2012). The errors in the satellite attitude determination were identified as a major contributor in the medium timescale details, equivalent to the frequency range from 3.125 to 12.5 mHz. This finding is consistent with the results presented by Inácio et al. (2015) and Bandikova et al. (2012).
Besides the previously known instrument error sources, longterm signatures due to eclipse transits of the satellites were identified. They appear as a bias term in the Kband rangerate observations. As this is a clearly deterministic effect, its influence can be reduced by coestimation of additional calibration parameters in the gravity field recovery process.
Analysis of the results from the implemented discrete wavelet transform brings new insights and a new understanding of the signals at the long timescale level. At this level, spectral analysis is unable to differentiate between the individual contributing sources, due to the nonstationary nature of the errors. Knowing that this scale level contains valuable information about the timevariable gravity field signal, we introduced nontidal mass variation and ocean tide models as the potential dominant sources. Comparing simulation results with the real data scenario, the EOT11a ocean tide errors are identified as the dominant error source within this scale. This means that using a more accurate ocean tide model can lower the residuals in this frequency band.
It has been shown that the waveletbased MRA approach can properly represent the major error sources in GRACE processing data. These error sources have the largest impact on the accuracy of gravity field solutions derived from observations by GRACE. Even if the purpose of this study is to find the degrading factors in monthly gravity field models, which are mainly affecting the observations in the millihertz (mHz) frequency band, the investigation will be further continued by looking for physical interpretations for features at the lower frequencies of the residuals. This can be achieved by using a wavelet base with higher vanishing moments and thus higher decomposition level.
Besides the rangerate observations, the presented framework is also beneficial for the data processing of the other sensors aboard GRACE or similar satellite missions. The results can potentially detect inconsistent time periods in each set of measurements and provide an initial interpretation of their possible origin.
The GRACE Level1B data (https://podaactools.jpl.nasa.gov/drive/files/GeodeticsGravity/grace/L1B, last access: 13 August 2019) as well as the ITSGGrace2016 gravity field solutions (https://doi.org/10.5880/icgem.2016.007; MayerGürr et al., 2016) are publicly available. The range rate residuals obtained from ITSGGrace2016 models are provided upon request (behzadpour@tugraz.at).
SB and TMG developed and carried out the analysis. JF, BK, and SG provided reviews and suggestions on the GRACE instrument characteristics. SB was the lead author and all coauthors contributed to drafting and editing of the paper.
The authors declare that they have no conflict of interest.
We would like to thank Srinivas Bettadpur from the Center for Space Research, The University of Texas at Austin for providing us the GRACE Level1A test datasets. We would also like to thank two anonymous reviewers for their suggestions and comments that contributed to improving the article.
This research has been supported by the German Research Foundation (DFG) (Relativistic Geodesy and Gravimetry with Quantum Sensors (geoQ), grant no. SFB 1128).
This paper was edited by Lev Eppelbaum and reviewed by two anonymous referees.
Bandikova, T. and Flury, J.: Improvement of the GRACE star camera data based on the revision of the combination method, Adv. Space Res., 54, 1818–1827, https://doi.org/10.1016/j.asr.2014.07.004, 2014. a
Bandikova, T., Flury, J., and Ko, U.D.: Characteristics and accuracies of the GRACE intersatellite pointing, Adv. Space Res., 50, 123–135, https://doi.org/10.1016/j.asr.2012.03.011, 2012. a, b, c
Bettadpur, S.: UTCSR Level2 Processing Standards Document for Level2 Product Release 0005, Tech. rep., Center for Space Research, The University of Texas at Austin, 2012. a
Carrere, L., Lyard, F., Cancet, M., and Guillot, A.: FES2014, a new tidal model on the global ocean with enhanced accuracy in shallow seas and in the Arctic region, Geophys. Res. Abstr., EGU20155481, EGU General Assembly 2015, Vienna, Austria, 2015. a
Dahle, C., Flechtner, F., Gruber, C., König, D., König, R., Michalak, G., and Neumayer, K.H.: GFZ GRACE Level2 Processing Standards Document for Level2 Product Release 0005, https://doi.org/10.2312/gfz.b10312020, Deutsches GeoForschungsZentrum (GFZ), 2012. a
Daubechies, I.: Ten Lectures on Wavelets, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, https://doi.org/10.1137/1.9781611970104, 1992. a
Ditmar, P., da Encarnação, J. T., and Farahani, H. H.: Understanding data noise in gravity field recovery on the basis of intersatellite ranging measurements acquired by the satellite gravimetry mission GRACE, J. Geodesy, 86, 441–465, https://doi.org/10.1007/s0019001105316, 2012. a, b
Dobslaw, H., Flechtner, F., Bergmann‐Wolf, I., Dahle, C., Dill, R., Esselborn, S., Sasgen, I., and Thomas, M.: Simulating high‐frequency atmosphere‐ocean mass variability for dealiasing of satellite gravity observations: AOD1B RL05, J. Geophys. Res.Ocean, 118, 3704–3711, https://doi.org/10.1002/jgrc.20271, 2013. a, b
Dobslaw, H., Bergmann‐Wolf, 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
Dransch, D., Köthur, P., Schulte, S., Klemann, V., and Dobslaw, H.: Assessing the quality of geoscientific simulation models with visual analytics methods – a design study, Int. J. Geogr. Inf. Sci., 24, 1459–1479, https://doi.org/10.1080/13658816.2010.510800, 2010. a
Ellmer, M. and MayerGürr, T.: High precision dynamic orbit integration for spaceborne gravimetry in view of GRACE Followon, Adv. Space Res., 60, 1–13, https://doi.org/10.1016/j.asr.2017.04.015, 2017. a
Folkner, W., Williams, J., Boggs, D., Park, R., and Kuchynka, P.: The Planetary and Lunar Ephemeris DE421, The Interplanetary Network Progress Report 42178, Jet Propulsion Laboratory, Pasadena, California, available at: http://ipnpr.jpl.nasa.gov/progress_report/42178/178C.pdf (last access: 9 August 2019), 2009. a
Harvey, N.: GRACE star camera noise, Adv. Space Res., 58, 408–414, https://doi.org/10.1016/j.asr.2016.04.025, 2016. a
Harvey, N., Dunn, C. E., Kruizinga, G. L., and Young, L. E.: Triggering Conditions for GRACE Ranging Measurement SignaltoNoise Ratio Dips, J. Spacecraft Rockets, 54, 327–330, https://doi.org/10.2514/1.a33578, 2017. a
Inácio, P., Ditmar, P., Klees, R., and Farahani, H. H.: Analysis of star camera errors in GRACE data and their impact on monthly gravity field models, J. Geodesy, 89, 551–571, https://doi.org/10.1007/s0019001507971, 2015. a, b
Keller, W.: Wavelets in Geodesy and Geodynamics, Walter de Gruyter, Berlin, https://doi.org/10.1515/9783110198188, 2004. a, b
Kim, J.: Simulation study of a lowlow satellitetosatellite tracking mission, PhD thesis, The University of Texas at Austin, available at: http://geodesy.geology.ohiostate.edu/course/refpapers/Kim_diss_GRACE_00.pdf (last access: 9 August 2019), 2000. a
Kim, J. and Tapley, B.: Error Analysis of a LowLow SatellitetoSatellite Tracking Mission, J. Guid. Control Dynam., 25, 1100–1106, https://doi.org/10.2514/2.4989, 2002. a
Klinger, B. and MayerGürr, T.: The role of accelerometer data calibration within GRACE gravity field recovery: Results from ITSGGrace2016, Adv. Space Res., 58, 1597–1609, https://doi.org/10.1016/j.asr.2016.08.007, 2016. a
Klinger, B., MayerGürr, T., Behzadpour, S., Ellmer, M., and Zehentner, A. K. N.: The new ITSGGrace2016 release, Geophys. Res. Abstr., EGU201611547, EGU General Assembly 2016, Vienna, Austria, https://doi.org/10.13140/rg.2.1.1856.7280, 2016. a
Ko, U.D., Tapley, B., Ries, J., and Bettadpur, S.: HighFrequency Noise in the Gravity Recovery and Climate Experiment Intersatellite Ranging System, J. Spacecraft Rockets, 49, 1163–1173, https://doi.org/10.2514/1.a32141, 2012. a, b
Koch, K.R.: Parameter Estimation and Hypothesis Testing in Linear Models, Springer, Berlin, Heidelberg, https://doi.org/10.1007/9783662039762, 1999. a
Mallat, S.: A theory for multiresolution signal decomposition: the wavelet representation, IEEE T. Pattern Anal., 11, 674–693, https://doi.org/10.1109/34.192463, 1989. a, b, c
MayerGürr, T.: Estimation of error covariance functions in satellite gravimetry, IAG General Assembly 2013, Potsdam, Germany, 2013. 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., Jäggi, A., Meyer, U., Prange, L., and Maier, A.: The combined satellite gravity field model GOCO05s, Geophys. Res. Abstr., EGU201512364, EGU General Assembly 2015, Vienna, Austria, https://doi.org/10.13140/rg.2.1.4688.6807, 2015. a
MayerGürr, T., Behzadpour, S., Ellmer, M., Kvas, A., Klinger, B., and Zehentner, N.: The new ITSGGrace2016 release, https://doi.org/10.5880/icgem.2016.007, GFZ Data Services, 2016. a, b, c
Meyer, Y.: Wavelets and Operators, Vol. 1, Cambridge University Press, Cambridge, https://doi.org/10.1017/CBO9780511623820, 1993. a, b
Montenbruck, O. and Gill, E.: Satellite Orbits, Springer Berlin Heidelberg, https://doi.org/10.1007/9783642583513, 2000. a
Petit, G. and Luzum, B. (Eds.): IERS Conventions (2010), Verlag des Bundesamts für Kartographie und Geodäsie, Frankfurt am Main, 2010. a
Savcenko, R. and Bosch, W.: EOT11A – Empirical Ocean Tide Model from MultiMission Satellite Altimetry, Deutsches Geodätisches Forschungsinstitut (DGFI), München, https://doi.org/10.1594/PANGAEA.834232, 2012. a, b
Tapley, B., Bettadpur, S., Watkins, M., and Reigber, C.: The gravity recovery and climate experiment: Mission overview and early results, Geophys. Res. Lett., 31, L09607, https://doi.org/10.1029/2004gl019920, 2004. a
van Dam, T. and Ray, R. D.: S1 and S2 Atmospheric Tide Loading Effects for Geodetic Applications, available at: http://geophy.uni.lu/ggfcatmosphere/tideloadingcalculator.html (last access: 9 August 2019), 2010. a, b, c
Vetterli, M. and Herley, C.: Wavelets and filter banks: theory and design, IEEE T. Signal Proces., 40, 2207–2232, https://doi.org/10.1109/78.157221, 1992. a
Zehentner, N. and MayerGürr, T.: Kinematic orbits for GRACE and GOCE based on raw GPS observations, IAG General Assembly, Potsdam, Germany, 2013. a