Analysis of Inter-satellite Ranging Precision for Gravity Recovery in a Satellite Gravimetry Mission

Abstract

In a satellite gravimetry mission similar to GRACE, the precision of inter-satellite ranging is one of the key factors affecting the quality of gravity field recovery. In this paper, the impact of ranging precision on the accuracy of recovered geopotential coefficients is analyzed. Simulated precise orbit determination (POD) data and inter-satellite range data of formation-flying satellites containing white noise were generated, and geopotential coefficients were recovered from these simulated data sets using the crude acceleration approach. The accuracy of the recovered coefficients was quantitatively compared between data sets encompassing different ranging precisions. From this analysis, a rough prediction of the accuracy of geopotential coefficients could be obtained from the hypothetical mission. For a given POD precision, a ranging measurement precision that matches the POD precision was determined. Since the purpose of adopting inter-satellite ranging in a gravimetry mission is to overcome the imprecision of determining orbits, ranging measurements should be more precise than POD. For that reason, it can be concluded that this critical ranging precision matching the POD precision can serve as the minimum precision requirement for an on-board ranging device. Although the result obtained herein is about a very particular case, this methodology can also be applied in cases where different parameters are used.

1. INTRODUCTION

Mapping the gravitational field of the Earth precisely is essential in many subfields of geoscience such as hydrology. For example, a map of water storage or glacier mass can be constructed from the knowledge of the precise temporal gravity field at that time. Thus, by tracking changes in the time-variable gravity field, we can identify changes in how mass is distributed on the Earth.

Satellite gravimetry, in which the Earth’s gravity field is determined from tracking satellites’ orbits, is a preferred method for gravity field recovery. Commonly used geopotential models such as EGM96 and JGM-3 employed tracking data of multiple satellites. There was a significant advancement in the resolution and precision of satellite gravimetry with the advent of missions designed exclusively for gravimetry, such as GOCE (Floberghagen et al. 2011) and GRACE (Tapley et al. 2004). GRACE sought to make an improvement by measuring timevarying range values between two identical satellites that orbit the same low-Earth orbit with a separation distance of about 200 km. The uneven mass distribution of the Earth causes the inter-satellite range to fluctuate as the twins fly over different parts of the Earth, and conversely the Earth’s gravity signal can be extracted from these time-varying inter-satellite range values, or their time derivatives (i.e., range rate and range acceleration). These range-related values, measured by K-band ranging, could be determined considerably more precisely than their orbits, and their high precision is the primary reason for the improvement in the quality of geopotential recovery. A number of geopotential models generated only from GRACE data, including GGM02S (Tapley et al. 2005) and EIGENGRACE02S (Reigber et al. 2005), have already been presented to the public since the public release of GRACE data. Most of these geopotential models mainly made use of GRACE range rate data in the functional model. The huge success of GRACE led to GRAIL (Hoffman 2009), a lunar gravimetry mission that was successfully accomplished with the same twin-satellite approach, and GRACE-FO (Flechtner et al. 2012), a successor mission to GRACE.

Autonomous Spacecraft Test Environment for Rendezvous In proximity (ASTERIX), composed of two spacecraft formation flying simulators and a testbed for the simulators, has been developed by Astrodynamics and Control Laboratory of Yonsei University (Eun et al. 2018). A femtosecond laser device that can be installed on the simulators to measure the range between the simulators is also being developed. Although laser ranging has been used extensively for satellite laser ranging (SLR) (Jo et al. 2011; Lim et al. 2011; Choi et al. 2014; Kim et al. 2015; Oh et al. 2017), the technique can also be employed for measuring an inter-satellite distance (Jung et al. 2012; Lee et al. 2018). In this study, it was assumed that the satellite navigation algorithms were properly performed to figure out satellite’s position and velocity vectors (Oh et al. 2016; Shin et al. 2016; Kim et al. 2017; Lee et al. 2017a, b). One possible scientific mission that could be achieved with such configuration is a satellite gravimetry mission akin to GRACE with formation flying. Although there are a lot of factors that could potentially affect the accuracy of geopotential models produced from this hypothetical mission’s data, a crude prediction can be made with the knowledge of two main factors: the laser ranging measurement precision and precise orbit determination (POD) precision.

Various methods to recover the gravitational field of the Earth from low-low satellite-to-satellite tracking data have been proposed and developed. These approaches can be classified into three major categories, and numerous variations exist in each category. The first category is the variational equations approach. The short-arc approach can also be considered a variation of this approach. Although these approaches demand high computational costs as they require a large number of orbit integrations, several high-quality geopotential models from GRACE measurements have been generated using these two approaches, including GGM05S (Ries et al. 2016) and ITSG-Grace2016 (Mayer-Gürr et al. 2016). The second category, the energy balance approach, is based on the law of conservation of energy. Han (2004) and Shang et al. (2015) were able to determine the gravity field using this approach. The third category is the acceleration approach, which is based on Newton’s second law of motion. One of its theoretical advantages is efficient computation owing to the linearity of its functional model. In practice, however, there are several drawbacks to this approach that make the method unfavorable for actual gravity recovery from real-life mission data. For example, some of the quantities required to utilize the method can be measured only with low precision or are not available at all for GRACE (Weigelt 2017). In spite of these disadvantages, we used the crude range acceleration approach suggested by Thomas (1999) in this paper. This approach is the simplest form among the variations within the acceleration approach. As stated, this approach is not to be used to compute precise geopotential coefficients from actual satellite data. Nevertheless, because it is not our intention to recover highly accurate geopotential coefficients based on actual satellite data, but to see the impact of ranging precision on coefficients, this crude approach sufficed for our purpose.

From this study, it is possible to make a rough estimate on the accuracy of geopotential models that can be generated with the data of this hypothetical mission. Additionally, we can estimate the minimum measurement precision that its onboard ranging device must satisfy for the range measurements to be meaningful. The accuracy of geopotential models derived from this approach relies heavily on both the POD precision and ranging precision, and therefore we can find the critical point at which the two precisions have the same order of influence on the accuracy of coefficients. We believe this critical ranging precision can be thought of as the minimum precision requirement. Detailed reasoning for this statement will be given in Sections 2 and 4.

In this paper, we will cover the whole process of the simulation, from generating simulated data obtainable from the hypothetical ranging device and POD to extracting gravity fields from the simulated data, and then we will determine the minimum precision requirement for the ranging system. In Section 2, the theory of the crude acceleration approach is explained in detail. In Section 3, the process of generating data to be used in the simulation, including satellite ephemerides and inter-satellite ranges, is illustrated. A filter to compute range accelerations from the range information is also discussed in the same section. In the last two sections, we perform geopotential recovery using the simulated data and discuss the results.

2. METHODS

Among numerous methods that have been developed to recover gravity field information from low-low satellite-tosatellite tracking and inter-satellite range data, we employed a simple and crude acceleration approach by Thomas (1999). One feature of this method is that it is heavily dependent on the velocity vector difference between the satellites in formation flying. The fact that this vector can be determined only by POD is what makes this approach less appropriate for real-life application compared to other approaches. In general, it is more desirable to use methods that rely more on precise ranging data than on less precise POD data. If we think the other way around, however, this method can be used to determine the minimum precision that a ranging device must satisfy: if the ranging device is imprecise to the point that the accuracy of geopotential models is limited by the imprecise range information even for this relatively highly POD-dependent approach, it would be useless to adopt other frequently used recovery approaches that rely more heavily on the range information. With that in mind, for a given POD precision, we can identify the lowest ranging precision where the accuracy of geopotential models is limited by the ranging precision rather than the POD precision, and we will call this the critical ranging precision. Installing a ranging device with a precision lower than the critical precision will result in gravitational field data highly hampered by its imprecise range information, which is not desirable. On the other hand, if we place a device with a precision higher than the critical precision, a more accurate recovery can be made by adopting other delicate recovery methods that are less dependent on the velocity vector difference or on other information obtainable from POD. For this reason, in a sense we can consider this critical precision as the minimum requirement for the ranging system.

2.1 The Crude Acceleration Approach

The crude acceleration approach (Thomas 1999) begins with Newton’s second law of motion. If we assume that the only force acting on the twin satellites is the Earth’s gravity, the acceleration of each satellite can be computed as a gradient of the Earth’s gravitational potential, as in Eq. (1).

\(\mathbf{a}=\nabla \mathrm{U}=\frac{\partial \mathrm{U}}{\partial \mathrm{r}} \widehat{\mathrm{u}}_{\mathrm{r}}+\frac{1}{\mathrm{r}} \frac{\partial \mathrm{U}}{\partial \phi} \widehat{\mathrm{u}}_{\phi}+\frac{1}{\mathrm{r} \cos \theta} \frac{\partial \mathrm{U}}{\partial \lambda} \widehat{\mathrm{u}_{\lambda}}\)       (1)

The Earth's gravitational potential U is given by Eq. (2) (Vallado 2013). In the equation, r, ϕ, and λ denote the distance between the satellite and the Earth’s center of mass, the geocentric latitude, and the longitude of the satellite, respectively. Here, μ is the Earth’s standard gravitational parameter, and R_⨁ is the equatorial radius of the Earth. C_l,m and S_l,m are the normalized geopotential coefficients of degree l and order m. P_l,m(x) is the associated Legendre polynomial of degree l and order m, and N_l,m is the Kaula normalization factor given by Eq. (3).

\(\mathrm{U}(\mathrm{r}, \phi, \lambda)=\frac{\mu}{\mathrm{R}_{\oplus}} \sum_{i=0}^{\infty}\left(\frac{\mathrm{R}_{\oplus}}{\mathrm{r}}\right)^{1+1} \sum_{\mathrm{m}=0}^{1}\left[\mathrm{N}_{1, \mathrm{m}} \mathrm{P}_{1, \mathrm{m}}(\sin \phi)\left\{\mathrm{C}_{1, \mathrm{m}} \cos (\mathrm{m} \lambda)+\mathrm{S}_{1, \mathrm{m}} \sin (\mathrm{m} \lambda)\right\}\right]\)       (2)

\(\mathrm{N}_{\mathrm{l}, \mathrm{m}}=\left\{\begin{array}{ll} \sqrt{\frac{(1-\mathrm{m}) !(2 \mathrm{l}+1)}{(\mathrm{l}+\mathrm{m}) !}}& (\mathrm{m} \neq 0) \\ \sqrt{\frac{2(1-\mathrm{m}) !(2 \mathrm{l}+1)}{(\mathrm{l}+\mathrm{m}) !}} & (\mathrm{m} \neq 0) \end{array}\right.\)       (3)

Denoting the leading satellite as p and the trailing satellite as q, the inter-satellite (scalar) range R is given by Eq. (4). For simplicity, we will often drop the ‘inter-satellite’ from now on. In the equation, R≡x_p-x_q is the position vector from the trailing satellite to the leading satellite.

\(\mathrm{R}=|\mathbf{R}|=\mathbf{R} \cdot \widehat{\mathbf{R}}\)       (4)

Differentiating R twice with respect to time, we can express the (scalar) range acceleration R̈ in terms of the satellites’ position and velocity vectors as in Eq. (5), where the definition of Ṙ_⊥ is given in Eq. (6).

\(\ddot{\mathrm{R}}=\ddot{\mathbf{R}} \cdot \hat{\mathbf{R}}+\dot{\mathbf{R}} \cdot \frac{\mathrm{d}}{\mathrm{dt}} \widehat{\mathbf{R}}=\ddot{\mathbf{R}} \cdot \hat{\mathbf{R}}+\frac{\dot{\mathbf{R}}_{\perp} \cdot \dot{\mathbf{R}}_{\perp}}{\mathrm{R}}\)       (5)

\(\dot{\mathbf{R}}_{\perp} \equiv \dot{\mathbf{R}}-(\dot{\mathbf{R}} \cdot \widehat{\mathbf{R}}) \widehat{\mathbf{R}}\)       (6)

Our objective was to determine the geopotential coefficients C_l,m and S_l,m from time-varying range accelerations R̈. To do so, we exploited the fact that the first term on the right side of Eq. (5), R̈·R̂, can easily be computed by Eq. (7) if we know both satellites’ positions and velocity vectors as well as range accelerations, which contain the geopotential coefficients. We would like to note that in Eq. (7), R̈ is influenced by the ranging precision, and Ṙ ⊥ by the POD precision. This is the reason this approach is heavily dependent on both precisions.

\(\ddot{\mathbf{R}} \cdot \widehat{\mathbf{R}}=\ddot{\mathrm{R}}-\frac{\dot{\mathbf{R}}_{\perp} \cdot \dot{\mathbf{R}}_{\perp}}{\mathrm{R}}\)       (7)

We can see how R̈·R̂ contains the geopotential coefficients by expressing it in terms of infinite sums involving partial accelerations of degree l and order m, as in Eq. (8). The partial acceleration differences (divided by the corresponding coefficients, to be precise) A_0,0, A_l,m,1, and A_l,m,2 are defined in Eq. (9). In Eq. (9), partial gravitational potential Ũ_0,0, Ũ_l,m,1, and Ũ_l,m,2 are defined in Eq. (10). It can easily be seen that the total gravitational potential U can be represented as the infinite sum of Ũ_l,m,1 (over l≥2, 0≤m≤l), Ũ_l,m,2 (over l≥2 and 1≤m≤l), and Ũ_0,0.

\(\begin{array}{l} \ddot{\mathbf{R}} \cdot \widehat{\mathbf{R}}=\mathbf{A}_{0,0} \cdot \widehat{\mathbf{R}} \\ +\left(C_{2,0} \mathbf{A}_{2,0,1}\right) \cdot \widehat{\mathbf{R}}+\left(\mathrm{C}_{2,1} \mathbf{A}_{2,1,1}\right) \cdot \widehat{\mathbf{R}}+\left(\mathrm{C}_{2,2} \mathbf{A}_{2,2,1}\right) \cdot \widehat{\mathbf{R}} \\ +\left(\mathrm{S}_{2,1} \mathbf{A}_{2,1,2}\right) \cdot \widehat{\mathbf{R}}+\left(S_{2,2} \mathbf{A}_{2,2,2}\right) \cdot \widehat{\mathbf{R}} \\ +\left(\mathrm{C}_{3,0} \mathbf{A}_{3,0,1}\right) \cdot \widehat{\mathbf{R}}+\left(\mathrm{C}_{3,1} \mathbf{A}_{3,1,1}\right) \cdot \widehat{\mathbf{R}} \\ +\cdots \end{array}\)       (8)

\(\begin{aligned} \mathbf{A}_{0,0} & \equiv\left(\nabla \widetilde{U}_{0,0}\right)_{\mathrm{q}}-\left(\nabla \widetilde{\mathrm{U}}_{0,0}\right)_{\mathrm{p}} \\ \mathbf{A}_{1, \mathrm{m}, 1} & \equiv \frac{\left(\nabla \widetilde{U}_{1, \mathrm{m}, 1}\right)_{\mathrm{q}}-\left(\nabla \widetilde{U}_{1, \mathrm{m}, 1}\right)_{\mathrm{p}}}{\mathrm{C}_{\mathrm{l}, \mathrm{m}}} \\ \mathbf{A}_{1, \mathrm{m}, 2} & \equiv \frac{\left(\nabla \widetilde{\mathrm{U}}_{1, \mathrm{m}, 2}\right)_{\mathrm{q}}-\left(\nabla \widetilde{\mathrm{U}}_{1, \mathrm{m}, 2}\right)_{\mathrm{p}}}{\mathrm{S}_{\mathrm{l}, \mathrm{m}}} \end{aligned}\)       (9)

\(\begin{aligned} \widetilde{\mathrm{U}}_{0,0} & \equiv \frac{\mu}{\mathrm{r}} \\ \widetilde{\mathrm{U}}_{1, \mathrm{m}, 1} & \equiv \frac{\mu}{\mathrm{R}_{\oplus}}\left(\frac{\mathrm{R}_{\oplus}}{\mathrm{r}}\right)^{1+1} \mathrm{N}_{1, \mathrm{m}} \mathrm{P}_{1, \mathrm{m}}(\sin \phi) \mathrm{C}_{\mathrm{l}, \mathrm{m}} \cos (\mathrm{m} \lambda) \\ \widetilde{\mathrm{U}}_{1, \mathrm{m}, 2} & \equiv \frac{\mu}{\mathrm{R}_{\oplus}}\left(\frac{\mathrm{R}_{\oplus}}{\mathrm{r}}\right)^{1+1} \mathrm{N}_{1, \mathrm{m}} \mathrm{P}_{1, \mathrm{m}}(\sin \phi) \mathrm{S}_{\mathrm{l}, \mathrm{m}} \sin (\mathrm{m} \lambda) \end{aligned}\)       (10)

It can be shown that A_0,0, A_I,m,1, A_I,m,2 and R̂ can be computed if we can determine both satellites’ positions. One thing to be noted here is that gradients computed in spherical coordinates result in vectors in position-specific coordinates; therefore, we must rotate them appropriately to express them in the same frame before the subtractions in Eq. (9) can be done. Furthermore, as we need to compute dot products of A_0,0, A_l,m,1, A_l,m,2, and R̂, all these vectors must be expressed in the same frame, such as the Earth-centeredEarth-fixed (ECEF) frame.

Going back to Eq. (7), we have already seen that R̈·R̂ can be computed easily. Accordingly, all terms on both sides of Eq. (8) are known with the exception of the geopotential coefficients, which can now be determined by the ordinary least squares method. We can write Eq. (8) for time 1 to n in matrix form

\(\mathbf{y}=\mathrm{Ap}\)       (11)

where

\(\mathbf{y}=\left[\begin{array}{l} \left\{\left(\ddot{\mathbf{R}}-\mathbf{A}_{0,0}\right) \cdot \widehat{\mathbf{R}}\right\}_{1} \\ \left\{\left(\ddot{\mathbf{R}}-\mathbf{A}_{0,0}\right) \cdot \widehat{\mathbf{R}}\right\}_{2} \\ \left\{\left(\ddot{\mathbf{R}}-\mathbf{A}_{0,0}\right) \cdot \widehat{\mathbf{R}}\right\}_{3} \\ \vdots \\ \left\{\left(\ddot{\mathbf{R}}-\mathbf{A}_{0,0}\right) \cdot \widehat{\mathbf{R}}\right\}_{n} \end{array}\right]\)       (12)

\(A=\left[\begin{array}{cccccc} \left(\mathbf{A}_{2,0,1} \cdot \widehat{\mathbf{R}}\right)_{1} & \left(\mathbf{A}_{2,1,1} \cdot \widehat{\mathbf{R}}\right)_{1} & \left(\mathbf{A}_{2,2,1} \cdot \widehat{\mathbf{R}}\right)_{1} & \cdots & \left(\mathbf{A}_{1, \mathrm{m}, 2} \cdot \widehat{\mathbf{R}}\right)_{1} & 1 \\ \left(\mathbf{A}_{2,0,1} \cdot \widehat{\mathbf{R}}\right)_{2} & \left(\mathbf{A}_{2,1,1} \cdot \widehat{\mathbf{R}}\right)_{2} & \left(\mathbf{A}_{2,2,1} \cdot \widehat{\mathbf{R}}\right)_{2} & \cdots & \left(\mathbf{A}_{1, \mathrm{m}, 2} \cdot \widehat{\mathbf{R}}\right)_{2} & 1 \\ \left(\mathbf{A}_{2,0,1} \cdot \widehat{\mathbf{R}}\right)_{3} & \left(\mathbf{A}_{2,1,1} \cdot \widehat{\mathbf{R}}\right)_{3} & \left(\mathbf{A}_{2,2,1} \cdot \widehat{\mathbf{R}}\right)_{3} & \cdots & \left(\mathbf{A}_{1, \mathrm{m}, 2} \cdot \widehat{\mathbf{R}}\right)_{3} & 1 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \left(\mathbf{A}_{2,0,1} \cdot \widehat{\mathbf{R}}\right)_{\mathbf{n}} & \left(\mathbf{A}_{2,1,1} \cdot \widehat{\mathbf{R}}\right)_{\mathrm{n}} & \left(\mathbf{A}_{2,2,1} \cdot \widehat{\mathbf{R}}\right)_{\mathrm{n}} & \cdots & \left(\mathbf{A}_{1, \mathrm{m}, 2} \cdot \widehat{\mathbf{R}}\right)_{\mathrm{n}} & 1 \end{array}\right]\)       (13)

\(\mathbf{p}=\left[\begin{array}{c} \mathrm{C}_{2,0} \\ \mathrm{C}_{2,1} \\ \mathrm{C}_{2,2} \\ \vdots \\ \mathrm{S}_{1, \mathrm{m}} \\ \mathrm{k} \end{array}\right]\)       (14)

A and p should be appropriately modified if we know the error model of range measurement, and then the error model parameters can also be estimated. For example, in Eq. (13) the last column of A, composed of 1’s, is added so that we can estimate the constant bias k in range data, although this is not an essential step for our study since we did not add a constant bias while making noisy range data.

As already stated, the values of R̈·R̂ and A_0,0·R̂ can be calculated at each time if we know the satellites’ positions and velocity vectors at each time, so y is known. As A_l,m,1·R̂ and A_l,m,2·R̂ can also be computed at each time if we know the satellites’ position vectors at each time, A is also known. The term p, which consists of the normalized geopotential coefficients up to degree l, is the only unknown. Therefore, we can estimate the coefficients p by the ordinary least squares, as in Eq. (15). The hat implies that it is an estimated value.

\(\widehat{\mathbf{p}}=\left(A^{\mathrm{T}} A\right)^{-1} A^{\mathrm{T}} \mathbf{y}\)       (15)

3. SIMULATION DATA

3.1 True Data and Noisy Data Generation

The General Mission Analysis Tool (GMAT) (Hughes 2016), open-source software developed by NASA, was used to generate simulation data. The satellites were deemed to be on coplanar, nearly polar, and nearly circular orbits, and their orbits were propagated for durations of about 15 days (1,296,120 sec to be precise). The orbital parameters of the twin satellites at epoch are listed in Table 1. The epoch time was chosen arbitrarily, but the position and velocity vectors at epoch are the actual values of the GRACE twin satellites at the time. At epoch, the initial distance between them is about 175 km. For propagation, it was assumed that the only force acting on the satellites was the Earth’s gravity. The EGM96 gravity model (Lemoine et al. 1998) truncated at degree 180 was used. In reality, other perturbations such as air drag and solar radiation pressure exist. Therefore, such perturbations should be directly measured or assumed by appropriate modeling, and then be taken into account during the recovery process. For instance, GRACE satellites have accelerometers that can measure all kinds of nongravitational accelerations (Flury et al. 2008).

Table 1. Orbital parameters of the satellites at epoch (February 21, 2011 23:58:45:000 UTC)

The output of the simulation, composed of the positions and velocities (in the Earth-centered J2000 equatorial frame) of both satellites at an interval of 0.01 sec, can be thought of as the ‘true’ orbits they are on. Also computable from this output are the ‘true’ ranges for the same time interval. Then we contaminated the true data to generate noisy data, which can be considered a proxy of the imprecise data we would have from POD and the on-board ranging system in a reallife mission. When GPS data is used in the POD process, it is natural that relative position and velocity between the satellites can be determined more precisely than with absolute values. With that in mind, the position and velocity values were contaminated with an appropriate amount of Gaussian noise. The true range values were also contaminated with Gaussian noise. The Gaussian noise of the position and velocity values corresponds to the precision of POD, and that of the range values to the precision of the laser ranging system. The Gaussian noise sigmas used for our simulation are given in Table 2. For each range noise sigma, we made 10 sets of random simulation data, so in total, 40 data sets were made.

Table 2. Gaussian noise sigmas used to simulate the noisy data

3.2 Filtering

The measurements we receive from the laser ranging system are inter-satellite ranges. If we are to use the crude acceleration approach, however, we need to know scalar range accelerations, denoted as r̈ in Eq. (7). To derive range acceleration values from range values, we used the CRN-class digital filter, which was also developed by Thomas (1999) and is currently used in processing GRACE mission inter-satellite range data (Wu 2006). By applying the filter to high-rate range data, we can obtain filtered high-rate range, range rate, and range acceleration data. There are several filter parameters, and we selected them so that our filtered result shares the same fundamental properties (e.g., the low-pass bandwidth and fit interval time) with GRACE Level 1B data. The filter parameters we used and the original parameters used to produce GRACE Level 1B data, are given in Table 3. The differences in the parameters are mainly due to their different high-rate raw data sampling frequencies: we assumed our laser ranging system measured ranges at 100 Hz, while the GRACE K-band ranging system operated at 10 Hz. More information about the filter and its parameters can be found in Thomas (1999) and Wu et al. (2006). Following the filtering process, low-rate (0.2 Hz) data was made by collecting the (filtered) high-rate values at intervals of 5 sec. Fig. 1 shows the low-rate range, range rate, and range acceleration values obtained by applying filters to noiseless high-rate range data.

Table 3. Parameters for the CRN-class digital filter used in GRACE data processing and our simulation

Fig. 1. Low-rate inter-satellite range (top), range rate (middle), and range acceleration (bottom) as a function of hours since epoch, obtained from true high-rate range values.​​​​​​​

A flowchart illustrating the process of making simulation data is shown in Fig. 2. First, the ‘true’ high-rate (100 Hz) position and velocity vectors (in the inertial J2000 frame) of each satellite and the ‘true’ high-rate range values were obtained by propagating their orbits with the help of GMAT. White noise was added to the true data to simulate noisy data (see Section 3.1). Then, the noisy positions and velocity vectors were rotated from the inertial frame to the ECEF frame in order that they could be used in the crude acceleration approach. Also, the noisy range data was filtered by the CRNclass digital filter to produce range acceleration data (see Section 3.2). Low-rate noisy data (0.2 Hz) was generated from the high-rate noisy data, and then geopotential coefficients were recovered from this low-rate noisy data.

Fig. 2. Flowchart for the simulation process.​​​​​​​

4. RESULTS

For each of the 40 noisy data sets, we recovered geopotential coefficients using the crude acceleration approach. It could be expected that the obtained coefficients would become more accurate as the precision of the ranging system improved, but there is a critical point at which a better ranging precision does not produce better coefficients because of the imprecision of POD. Below that ranging precision, the accuracy of geopotential coefficients is limited by imprecise ranging; above that precision, the accuracy is limited by imprecise POD. In other words, the accuracy of the coefficients can be limited by either of the precisions depending on their values. We can think of that critical ranging precision as the minimum requirement of the ranging device. A ranging precision lower than this requirement would make the whole idea of using inter-satellite ranges pointless, but we could make use of more precise range data if we adopted other advanced recovery approaches.

Since the number of geopotential coefficients up to degree and order 100 equals 10,197, it is impractical to inspect their accuracies one by one. Instead, we used three different ways to quantify or visualize the differences between the ‘true’ gravitational coefficients and the obtained ones. Fig. 3 shows the root mean square (RMS) coefficient differences per degree. The RMS coefficient difference for degree j, where j ranges from 2 to 100, can be computed by Eq. (16). In the equation, hatted quantities refer to the recovered (estimated) coefficients, and barred quantities to the true coefficients.

Fig. 3. Averaged RMS coefficient differences per degree. Light-colored lines are drawn for each data set. Bright-colored lines are drawn by averaging the ten same light-colored lines. Smaller values imply more accurate coefficients. The RMS coefficients per degree of EGM96 are plotted as a black dotted line for comparison. Note that all values are unitless.​​​​​​​

\(\mathrm{e}_{\mathrm{j}}=\sqrt{\frac{\Sigma_{\mathrm{m}=0}^{\mathrm{j}}\left(\hat{\mathrm{C}}_{\mathrm{j}, \mathrm{m}}-\overline{\mathrm{C}}_{\mathrm{j}, \mathrm{m}}\right)^{2}+\sum_{\mathrm{m}=1}^{\mathrm{j}}\left(\hat{\mathrm{S}}_{\mathrm{j}, \mathrm{m}}-\overline{\mathrm{S}}_{\mathrm{j}, \mathrm{m}}\right)^{2}}{2 \mathrm{j}+1}}\)       (16)

Then, ten lists of RMS values obtained from ten different data sets that had the same ranging precision were averaged. In Fig. 3, light-colored lines represent each list of RMS values from a single data set, and bright-colored lines represent the averaged RMS values. Except for the bright red line, all three other bright lines are indistinguishable from one another. This implies that more precise ranging does not yield more accurate coefficients when the one-sigma ranging precision is below a few millimeters.

The colors in Fig. 4 show the RMS coefficient differences per each coefficient. Ten coefficient differences for a single coefficient, each obtained from ten data sets with the same ranging precision, were root-mean-squared. As we have seen in Fig. 3, the coefficients’ accuracy starts deteriorating when the ranging precision sigma becomes bigger than a few millimeters. One feature that was not found in Fig. 3, however, is that as the one-sigma ranging precision is improved from 1 millimeter (Fig. 4(b)) to 100 micrometers (Fig. 4(c)), there is a slight but definite improvement in coefficients of higher degree and lower order. The reason this improvement was not detectable in Fig. 3 is that the coefficients of higher degree and lower order are closer to their true values when compared with the coefficients of the same degree but higher order. Thus, the improvement in those coefficients only had a minute impact on the RMS differences per degree. Still, as in Fig. 3, not much difference can be seen between the 100 micrometer case (Fig. 4(c)) and the 10 micrometer case (Fig. 4(d)).

Fig. 4. RMS coefficient differences per each coefficient for different one-sigma ranging precisions: (a) 10 mm (1e-5 km), (b) 1 mm (1e-6 km), (c) 100 μm (1e-7 km), and (d) 10 μm (1e-8 km). Smaller values imply more accurate coefficients. Note that all values are unitless.​​​​​​​

Finally, we can visualize the results in the form of geoid height differences. Geoid heights are slightly different depending on the coefficients used to compute them, and if the differences between the used coefficients and the true coefficients are small, their geoid height differences would also be small. Geoid height differences were computed with Eq. (17) (Heiskanen & Moritz 1967). In the equation, r(ϕ) and γ₀ (ϕ) (respectively) refer to the local ellipsoidal height and local normal gravity at geocentric latitude ϕ, computed with the constants from WGS84. Note that the maximum degree and order used to compute geoid height differences were truncated at 80, not 100.

\(\mathrm{N}_{\mathrm{e}}(\lambda, \phi)=\frac{\mu}{\mathrm{r}(\phi) \gamma_{0}(\phi)} \sum_{1=2}^{80}\left(\frac{\mathrm{a}}{\mathrm{r}(\phi)}\right)^{1} \sum_{\mathrm{m}=0}^{1}\left[\mathrm{N}_{\mathrm{l}, \mathrm{m}} \mathrm{P}_{1, \mathrm{m}}(\sin \phi)\left\{\left(\hat{\mathrm{C}}_{1, \mathrm{m}}-\overline{\mathrm{C}}_{\mathrm{l}, \mathrm{m}}\right) \cos (\mathrm{m} \lambda)\right.\right.\\ \left.\left.+\left(\widehat{\mathrm{S}}_{1, \mathrm{m}}-\overline{\mathrm{S}}_{\mathrm{l}, \mathrm{m}}\right) \sin (\mathrm{m} \lambda)\right\}\right]\)​​​​​       (17)

For each data set, using the corresponding recovered coefficients, geoid height differences were computed at 64,080 points over the globe: from 0 to 360 degrees with 1 degree intervals in longitude, and -89 degrees to 89 degrees with 1 degree intervals in latitude. These 64,080 geoid height differences of that specific data set were root-mean-squared, so that we can easily compare the RMS values of different data sets. The RMS geoid height differences for all data sets are listed in Table 4. In Table 4, the RMS values for ten sets sharing the same ranging precision are listed in the same row, and then the average RMS value for that ranging precision is listed in the rightmost column. Fig. 5 shows on a world map the geoid height differences, using the worst data set (i.e., the data set with the largest RMS geoid height difference) among each range noise sigma. In Fig. 5, regions with higher geoid heights compared to the ‘true’ geoid height are shown in red, and regions with lower geoid heights are in blue. If the geopotential recovery process and data were perfectly ideal, every part of the Earth must have the same geoid height as the true geoid height, making the map green everywhere. From Table 4 and Fig. 5, we can easily conclude that the geopotential coefficients are undeniably poor in accuracy when the range noise sigma is 10 millimeters (Fig. 5 (a)), compared to the other three cases with higher precision (Fig. 5(b)-5(d)). This outcome is consistent with the conclusions made from Figs. 3 and 4. In Table 4, a slight difference in the averaged RMS values can be found between the 1 millimeter case and the other two cases with higher precision. We think that this slight difference is concordant with the slight but definite improvement in high-degree and low-order coefficients that was discussed while we were interpreting Fig. 4.

Table 4. RMS geoid height differences for all data sets

¹ Smaller values imply a more accurate geopotential model.

Fig. 5. Map of geoid height differences for the worst data set among the ten sets sharing the same range noise sigma of (a) 10 mm (1e-5 km), (b) 1 mm (1e-6 km), (c) 100 μm (1e-7 km), and (d) 10 μm (1e-8 km).​​​​​​​

From what we have discussed, we conclude that the critical ranging precision (one-sigma) lies between several hundred micrometers and a few millimeters, and that this critical precision can be assumed to be the minimum requirement of the ranging system. One could expect that if the POD precision were improved, this minimum requirement would become smaller (more stringent) so that the critical ranging precision could be matched with the improved POD precision. To check if this assumption is correct, the same procedures were performed again with 40 totally different data sets in which the data were perturbed with ten times more precise POD than the original case (e.g., the absolute position 3D noise sigma was 1e-5 km instead of 1e-4 km). The result, shown in the same format as in Fig. 3, is displayed in Fig. 6. There are two points to be mentioned about Fig. 6 in comparison with Fig. 3. One is the overall improvement in the accuracy of the geopotential models obtained. The result itself it not something to be surprised about, because there was an improvement in the POD precision while the ranging precision was kept at the same level. However, one specific aspect to be noted is the different amount of improvement in the accuracy of geopotential models. Only a slight improvement was found for the 10 millimeter case (red lines), while a major improvement of about one order of magnitude can be seen in the other three cases. In the 10 millimeter case, the accuracy was chiefly limited by the imprecise range information, so the introduction of more precise POD could not bring about a considerable amount of change. In contrast, in the other three cases with higher ranging precision, the accuracy was primarily limited by the imprecise POD, and consequently the enhanced POD resulted in a significant improvement. The other point to be mentioned is the stricter minimum precision requirement, as we predicted. In Fig. 6, unlike in Fig. 3, the bright green line that corresponds to the 1 millimeter case can be distinguished from the blue and purple lines, implying that the critical ranging measurement precision (one-sigma) is a few hundred micrometers, which is smaller than that of the original (less precise POD) cases.

Fig. 6. Averaged RMS coefficient differences per degree, for the improved POD case. See Fig. 3 for a more detailed description.​​​​​​​

5. CONCLUSIONS

In this study, geopotential coefficient recovery was performed on data sets generated and processed by assuming a hypothetical GRACE-like mission. The crude acceleration approach was used in recovering coefficients. We were able to make a rough guess of the accuracy of the coefficients obtainable from the hypothetical mission. During the process, we were also able to estimate the minimum requirement for the on-board ranging system to be adopted in such mission.

Assuming the POD precision given in Table 2, we conclude that the minimum one-sigma precision requirement for the ranging device is between several hundred micrometers and a few millimeters. It should be emphasized that this is the minimum requirement, not the ultimate one. We also conclude that the expected accuracy of the normalized geopotential coefficients is approximately between the order of 10^-9 and 10^-10, if the ranging precision satisfies the minimum requirement we just mentioned. Nevertheless, this result is a rough figure and there is room for improvement if we employ more delicate techniques that can take advantage of highly precise range information. On the other hand, there are some factors we have not considered in our simulation that negatively affect the recovery process. For instance, if it is decided to install an on-board accelerometer to directly measure the non-gravitational accelerations, errors in the accelerometer measurements would have an adverse effect on the accuracy of geopotential recovery. It should also be noted that there are numerous other parameters that could also affect the quality of the resulting coefficients, such as the initial distance between the satellites (Wei et al. 2012), the amount of data used in recovery, the maximum degree and order of coefficients to be recovered (Gunter et al. 2006), the sampling frequency of the ranging system, and the type of filter used in processing the raw data. Thus, the impacts of all such factors should be thoroughly investigated before an actual mission is attempted. Moreover, since the parameters used in this study are fixed at specific values, the result obtained herein is applicable to a very limited number of cases. Still, the methods we have followed can be applied to other cases in which different parameters are used, or possibly to parametric studies of different parameters.