Influence of the Low-Frequency Error of the Residual Orbit on Recovering Time-Variable Gravity Field from the Satellite-To-Satellite Tracking Mission

When using the dynamic approach to recover the time-variable gravity field, the reference orbit generated by the perturbation model and the non-conservative force observed from the accelerometer should be introduced at first, and then the observation equations of the residual orbit and the residual range rate are established. This introduces a perturbation model error and instrument noise. Thus, there are low-frequency errors in the residual orbit and the residual range rate. Currently, most studies only focus on the low-frequency error of the residual range rate, neglecting the influence of the low-frequency error in the residual orbit. Therefore, under the condition of the perturbation model error and instrument noise including the constant term and 1CPR term, the low-frequency error formulas of the residual orbit and residual range rate are derived according to the characteristics of the solution of the Hill equation. Then, the influence of the low-frequency error on the residuals is analyzed by using the simulation and the real data processing respectively. In the simulation and real data processing, the accuracy of the recovered gravity field can maintain a good consistency for different arc lengths by removing the low-frequency error in the residual orbit. Finally, the time-variable gravity field model UCAS-IGG (University of Chinese Academy of Sciences-Institute of Geodesy and Geophysics) was solved from January 2005 to February 2010 by removing the low-frequency error of the residual orbit and residual range rate. Compared with the official institutions, the UCAS-IGG presents a good consistency in the estimating time-variable gravity field signal. This study demonstrates how the effect of the low-frequency error of the residual orbit should be taken into consideration when the longer arc length is used to recover a time-variable gravity field. Using a long arc length can reduce the variables of the initial state and recover the influence of the small force.

The essence of the dynamic approach is to generate the reference orbit and reference range rate by using the conservative force model (N-body perturbation, Solid Earth tides, Ocean tide, and so on) and the non-conservative force data observed by the accelerometer, to then obtain the residual orbit and residual range rate by the differences between the GRACE data and the reference data, and build the observation equations of the residual orbit and residual range rate to recover the time-variable gravity field. In fact, the used perturbation model and the non-conservative force are not consistent with the real forces on the satellite. This leads to the perturbation model error and instrument noise. Thus, the residual range rate includes the low-frequency error [29]. The low-frequency error can propagate globally during orbit integration and accelerometer calibration procedure, and could cause errors at low degrees [30]. For the low-frequency error, Tapley et al. [3], Visser [31], Bruinsma et al. [32], Liu et al. [18,19], Meyer et al. [33], Ditmar et al. [34], Hashemi Farahani et al. [35], Wang et al. [36], Luo et al. [37], and Zhou et al [8]. used an empirical formula to absorb the low-frequency error of the residual range rate. For utilizing this method, Zhao et al. [30] pointed out that the co-estimated empirical parameters and the coefficients of the gravity field have the minimum influence on the time-variable gravity field; Zhou et al. [8] pointed out that the low-frequency error is simultaneously removed on both sides of the observation equations by using empirical formulas, which also has a minimal effect on the time-variable gravity field. In fact, the residual orbit also includes the low-frequency error and affects the accuracy of the recovered gravity field. However, most studies mainly focus on the low-frequency error of the residual range rate and ignore the low-frequency error of the residual orbit.
The forms of the low-frequency error of the residual orbit and residual range rate are related to the solutions of the Hill Equation [29]. Colombo [38] solved the Hill equation based on the Laplace Transform. There are characteristic roots 0 and ±ωi (i = √ −1 is a complex unit, ω is the satellite's mean angular velocity) in the homogeneous Hill equation. If the inhomogeneous term (perturbation model error and instrument noise) of the Hill Equation has a constant term and 1CPR (Cycle Per Revolution) term, there are t, t 2 , t cos ωt and t sin ωt (where t is time), which are resonance terms and belong to the low-frequency form, in the special solution of the Hill Equation according to the differential Equation solution rule. Besides, Seo et al. [39] pointed out that the order 15 and multiples of 15 of the design matrix of the coefficients of the disturbing potential also have the same resonance terms according to Kaula's [40] orbital resonance theory. For a longer arc length, the resonance terms will increase rapidly and the accuracy of the recovered time-variable gravity field becomes worse, especially the accuracy of the multiples of the order 15, when least-squares is used to recover the time-variable gravity field. Colombo [38] gives the general solution of the homogeneous Hill Equation based on the Laplace Transform, but this method is in an algebraic form to solve the differential equation, and it is difficult to find the physical meaning, thus it is not convenient for practical applications. For this purpose, we re-derive the solution of the Hill Equation in the satellite moving frame, which gives the low-frequency error formulas of the residual orbit and the residual range rate.
The outline of this study is as follows. In Section 2, we establish the observation equations of the residual orbit and the residual range rate with the non-linear corrections based on the orbital perturbation differential equations. In addition, the solutions of the Hill Equation are re-derived based on the coordinate transformation. In Section 3, we do the simulation and real data processing to illustrate that the low-frequency error of the residual orbit should be considered. In Section 4, the time-variable gravity field model entitled UCAS_IGG is solved from January 2005 to February 2010 by removing the low-frequency error of the residual orbit and residual range rate. Finally, the conclusions are summarized in Section 5.

Non-Linear Correction
The dynamic approach is derived based on Taylor's expansion, where the linear form is kept and the non-linear term is ignored [41][42][43]. Due to using the longer arc length (i.e., 48 h) to recover the time-variable gravity field model, the impact of the non-linear term should be considered. For using a longer arc length to recover the time-variable gravity field, Yu et al [44]. and Xu [45] gave the analysis. Xu [45] pointed out that the error caused by linearization increases rapidly as the arc length increases for the dynamic approach, and proposed that the satellite's observation position (Satellite to Satellite Trace) is used as a reference orbit to make sure the linearization errors do not accumulate largely and prolong the arc length. However, at this time the Volterra's integral equation is introduced to establish the observation equations. Yu et al. [44] introduces the non-linear corrections to compensate for the rapid increase of the residual, which greatly relaxes the application limits of the dynamic approach and prolongs the arc length of the recovered time-variable gravity field. This section chooses the dynamic approach with the non-linear correction to build the observation equations of the residual orbit and residual range rate since the Volterra's integral equation is more difficult to calculate.
For the convenience of description, without special declarations, mathematical notations are defined as follows: r is the position vector of the satellite (r is module of r), v or . r is the velocity vector of the satellite, the corresponding subscript "0" indicates the initial state of the satellite, the subscript "1" denotes the real position of the satellite, the subscript "2" is the reference position of the satellite, the superscript "(A)" and "(B)" indicate the GRACE-A and GRACE-B respectively. For example: represents the velocity vector of the reference orbit of GRACE-A.
The orbital perturbation differential equations of the single satellite with the non-linear corrections can be written [44]: where β = r 1 − r 2 is the vector of the residual orbit, H 3×3 (t) is the Hessian matrix of the reference gravity field on the reference orbit, b 1 (t) is the non-linear correction term, T k is the coefficient of the disturbing potential, and b k (t)(k ≥ 2) is the value of the spherical harmonic function of the different degree and order along the orbit. In Equation (1), the T k needs to be estimated.
The system of differential Equation (1) can be solved by using the superposition principle of the solution because Equation (1) is linear with respect to β. Then, it can be decomposed into the following three equations: where O 3 and I 3 are the three-order zero matrix and the identity matrix, respectively. The solution is the non-linear correction term when k = 1 in the system of Equation (4). The system of differential Equations (2) and (3) are usually called state transition equations, and the system of the differential Equations (4) are usually called variation equations in the satellite precise orbit theory. Of course, b k (t) not only indicates the spherical harmonic function and non-linear correction term in the inertial coordinate system, but also indicates the perturbation model error and the instrument noise in the inertial coordinate system. Therefore, the solution of the system of differential equations (4) is the low-frequency error of residual orbit when b k (t) is the perturbation model error and instrument noise.
Because the differential equations are linear about their initial problem, their solutions exist and are unique. If their solutions are denoted by Φ 3×3 (t), Φ 3×3 (t),Ψ 3×3 (t) and B k (t)(k = 1, 2...), the solution of the orbital perturbation differential equations (1) can be written as Since Equations (2)-(4) are solved based on numerical integration, we can simultaneously obtain: .
The above is the observation equation of the residual orbit with a non-linear correction term for a single satellite. Next, we will establish the observation equation of the residual range rate with a non-linear correction term. In fact, the range rate of GRACE can be written .
1 is the position difference between GRACE-A and GRACE-B, α 1 is the corresponding mode, and is the velocity difference between GRACE-A and GRACE-B. About the range rate . ρ 1 (t), the corresponding observation equation can be established as follows: . where In the observation Equation (7), the variables with the subscript "2" are obtained from the known reference orbit, . ρ 1 (t) is the observation value of the range rate, T k is the coefficient of the disturbing potential that needs to be estimated, where the variables contained in ε(t) and . ε(t) are from the solution of the system Equation (1).

The Low-Frequency Error
The system Equation (4) is the Hill Equation when the b k (t) represents the perturbation model error and instrument noise. To solve the Hill Equation, we first need to solve the general solutions of the differential equation For GRACE, its orbit is a near-circular polar orbit. The zero-order term GM r (where GM is the gravitational constant of the Earth, r = x 2 + y 2 + z 2 , (x, y, z) is the components of the coordinate system of the satellite orbit plane) of the gravity field is the main term for generating the reference orbit [46,47]. Therefore, under the condition of neglecting the influence of the term J 2 (J 2 is the oblateness of the Earth and the order of magnitude is about 10 −3 ), the Hessian matrix can be written as [44] and the satellite's position in the coordinate system of the satellite orbit plane can be written as where a is the mean orbit radius of satellite, ω = GM a 3 is the average angular velocity of satellite. From Equation (11), the O − xz plane is the orbit plane of the satellite and the y-axis is the normal direction of the orbit plane. Substituting (10) and (11) into the state transition Equation (2) or (3), the state transition equation can be written in the following simplified form: ..
Obviously, the general solution to Equation (14) is This solution reflects the characteristics of the residual orbit in the normal direction of the orbit plane. Using the matrix diagonalization method for Equation (13), we can obtain 1 − 3 cos 2 ωt cos ωt sin ωt cos ωt sin ωt 1 − 3 sin 2 ωt = cos ωt − sin ωt sin ωt cos ωt Let Thus, the system of Equation (13) is reduced to ..
By a further simplification, we can obtain In the system of Equation (19), there are characteristic roots 0 and ±ωi for η 1 and there are characteristic roots ±ωi and double roots 0 for η 3 . Therefore, the general solution can be found as In the transformation Formula (17), η 1 is the radial component of the residual orbit (outer normal direction r 0 ), and η 3 is the moving direction component of the residual orbit (tangential direction τ ). Thus, according to Equation (15) and Equation (17), we know that (η 3 , β 2 , η 1 ) is the component of the residual orbit in the satellite tangential coordinate system (because GRACE is a near-circular polar orbit, the tangential coordinate system is basically the same as the radial coordinate system). Let n represent the normal direction of the orbital plane, then {τ, n,r 0 } constitutes the satellite's moving coordinate system, and τ ≈ −r 0 . So far, we have solved the homogeneous Hill Equation, and (η 3 , β 2 , η 1 ) is the component of the residual orbit in the satellite moving coordinate system {τ, n, r 0 }. Now, we will discuss the structure of the Hill Equation solution. Assuming that the components of the perturbation model error and instrument noise are {F 1 , F 2 , F 3 } in the satellite moving coordinate system {τ, n, r 0 }, the Hill Equation in the moving coordinate system can be written as [29,38,48] follows: If the perturbation model error and the instrument noise include a constant term and 1CPR term in the satellite moving coordinate system {τ, n, r 0 }, the special solution of the Hill equations can be written as follows: where c k , D k , E k , F k , G k (k = 1, 2, 3), b 1 , b 3 , a 3 are constants and need to be estimated, η * 1 , β * 2 , η * 3 are residual orbit components caused by the perturbation model error and instrument noise. Formula (22) is defined as the low-frequency error of the residual orbit. Actually, the design matrix of the coefficients of the disturbing potential contains the general solution of the Hill Equation, and the design matrix also includes the same resonance terms as Formula (22) in the multiples of the order 15 [39,49]. If the low-frequency error of the residual orbit is not processed, the low-frequency error will accumulate faster when the longer arc is used to recover the time-variable gravity field. This leads to a decrease in the accuracy of the recovered gravity field in using least-squares to solve the observation equations. Therefore, we suggest that the low-frequency error of the residual orbit need to be processed using the orbit data to recover the time-variable gravity field. At the same time, the low-frequency error of the residual range rate can be derived. For the range rate, we can get For the GRACE, satellite A and B are close. Therefore, the unit vector and tangent τ are approximately equal, the Equation (23) can be simplified as Substituting Formula (22) into Formula (24), we can derive the low-frequency error of the residual range rate: . ρ bt = a 0 + a 1 t + a 2 t 2 + a 3 cos ωt + a 4 sin ωt + a 5 t cos ωt + a 6 t sin ωt (25) where . ρ bt is the low-frequency error of the residual range rate caused by perturbation model error and instrument noise, a k (k = 1, 2...6) are constant terms that can be estimated by using least-squares.
Compared with the solutions of the Hill Equation with Colombo [38], our solutions are more convenient for the practical application. For the low-frequency error processing of the residual orbit, the observation equations of the residual orbit are first converted into the satellite moving coordinate system according to Formula (17), and then the low-frequency error of the residual orbit is removed according to Formula (22). The low-frequency error of the residual range rate is directly processed according to Formula (25). For processing the low-frequency error method, Zhao et al. [30] and Zhou et al. [8] both gave the empirical parameter processing method and the results were basically the same. In this study, we chose to remove the low-frequency error on both sides of the observation equations for every 1.5 h, and then recover the time-variable gravity field.

The Low-Frequency Error Processing
We will illustrate whether the low-frequency error of the residual orbit needs to be processed on the recovering gravity field by using simulation and real data processing.

Simulation
In the simulation, the real position and velocity of GRACE at 00:00 h on January 1 2010 are chosen as the initial position and velocity of the simulation. The first degree and order 60 of EGM08 and GIF48 are selected as the real gravity field and reference gravity fields, respectively. The difference between the EOT11A and FES2004 tide model is used as the perturbation model error. In addition, we also assume that there is a deviation of 3 cm and 0.05 cm/s in the initial position and velocity, respectively. Furthermore, there are 2 cm and 1.0 × 10 −7 m/s random errors in the orbit and range rate, respectively. The arc lengths of the recovered gravity field are 6, 12, 24, and 48 h respectively. The running time of the satellite is 30 days. The green line expresses the low-frequency error is removed and the magenta line expresses the low-frequency error is not removed. The degree variance and the logarithm of the formal errors are calculated according to Appendices A and B, respectively. In the inverted triangles, the first row indicates that the low-frequency error of the residual orbit is not removed, while the second row is removed.
The orbit data are used to recover the gravity field. The degree variance of the recovered gravity field for different arc lengths is shown in Figure 1. The accuracy of the gravity field model can be improved by removing the low-frequency error of the residual orbit, that is, the green line is slightly lower than the magenta line. For different arc lengths, the degree and order of the recovered gravity field model orbit is about 45, without removing the low-frequency error of the residual, while the degree and order is about 50 by processing the low-frequency error of the residual orbit. The logarithm of the formal errors of the recovered spherical harmonic coefficient in the form of inverted triangle plots are shown in Figure 2 for different arc lengths. It can be seen that the accuracy can be improved in the second row compared with the first row.  The logarithm of formal errors (comparing the EGM08) in the spherical harmonic coefficients of the recovered gravity field by using 6, 12, 24, and 48 h arcs, when using the orbit data to recover the gravity field.
We present the results for the different arc lengths by using the range rate data to recover the gravity field. In Figure 3, the accuracy of the recovered gravity field by removing the low-frequency error is better than the one that does not remove the low-frequency error. The gravity field models are presented by inverted triangles as shown in Figure 4. Compared with the first row, the accuracy of the low-order part of each degree of the recovered gravity field, especially within the order 15, is mainly improved in the second row, and there are obvious blue band areas. Figure 3. The degree variance of the recovery gravity field by using 6, 12, 24 and 48 h arcs, when using the range rate data to recover the gravity field. Figure 4. The logarithm of formal errors (comparing the EGM08) in the spherical harmonic coefficients of the recovered gravity field by using 6, 12, 24, and 48 h arcs, when using the range rate data to recover gravity field.
Combine the orbit data and the range rate data to recover the gravity field, where the low-frequency error of the residual range rate is removed. The degree of variance of the recovered gravity field for different arc lengths is shown in Figure 5. Whether or not the low-frequency error of the residual orbit is removed, the accuracy of the first degree and order 50 of the recovered gravity field is basically the same, at 6, 12, and 24 h arcs. But after the degree and order 50, the accuracy of the recovered gravity field by removing the low-frequency error of the residual orbit is better. For 48 h arcs, we can see that the large jump appears around the 30th degree when the low-frequency of the residual orbit is not removed. However, the large jump does not appear when removing the low-frequency error of the residual orbit. Figure 5. The degree variance of the recovered gravity field by using 6, 12, 24 and 48 h arcs, when combining the orbit data and the range rate data to recover the gravity field.
The degree variance of the recovered gravity field for different arc lengths in Figure 5 is plotted in Figure 6 according to whether or not we are processing the low-frequency error of the residual orbit. It can be seen that the degree of variance of the gravity field obtained by removing the low-frequency error of the residual orbit has a good consistency. However, the recovered gravity field without processing the low-frequency error of residual orbit has a good consistency at 6, 12 and 24 h arcs, while the degree variance has a large jump at 48 h arcs, which cannot maintain a good consistency.
The logarithm of the formal errors of the recovered spherical harmonic coefficient for different arc lengths in the form of inverted triangle plots are shown in Figure 7. It can be seen that the accuracy of the multiples of the order 15 of the recovered gravity field and the low order (within the order 15) can be improved. There are obvious blue band areas in the second row of inverted triangles, while the first row does not show such blue bands. In addition, the accuracy of higher-order parts is also improved, just as the red areas on both sides of the inverted triangle are reduced.  The logarithm of formal errors (comparing the EGM08) in the spherical harmonic coefficients of the recovered gravity field by using 6, 12, 24, and 48 h arcs, when combining the orbit and range rate to recover the gravity field.
In the case of combining the orbit and the range rate to recover gravity field, when removing the low-frequency error of the residual orbit, the advantage of the low-frequency error processing of the residual range rate in Figure 4 is retained in Figure 7, that is, the same blue band region is displayed in the second row of Figure 7. However, this advantage cannot be retained when only the low-frequency error of the residual rate is removed.
From the simulation, we proved that this method can be used to recover the gravity field, the accuracy of the recovered gravity field can be improved, especially the accuracy of the multiples of the order 15, and the accuracy of the recovered gravity field has a good consistency for different arc lengths.

Background Model
According to Formulas (5) and (7), the observation equations of the residual orbit and the residual range rate are established to recover the time-variable gravity field. The reference orbit is required to remove the time-variable signal caused by the N-body perturbations (International Earth Rotation Service (IERS) 2010 conventions), the solid tides (IERS 2010 conventions), the ocean tides (EOT11a, Nmax = 100), the solid earth pole tides (IERS 2010 conventions, C 21 &S 21 ), the ocean pole tides (Desai model, Nmax = 100), the atmosphere and oceanic variability (AOD1B RL06, Nmax = 180), and the general relativistic perturbations (IERS 2010 conventions). The background field model selected for this purpose is shown in Table 1. In order to reduce the high-frequency signal noise, the first degree and order 180 of the GGM05C model [50] is selected as the static gravity field model. The degree and order of the recovered time-variable gravity field model is 60. GRACE Level-1b data is processed and released by JPL, and the released version is from RL01 to RL03, where every new inversion of GRACE Level-1b data consistently brought clear improvements to gravity field estimates [51]. The latest version currently released is RL03, which improves the quality of the star camera and range rate data comparing with RL02 [52]. Accelerometer calibration is involved in the GRACE data processing. For the accelerometer calibration, we chose to calibrate the hourly biases [11,53]. The process of the recovered time-variable gravity field model is as follows. Firstly, the accelerometer parameters are calibrated according to the existing static gravity field model. Then, we recalculated the reference orbit based on the calibrated accelerometer parameters, and constructed the observation equations of the residual orbit and residual range rate. Finally, we recovered the time-variable gravity field by using the least-squares [20].
In order to ensure that there was enough data to establish the observation equations when the longer arc length is used to recover the time-variable gravity field model, we choose the GRACE data of January 2010 which has fewer discontinuities to recover timevariable gravity fields by using 6, 12, 24, 36, and 48 h arcs, respectively. We present the results according to processing the low-frequency error of the residual range rate or not, where the low-frequency error of the residual range rate is removed. In addition, we select the RL06 model released by CSR, GFZ and JPL, and ITSG-2018 model for comparison. The degree variance of the recovered time-variable gravity field is shown in Figure 8.
In Figure 8, whether or not the low-frequency error of the residual orbit is removed, the accuracy of the recovered time-variable gravity field by using 6, 12 and 24 h arcs is basically the same. However, the difference of the degree variance of the recovered time-variable gravity field gradually increases after 36 h arcs. When removing the low-frequency error of the residual orbit, it can be seen that the degree variance of the recovered time-variable gravity field for different arc lengths has a good consistency with CSR, JPL, GFZ and ITSG-2018.

The Monthly Time-Variable Gravity Field Analysis
According to the above simulation, the accuracy of the recovered gravity field can be improved, especially the accuracy of the multiples of the order 15. From the real data processing, the accuracy of the recovered time-variable gravity field models can maintain a good consistency for different arc lengths. In order to further evaluate the performance of solutions derived utilizing this method to consider low-frequency errors, we solved the monthly gravity field model from the periods January 2005 to February 2010, and compared them with the official agencies in terms of the annual amplitude of the mass change and the time series of the characteristics of typical basins. Due to the gaps in the GRACE data and in order to ensure that there are enough data to recover time-variable gravity fields [14], it is not guaranteed to use 48 h arcs to recover time-variable gravity field. Therefore, we used a shorter arc length to recover the time-variable gravity field in some months. According to the statistics, about 35% of the time-variable gravity field is used over 24 h arcs (where 30% is 48 h arcs and 5% is 36 h arcs), about 46% uses 24 h arcs, and about 19% uses less than 24 h arcs (where 6% is 6 h arcs and 13% is 12 h arcs). Because the time-variable gravity field models have exceeded the degree and order 60 provided by CSR, JPL, and GFZ, they are truncated to the degree and order 60 in order to ensure a fair comparison. To get the global mass changes, the post-process of time-variable gravity, including: (1) the C 20 is replaced by the result of SLR [54]; (2) the effect of north-south stripes is processed proposed by Swenson and Wahr [55]; (3) the high-frequency noise is filtered by Gaussian smoothing with a radius of 300km; (4) the spherical harmonic coefficients are converted into a global mass changes equivalent water height [56].
In Figure 9, we show the annual amplitudes of the global mass changes from UCAS_IGG, CSR, GFZ, and JPL, respectively. From the spatial distribution, similar large amplitudes are shown, such as in the Amazon basin, central Africa, Tibet Plateau, and north Australia. The correlation coefficients of the mean annual amplitudes of the global mass change among UCAS_IGG and CSR, GFZ, JPL is 0.9942, 0.9933, and 0.9940, respectively. For further comparisons, we chose the Amazon Basin, Yangtze Basin, Greenland, and Sahara Desert as the study areas [57,58]. The time series are shown in Figure 10, the UCAS_IGG has a good agreement with other models. The correlation coefficients between UCAS_IGG and CSR, GFZ, JPL are shown in Table 2. The correlation coefficients are roughly above 0.9 in the Yangtze River Basin, Greenland, and the Amazon Basin. In Sahara Desert, the standard deviations are 1.58 cm, 1.48 cm, 1.51 cm, and 1.54 cm, respectively. The four institutions are basically at the same level.

Discussion
In this study, we analyzed the influence of the low-frequency error of the residual orbit on recovering the time-varying gravity field model by the dynamic approach. According to the results of simulation and real data processing, removing the low-frequency error of the residual orbit or not is related to the arc length used. For the long arc length (such as 48 h arcs), there would be less variables in the initial state when establishing the observation equations, which made the solutions of the observation equations easier. Besides, the method of using the longer arc length could minimize the influence of random errors and means that some small force can be captured [59]. Due to the errors existing in the background model, there are low-frequency errors in the residual orbit and the residual range rate [30,39], which are usually neglected for the shorter arc length. However, this should be modeled for the longer arc length, such as an arc length of more than 24 h. In fact, the low-frequency terms of the residual orbit generated by the time-varying signal and perturbation model error are the same according to formulas (2)-(3), thus it is difficult to separate the perturbation model error and the time-varying signal from the low-frequency part of the residual orbit (within 2CPR) [30,60]. The low-frequency error removing method proposed in this paper not only removes the influence of the perturbation model errors, but also removes the influence of the time-varying signals. However, this method is feasible in terms of the accuracy of the existing background field model, and the arc length of the recovered time-varying gravity field can be prolonged and the time-varying signal can be accurately obtained. In the future, we will need to focus on evaluating the relationship between the accuracy of the background model and the low-frequency error removal method.

Conclusions
Based on the orbital perturbation differential equation, we deduced the observation equations of the residual orbit and the residual range rate with the non-linear corrections, respectively. The general solution of the Hill Equation is solved by neglecting the influence of the magnitude of the J 2 term in the satellite moving coordinate system. Under the condition of the perturbation model error and the instrument noise that includes a constant term and 1CPR term, the low-frequency error of the residual orbit is derived in the satellite moving coordinate system. Then, the low-frequency error of the residual range rate is also derived according to the low-frequency error of the residual orbit. There are resonance terms in the low-frequency error formulas of the residual orbit and the residual range rate. As the arc length increases, these resonance terms rapidly enlarge, and thus affect the accuracy of the recovered time-variable gravity field, especially the accuracy of the multiples of the order 15. We recommend that the low-frequency error processing of the residual orbit and residual range rate should both be taken into consideration on recovering time-variable gravity field, especially for the longer arc length.
We use the simulation and real data processing to explain the influence of the lowfrequency error of the residual orbit. In the simulation, the accuracy of the recovered gravity field can be improved by removing the low-frequency error, whether only the orbit data or combined orbit data and range rate data are used. In the real data processing, the time-variable gravity field models are recovered by using 6, 12, 24, 36, and 48 h arcs respectively based on January 2010 GRACE data. The accuracy of the recovered timevariable gravity field models by removing the low-frequency error of the residual orbit can always maintain a good consistency for different arc lengths.
We developed the time-variable gravity field model UCAS_IGG from January 2005 to February 2010 by removing the low-frequency error of the residual orbit and the residual range rate. By comparing the mean amplitudes of the global mass changes and the time series of the characteristics of the basin, our recovered time-variable signal levels have a good consistency with the international agencies.