Whole cycle ambiguity fixed orbit and decorrelation random model gravitational field inversion method

Abstract

The invention relates to the cross technical field of satellite gravity mechanics, satellite navigation science and the like, in particular to a whole-cycle ambiguity fixed orbit and decorrelation random model gravity field inversion method capable of remarkably improving the overall precision of gravity field recovery, and adopts a celestial mechanics method CMA to recover the earth gravity field. In the ambiguity solution of the single receiver, a single-difference IAR method is adopted, and a high-precision solution is obtained through ambiguity interpolation between fixed arc segments; using a variance component estimation (VCE) frame to iteratively update the observation value weight; according to the method, the integer ambiguity is fixedly introduced into the orbit reconstruction process, the precision of the GFO satellite orbit in the along-orbit direction and the transverse-orbit direction can be remarkably improved, meanwhile, empirical covariance estimation is carried out based on the residual error after fitting, the uncertainty of a background model is absorbed more reasonably, and therefore the overall precision of gravity field recovery is improved.

Description

Technical Field

[0001] This invention relates to the interdisciplinary fields of satellite gravimetry and satellite navigation, and more specifically to a gravity field inversion method using integer ambiguity fixed orbit and decorrelation stochastic model that can significantly improve the overall accuracy of gravity field recovery. Background Technology

[0002] The development of Global Positioning System (GPS) technology has significantly improved the orbit determination accuracy of Low Earth Orbit (LEO) satellites. Using geodetic-grade GPS receivers, centimeter-level precision orbit determination (POD) can be achieved. This accuracy not only supports various applications such as remote sensing, communication, and sea surface altimetry, but also plays a crucial role in long-wave gravity field signal recovery. Missions such as CHAMP, GRACE, GOCE, and GRACEFollow-On (GFO) have all verified the core role of GPS in satellite gravity field recovery.

[0003] LEO precise kinematic orbits are a crucial foundation for Earth's gravity field inversion. Unlike simplified or fully dynamic methods, kinematic orbits are unaffected by dynamic model biases and fully preserve variance-covariance information. In the Celestial Mechanics Approach (CMA), these orbits are treated as quasi-observations, and their formalized covariance information can be derived from variance propagation, enabling robust simultaneous estimation of satellite orbital parameters and spherically harmonic gravity field coefficients. However, the accuracy of kinematic orbits used for gravity field recovery is often limited. Because the orbits originate from unambiguous GPS pseudorange and carrier phase observations, their accuracy is significantly lower than that of integer ambiguity resolution (IAR) solutions. Nevertheless, floating-point solution orbits remain the standard input for many operational gravity field products.

[0004] Besides orbital accuracy, the accurate construction of the stochastic model is equally crucial. In gravity field recovery, the orbital covariance matrix can serve as prior information. However, due to the enormous storage and computational overhead of the full-moment covariance, simplified models are typically used in research, considering only intra-epoch correlations while ignoring inter-epoch correlations. To better describe the inter-epoch correlation structure, researchers often supplement this with empirical covariance estimation, i.e., constructing an empirical model based on the a priori residual autocovariance of the dynamic solution. This method is also common in KBR-based gravity field recovery, used to compensate for uncertainties in the background model. In practice, these empirical covariances are also simplified, typically considering only correlations within approximately 50 minutes. Although combining inter-epoch covariance with empirical stochastic model estimation can improve low-order gravity field recovery, prior information has not yet been fully utilized.

[0005] Some studies have explored the role of IAR in gravity field recovery, but the conclusions show that its potential and limitations coexist. Lasser et al. (2020) found that IAR can reduce long-period noise in kinematic orbits, but considering only short-term correlations is insufficient to fully characterize the noise characteristics, resulting in gravity field results comparable to floating-point solutions. At the 28th International Union of Geodesy and Geophysics (IUGG) General Assembly, Lasser et al. (2023) reported that IAR only marginally improves the accuracy of low-order gravity fields. Meanwhile, Gao et al. (2023) pointed out that IAR in dynamic methods may amplify resonance errors in gravity field coefficients. Lasser et al. (2024) further showed that the gravity field signal of the IAR solution is similar to that of the floating-point solution. It is evident that although IAR continues to improve orbit accuracy and reduce correlations in covariance, its negative impact on gravity field inversion results still needs attention. In particular, the potential for high-precision gravity field reconstruction using a combination of kinematic orbit and inter-satellite K-band range-rate (KRR) data has not yet been systematically evaluated. Summary of the Invention

[0006] This invention addresses the shortcomings and deficiencies of existing technologies. Within the CMA framework, it examines the impact of kinematic trajectory accuracy and covariance information integrity on gravity field inversion, and proposes a gravity field inversion method using integer ambiguity fixed trajectories and decorrelation stochastic models. This invention achieves its purpose through the following measures: A gravity field inversion method using integer ambiguity fixed orbit and decorrelation stochastic model employs the celestial mechanical method (CMA) to recover the Earth's gravity field. Its key feature is the use of a single-difference interpolation (IAR) method in single-receiver ambiguity resolution, employing ambiguity interpolation between fixed arc segments to obtain a high-precision solution; and the use of a variance component estimation (VCE) framework to iteratively update observation weights. Specifically, the celestial mechanical method (CMA) for recovering the Earth's gravity field involves: first, processing GPS observations to derive the kinematic low-Earth orbit (LEO) satellite position and its covariance information at the epoch; second, using the LEO satellite position as a virtual observation and weighting it according to its covariance information, and then, within the generalized dynamic orbit determination framework, utilizing these virtual observations... Daily normal equations are constructed to solve for unknown gravity field coefficients. Based on the timeliness of the parameters, the estimated parameters are divided into local parameters and global parameters. Local parameters include the kinematic empirical parameters of satellite position and velocity, K-band range variability (KRR), and accelerometer bias parameters. These local parameters are eliminated before estimating the global parameters. Finally, the daily normal equations are corrected and accumulated to form a monthly normal equation system to solve for the global parameters, including accelerometer scale parameters and gravity field spherical harmonic coefficients. By solving the combined normal equations and substituting back, both global and local parameters can be recovered. After parameter estimation, the residuals are calculated using the observation equations, and the variance component (VCE) is estimated. The weighting factors are then updated, and the next iteration begins.

[0007] In this invention, the monthly full-scale factor matrix was estimated, which includes three diagonal terms, three symmetric shear terms, and three antisymmetric rotation terms. All three axes were modeled using quadratic polynomials for bias, as shown in the following formula:

[0008] Where, represents the deviation value of the final calibration, , and These are the coefficients of the constant term, linear term, and quadratic term corresponding to the axis. It is the current time. It is the start time of the parameter estimation arc segment. Different estimation intervals are set for different axes: the most sensitive X-axis deviation factor is estimated every 24 hours, the second most sensitive Z-axis is estimated every 12 hours, and the least sensitive Y-axis is estimated every 6 hours.

[0009] The present invention employs a single-difference IAR method in single-receiver ambiguity resolution to obtain a high-precision solution by fixing the ambiguity difference between arc segments. The specific process is as follows: First, the observation data is divided into continuous carrier phase tracking arc segments, and the average HMW combination value of each arc segment is calculated using dual-frequency pseudorange and carrier phase observations. Then, the inter-arc difference is corrected by combining the specific wide-lane deviation of GPS satellites, so that the wide-lane ambiguity is fixed to an integer value. On this basis, the inter-arc difference is determined and constrained to an integer value using these fixed wide-lane ambiguities and the floating-point ambiguity estimate of the non-ionospheric IF combination carrier phase obtained by least squares adjustment. Once the HMW and narrow-lane NL ambiguities are successfully fixed, the IF ambiguity can be reconstructed, and the normal equation is introduced with high weight as a constraint.

[0010] This invention scales the empirical observation error covariance based on the residual statistics obtained from data processing, and uses a variance component estimation (VCE) framework to iteratively update the observation weights. For observation type i, the initial prior weighting matrix... Assuming it is a diagonal matrix ， It is given by the following formula:

[0011] in, Representation type empirical standard deviation This represents the total number of observation types. For kinematic orbit data, the daily covariance matrix of each satellite is considered as an independent observation type. Within a one-month arc, 60 kinematic datasets for the binary satellites (GFO-C and GFO-D) are considered, with their initial weights derived from the inverses of their respective covariance matrices. For KRR observations, epoch-independent observations are assumed, and an identity matrix is ​​assigned as the initial weighting matrix (i.e., unit weights). The weighting factor after one iteration is defined as:

[0012] in, Represents the weighting factor. It is an observation-subtraction computation (OC) vector. ， Indicates the parameter to be estimated. It is the corresponding design matrix ， It represents the number of observations for each observation type. It is the total number of types, the observation weight matrix. The weighting factors are updated as follows:

[0013] in, It requires multiple iterations to converge.

[0014] This invention demonstrates that introducing integer ambiguity fixing into the orbit reconstruction process can significantly improve the accuracy of GFO satellite orbits in the along-orbit and transverse directions. However, this improvement does not automatically translate into a higher-quality gravity field solution in GPS-only gravity field recovery. The ambiguity fixing process does not completely eliminate time correlation, but only slightly weakens it, while limiting the ability to model these correlations. Although IAR significantly improves the solution accuracy, it does not necessarily lead to a better characterization of the redistribution of Earth's mass. This invention uses empirical covariance estimation based on the fitted residuals to more reasonably absorb the uncertainty of the background model, thereby improving the overall accuracy of gravity field recovery. Attached Figure Description Figure 1 It is the radial, along-track, and transverse RMS distribution (cm) of the predicted orbits of the GFO-C satellite floating-point solution and IAR solution relative to the JPL Precision Scientific Orbit (PSO).

[0015] Figure 2 It is the KBR residual amplitude spectral density (ASD) (cm) of the floating-point solution and IAR solution (30 s orbit) of the GFO satellite on November 1, 2019, with the frequency expressed in cycles per orbit (CPR).

[0016] Figure 3 The kinematic orbital covariance structure and epoch correlation of GPS data sampled by GFO-C satellite for 30 s on November 16, 2019 are shown: (a) is the floating-point solution covariance matrix (log10); (b) is the IAR solution covariance matrix (log10). The horizontal and vertical axes are both orbital epoch numbers (0–8636). The time correlation of the floating-point solution can be extended to about 90 min, while the IAR solution is limited to about 3 min.

[0017] Figure 4 It is the radial, along-orbit, and transverse RMS distribution of the orbit relative to the JPL Precision Science Orbit (PSO) after fitting the GFO-C satellite floating-point solution with three covariance strategies (epoch independent, 1-hour epoch, and full moment) and the IAR solution.

[0018] Figure 5 This is a comparison of the geoid height difference (mm) of the GPS-only kinematic orbit gravity field solution of the GFO-C satellite in November 2019 under different configurations with the static reference model GGM05C. The floating-point solution includes three covariance strategies (epoch independent, 1-hour epoch, full matrix) and is compared with the IAR solution.

[0019] Figure 6It is the radial, along-orbit, and transverse RMS distribution (cm) of the orbits of the GFO-C satellite after fitting three schemes (Float-OW, IAR-OW, IAR-FW) in November 2019 relative to the JPL Precision Science Orbit (PSO).

[0020] Figure 7 It is the residual distribution (μm / s) after KRR fitting for different processing schemes (Float-OW, IAR-OW, IAR-FW) of GFO-C satellite.

[0021] Figure 8 This is a comparison of the successive geoid height differences (mm) of the GFO satellite monthly gravity field solution in November 2019: JPL, CSR, and GFZ RL06.1 models, as well as the three schemes of this invention (Float-OW, IAR-OW, IAR-FW), with the static reference model GGM05C.

[0022] Figure 9 This is a comparison of the successive geoid height differences (mm) of the monthly gravity field solution from the GFO satellite in November 2019 with three schemes (Float-OW, IAR-OW, IAR-FW) relative to the JPL RL06.1 model.

[0023] Figure 10 The differences (cm) in equivalent water height (EWH) of the GFO satellite monthly gravity field solution in November 2019 are as follows: (a–f) differences between the three official RL06.1 products (JPL, CSR, GFZ) and the three schemes (Float-OW, IAR-OW, IAR-FW) and the static reference model GGM05C; (g–i) differences between the three schemes and the JPL RL06.1 model. All results are filtered by a 350km Gaussian filter.

[0024] Figure 11 This is the difference (mm / s) between the inter-satellite relative velocity derived from the floating-point solution and IAR solution dynamic velocity in the GPS-only case of the GFO satellite in November 2019 and the observed K-band distance variability (KRR). Detailed Implementation

[0025] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0026] This study investigates the impact of kinematic orbit accuracy and covariance information completeness on gravity field inversion within the CMA framework. Using GRACE-FO GPS and KRR data from November 2019, monthly gravity field solutions of orders up to 60 (GPS-only) and 96 (GPS+KRR) were constructed. First, the importance of covariance information was assessed, and a covariance selection method based on inter-epoch correlation analysis was proposed to balance completeness and computational complexity. Then, the performance of floating-point and IAR orbits and their covariances within the CMA framework was compared, and the results were evaluated using orbit accuracy, apologetic KBR and KRR residuals, and the gravity field.

[0027] In the CMA method, the core of gravity field recovery is solving the generalized orbit determination problem for each GRACE satellite. This example uses a three-step method to recover the Earth's gravity field from GRACE Follow-On (GFO) GPS data: First, the GPS observations are processed to derive the kinematic low Earth orbit (LEO) satellite positions and their covariance information at each epoch. Second, these kinematic positions are used as virtual observations and weighted according to their covariance information. Within the framework of generalized dynamic orbit determination, daily normal equations are constructed using these virtual observations to solve for the unknown gravity field coefficients. Based on the timeliness of the parameters, the estimated parameters are divided into local parameters and global parameters. The local parameters include the satellite position and velocity, the kinematic empirical parameters of the K-band range-rate (KRR), and the accelerometer bias parameters. To reduce the size of the normal equations, these local parameters are eliminated before estimating the global parameters. Finally, the daily normal equations are corrected and accumulated to form a monthly normal equation system to solve for the global parameters, including accelerometer scale parameters and gravitational field spherical harmonic coefficients. By solving the combined normal equations and substituting back, both global and local parameters can be recovered. After parameter estimation, the residuals are calculated using the observation equations, and variance component estimation (VCE) is performed. The weighting factors are then updated accordingly, and the next iteration begins. Typically, 4–6 iterations are required to ensure convergence. It should be noted that no constraints are imposed on the gravitational field parameters throughout the process to maintain the physical fidelity of the solution.

[0028] This example estimates the monthly full-scale factor matrix, which includes three diagonal terms, three symmetric shear terms, and three antisymmetric rotation terms. For bias modeling, an empirical strategy was developed based on extensive experiments. Specifically, all three axes are modeled using quadratic polynomials, as follows:

[0029] in, This represents the deviation value of the final calibration. , and These are the coefficients of the constant term, linear term, and quadratic term corresponding to the axis. It is the current time. This is the start time of the parameter estimation arc. Different estimation intervals were set for different axes: the most sensitive X-axis deviation factor was estimated every 24 hours, the second most sensitive Z-axis every 12 hours, and the least sensitive Y-axis every 6 hours. To account for the uncertainty of the background model, empirical KBR kinematic parameters were applied—including the KRR constant deviation and linear drift estimated every 45 minutes, and the 1-CPR sine term estimated every 90 minutes—consistent with Kim (2000).

[0030] Integer Ambiguity Resolution (IAR) is a crucial step in this example. In single-receiver ambiguity resolution, a single-difference IAR method is employed to obtain a high-precision solution by fixing the ambiguity differences between arc segments. The specific process is as follows: First, the observation data is divided into continuous carrier phase tracking arc segments, and the average Hatch-Melbourne-Wübbena (HMW) combination value for each arc segment is calculated using dual-frequency pseudorange and carrier phase observations. Then, the inter-segment differences are corrected using GPS satellite-specific wide-lane bias corrections provided by Wuhan University (WHU) products, fixing the wide-lane ambiguity to integer values. Based on this, the inter-segment differences are determined and constrained to integers using these fixed wide-lane ambiguities and the ionosphere-free (IF) combined carrier phase floating-point ambiguity estimates obtained through least-squares adjustment. Once the HMW and Narrow-Lane (NL) ambiguities are successfully fixed, the IF ambiguities can be reconstructed and introduced into the normal equations with high weights as constraints.

[0031] Unlike simplified dynamic PODs, ambiguity fixation in kinematic PODs faces significant challenges, as its stability is affected by the quantity and quality of spaceborne GNSS observations and the accuracy of modeling. To overcome these issues, an improved strategy proposed by Gao et al. (2023) is adopted, incorporating the clock bias modeling method of the GRACE satellite ultra-stable oscillator (USO). This method reduces the number of receiver clock bias parameters that need to be estimated in the kinematic observation equations, thereby improving the accuracy of ambiguity resolution and the fixation success rate. Once the ambiguity is reliably resolved, it is substituted into the observation equations, while simultaneously removing the constraints related to random walk clock bias in the normal equations.

[0032] In addition to the methods mentioned above, some scholars have proposed another approach to ambiguity resolution. Arnold et al. (2019) based their approach on the floating-point ambiguity of the simplified dynamic POD solution (which is usually more accurate than the kinematic POD solution). They first fixed the ambiguity and then introduced the successfully resolved ambiguity into the kinematic processing to reconstruct the unbiased carrier phase observation. This method is particularly suitable when it is necessary to preserve kinematic properties.

[0033] In the presence of heterogeneous observation types, optimizing weighting is crucial to ensuring that each dataset contributes appropriately to the estimation process. The empirical observation error covariance is scaled based on the residual statistics obtained from data processing. The observation weights are iteratively updated using a variance component estimation (VCE) framework. For observation type i, the initial prior weighting matrix... Assuming it is a diagonal matrix ， It is given by the following formula:

[0034] in, Representation type empirical standard deviation This represents the total number of observation types. For kinematic orbit data, the daily covariance matrix of each satellite is considered as an independent observation type. Within a one-month arc, 60 kinematic datasets for the binary system (GFO-C and GFO-D) are considered, with initial weights derived from the inverses of their respective covariance matrices. For KRR observations, epoch-independent observations are assumed, and an identity matrix is ​​assigned as the initial weighting matrix (i.e., unit weights). Observation Type The weighting factor after one iteration is defined as:

[0035] in, Represents the weighting factor. It is an observation-subtraction computation (OC) vector. ， Indicates the parameter to be estimated. It is the corresponding design matrix ， It represents the number of observations for each observation type. This represents the total number of types. Observation weight matrix. The weighting factors are updated as follows:

[0036] in, It requires multiple iterations to converge.

[0037] The accuracy of the kinematic trajectory and its covariance information is an important prior condition for trajectory reconstruction and gravity field recovery. Therefore, before experimental analysis, it is necessary to examine the characteristics of the kinematic trajectory and its covariance matrix for both the floating solution and integer ambiguity resolution (IAR). Figure 1 A comparison of IAR and non-IAR orbits with respect to Precise Science Orbits (PSOs) is presented. PSOs are generated by the Jet Propulsion Laboratory (JPL) using a simplified dynamical POD method that incorporates double-difference IARs and is validated using satellite laser ranging (SLR) and K-band ranging (KBR) residuals, achieving accuracies of approximately 1 cm and 2 m, respectively. Given the extremely high accuracy of PSOs, they are used as a reference benchmark. Figure 1 It is evident that IAR can significantly improve the accuracy of kinematic orbits. Taking the GFO-C satellite as an example, IAR reduced the radial error from 2.3 cm to approximately 1.5 cm, the along-orbit error from 2.4 cm to approximately 1.0 cm, and the transverse error from 1.5 cm to approximately 0.8 cm.

[0038] The KBR ranging system of the GFO satellite can measure inter-satellite offset distances independently of GPS with micrometer-level accuracy. Due to its high sensitivity along the orbital direction, KBR residuals are often used to evaluate orbital accuracy in that direction. After fixing integer ambiguity, the standard deviation of the KBR residuals decreased from 2.0 cm to 1.0 cm, further validating the accuracy improvement effect of the IAR orbit. Furthermore, Figure 2 The residual spectrum analysis shown indicates that, compared with the floating-point solution, IAR significantly suppresses low-frequency noise (less than 2 cycles per revolution, CPR), while having little effect on high-frequency components, which is consistent with the conclusions of Lasser et al. (2020).

[0039] To further analyze the covariance characteristics, the normal equations after eliminating receiver clock bias and ambiguity parameters from the normal equation matrix are as follows:

[0040] in, The normal equation matrix is ​​obtained after eliminating the GFO receiver clock bias and ambiguity parameters. The remaining parameters to be estimated are the orbital corrections for each epoch. Covariance matrix. Initially represented in the Earth-fixed system, the orbital corrections were then transformed into radial (r), along-orbit (a), and transverse (c) components by rotating the unit vector derived from the simplified dynamic orbit of JPL to the Local Orbit Frame (LOF).

[0041] Figure 3 a presents the floating-point covariance structure of the GFO-C satellite on November 16, 2019. Figure 3 b represents the corresponding Integer Ambiguity Fixed (IAR) solution. Both are expressed with absolute values ​​of covariance and logarithmic scale (log10). The observation quality for this day was high, with 2879 epochs successfully processed, and the narrow-lane ambiguity fixing success rate reached 95%. From Figure 3 As can be seen, the floating-point solution exhibits significant inter-epoch correlation: the values ​​are higher near the diagonal and gradually decrease with increasing time interval, reflecting the correlation caused by the ambiguity of the unsolved arc segment. In contrast, the IAR solution ( Figure 3 b) shows a significant decrease in the off-diagonal terms, indicating that fixing the ambiguity effectively weakens the temporal correlation. However, its covariance matrix still exhibits a local blocky structure near the diagonal, because some ambiguities were not fixed. When there are fewer than 4 fixed ambiguities in a single epoch, a certain degree of correlation will still be retained between adjacent epochs.

[0042] To quantitatively characterize the temporal evolution of correlation between epochs, the correlation curves between the floating-point solution and the IAR solution were calculated in three directions: radial (r), along the orbit (a), and transverse (c). The correlation function is defined as follows:

[0043] in, Indicates the interval The number of covariances in each epoch. Let be the covariance of the u and v components between epochs. Figure 3 c and Figure 3 Figure d presents the six covariance function components of the GFO-C satellite in the local orbit coordinate system, corresponding to the floating-point solution and the IAR solution, respectively. The results show that the correlation of the floating-point solution gradually decreases with time lag. At minute 1 minute, the amplitudes of all covariance components decreased to below 1 mm², and at that time... The correlation approaches zero at the minute mark. The correlation is most significant in the horizontal direction, which is related to... Figure 3 The lattice structure in solution 'a' is consistent; however, the correlations in the radial and along-track directions are weak, and the cross-correlation with the transverse track is almost negligible. In contrast, almost all six components of the IAR solution converge to zero within 30 s, indicating that fixed ambiguity can effectively eliminate inter-epoch correlations in floating-point solutions, making the system more statistically rigid. Figure 3 c and Figure 3The comparison of d shows that the floating-point solution preserves long-range correlations on the tens of minutes scale, while the IAR solution almost completely suppresses such effects.

[0044] This example evaluates the gravity field recovery performance of a 30-second sampled kinematic orbit in November 2019, focusing on comparing the differences between floating-point solutions and integer ambiguity fixed-AR solutions in orbit reconstruction and gravity field estimation. Since the results of GFO-D are basically consistent with those of GFO-C, this invention only presents the solution results of GFO-C. For floating-point solutions, the full-moment covariance matrix can fully characterize the inter-epoch correlation, but the computational cost is extremely high. For example, for a day's data sampled for 30 seconds, the covariance matrix contains about 37 million elements, requiring about 280 MB for storage alone. To alleviate this problem, a truncated stochastic model is used, retaining only the correlation within 1 hour. Three configurations were tested for floating-point solutions: (1) independent covariance between epochs; (2) 1-hour inter-epoch covariance; and (3) full-moment covariance. For IAR solutions, the independent covariance model between epochs is used. Based on these four configurations, fitted orbits and unconstrained spherical harmonic gravity field models up to order 60 were estimated within the CMA framework.

[0045] Figure 4 The RMS distribution of orbital discrepancies after fitting the GFO-C satellite is presented. The results show that, within the floating-point solution framework, the independent covariance configuration between epochs yields the worst accuracy in the radial and along-orbit directions; in contrast, the 1-hour covariance and full-moment covariance significantly improve orbital accuracy in these two directions. However, the improvement is limited in the transverse direction, indicating that this direction is insensitive to time-dependent modeling. Notably, the accuracy of the 1-hour covariance solution is almost equivalent to that of the full-moment covariance solution, demonstrating that this approximation scheme achieves comparable accuracy to the full-moment method while significantly reducing computational and storage requirements.

[0046] Figure 5This paper presents the order-wise geoid height differences of the gravity field model based on the GPS-only solution from the GFO-C satellite under different configurations. The results show that for the floating-point solution, the three covariance configurations maintain good consistency at orders 20 and below; after this, the error of the epoch-independent covariance solution increases significantly, while the 1-hour epoch-interval covariance solution shows a significant improvement. The 1-hour covariance and full-moment covariance solutions are almost indistinguishable at orders 60 and below, verifying the effectiveness of the 1-hour approximation method in maintaining the accuracy of the gravity field model. The IAR-based model maintains high consistency with the 1-hour and full-moment covariance floating-point solutions at orders 45 and below, but shows significant deviations at higher orders. This may stem from the systematic error introduced by IAR during the kinematic POD process: although it effectively reduces orbital noise, it introduces a systematic impact on high-frequency gravity signals. This issue will be further analyzed in the discussion section. It should be noted that the effective resolution of gravity field recovery under GPS-only observation conditions is typically limited to the 20th–30th order (corresponding to a spatial resolution of approximately 1300–2000 km), a conclusion confirmed by missions such as CHAMP, GRACE, GOCE, and Swarm. Therefore, different stochastic model configurations have a limited impact on the overall accuracy of the gravity field. Since the estimation of lower-order terms is highly sensitive to other systematic errors (such as insufficient modeling of antenna phase center changes), Figure 5 The differences shown to some extent obscure the role of covariance modeling, making it difficult to assess its impact on low-order gravitational fields independently.

[0047] This invention evaluates the accuracy of gravity field reconstruction after fusing GFO satellite orbit and KBR ranging observations within a joint estimation framework. Compared to methods relying solely on GPS, introducing K-band distance variability observations with micrometer / second precision significantly improves sensitivity to high-frequency gravity signals. Using GRACE-FO L1B KRR data (5-second sampling) from November 2019, this invention combines orbit and covariance information from floating-point and integer ambiguity fixed (IAR) solutions to jointly invert to a 96th-order monthly gravity field model within a CMA framework.

[0048] In the floating-point solution, a 1-hour epoch covariance structure is used because this method has been proven to approximate the full-moment covariance in terms of recovery accuracy. The IAR solution uses independent epoch covariance to fully utilize the decorrelation properties of the orbit after ambiguity fixation. Since the IAR orbit has higher accuracy, kinematic position is given higher weight in variance component estimation, thus relatively weakening the contribution of KBR. To balance the effects of the two types of observations, this invention introduces a fixed-weight (FW) strategy: for the IAR-FW solution, the standard deviations of kinematic position and KBR observations are fixed at 5 cm and 0.05 μm / s, respectively, to enhance the influence of KBR in the solution. This scheme serves as a comparison with the optimal weight (OW) strategy for KBR.

[0049] First, the quality of the orbital solution was evaluated by comparison with the JPL Precision Scientific Orbit (PSO) and by analyzing the KBR and KRR residuals. Then, the lunar gravity field solution was compared with the official solution published by the GRACE-FO Science Data System (SDS), including geoid height differences and spatial domain consistency. SDS products include CSR RL06.1 (University of Texas Space Research Center), GFZ RL06.1 (Germany Geosciences Centre), and JPL RL06.1.

[0050] Figure 6 The orbital accuracy of the GFO-C satellite after fitting under the GPS+KBR joint processing framework is demonstrated (due to the consistent performance of GFO-D). The results show that, compared to Float-OW, IAR-OW significantly improves orbital reconstruction accuracy in both the along-orbit and transverse directions. The RMS in the along-orbit direction decreases from 1.4 cm to 0.7 cm (a 50% improvement), and in the transverse direction from 1.9 cm to 0.6 cm (a 68% improvement). In the radial direction, IAR-OW has a negligible effect, showing only a slight improvement. The results of IAR-FW are almost identical to those of IAR-OW, validating the rationality of the fixed weights (kinematic position 5 cm, KRR 0.05 μm / s) set in this invention.

[0051] Figure 8 This paper presents the order-wise geoid high error of the gravity field model based on the GFO GPS+KBR joint solution from November 2019, compared to the static reference model GGM05C. The results show that the Float-OW solution maintains high consistency with the GRACE Scientific Data System (SDS) models (JPL, CSR, GFZ RL06.1) across the entire spectral range (up to order 96). The IAR-OW solution exhibits comparable accuracy to the SDS model at orders 25 and below, but shows significant deviations at higher orders, indicating limitations of the IAR method in recovering high-frequency gravity signals. Combined with... Figure 7 The KRR residual analysis in the model shows that the optimal weighting (OW) strategy failed to fully utilize the role of KRR observations when using IAR orbits. This may be because the orbits are assigned relatively high weights in the adjustment, thus weakening the contribution of KRR. In contrast, the IAR-FW solution is highly consistent with the SDS model at almost all orders, highlighting its advantage in enhancing the relative contribution of KRR. To further clarify the differences between the three configurations, Figure 9The successive geoid height differences between the Float-OW, IAR-OW, and IAR-FW solutions and the JPL RL06.1 model are presented. The results show that the IAR-FW configuration has the best agreement with the official JPL model. Notably, both the IAR-OW and IAR-FW solutions exhibit significant spectral anomalies around the 46th order.

[0052] In spatial domain analysis, Figure 10 (a)–(f) illustrate the difference in equivalent water height (EWH) between the GFO satellite monthly gravity field solution in November 2019 and the static gravity model GGM05C. The comparison schemes include three official RL06.1 products (JPL, CSR, GFZ) and the solutions obtained in this invention based on three processing strategies (Float-OW, IAR-OW, IAR-FW). These distribution plots are used to assess the consistency of different solutions in characterizing global-scale mass redistribution.

[0053] from Figure 10 As can be seen from (a)–(f), the six solutions are highly consistent in characterizing the dominant features of large-scale mass changes, verifying their ability to capture the main gravitational signals. Statistically, the results of this invention are largely consistent with the extreme values, mean, and RMS of the official model, indicating the reliability of the obtained gravitational field. However, differences still exist. For example, the mean EWH of the Float-OW and IAR-OW solutions is approximately 0.017 cm, about 30% smaller than the 0.025 cm of the RL06.1 product. This difference may stem from different accelerometer parameter modeling strategies, while the SDS product employs a different method. Furthermore, the IAR-OW solution ( Figure 10 e) It exhibits more pronounced north-south strip noise than Float-OW and IAR-FW. The EWH RMS is slightly higher (38.724 cm, compared to 37.88 cm for Float-OW and 38.05 cm for IAR-FW), indicating that IAR-OW weakens the contribution of KRR observations in least squares adjustment, thereby exacerbating strip noise.

[0054] To directly compare the noise characteristics of the three schemes, this invention further calculated the EWH difference between the monthly solutions and the JPL RL06.1 model, and the results are as follows: Figure 10As shown in (g)–(i). The results indicate that the north-south striping of IAR-OW is the most prominent, followed by Float-OW, while the striping effect of IAR-FW is the weakest. The mean EWH values ​​of the three are 0.008 cm for Float-OW and IAR-OW, and 0.005 cm for IAR-FW, showing that IAR-FW reduces the bias by 38%. The corresponding RMS values ​​are 6.92 cm, 9.48 cm, and 5.40 cm, respectively, indicating that IAR-OW increases noise by 37% compared to Float-OW, while IAR-FW reduces the noise level by 22% compared to IAR-OW. In summary, these results highlight the importance of a reasonable weighting strategy when applying IAR orbits, with fixed weights (IAR-FW) achieving a more balanced fusion of GPS and KBR data, thereby significantly improving the stability of the solution.

[0055] This invention demonstrates that introducing integer ambiguity fixation into the orbit reconstruction process can significantly improve the accuracy of GFO satellite orbits in both along-orbit and transverse directions. However, this improvement does not automatically translate into a higher-quality gravity field solution in GPS-only gravity field recovery. Figure 5 As shown, compared to floating-point solutions, IAR orbits do not exhibit an advantage in high-frequency gravity signal recovery, a problem inevitably reflected in the orbit reconstruction results. The ambiguity fixing process does not completely eliminate temporal correlations, only slightly weakens them, and limits the ability to model these correlations. If the failure to improve gravity field recovery accuracy truly stems from IAR weakening key stochastic information (such as temporal correlations), then the orbit accuracy after fitting the IAR solution should be comparable to that of the floating-point solution. However, this is not the case; IAR orbits consistently demonstrate higher positioning accuracy.

[0056] While IAR significantly improves computational accuracy, it doesn't necessarily provide a better characterization of Earth's mass redistribution. IAR reduces the number of unknown parameters by fixing ambiguities, but simultaneously introduces additional background model dependencies and potential uncertainties (such as centroid correction and hardware bias). Therefore, after fixing the ambiguities, the GPS carrier phase residuals may actually increase. IAR makes the estimation system more statistically "rigid" and robust, which is particularly advantageous for orbit determination primarily targeting long-wavelength gravity signals. In contrast, floating-point solutions retain more degrees of freedom and flexibility. Although less robust and noise-resistant than IAR, they can fully propagate and characterize the correlations between parameters and epochs. In summary, IAR constructs a more robust system for orbit determination by reducing the number of parameters; while floating-point solutions maintain higher flexibility but are less robust. The differences in orbit accuracy and gravity field recovery performance reflect the methodological complementarity between IAR and floating-point solutions.

[0057] In GPS-only gravity field recovery, the standard deviation of the relative velocity residuals increases from 0.023 mm / s in the floating-point solution to 0.032 mm / s in the IAR solution. Figure 11 Time series of KRR observations and relative velocity differences from GFO-C / D satellites are presented. Since KRR observations are independent and extremely accurate, their residuals can serve as a reliable reference for assessing velocity accuracy. The results indicate that IAR may unintentionally reduce the accuracy of the velocity solution, possibly due to the propagation of systematic errors or the absorption of additional background model uncertainties into the solution process. To address this issue, this invention employs empirical covariance estimation based on the fitted residuals to more reasonably absorb background model uncertainties, thereby improving the overall accuracy of gravity field recovery.

Claims

1. A method for inverting the gravity field of an integer ambiguity fixed orbit and a decorrelation stochastic model, employing the celestial mechanical method CMA to recover the Earth's gravity field, characterized in that... In single-receiver ambiguity resolution, the single-difference IAR method is employed, using ambiguity interpolation between fixed arc segments to obtain a high-precision solution. The variance component estimation (VCE) framework is used to iteratively update the observation weights. Specifically, the celestial mechanical method (CMA) for Earth's gravity field recovery involves: first, processing GPS observations to derive the kinematic LEO satellite position and its covariance information at the epoch; second, using the kinematic LEO satellite position as a virtual observation and weighting it based on its covariance information, and within the generalized dynamic orbit determination framework, constructing daily normal equations using these virtual observations to solve for the unknown gravity field coefficients, based on the parameters... To ensure timeliness, the estimated parameters are divided into local and global parameters. Local parameters include the empirical kinematic parameters of satellite position and velocity, K-band range variation (KRR), and accelerometer bias parameters. These local parameters are eliminated before estimating the global parameters. Finally, the daily normal equations are corrected and accumulated to form a monthly normal equation system to solve for the global parameters, including accelerometer scale parameters and gravity field spherical harmonic coefficients. By solving the combined normal equations and substituting back, both global and local parameters can be recovered. After parameter estimation, the residuals are calculated using the observation equations, and the variance component (VCE) is estimated. The weighting factors are then updated accordingly, and the next iteration begins.

2. The gravity field inversion method for integer ambiguity fixed orbit and decorrelation stochastic model according to claim 1, characterized in that, The monthly full-scale factor matrix was estimated, containing three diagonal terms, three symmetric shear terms, and three antisymmetric rotation terms. Bias modeling was performed on all three axes using quadratic polynomials, as shown in the following formula: ， in, This represents the deviation value of the final calibration. , and These are the coefficients of the constant term, linear term, and quadratic term corresponding to the axis. It is the current time. It is the start time of the parameter estimation arc segment. Different estimation intervals are set for different axes: the most sensitive X-axis deviation factor is estimated every 24 hours, the second most sensitive Z-axis is estimated every 12 hours, and the least sensitive Y-axis is estimated every 6 hours.

3. The gravity field inversion method for integer ambiguity fixed orbit and decorrelation stochastic model according to claim 1, characterized in that, In the single-receiver ambiguity resolution, the single-difference IAR method is adopted to obtain a high-precision solution by fixing the ambiguity difference between arc segments. The specific process is as follows: First, the observation data is divided into continuous carrier phase tracking arc segments, and the average HMW combination value of each arc segment is calculated using dual-frequency pseudorange and carrier phase observations. Then, the difference between arc segments is corrected by combining the specific wide-lane deviation of GPS satellites, so that the wide-lane ambiguity is fixed to an integer value. On this basis, the difference between arc segments is determined and constrained to an integer value using these fixed wide-lane ambiguities and the floating-point ambiguity estimate of the non-ionospheric IF combination carrier phase obtained by least squares adjustment. Once the HMW and narrow-lane NL ambiguities are successfully fixed, the IF ambiguity can be reconstructed and introduced into the normal equation with high weight as a constraint.

4. The gravity field inversion method for integer ambiguity fixed orbit and decorrelation stochastic model according to claim 1, characterized in that, The empirical observation error covariance is scaled based on the residual statistics obtained from data processing, and the observation weights are iteratively updated using the variance component estimation (VCE) framework. For observation type i, the initial prior weighting matrix... Assuming it is a diagonal matrix ， It is given by the following formula: （2）， in, Representation type empirical standard deviation This represents the total number of observation types. For kinematic orbit data, the daily covariance matrix of each satellite is considered as an independent observation type. Within a one-month arc, 60 kinematic datasets for the binary satellites (GFO-C and GFO-D) are considered, with their initial weights derived from the inverses of their respective covariance matrices. For KRR observations, epoch-independent observations are assumed, and an identity matrix is ​​assigned as the initial weighting matrix (i.e., unit weights). The weighting factor after one iteration is defined as: (3), where represents the weighting factor. It is an observation-subtraction computation (OC) vector. ， Indicates the parameter to be estimated. It is the corresponding design matrix ， It represents the number of observations for each observation type. It is the total number of types, the observation weight matrix. The weighting factors are updated as follows: (4), of which, It requires multiple iterations to converge.

Images