Simulation analysis of recovering time-varying gravity fields based on Starlink-like constellation

SUMMARY

Time-varying gravity fields play a crucial role in understanding and analysing geodynamic processes, particularly the migration of matter across the Earth's surface. However, the current limitations in spatiotemporal resolution hinder their accurate representation. In this context, the use of a giant constellation of low-orbit satellites holds great potential for accurately recovering time-varying gravity fields with high spatiotemporal resolution. Based on the orbital parameters of 5199 satellites in 123 different orbital planes in the first phase configuration of the Starlink constellation and the orbital parameters of the Bender constellation in the next generation gravity mission, we conducted a closed-loop simulation to analyse the recovery ability of time-varying gravity field in 9 d using the short-arc integral method. The errors of aliasing AOHIS signal (Atmosphere, Ocean, Hydrology, Ice and Solid Earth), ocean tide models, orbit positions, intersatellite range rates and accelerometer observations were considered in the numerical simulation. Compared with the Bender constellation, the Starlink-like constellation can effectively decrease the aliasing errors in the spatial- and frequency domain when the observation noise is not considered. The Starlink-like constellation can also effectively improve the reliability of low-degree coefficients (below degree 15) of retrieved time-varying gravity field models and present higher time resolution (within 9 d) for the full-degree spherical harmonic solutions than the Bender constellation when the observation noise is considered. The aliasing effect on the low-degree part of the Bender constellation can be significantly decreased by combining the Starlink-like and Bender constellations, and the accuracy of the recovered time-varying gravity field within degree 30 can be improved by about 0.5–1 order of magnitude. Our results can provide a technical reference for the design of future gravity satellite mission.

1 INTRODUCTION

The successful implementation of the Gravity Recovery and Climate Experiment and its Follow-on satellite mission (GRACE/GRACE-FO) has greatly promoted the understanding of surface mass changes in Earth system (Tapley et al. 2004a,b, 2019; Kornfeld et al. 2019). For the first time, the low-low Satellite-to-satellite Tracking (ll-SST) measurement mode was adopted by GRACE/GRACE-FO, which has accumulated a large amount of observations during their 20 yr of operation and recovered several high-precision static gravity field models (Akyilmaz et al. 2016; Zhou et al. 2017; Chen et al. 2018; Mayer-Gürr et al. 2018). In addition, the monthly observations from the GRACE/GRACE-FO satellites was used to obtain time-variable gravity field products with spatial resolution of approximately 300 km, playing an irreplaceable role in studying changes in terrestrial water storage (Feng et al. 2018; Xu et al. 2023), glacier mass (Luthcke et al. 2008; Chen et al. 2022), post-seismic deformation (Liang et al. 2018; Ghobadi-Far et al. 2020) and other aspects (Chao 2015; Li et al. 2023).

However, the gravity field models retrieved from GRACE/GRACE-FO observations are still subject to limitations due to measurement accuracy, anisotropy of observations and aliasing errors (Wiese et al. 2011; Daras & Pail 2017; Purkhauser & Pail 2019). The intersatellite observations of GRACE/GRACE-FO reflect the orbital differences between the two satellites in the line-of-sight direction, making them more sensitive in the along-track direction. Consequently, this introduces an anisotropic error spectrum and north–south striping characteristics (Daras & Pail 2017). Furthermore, the satellite's inability to revisit the same area of Earth in a short period results in insufficient time sampling of time-varying signal, leading to aliasing effects. De-aliasing and ocean tide model errors can similarly impact the accuracy of the gravity field model in the form of aliasing errors (Zenner et al. 2012). As measurement technology continues to advance, future satellite gravity missions are anticipated to utilize higher precision instruments and improved technologies, such as drag-free control systems. Under these advancements, it is expected that the measurement accuracy would no longer be the primary limiting factor (Wiese et al. 2011), while the research focus of future satellite gravity missions has shifted towards addressing the aliasing problem and the observed anisotropy.

Previous publications have extensively studied the observed anisotropy and aliasing errors of GRACE/GRACE-FO. For instance, Bender et al. (2008) proposed a new constellation configuration (called Bender configuration), which adds a satellite pair on an inclined orbit to a pair of polar-orbiting satellites to improve the spatial sampling of geophysical signal. This observation method not only ensures nearly global coverage of satellite observations but also increases the observation of the cross-orbit direction component of the Earth's gravity field, as well as significantly reduces the north–south striping error existing in GRACE/GRACE-FO. The Bender constellation has been widely discussed and studied in the Next Generation Gravity Mission (NGGM), which aims to improve the monitoring of changes in Earth's gravity field by improving spatiotemporal sampling capabilities and using a new generation of higher precision instruments (Cesare & Sechi 2013). According to the requirements of the NGGM, the Additional Constellation and Scientific Analysis of the Next Generation Gravity Mission Concept (ADDCON) project extensively studied the inversion results of the gravity field using the Bender constellation under different orbital parameters and conducted a detailed analysis of their performance (Pail et al. 2018).

At present, there are three methods used to address the aliasing problem (Wiese et al. 2011). First, increasing the number of satellites to improve the spatiotemporal sampling of time-varying signal in the gravity field. Second, improving the accuracy of atmospheric, oceanic, and tidal models to mitigate the impact of model errors. The third approach involves reducing the time-aliasing errors by co-estimating the parameters. Wiese et al. (2011) proposed a method to estimate the gravity field at high frequency and low resolution, which has been widely used in NGGM research as a simple and effective de-aliasing method. However, because the periods of some ocean tide models are shorter than one day, daily estimation alone is not adequate to reduce the aliasing effect caused by ocean tide model errors. These errors remain significant contributors to aliasing errors even after the application of Wiese's method (Daras & Pail 2017). To address the aliasing problem attributed to ocean tide models, Zhao et al. (2015) investigated the performance of different satellite formations in reconstructing the gravity field and evaluated the feasibility of various formations in mitigating aliasing errors in ocean tide models. Hauk & Pail (2018) conducted a co-estimation of the parameters of ocean tides and utilized the estimated ocean tide model for de-aliasing, leading to a significant reduction in the aliasing errors of the tidal model. Additionally, Liu & Sneeuw (2021) elaborated the two-step aliasing mechanism of ocean tides and demonstrated it through closed-loop simulation. These studies have indeed reduced the impact of aliasing effects to some extent. However, the limitation of insufficient time sampling in the GRACE/Bender constellation persists, and the root cause of the aliasing effect remains. Therefore, there is a need to enhance the time sampling of time-varying signal and improve the time resolution of the gravity field. Pfaffenzeller & Pail (2023) used up to 18 pairs of GRACE/Bender-type satellites to carry out simulation experiments of different configurations to retrieve time-varying gravity field, and the results showed the advantages of a network of miniaturized satellites in reducing aliasing problems.

Giant constellations, widely used in communication, navigation, remote sensing and other fields (Xue et al. 2022; Zhang et al. 2022), also hold significant potential for retrieving time-varying gravity fields. In recent years, numerous countries have been competing to develop giant constellations. Among them, the Starlink constellation established by the United States Space Exploration Technology Company (SpaceX) has emerged as the most mature. The first phase configuration of the Starlink constellation comprises 12 196 satellites, with the majority positioned in low Earth orbit (LEO) and very-low Earth orbit (VLEO). Given the Starlink constellation's ability to rapidly cover vast ground areas, there is theoretical potential to retrieve the time-varying gravity field using orbital data within a short period, thus allowing for higher time resolution in gravity field solutions. Furthermore, increasing spatiotemporal sampling holds the promise of reducing the aliasing effect caused by time-varying signal undersampling and ocean tide model errors. To this end, this study conducted a closed-loop simulation to retrieve the Earth's nontidal geophysical processes ‘AOHIS’ [Atmosphere (A), Ocean (O), Hydrology (H), Ice (I), and Solid Earth (S)] within a 9-d period, using simulated observations from the Starlink phase I constellation configuration and the ADDCON's Bender constellation. We conducted analyses under various error sources and recovery periods (1, 3, 6, and 9 d) through simulation to evaluate the performance of different constellations in recovering time-varying signal and reducing aliasing effect. Our results would allow us to explore the respective strengths and weaknesses of the Starlink-like and Bender constellations.

2 METHOD

The gradient-corrected short-arc integral method was proposed by Mayer-Gürr (2006) to retrieve the gravity field based on high-low satellite-to-satellite tracking (hl-SST) and ll-SST observation equations. It has been successfully applied to the solution of ITG-Grace2010 (Mayer-Gürr et al. 2010) and other gravity field models. In this simulation experiment, Starlink-like constellation utilized hl-SST observation mode, while GRACE/Bender constellation used hl-SST and ll-SST observation modes. The observation equations of hl-SST and ll-SST were established based on the data of orbital position and intersatellite range rate, respectively.

By integrating the satellite acceleration based on Newton's equation of motion and using the gradient correction, we can obtain the hl-SST observation eq. (1), which utilizes the satellite's orbital position as the observation:

where

In eqs (1) and (2), r_g is the observation of satellite's orbital position, r₀ is the satellite's position after gradient correction, K is the numerical integral matrix, T is the gradient matrix of the force models, |${\bf\bar {B}}$| is the coefficient matrix of the orbital boundary vector, G is the design matrix linking the force f to the spherical harmonic coefficients, x is the spherical harmonic coefficients correction and b is the orbital boundary correction.

When the observation is the intersatellite range rate, the ll-SST observation equation can be established:

where

In eqs (3) and (4), the subscripts ‘A’ and ‘B’ distinguish the preceding and trailing satellites. |${{{{\bf r}}}_{{\rm{AB}}}}$| and |${{{{\bf \dot{r}}}}_{{\rm{AB}}}}$| represent the relative position and velocity of the two satellites. |$\dot{\rho }$| is the observation of intersatellite range rate and |${{\dot{\rho }}_0}$| is the intersatellite range rate after gradient correction.

Using the least-squares criterion and eliminating local variables in the observation equation, we can get the normal equation for the global variable x, and the potential coefficients of the gravity field can be further obtained.

When calculating multiday solutions, Wiese's method can be used to parametrize low-degree terms of the gravitational time-varying signal daily to mitigate the time-aliasing effect of the AOHIS signal (Wiese et al. 2011). This method maintains the separation of low-degree terms from high-degree terms through block inversion (Ullman 1997), makes overall estimation of high-degree terms, and substitutes them into daily normal equations to obtain low-degree terms, and finally obtains time-varying gravity fields with improved spatial resolution.

3 NUMERICAL SIMULATION

3.1 Constellation configuration

3.1.1 Starlink-like constellation

Presently, satellites in the Starlink constellation are mainly positioned in LEO and VLEO. However, their orbital data are not publicly accessible. As a result, this study utilized the Starlink phase I configuration provided by SpaceX to the US Federal Communications Commission in 2020 April (Table 1), as well as the Satellite Two Line Element (TLE) available on the CelesTrak website (http://celestrak.org), to simulate Starlink satellite orbit parameters.

The phase I configuration consists of LEO-type orbits at an altitude of approximately 550 km and VLEO-type orbits at an altitude of approximately 340 km. Currently, information about the number of orbital planes and satellites per plane for the VLEO-type orbits is unavailable, represented by the symbol ‘/’ in the table. The orbital inclination of the LEO satellite is concentrated around 97.6°, 70° and 53°. In contrast, the VLEO-type orbits exhibit a small inclination difference between the three parts, with a maximum inclination of 53°.

Because of the large number of satellites and the significant computational workload, the representative parts of the constellation were selected for simulation and calculation, and different orbital inclinations and altitudes of the Starlink constellation were considered. The constellation utilized in this study includes four parts (Table 2), consisting of 123 orbital planes and a total of 5199 satellites. Parts a, b and c are from LEO-type orbits, and part d is from VLEO-type orbits. Because of the lack of precise information, only the number of satellites is provided for part d. Given the even distribution of satellites in space and the modest number of orbital planes (1, 3, 9, 283, 849 and 2547), this study posits that the most probable scenario involves nine orbital planes, each accommodating 283 satellites.

Based on the constellation configuration presented in Table 2, the orbital parameters of satellite j on the orbital plane i in each part of the Starlink-like constellation can be calculated, as shown in Table 3. The semimajor axis is determined by the orbital altitude and Earth radius of 6 378 137 m. Analysis of the TLE data from the operational Starlink satellites reverals that the eccentricity of most satellites is approximately 0.00015, and the arguments of the perigee of parts a, b, c and d are concentrated around 90°, 265°, 90° and 90°, respectively, for the establishment of the corresponding parameters. Here, N_S denotes the number of satellites in a single orbital plane, while N_O represents the number of orbital planes.

Assuming that the satellites are uniformly distributed along the orbital plane, and that the ascending node of the satellites on the same orbital plane is the same, it is taken into account that the difference in mean anomaly between two adjacent satellites on the same orbital plane is 360/N_S, and the ascending node difference between two adjacent orbital planes is 360/N_O. These parameters are utilized to calculate the initial orbital Kepler parameters of each satellite, as illustrated in Table 3. Fig. 1 displays the initial positions of the simulated Starlink satellites and the spatial distribution of the ground track within 10 min, demonstrating that the Starlink-like constellation can offer high-density coverage of the global area, except for the poles, in a short period. The blue, red, green and magenta circles represent the initial positions of satellites in parts a, b, c, and d, respectively. The corresponding curves represent the satellites' ground track.

3.1.2 Bender constellation

The Bender constellation, which is the preferred configuration for the next generation of gravity satellite missions, consists of two pairs of satellites. One pair operates in a near-polar orbit with an inclination of approximately 90°, while the other pair operates in an inclined orbit with an inclination of approximately 70° or 110° (Pail et al. 2018). The two satellites forming a pair fly in the same orbital plane, with the distance between the satellites typically in the range 100–200 km. The primary observation is focused on the range or range rate between the satellites. In this study, the simulation of the Bender constellation was conducted based on scenario 4 in the ADDCON project, following the NGGM simulation process (Pail et al. 2019). The initial orbit parameters are shown in Table 4, with A1 and A2 representing polar orbit satellites, and B1 and B2 representing inclined orbit satellites. Additionally, for comparative analysis, a pair of GRACE-type observations composed of A1 and A2 were also used separately in the calculation.

3.2 Numerical simulation process for Starlink-like and Bender constellations

To analyse the performance of the Starlink-like and Bender constellations in recovering the time-varying gravity field, simulation experiments were conducted following the process shown in Fig. 2. First, the true models were chosen to simulate the satellite observations, and subsequently, reference models with deviations from true models were introduced. The impact of the reference models was deducted in the calculations, resulting in the derivation of residual observations to recover the time-varying signal of the gravity field. When constructing the observation equation, the gradient correction is carried out using the GOCO05s model (Mayer-Gürr 2015) within degree 5. These simulation and calculation processes are mainly based on the GROOPS software (Mayer-Gürr et al. 2021).

3.2.1 Background models

The aliasing effect primarily results from the undersampling of recovered signal (e.g. hydrological variations), de-aliasing model errors (related to non-tidal time-varying signal in the atmosphere and ocean) and errors in ocean tide models (Wiese et al. 2011; Daras & Pail 2017). To study the impact of aliasing on the recovery of time-varying signal in different constellations, the true models used in the orbit simulation process and the reference models employed in the inversion process of the time-varying gravity field were established according to Table 5.

During the simulation processes of this study, the Earth's non-tidal time-varying signal (i.e. ‘AOHIS’) was directly recovered, introducing the aliasing effect due to inadequate time sampling of AOHIS signal. To introduce ocean tide model errors, this study utilized ocean tide models from two different institutions to simulate the aliasing effect caused by such errors. Specifically, the EOT11a model (Savcenko & Bosch 2012) was chosen as the true model, while the GOT4.7 model (Ray 1999) was selected as the reference model. All of these models have a maximum degree and order (d/o) of 120.

3.2.2 Simulation of satellite observations

First, we need to set the time period for the simulation. Assessments by research institutions of the Bender constellation's capability to recover the time-varying gravity field were typically conducted over the period of a repeat orbit to provide the densest possible ground track pattern (Pail et al. 2019). In this experiment, the orbital repetition period of Bender constellation is 9 d. In order to make a unified comparison with the results of Starlink-like configuration, the simulation period was set to a duration of 9 d commencing from 0:00 on 2002 January 1, with a sampling interval of 5 s. AOHIS data (Dobslaw et al. 2015) during the corresponding period were selected and linearly interpolated for each epoch. The initial positions of the satellites were determined according to the orbital parameters of the Starlink-like and GRACE/Bender constellations in Section 3.1. Utilizing the forces acting on the satellites in the true models, the orbit was recursively integrated using the Runge–Kutta 4 and Adams–Bashforth–Moulton methods. Ultimately, orbit positions and velocities of all satellites in each epoch of the Celestial Reference Frame (CRF) were obtained. For the Bender satellites, the intersatellite range rate was calculated based on the simulated error-free orbit positions of the two satellites in the same orbit, serving as the observation for subsequent inversion.

In this paper, we assumed that all satellites were equipped with accelerometers to measure non-conservative forces. During the simulation experiment, the accelerometer data were not simulated and involved in the calculation, but the influence of its measurement error was directly considered. Consequently, the satellite attitude data were needed to convert the accelerometer measurement error into data in CRF coordinate system. Therefore, based on the simulated satellite velocity, we obtained its motion in CRF, calculated the attitude data corresponding to each epoch, and expressed it in quaternion form.

3.2.3 Noise simulation

Since the accuracy level of orbit determination of low-orbit satellites can reach 1–2 cm at present (Bock et al. 2014), Gaussian white noise with a standard deviation of 1 cm was applied to the three-directional components of the orbit positions of Starlink and Bender satellites in the CRF coordinate axes. The intersatellite range rate of the Bender satellite was obtained using the Laser Ranging Interferometer, with coloured noise added according to eq. (5). Coloured noise for the accelerometer's sensitive and non-sensitive axes was added according to eqs (6) and (7), respectively. Fig. 3 shows the amplitude spectral density (ASD) of coloured noise added in the experiment:

3.2.4 Calculating the time-varying gravity field from simulated Starlink-like and GRACE/Bender constellations

1. Decorrelation of coloured noise

   Due to the presence of coloured noise in the accelerometer and intersatellite range rate observations, correlations exist between the observations of different epochs. To properly determine the weight of the observations, a covariance matrix was established based on the residuals of the observations for de-correlation. The concept is as follows: initially, equal weights of the observations in different arcs and epochs are used to construct the normal equations. The parameters to be obtained are then calculated after combining all the arcs. Subsequently, the calculated results are substituted into the design matrix, and the residual is obtained through comparison with the observations. Following this, the PSD of the residual is estimated and transformed into a covariance matrix using a cosine transform. The new covariance matrix is then utilized to decorrelate the observations and construct a new normal equation. The covariance matrix, refined through multiple iterations, was ultimately employed for the decorrelation of observations.

2. Maximum degree of normal equation

   In scenarios where the recovery period was curtailed to short durations, the sample size of GRACE/Bender constellation was limited, and insufficient to support the restoration of the time-varying gravity field to a higher degree. Consequently, it is imperative to set the maximum degree of the normal equation reasonably. According to the Colombo-Nyquist rule (Colombo 1984), the maximum resolvable degree of the gravity field equals half the number of satellite revolutions at even ground track spacing. In the context of experiments conducted over short periods (spanning 1, 3 and 6 d), although there was no complete repetition period and the orbital spacing varied, the Colombo-Nyquist rule can still be considered according to the number of samples on the equator, and the maximum degree of the normal equation was set to half the number of satellite revolutions, as shown in Table 6, while the influence of uneven sampling was outside the scope of this paper. Because the Colombo-Nyquist rule is too conservative and theoretically possible to obtain better spatial resolution (Visser et al. 2012; Weigelt et al. 2013), the application of this rule will not have a significant impact on the comparison with the Starlink-like constellation. For comparison with the results of other institutions (Pail et al. 2019), the maximum degree of the normal equation for each constellation was not exceed 70. When using the daily parameterization method, due to the difference in the number of observations within a single day and the constraints on the degree of the normal equation, the maximum degrees of the daily estimation of the normal equations for the Starlink-like, Bender, and GRACE configurations were set to 15, 15 and 7, respectively.

4 RESULTS AND DISCUSSION

The root mean square (RMS) of spherical harmonic coefficients degree errors is typically calculated using eq. (8) to demonstrate the variance between the calculated results and the true signal. The superscript ‘est’ represents the recovery result of the gravitational time-varying signal, while ‘ref’ denotes the average value of the AOHIS signal in the true model during the recovery period.

When comparing the degree error RMS of the Starlink-like constellation with those of other constellations, it is important to note that the lower order terms may be affected by the polar gap (PG) problem, leading to a larger degree error RMS. Sneeuw & van Gelderen (1997) gave a rule of thumb for estimating the maximum order affected by the PG problem:

where I is the satellite's orbital inclination and n is the degree. When comparing the degree error RMS of the solution results from different constellation configurations, the low-order coefficients affected by the PG problem, caused by the orbital inclination of 97.6°, should be removed first. Subsequently, the degree error RMS can be calculated using eq. (10), allowing the effect of the PG problem on the degree error RMS to be disregarded:

4.1 Arc length selection

Before using the short-arc integral method to solve the time-varying gravity field, it is necessary to consider the arc length of the integral. In the traditional short-arc integral method, the arc length is generally less than 1 hr, requiring a gradient correction of the satellite orbit. Daras & Pail (2017) set the short-arc length of integration to 30 min in their simulation, while Chen et al. (2015) proposed the modified short-arc approach, which set the short-arc length of integration to 2 hr after comparing the calculation effect of different arc lengths. The impact of arc length on the results is multifactorial. A longer arc length can lead to a greater accumulation error of the force and integral models (You et al. 2011). At the same time, however, a longer arc may amplify the small force acting on the spacecraft, thereby aiding in better estimation of the tesseral harmonic coefficients (Xu 2008; Chen et al. 2019). On the other hand, in terms of decorrelating observations, the covariance matrix corresponding to the long arc has more elements, allowing for a more accurate reflection of the correlation between observations, whereas the covariance matrix corresponding to the short arc has a more limited reflection of the correlation between observations.

To determine the appropriate arc length for the short-arc integral method, the 1-d solution of 2002 January 1, was resolved with short-arc lengths of 0.5, 1 and 2 hr, respectively. The results were then compared with the true AOHIS signal. The degree errors of the 1-d solutions for the GRACE/Bender and Starlink constellations under different arc lengths were calculated according to eq. (8), as shown in Fig. 4. As can be seen from these figures, for the GRACE/Bender constellation, the 2 hr arc length yields the best results, with better accuracy at lower degrees. Consequently, an integral arc length of 2 hr was selected for the recovery of the time-varying gravity field for the GRACE/Bender constellation. For the Starlink-like constellation, the accuracy of the results obtained using different arc lengths is slightly different. However, due to the inability of the Starlink satellite orbits to fully cover the polar regions, the accuracy of the lower order terms of the potential coefficients is compromised (Sneeuw & van Gelderen 1997). The results show that the 0.5 hr arc length is less affected by the polar gaps (Fig. 5), leading to the determination of the arc length for the time-varying gravity field of the Starlink-like constellation.

4.2 Error analysis

In the simulation of the Starlink-like and GRACE/Bender constellations, errors in the recovered time-varying gravity field model may primarily stem from the following sources: the aliasing effect caused by insufficient time sampling of the AOHIS signal, the aliasing effect caused by ocean tide model errors, and observation noise. In addition, due to the recovery of time-varying signal to degree 70, time-varying signal higher than degree 70 in real models will be mapped to low-degree terms, which will also cause aliasing effects. To analyse the impact of aliasing errors on the result of recovering the time-varying gravity field in the Starlink-like and GRACE/Bender constellations, the influence of each error source was assessed individually. To this end, five schemes were devised for recovering 9-d time-varying gravity field, with specific information shown in Table 7.

Scheme 1 (no errors) is the most ideal error-free simulation scenario without considering any aliasing effects and observation noise. In scheme 2 (high-degree signal aliasing), we only consider the impact of the part of the time-varying signal higher than degree 70 on the solution result. In scheme 3 (AOHIS aliasing), we only consider the aliasing effect caused by insufficient time sampling of AOHIS signal. In scheme 4 (tidal aliasing), we only consider the aliasing effect caused by ocean tide model errors. In scheme 5 (sensor noise), we only consider the observation noise. Without considering the aliasing effect caused by the mapping of higher to lower degree terms, the maximum degree of time-varying signal in the true model is 70. When the aliasing effect caused by ocean tide model errors is not considered, EOT11a is used in both the true and the reference models. When the aliasing effect caused by undersampling of the non-tidal time-varying signal in time is not considered, the time-varying field in the true model is set to the mean value of the AOHIS over 9 d.

Fig. 6 depicts the degree error RMS of the 9-d solution compared with the mean value of AOHIS signal in the true model in schemes 1–5, without daily parametrization. In Fig. 6(a), it is evident that when all the aforementioned error sources are disregarded, all constellations yield favourable recovery outcomes, with GRACE/Bender yielding superior results. Fig. 6(b) illustrates that the degree error increases to some extent compared to the error-free result, due to the mapping of background field signals beyond 70° to those within 70°. In Figs 6(c) and (d), when considering the aliasing effect caused by undersampling of AOHIS signals or ocean tide model errors, compared to the GRACE/Bender constellation, the Starlink-like constellation showed a significant improvement in results, indicating that observations from giant constellations can effectively mitigate the impact of aliasing effects. Fig. 6(e) reveals that when considering the influence of observation noise, the recovery model error of the Starlink-like constellation increases rapidly with degree, and it is only marginally better than GRACE up to the degree 10 and significantly worse than the GRACE/Bender constellation at higher degrees.

At the same time, we note that in Figs 6(b)–(d), the Starlink abcd configuration (more satellites) exhibits inferior performance to Starlink abc (fewer satellites). The Starlink abcd configuration, which has a larger number of satellites with an orbital inclination of 53°, will be subject to more severe PG problem than the abc configuration. As can be seen from Fig. 7, the maximum order affected by PG problem at the degree 70 is about 45, which accords well with the rule of thumb given by Sneeuw and van Gelderen (1997). Compared with the Starlink abc configuration, abcd configuration has a greater error in the wedge region affected by PG problem in the spectrum and performs slightly worse in Figs 6(b)–(d).

As can be seen from Fig. 8, for Starlink-like constellation, the influence of observation noise is greater than that of other error sources, and the aliasing error caused by the mapping of high-degree signals is only more significant in the part above degree 40 compared to other error sources. For GRACE/Bender constellation, AOHIS aliasing error is the main influencing factor, far greater than the influence of observation noise. The above analysis indicates that under the current numerical simulation conditions, the recovery results of the Starlink-like constellation are primarily constrained by the accuracy of the observation. In contrast, those of the GRACE/Bender constellation are mainly restricted by the aliasing effect caused by the undersampling of the AOHIS signal and errors of the ocean tide models.

4.3 Analysis of short-period gravity fields

We also investigated the recovery of time-varying gravity fields for various periods (1, 3, 6 and 9 d) used the Starlink-like and GRACE/Bender constellations, while considering all error sources. Figs 9 –12 display the degree error RMS between the calculated model and the true model's AOHIS signal, both before (left) and after (right) utilizing Wiese's daily parametrization. These figures reveal that the Starlink abcd configuration can effectively recover signal up to degree 20 with just one day of observations, demonstrating the capacity to retrieve low-degree time-varying gravity field signal in a short time frame. This performance is notably superior to the results obtained from the GRACE/Bender constellation. For multiday solutions, the Wiese's method can partly improve the solution accuracy of the GRACE/Bender constellation. In contrast, the Starlink-like constellation effectively mitigates aliasing errors by increasing the sampling, thus resulting in no significant improvement after applying Wiese's method. Without implementing Wiese's method, the solutions using the Starlink abcd configuration are superior to those of the Bender constellation up to degree 20. With Wiese's method, the solutions using the Starlink abcd configuration surpass those of the Bender constellation up to degree 15.

Compared with Fig. 6, it can be seen that for the 9-d time-varying gravity field recovery results, the Starlink-like constellation mitigates the aliasing effect caused by insufficient time sampling of the AOHIS signal and the errors in ocean tide models, particularly for degrees less than 25, thereby improving the accuracy of low-degree terms. Above degree 25, the error is greater than the time-varying gravity field signal because of the limitation of observation accuracy.

For a more comprehensive analysis of the results, the coefficients errors of the 9-d solutions of the Starlink-like and GRACE/Bender constellations are presented in Fig. 13, which indicates that the solution results of the Starlink-like constellation outperform in terms of sectorial harmonics. However, the zonal harmonic coefficients exhibit inferior performance due to polar gap regions in the orbit coverage. Additionally, the polar blank effect becomes increasingly pronounced with the inclusion of low-inclination satellites (parts c and d). Furthermore, the addition of parts b, c and d in part a of Starlink significantly enhances the accuracy of the sectorial harmonic terms and partial tesseral harmonic terms. This improvement is attributed to adding low-inclination satellites to near-polar orbit satellites, which increases the amount of observation in other directions. This indicates that multidirectional gravity field observation is conducive to improving the accuracy of the model solution. In comparison to the GRACE configuration, the Bender configuration demonstrates improved accuracy in the sectorial harmonic coefficients. Although the Starlink abcd configuration outperforms the Bender configuration at lower degrees, the latter yields superior results at higher degrees, owing to its high-precision ll-SST observation mode. This is consistent with the degree variance analysis result.

The equivalent water height (EWH) degree error RMS between the 9-d solution result of the Bender constellation and the AOHIS signal according to eq. (11) is shown in Fig. 14. By comparing these results with those of other institutions using the same background model (Pail et al. 2019), it can be found that the simulation results presented in this paper are generally consistent with them. Notably, the accuracy closely resembles that reported by Tongji University, confirming the reliability of the results in this study:

Figs 15 and 16 show the spatial distribution of EWH errors between the recovery results of different constellation configurations and the true model's AOHIS signal within the first 15° under varying recovery periods. The inversion results of the Starlink-like constellation exhibit an east–west striping pattern in the error distribution, while the GRACE/Bender constellation shows a north–south striping pattern. Additionally, the addition of low-inclination satellites effectively improves the calculation results in middle and low latitudes. Given its denser coverage in the middle and low latitudes, along with sparse coverage in the high latitudes and no coverage in the polar regions, the Starlink-like constellation exhibits superiority over the GRACE/Bender constellation in recovering time-varying signal (up to degree 15) at middle and low latitudes. Conversely, at high latitudes, the results from the Starlink-like constellation are frequently poor due to the absence of polar observations, and the polar error becomes more serious with an extended recovery period.

Based on the preceding analysis, it is evident that in instances where ocean tide models exhibit errors and the GRACE/Bender satellite constellation inadequately samples the time-varying signal of AOHIS, reducing low-degree aliasing errors can be achieved by increasing the sampling and multidirectional observation capabilities of the Starlink-like constellation. Below degree 15, the solution results of the Starlink abcd configuration outperform those of the Bender constellation. In spatial distribution, the inversion results of the Starlink-like constellation are superior in the middle and low latitudes, while the inversion results of the Bender constellation are superior in the high latitudes.

4.4 The combination of Starlink-like and GRACE/Bender constellations

Given that the Starlink-like and GRACE/Bender constellations each have distinct advantages across different degrees, orders and spatial regions, the normal equations of these constellations are considered to be jointly solved. The performance of the combined Starlink + GRACE/Bender configuration in recovering the time-varying signal (AOHIS) was assessed by comparing the outcomes before and after the integration.

After implementing the daily parametrization method, Figs 17 and 18 display the degree error RMS and coefficients errors of the Starlink abcd configuration, GRACE/Bender constellation, and the combined solution results. As can be seen from these figures, the combined Starlink and Bender solution provides the best performance. The combined results exhibit a marked enhancement across all degrees, and the recovered time-varying gravity field within degree 30 can be improved by 0.5–1 order of magnitude compared with the individual solutions of GRACE/Bender. The error spectrum shows that the combined results have improved the precision in the zonal and sectorial harmonic terms. The addition of Starlink satellites enhances observation in the cross-orbit direction, addressing the sampling limitations of the GRACE/Bender satellites and greatly improving spatial coverage. Conversely, the inclusion of GRACE/Bender satellites fills in observation gaps in the polar region, ameliorating the accuracy of the zonal harmonic terms where Starlink satellites are lacking.

Fig. 19 shows the spatial distribution of the EWH error calculated using spherical harmonic coefficients in the first 30 degrees. As can be seen from the figure, the combined solution results mitigate the striping effect in the high-latitude, and middle- and low-latitude regions compared to the Starlink-like and GRACE/Bender configurations, respectively. Contrastingly, when compared to the results of the Bender configuration, the increased sampling significantly reduces the aliasing effect caused by insufficient sampling of the AOHIS signal and errors in ocean tide models in the low-degree terms, ultimately enhancing calculation accuracy.

5 CONCLUSIONS

In this study, we conducted simulations using 9 d of observations, employing the first phase configuration of Starlink and the Kepler parameters of the Bender constellation from the ADDCON project. In the context of AOHIS signal aliasing errors, ocean tide model errors, orbit position errors, intersatellite range rate errors and accelerometer observation errors, a closed-loop simulation was performed to analyse the recovery ability of the complete spectrum of the Earth's nontidal geophysical processes (AOHIS), using the short-arc integral method. Our analysis focused on comparing the performance of the Starlink-like, GRACE and Bender constellations, with particular emphasis on evaluating the influence of different error sources, the recovery period of time-varying models, and different constellation configurations. Additionally, we combined different constellations to assess their collective capability in recovering the time-varying gravity field. The key conclusions derived from our study are as follows:

1. Without considering observation noise, the Starlink-like constellation significantly increases the number of observed samples within the same period compared to the GRACE/Bender constellation. As a result, it effectively mitigates the aliasing effect caused by insufficient time sampling of the AOHIS signal and the errors of ocean tide models.

2. When considering all error sources, the Starlink-like constellation outperforms the GRACE/Bender constellation in recovering low-degree signal of the short-period time-varying gravity field. The recovered 1, 3, 6 and 9-d time-varying gravity field models demonstrate superior performance within degree 15. Notably, the signal-to-noise ratio (SNR) of the 1-d solution reached approximately 1 at degree 20, and the SNR of the 9-d solution reached approximately 1 at degree 25.

3. In the case of observation error simulation, the aliasing error is the main limiting factor for the accuracy of solving the time-varying gravity field by the GRACE/Bender constellation. To a certain extent, the daily parametrization method can improve the accuracy of recovering time-varying signal using the GRACE/Bender constellation. However, for the Starlink-like constellation, the improved spatiotemporal sampling effectively reduces the aliasing error, making the daily parametrization method is not obvious.

4. Jointly solving the time-varying gravity field model with the Starlink-like constellation and the GRACE/Bender constellation provides superior results compared with only using individual constellation. This is primarily attributed to the enhanced accuracy of the zonal harmonic terms and sectorial harmonic terms. The effective degree of the 9-d gravity field model obtained by the joint solution can reach 36 (SNR reaches 1). The combination of the Bender constellation with the Starlink-like constellation reduces the aliasing effect and improves the accuracy of solving low-degree terms.

It is important to note that this study simulated and analysed the ability of the Starlink-like constellation to invert the time-varying gravity field under ideal conditions. The satellite orbit accuracy reached 1 cm, whereas the real Starlink satellite orbit accuracy used for communication is far lower. At the same time, the satellite loads, observation accuracy, orbit parameters, flight control and other aspects under real conditions were not discussed in detail. In addition, because of the limitations in computing and storage space, only a subset of the Starlink constellation configuration was investigated in this study. If the entire (12 196) satellites of the Starlink phase I configuration were used, the recovery results for the time-varying gravity field model would likely be improved at the current observation accuracy level. Our results can provide a technical reference for the requirement analysis of star network missions and the design of future satellite gravity missions.

DATA AVAILABILITY

The data sets used in this study are as follows: the complete data set of the updated ESM (Dobslaw et al. 2015) is publicly available at doi:10.5880/GFZ.1.3.2014.001, the code for the GROOPS software is available as an open source on github (https://github.com/groops-devs/groops).

Acknowledgement

This study was funded by the National Natural Science Foundation of China (nos 42192533, 42388102 and 42074019) and the Fundamental Research Funds for the Central Universities (no. 2042022dx0001).

References