Method for enhancing orbit determination stability based on system bias detection and receiver clock error constraint

Abstract

The present application relates to the cross technical field of satellite navigation, satellite gravitation, and specifically to a method for enhancing orbit determination stability based on system bias detection and receiver clock bias constraint, a two-stage quality control method is constructed by combining arc segment system bias detection and clock jump adaptive receiver clock bias constraint, the arc segment level processing uses the consistency test of the posterior residual to iteratively identify and eliminate the carrier phase arc segment containing system bias; the epoch level processing detects clock jump according to the clock bias increment of adjacent epochs, and reduces the weight of the cross clock jump epoch constraint; the present application can weaken the influence of system bias and clock jump on orbit determination from the arc segment scale and the epoch scale respectively, improve the ambiguity fixing effect, suppress the abnormal jump of the orbit, thereby improving the stability of the kinematic precise orbit determination result and the consistency with the reference orbit.

Description

Technical Field

[0001] This invention relates to the interdisciplinary fields of satellite navigation and satellite gravimetry, specifically a method for enhancing orbit determination stability based on system deviation detection and receiver clock bias constraints. Background Technology

[0002] High-precision orbits are fundamental to geodesy and Earth observation data processing for low Earth orbit (LEO) satellites. Kinematic precise orbit determination (KPOD) uses observations from the satellite-borne Global Positioning System (GPS) to estimate the satellite's position epoch by epoch, independent of satellite force models. Therefore, it can generate independent orbit products and be used for subsequent gravity field calculations. Because each epoch is calculated independently, the kinematic orbit is sensitive to satellite geometry, observation system errors, and the correlation between receiver clock errors and radial position. Yunck et al. (1990) were among the first to verify the feasibility of using GPS for precise tracking of LEO satellites, and subsequently, related methods were applied to real-time navigation, kinematic orbit determination, and simplified dynamic orbit determination. For the Gravity Recovery and Climate Experiment (GRACE) and its follow-up mission GRACE Follow-On (GRACE-FO), orbital errors also affect the processing of inter-satellite ranging observations and gravity field inversion.

[0003] Fixed carrier phase integer ambiguity can fully leverage the high precision of carrier phase observations and is an effective way to improve the quality of KPOD solution. For LEO satellite orbit determination, existing studies have proposed ambiguity fixing strategies such as non-differential, inter-satellite single-differential, and inter-track differential. Related comparative results show that the effectiveness of ambiguity fixing depends not only on the processing model and satellite configuration but also on the estimation quality of floating-point ambiguity. When observational anomalies or model residuals exhibit persistent deviations or slow drifts within a carrier phase arc, both floating-point ambiguity and its reference selection may be affected. The GPS tracking arc of LEO satellites is relatively short, and ambiguity parameters are closely related to the epoch-by-epoch receiver clock bias; arc-segment system biases may further propagate to clock bias and orbital parameters. Although current processing methods consider factors such as carrier phase entanglement, antenna phase center changes, satellite attitude, ionospheric effects, and cycle slip detection, unmodeled error characteristics are still visible in the GRACE-FO observation residuals. Therefore, it is necessary to perform residual consistency checks on the carrier phase arc before proceeding to integer ambiguity fixing.

[0004] Constraining receiver clock bias variations between adjacent epochs can also enhance the stability of KPOD solutions. The ultra-stable oscillator (USO) onboard GRACE exhibits good frequency stability over short periods. Previous studies have shown that constraining receiver clock bias increments can improve GRACE kinematic orbit determination results, particularly for periods with poor GPS satellite geometry. This method has also been extended to clock bias modeling and joint clock bias-orbit estimation for other LEO satellites. However, when a receiver clock jump occurs, the assumption of continuous clock bias between adjacent epochs no longer holds. If strong constraints are still applied across clock jump epochs, discontinuities not absorbed by the clock bias parameters will be transferred to the orbital parameters, primarily affecting orbital components strongly correlated with the receiver clock bias. Summary of the Invention

[0005] This invention addresses the issue of how to reduce the impact of arc segment system bias and receiver clock skipping on the ambiguity fixed KPOD by proposing a method to enhance orbit determination stability based on system bias detection and receiver clock bias constraint.

[0006] This invention achieves its purpose through the following measures: A method for enhancing orbit determination stability based on system bias detection and receiver clock bias constraints is characterized by constructing a two-level quality control method that integrates arc-segment system bias detection and clock jump adaptive receiver clock bias constraints. The arc-segment level processing utilizes post-hoc residual consistency checks to iteratively identify and eliminate carrier phase arc segments containing system biases. The epoch-level processing detects clock jumps based on adjacent epoch clock bias increments and reduces the weight of cross-clock jump epoch constraints. Specifically, it includes the following steps: Step 1: Arc segment system bias screening: Before the integer ambiguity is fixed, based on the post-hoc residual consistency test, the carrier phase arc segments containing system bias are identified and iteratively eliminated segment by segment to improve the estimation quality of floating-point ambiguity and ambiguity reference. Step 2: Ambiguity Fixation and Reliability Verification: Construct and fix the wide-lane and narrow-lane integer ambiguities of the observation data after arc segment screening, and perform carrier phase distance consistency verification to remove satellite observations that do not meet the consistency requirements; Step 3: Epoch-level quality control: The receiver clock error increments of adjacent epochs are checked to detect clock jumps. After a clock jump is detected, the weight of the clock error constraint between clock jump epochs is adaptively reduced to avoid the transmission of clock error discontinuity to the orbital parameters and improve the stability of the kinematic precision orbit determination results.

[0007] In step 1 of this invention, it is set and These are the post-hoc residual vectors and their weight matrices for all observations. and Let the number of observations and the number of parameters to be estimated be respectively. Then the estimated unit weighted variance after verification is: , For ambiguity arc segment , and These are the residual sub-vectors and weight matrix sub-blocks corresponding to the arc segment, respectively, with the null hypothesis... Assuming the residual of this arc segment follows a globally stochastic model, the alternative hypothesis is... The test statistic is defined as follows: If the weighted sum of squared residuals for this arc segment is significantly larger than expected, then the test statistic is defined as follows: , Where, in the formula, Describing the degrees of freedom as A chi-squared distributed random variable, For arc segment The corresponding degrees of freedom, at the significance level Under the following conditions: , Refusal ,in The upper quantile of the chi-square distribution is given. Since the common parameters and unit weight variances are estimated jointly by all observations, the arc segment residuals are not strictly independent. Therefore, equation (8) serves as an approximate statistical criterion for arc segment selection. The selection is carried out iteratively, and the test statistics of all remaining arc segments are recalculated in each round. If all arc segments are accepted... If the condition is met, the iteration ends; otherwise, the over-limit ratio is calculated using the following formula: Arcade segments that fail the test are sorted and then removed. The largest arc segment is identified, and then the floating-point solution, unit weight variance, and test statistics for the remaining arc segments are re-estimated. Only one arc segment is removed in each round to avoid misjudging other normal arc segments after strong anomalous arc segments disturb the clock bias of the common receiver. The iteration terminates when all remaining arc segments pass the test, or when continuing to remove them would make the adjustment model unsolvable.

[0008] After arc segment selection in step 2 of this invention, wide-lane floating-point ambiguity is first constructed, and its fractional part is used to jointly estimate the fractional-cycle bias (FCB) between the receiver and satellite ends. After subtracting the FCB, candidate values ​​for wide-lane integer-cycle ambiguity are determined according to the nearest integer criterion. , In equation (10), The wide-lane integer ambiguity candidate values ​​for GPS satellites are given in weeks. The corresponding wide-lane fractional-cycle bias correction is given, in weeks. A reference satellite is then selected to establish an ambiguity baseline, with all other ambiguities expressed relative to this baseline. After fixing the wide-lane ambiguity, the short-time stability of the receiver's FCB is used to connect the narrow-lane ambiguities of adjacent epochs. The baseline narrow alley ambiguity candidate value is In equation (11), For the corresponding The narrow-lane integer ambiguity reference of GPS satellites in the arc segment, in weeks; For the corresponding The fractional-cycle deviation correction for the width of the arc segment, in cycles, is defined as follows: the fractional part of a single floating-point ambiguity is defined as... The fractional parts of the ambiguity of multiple satellites are used together to estimate the FCB at the receiver and satellite ends. Equations (10) and (11) are used to summarize the conventional ambiguity fixing process and to explain the impact of abnormal floating-point ambiguity on FCB estimation and ambiguity reference selection. After fixing the ambiguity, a carrier phase distance consistency check is further performed. The fixed integer ambiguity is substituted back into the carrier phase observation to obtain the ambiguity-free carrier phase distance. In this check, the prior simplified dynamic trajectory is fixed, and only the receiver clock error and carrier phase residual are re-estimated. The test statistic is: In the formula, and These are the carrier phase distance residual vector and its weight matrix, respectively. For reference variance, For the corresponding degrees of freedom, at the significance level Under the following conditions: , Accept the ambiguity fixation result for this epoch. If the test fails, remove the satellite that contributes the most to the weighted sum of squared residuals and re-perform the ambiguity fixation and consistency test. Only one satellite is removed in each round. Stop the iteration after the test passes. If there are fewer than 4 satellites remaining, do not perform ambiguity fixation for this epoch.

[0009] In step 3 of this invention, adjacent epochs are... and The receiver clock bias increment is used to construct a pseudo-observation with a mean of zero: , In the formula, For the corresponding arc segment The receiver clock error of an epoch, measured in seconds; For the corresponding arc segment The receiver clock error of an epoch, measured in seconds; Let be the random clock bias increment expressed in terms of distance; under the condition of continuous clock bias, The expected value is 0, and its variance varies with the epoch interval. change; The variance of the nominal clock bias increment is calculated using the following formula according to the original procedure: , , In the formula, Pick empirical coefficient Take 3, based on the estimated clock bias increment. Perform clock jump detection when: , At that time, the epoch was determined. When a clock jump occurs, the detection threshold is directly taken as the nominal standard deviation given by equation (15), without setting a sliding time window or local scale estimation. After detecting a clock jump, only the standard deviation is reduced. and Weights of clock bias constraints between: , The nominal constraint variance is used for normal epochs; the larger constraint variance is used for clock jump epochs, so that the receiver clock error parameters can absorb discontinuities and avoid clock jump errors from being transmitted to the orbit coordinates. The original GNSS observations of the clock jump epochs still participate in the solution, and the nominal clock error constraints are used for the remaining epochs.

[0010] This invention proposes a two-stage quality control method that integrates arc segment system deviation detection and clock jump adaptive receiver clock bias constraint. Before ambiguity fixation, based on post-hoc residual consistency test, carrier phase arc segments containing system deviations are identified and iteratively eliminated segment by segment to improve floating-point ambiguity and ambiguity reference estimation. Simultaneously, a 3σ test is performed on the clock bias increment of adjacent epochs. After a clock jump is detected, the weight of the cross-clock jump epoch constraint is adaptively reduced to avoid the transmission of clock bias discontinuity error to orbital parameters. Experiments show that the arc segment system deviation screening of this invention increases the ambiguity fixation rate from 85.6% to 89.6%, a relative improvement of approximately 4.7%. At the representative clock jump epoch of SOD 53820 s on August 3, 2019, conventional clock bias constraints produced a radial orbital anomaly of approximately 0.25 m. This anomaly was effectively suppressed after adopting clock jump adaptive constraints. After comprehensively employing arc-segment systematic bias screening, ambiguity fixing, and clock jump adaptive constraints, the root mean square of the radial orbital differences of the four satellites decreased by 34.2%, 40.5%, 41.7%, and 63.0% respectively compared to the original floating-point solution. The results show that the present invention can reduce the influence of systematic bias and clock jump on orbit determination solutions at both the arc-segment and epoch scales, improve the ambiguity fixing effect, suppress orbital anomaly jumps, and thus improve the stability of the kinematic precision orbit determination results and their consistency with the reference orbit. Attached Figure Description

[0011] Figure 1 This is a flowchart of the present invention.

[0012] Figure 2 It is the IF carrier phase prefit residual of GRACE-A within the selected time period on November 1, 2009, in this embodiment of the invention.

[0013] Figure 3 This refers to the carrier phase aberration segment identified by GRACE-A in this embodiment of the invention on November 1, 2009, through iterative verification.

[0014] Figure 4 This refers to the narrow-lane floating-point ambiguity during a certain period on November 1, 2009, in an embodiment of the present invention, which was not filtered by the arc-segment system deviation.

[0015] Figure 5 This refers to the narrow-lane floating-point ambiguity after arc-segment system deviation filtering for a portion of the time period on November 1, 2009, in this embodiment of the invention.

[0016] Figure 6 This refers to the epoch-by-epoch ambiguity fixation rate of GRACE-A before and after the arc segment system deviation screening during a certain period on November 1, 2009, in an embodiment of the present invention.

[0017] Figure 7This refers to the receiver clock error of GRACE-A before and after the arc segment system deviation screening during a certain period on November 1, 2009, in an embodiment of the present invention.

[0018] Figure 8 In this embodiment of the invention, GRACE-A has radial, along-track, and normal orbital differences relative to the simplified dynamic orbit of JPL before and after the arc segment system deviation screening on November 1, 2009.

[0019] Figure 9 In this embodiment of the invention, the radial and along-track differences of GRACE-C relative to the simplified dynamics orbit of JPL before and after the application of conventional receiver clock bias constraints on August 3, 2019 are described.

[0020] Figure 10 In this embodiment of the invention, GRACE-C's radial, along-track, and normal orbital differences relative to the JPL simplified dynamic orbit near the SOD 53820 s clock jump on August 3, 2019 are described.

[0021] Figure 11 In this embodiment of the invention, GRACE-C has the IF combined carrier phase prefit residual, GDOP, and receiver clock error near the SOD 53820 s clock jump on August 3, 2019. Figure 11 Image (a) shows a schematic diagram of the residuals from the prefitting of the IF combined carrier phase. Figure 11 (b) is a schematic diagram of GDOP. Figure 11 (c) is the receiver clock bias diagram.

[0022] Figure 12 This refers to the receiver clock jump distributions of GRACE-A / B in November 2009 and GRACE-C / D in August 2019 in this embodiment of the invention.

[0023] Figure 13 The radial, along-track, and normal orbital differences of GRACE-C in this embodiment of the invention when using conventional and clock-jump adaptive constraints near the SOD 53820 s clock jump on August 3, 2019.

[0024] Figure 14 It is the orbital difference between the original floating-point solution and the final fixed ambiguity solution of GRACE-A on November 1, 2009, in the embodiment of the present invention.

[0025] Figure 15 It is the orbital difference between the original floating-point solution and the final fixed ambiguity solution of GRACE-C on August 3, 2019, in the embodiment of the present invention. Detailed Implementation

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

[0027] This invention constructs a two-stage quality control method that integrates arc-segment system bias detection and clock-jump adaptive receiver clock bias constraint. Arc-segment level processing utilizes post-hoc residual consistency checks to iteratively identify and eliminate carrier phase arcs containing system biases. Epoch-level processing detects clock jumps based on adjacent epoch clock bias increments and reduces the weight of cross-clock-jump epoch constraints. Experiments were conducted using one month each of GRACE-A / B and GRACE-C / D satellite-borne GPS observations to analyze the impact of arc-segment selection on ambiguity fixation and orbit continuity, the role of receiver clock bias constraints in normal and clock-jump epochs, and the improvement of the complete processing scheme on the orbit determination results of four satellites.

[0028] In this invention, the spaceborne receiver GPS satellites At frequency The equations for the non-differenced pseudorange and carrier phase observations are Equations (1) and (2), respectively. Where, This represents the geometric distance from the satellite to the receiver; It is the speed of light in a vacuum. and These are receiver clock bias and GPS satellite clock bias, respectively. Reference frequency The first-order ionosphere delay, and according to Mapping to frequency ; and These are the carrier wavelength and the integer ambiguity in cycles, respectively. This indicates the hardware delay at the receiver or satellite end, measured in meters (m). and These are the pseudorange hardware delay deviations for low-Earth orbit satellites and navigation satellites, respectively. and These are the carrier phase hardware delay deviations for low-Earth orbit satellites and navigation satellites, respectively. Includes observation noise and residual model error. and These are pseudorange and carrier observation errors, respectively.

[0029] Since the signal propagation path from the GPS satellite to the LEO receiver does not pass through the troposphere, the observation equation does not contain a tropospheric delay term. The data processing also takes into account conventional model corrections such as GPS satellite and LEO receiver antenna phase center correction, attitude-related phase center change, carrier phase entanglement, and relativistic effects, which are not listed item by item in equations (1) and (2).

[0030] , , After eliminating first-order ionospheric delay using ionosphere-free (IF) linear combination and correcting GPS satellite clock errors, the IF combination pseudorange and carrier phase observation equations can be written as follows: , , In equations (3) and (4), This is pseudo-distance clock error; and These represent the noise levels of combined observations without an ionosphere.

[0031] According to the notation conventions of equations (1) and (2), the IF combined ambiguity can be expressed as: , In equation (5), The carrier hardware delay after low-Earth orbit satellites are combined without an ionosphere; The pseudorange hardware delay after combining low-Earth orbit satellites without an ionosphere; The carrier hardware delay after the navigation satellites are combined without an ionosphere; The pseudorange hardware delay after the combination of navigation satellites without an ionosphere.

[0032] The IF combined ambiguity includes a combination of dual-frequency integer ambiguity and pseudorange and carrier phase hardware delay, and therefore does not possess integer characteristics. In practical processing, wide-lane and narrow-lane combinations are constructed separately to recover and fix the corresponding integer ambiguity. This invention follows the ambiguity fixing methods of Li et al. (2016), Zhang et al. (2021), and Gao et al. (2023), and focuses on arc segment system deviation verification and clock jump adaptive clock error constraint.

[0033] This invention employs a two-level quality control system: arc-level and epoch-level. Arc-level processing is performed before the integer ambiguity is fixed; the process is as follows. Figure 1 First, it is verified whether the weighted post-verification residuals of each carrier phase arc segment are consistent with the stochastic model obtained from the global adjustment. Arc segments that fail the verification are removed one by one, and the floating-point solution is recalculated after each removal until all remaining arc segments pass the verification. This step mainly identifies systematic biases or slow drifts that persist within the complete arc segment, so as to reduce their impact on floating-point ambiguity estimation and ambiguity reference selection.

[0034] Epoch-level processing addresses the clock bias increments of receivers at adjacent epochs. When clock bias changes satisfy the continuous model, inter-epoch clock bias constraints are applied; when a clock jump is detected, only the weight of the cross-clock jump epoch constraint is reduced, while nominal constraints remain in effect for other epochs. These two levels of processing correspond to arc-segment systematic bias and inter-epoch clock bias discontinuities, respectively: the former cannot identify true clock jumps, and the latter cannot eliminate persistent systematic biases within an arc. The combination of these two processes ensures both the quality of floating-point ambiguity estimation and the applicability of the receiver clock bias continuous model.

[0035] Arc segment system deviation screening: Let and These are the post-hoc residual vectors and their weight matrices for all observations. and Let the number of observations and the number of parameters to be estimated be respectively. Then the post-hoc unit weighted variance estimate is: , For ambiguity arc segment , and These are the residual sub-vectors and weight matrix sub-blocks corresponding to the arc segment, respectively. Null hypothesis. Assuming the residual of this arc segment follows a globally stochastic model, the alternative hypothesis is... The weighted sum of squared residuals for this arc segment is significantly larger than expected. The test statistic is defined as follows: , In the formula, Describing the degrees of freedom as A chi-squared distributed random variable, For arc segment The corresponding degrees of freedom, at the significance level Next, when satisfied , sometimes refuse ,in is the upper quantile of the chi-square distribution. Since the common parameters and unit weight variances are estimated jointly by all observations, the arc segment residuals are not strictly independent. Therefore, equation (8) serves as an approximate statistical criterion for arc segment selection.

[0036] The screening is conducted iteratively. In each round, the test statistic for all remaining arc segments is recalculated; if all arc segments accept... If the condition is met, the iteration ends. Otherwise, calculate the over-limit ratio using the following formula. Arcade segments that fail the test are sorted and then removed. The largest arc segment is selected, and then the floating-point solution, unit-weighted variance, and test statistics for the remaining arc segments are re-estimated. Only one arc segment is removed in each round to avoid misclassification of other normal arc segments due to strong anomalous arc segments disturbing the common receiver clock bias. The iteration terminates when all remaining arc segments pass the test, or when further removal would render the adjustment model unsolvable.

[0037] The above tests are only used to identify arc segments with statistically abnormal characteristics, not to determine the physical source of the bias. Inconsistent hardware delays, antenna model residuals, missed detections or incorrectly repaired cycle slips, and other observational anomalies may all lead to the exclusion of arc segments. Passing the test does not mean that the arc segment is completely free of systematic bias; it only indicates that its post-test residuals are consistent with the used stochastic model at a given significance level.

[0038] Ambiguity Fixation and Reliability Verification: After arc segment selection, wide-lane floating-point ambiguities are first constructed, and their fractional parts are used to jointly estimate the fractional-cycle bias (FCB) at the receiver and satellite ends. After subtracting the FCB, candidate values ​​for wide-lane integer ambiguities are determined according to the nearest integer criterion. In equation (10), The wide-lane integer ambiguity candidate values ​​for GPS satellites are given in weeks. This is the corresponding decimal-cycle deviation correction for the wide lane, in cycles. Subsequently, a reference satellite is selected to establish an ambiguity baseline, and all other ambiguities are represented relative to this baseline. After the wide-lane ambiguities are fixed, the short-time stability of the receiver's FCB is used to connect the narrow-lane ambiguities of adjacent epochs. At epoch... The baseline narrow alley ambiguity candidate value is , In equation (11), For the corresponding The narrow-lane integer ambiguity reference of GPS satellites in the arc segment, in weeks; For the corresponding The fractional-cycle deviation correction for the width of the arc segment, in cycles, is defined as follows: the fractional part of a single floating-point ambiguity is defined as... The fractional parts of the ambiguities from multiple satellites are used together to estimate the FCB at the receiver and satellite ends. Equations (10) and (11) are used to summarize the conventional ambiguity fixing process adopted in this invention and to explain the impact of abnormal floating-point ambiguity on FCB estimation and ambiguity benchmark selection; this invention does not involve new integer search algorithms.

[0039] After fixing the ambiguity, a carrier phase distance consistency check is performed. The fixed integer ambiguity is substituted back into the carrier phase observation to obtain the ambiguity-free carrier phase distance. In this check, the prior simplified dynamic trajectory is fixed, and only the receiver clock error and carrier phase residual are re-estimated. Epoch The test statistic is In the formula, and These are the carrier phase distance residual vector and its weight matrix, respectively. For reference variance, For the corresponding degrees of freedom. At the significance level. Next, when satisfied , The ambiguity fixing result for that epoch is accepted. If the test fails, the satellite that contributes the most to the weighted sum of squared residuals is removed, and the ambiguity fixing and consistency test is repeated, removing only one satellite per round. Iteration stops after the test passes; if fewer than four satellites remain, ambiguity fixing is not performed for that epoch. This test reflects the degree of consistency between the fixed ambiguity and the simplified dynamic reference orbit and observation model, and cannot be used as an independent assessment of orbit accuracy.

[0040] Receiver clock bias model and clock jump adaptive constraints: For adjacent epochs and The receiver clock bias increment is used to construct a pseudo-observation with a mean of zero: , In the formula, For the corresponding arc segment The receiver clock error of an epoch, measured in seconds; For the corresponding arc segment The receiver clock error of an epoch, measured in seconds; Let be the random clock bias increment expressed in terms of distance; under the condition of continuous clock bias, The expected value is 0, and its variance varies with the epoch interval. When GPS satellite geometry is poor, this constraint can enhance the stability of the epoch-by-epoch coordinate solution and reduce the parametric correlation between receiver clock bias and radial position. Similar methods have been used for GRACE KPOD and other LEO satellite clock bias estimation.

[0041] The variance of the nominal clock bias increment is calculated using the following formula according to the original procedure: , In the formula, Pick empirical coefficient Take 3. Based on the estimated clock bias increment. Perform clock jump detection when: , At that time, the epoch was determined. A clock jump occurs. The detection threshold is directly the nominal standard deviation given by equation (15), without setting a sliding time window or local scale estimation. After detecting a clock jump, only the value is reduced. and Weights of clock bias constraints between: , Nominal constraint variance is used for normal epochs; a larger constraint variance is used for clock-jump epochs, allowing the receiver clock error parameters to absorb discontinuities and preventing clock-jump errors from being propagated to the orbit coordinates. The original GNSS observations for clock-jump epochs are still used in the solution, while nominal clock error constraints continue to be used for the remaining epochs.

[0042] Experimental verification: The experimental data include satellite-borne GPS observations from GRACE-A / B in November 2009 and GRACE-C / D in August 2019. Arc-segment system bias analysis was performed on November 1, 2009 (GRACE-A), and receiver clock bias constraint analysis was performed on August 3, 2019 (GRACE-C).

[0043] To facilitate result comparison, this invention defines the following processing schemes: "Original floating-point solution" is the KPOD result without arc segment filtering and ambiguity fixing; "Filtered floating-point solution" is the floating-point solution recalculated after removing abnormal arc segments; "Solution without clock bias constraints" is the epoch-by-epoch solution without applying inter-epoch receiver clock bias constraints; "Conventional clock bias constraint solution" is the result with nominal clock bias constraints applied to all epochs; "Clock jump adaptive constraint solution" is the result after reducing the constraint weight at clock jump epochs; "Final fixed solution" is the solution that comprehensively applies arc segment filtering, ambiguity fixing, and clock jump adaptive clock bias constraints. Not all of the above schemes are combined; comparisons between groups are only used to evaluate the corresponding processing steps.

[0044] The orbit determination results were evaluated by comparing them with the simplified dynamic orbit of JPL GNV, and the radial, along-track, and normal orbit differences and their root-mean-square (RMS) were statistically analyzed. The ambiguity fixation rate was defined as the ratio of the number of successfully fixed ambiguities to the total number of ambiguities involved in the fixation; the receiver clock jump amplitude was used... Converted to distance, unit is meters. Since the JPL reference orbit uses relevant tracking data and is not completely independent, the orbit difference in this example is only used to measure the degree of consistency between different schemes and the reference orbit, and is not used as absolute orbit determination accuracy.

[0045] Figure 2 — Figure 8The GRACE-A arc segment systematic bias screening method and its impact on ambiguity fixation and trajectory results are presented. Figure 9 — Figure 13 Analyze the performance of receiver clock bias constraints in normal epochs and clock jump epochs, and statistically analyze clock jumps throughout the month; Figure 14 and Figure 15 Compare the original floating-point solution with the final fixed solution. Since the final fixed solution includes three processes, its difference from the original floating-point solution reflects the combined effect of the complete processing flow and cannot be decomposed into the independent contribution of individual measures.

[0046] Results analysis: Figure 2 The prefit residuals of the IF combined carrier phase for a portion of GRACE-A on November 1, 2009, are presented. Most synchronous observation arcs show similar trends, mainly reflecting common variations in receiver clock bias; however, arc G08 deviates significantly from this trend, with a residual variation of approximately 0.5 m within a single arc. This indicates that arc G08 is inconsistent with other observations from the same period, but the physical cause cannot yet be determined from this figure. The arc selection results after iterative convergence are shown in [figure missing]. Figure 3 The final abnormal arc segments identified included G08, G27, and other arc segments whose weighted residual statistics exceeded the critical value. Because parameter estimation and arc segment testing were performed again after each removal, Figure 3 The results shown are the complete iterative results, not the set of anomalous arcs obtained from a single test. These results indicate that the residual statistical characteristics of some carrier phase arcs do not conform to the global stochastic model. Current experiments can only identify statistical anomalies and cannot distinguish specific causes such as hardware delay variations, antenna model residuals, cycle slip processing errors, or other observational anomalies.

[0047] Impact of arc segment filtering on fixed solution trajectories: Without arc segment filtering, the floating-point ambiguity of some narrow lanes deviates significantly from the cluster centers near integers. Figure 4 After removing outlier arcs and recalculating the floating-point solution, the remaining narrow alley ambiguity is closer to an integer. Figure 5 Because the screening process directly eliminated complete observation arc segments, Figure 4 and Figure 5 The sets of ambiguities in the data are not entirely identical. This result indicates that arc filtering can improve the estimation quality of floating-point ambiguities without changing the subsequent integer ambiguity search method.

[0048] During this period, the ambiguity fixation rate increased from 85.6% before screening to 89.6% after screening, a relative increase of approximately 4.7%. Figure 6 ). Figure 7The receiver clock bias before and after arc segment selection is further presented. The unselected solution experiences data interruptions during periods with insufficient fixed ambiguities; at approximately 81,000 seconds of the second of day (SOD), ambiguity reference switching and anomalous jumps in the clock bias sequence occur simultaneously. After arc segment selection, the continuity of the clock bias sequence improves and is more consistent with the clock bias changes derived from the simplified dynamic orbit determination solution. This indicates that arc residuals may abnormally affect receiver clock bias estimation through ambiguity reference selection, but the existing results are insufficient to determine a unique error propagation path.

[0049] Figure 8 The radial, along-track, and normal orbital differences were compared before and after arc segment screening. Screening increased the number of effective orbital epochs, and some obvious anomalous segments were suppressed. Since this example does not define a quantitative index for orbital continuity, only a qualitative comparison is made, and the percentage improvement in continuity is not given. For this GRACE-A example, the orbital difference sequence after screening is more continuous.

[0050] Receiver clock bias constraint: For GRACE-C data from August 3, 2019, after applying the conventional inter-epoch clock bias constraint, the radial orbital error RMS decreased from 1.6 cm to 1.4 cm, a reduction of 12.5%; the along-orbit and normal RMS remained at 0.7 cm and 0.6 cm, respectively. Figure 9 As can be seen, clock bias constraints primarily reduce the dispersion of radial orbital error, which is consistent with the strong parametric correlation between receiver clock bias and radial position. The above results only correspond to single-day GRACE-C data and do not represent the same level of improvement that can be achieved for the entire month or the entire mission period.

[0051] When a receiver clock jump occurs, conventional clock bias constraints may lead to orbital anomalies. At SOD 53820 s, the conventional constraint solution exhibits a radial orbital anomaly of approximately 0.25 m and a normal orbital anomaly of approximately 0.04 m, while the along-track component shows relatively small variations. Figure 10 The prefit residuals of the IF combined carrier phase at the same epoch exhibit a common jump of approximately 0.17 m. Figure 11 In (a) of this case, the geometric dilution of precision (GDOP) is approximately 3.5. Figure 11 (b)). The receiver clock errors obtained from the simplified dynamic orbit determination solution and the clock error-free solution both exhibit synchronous jumps, while the clock error of the conventionally constrained solution remains continuous. Figure 11 (c) The synchronous changes in the observed residuals and clock bias sequences under different schemes, along with the satellite geometry which did not significantly deteriorate at this epoch, support the view that the anomaly was caused by receiver clock skipping. Conventional constraints force clock bias continuity, causing some minute skipping errors to propagate to the orbital parameters. Figure 10The 0.25 m in the figure represents the radial orbital anomalous amplitude of that epoch, not the clock jump amplitude.

[0052] When the clock jump detection method was applied to the clock difference sequence of a full-month receiver, GRACE-A, GRACE-B, GRACE-C, and GRACE-D detected 5, 43, 93, and 67 clock jumps, respectively. Figure 12 The number of times the clock jump amplitude equivalent distance exceeded 0.1 m were 0, 9, 55, and 28, respectively. The largest clock jump amplitude equivalent distance was 0.70 m, occurring during the GRACE-FO data period. This value indicates... The receiver clock jump variable is not an orbital difference caused by clock jumps. Monthly statistics show that clock jumps are not isolated phenomena, but since each mission only analyzes one month, and the detection results are affected by the clock jump stochastic model and... The threshold effect prevents a direct comparison of clock jump rates between GRACE and GRACE-FO throughout the entire mission.

[0053] Clock jump adaptive constraint: For the clock jump at SOD 53820 s, by only reducing the weight of the cross-clock jump epoch constraint, the radial and normal orbital differences in the conventional clock error constraint solution are often significantly suppressed, while the adjacent normal epochs and the along-track components remain basically unchanged. Figure 13 Therefore, it can be seen that the clock jump adaptive constraint only adjusts the clock difference relationship between abnormal epochs, while the clock difference constraint between epochs is still retained in normal epochs.

[0054] To evaluate the overall effect of the complete processing flow, the original floating-point solution and the final fixed solution were further compared. The trajectory differences in all three directions of the GRACE-A final fixed solution converged significantly, and outliers were significantly reduced. Figure 14 The changes in GRACE-B are similar to those in GRACE-A, therefore no further plotting is needed. See the results for GRACE-C. Figure 15 GRACE-D and GRACE-C showed similar performance, and their time series plots are also omitted. The radial orbital difference RMS of GRACE-A and GRACE-B decreased from 3.8 cm to 2.5 cm and from 3.7 cm to 2.2 cm, respectively, representing decreases of 34.2% and 40.5%; while those of GRACE-C and GRACE-D decreased from 2.4 cm to 1.4 cm and from 2.7 cm to 1.0 cm, respectively, representing decreases of 41.7% and 63.0%. The radial orbital difference RMS of all four satellites decreased to varying degrees.

Claims

1. A method for enhancing orbit determination stability based on system deviation detection and receiver clock bias constraints, characterized in that, A two-level quality control method integrating arc segment system deviation detection and clock jump adaptive receiver clock error constraint was constructed. The arc segment level processing utilizes post-hoc residual consistency test to iteratively identify and eliminate carrier phase arc segments containing system deviation. The epoch-level processing detects clock jumps based on the clock difference increment between adjacent epochs and reduces the weight of cross-clock jump epoch constraints. Specifically, the following steps are included: Step 1: Arc segment system bias screening: Before the integer ambiguity is fixed, based on the post-hoc residual consistency test, the carrier phase arc segments containing system bias are identified and iteratively eliminated segment by segment to improve the estimation quality of floating-point ambiguity and ambiguity reference. Step 2: Ambiguity Fixation and Reliability Verification: Construct and fix the wide-lane and narrow-lane integer ambiguities of the observation data after arc segment screening, and perform carrier phase distance consistency verification to remove satellite observations that do not meet the consistency requirements; Step 3: Epoch-level quality control: The receiver clock error increments of adjacent epochs are checked to detect clock jumps. After a clock jump is detected, the weight of the clock error constraint between clock jump epochs is adaptively reduced to avoid the transmission of clock error discontinuity to the orbital parameters and improve the stability of the kinematic precision orbit determination results.

2. The method for enhancing orbit determination stability based on system deviation detection and receiver clock bias constraint as described in claim 1, characterized in that, In step 1, let and These are the post-hoc residual vectors and their weight matrices for all observations. and Let the number of observations and the number of parameters to be estimated be respectively. Then the estimated unit weighted variance after verification is: ， For ambiguity arc segment , and These are the residual sub-vectors and weight matrix sub-blocks corresponding to the arc segment, respectively, with the null hypothesis... Assuming the residual of this arc segment follows a globally stochastic model, the alternative hypothesis is... The test statistic is defined as follows: If the weighted sum of squared residuals for this arc segment is significantly larger than expected, then the test statistic is defined as follows: ， Where, in the formula, Describing the degrees of freedom as A chi-squared distributed random variable, For arc segment The corresponding degrees of freedom, at the significance level Under the following conditions: ， Refusal ,in The upper quantile of the chi-square distribution is given. Since the common parameters and unit weight variances are estimated jointly by all observations, the arc segment residuals are not strictly independent. Therefore, equation (8) serves as an approximate statistical criterion for arc segment selection. The selection is carried out iteratively, and the test statistics of all remaining arc segments are recalculated in each round. If all arc segments are accepted... If the condition is met, the iteration ends; otherwise, the over-limit ratio is calculated using the following formula: ， Arcade segments that fail the test are sorted and then removed. The largest arc segment is identified, and then the floating-point solution, unit weight variance, and test statistics for the remaining arc segments are re-estimated. Only one arc segment is removed in each round to avoid misjudging other normal arc segments after strong anomalous arc segments disturb the clock bias of the common receiver. The iteration terminates when all remaining arc segments pass the test, or when continuing to remove them would make the adjustment model unsolvable.

3. The method for enhancing orbit determination stability based on system deviation detection and receiver clock bias constraint as described in claim 2, characterized in that, After the arc segment selection in step 2 is completed, the wide-lane floating-point ambiguity is first constructed, and the fractional-cycle bias (FCB) at the receiver and satellite ends is jointly estimated from its fractional part. After subtracting the FCB, the candidate values ​​of the wide-lane integer ambiguity are determined according to the nearest integer criterion. ， In equation (10), The wide-lane integer ambiguity candidate values ​​for GPS satellites are given in weeks. The corresponding wide-lane fractional-cycle bias correction is given, in weeks. A reference satellite is then selected to establish an ambiguity baseline, with all other ambiguities expressed relative to this baseline. After fixing the wide-lane ambiguity, the short-time stability of the receiver's FCB is used to connect the narrow-lane ambiguities of adjacent epochs. The baseline narrow alley ambiguity candidate value is ， In equation (11), For the corresponding The narrow-lane integer ambiguity reference of GPS satellites in the arc segment, in weeks; For the corresponding The fractional-cycle deviation correction for the width of the arc segment, in cycles, is defined as follows: the fractional part of a single floating-point ambiguity is defined as... The fractional parts of the ambiguity of multiple satellites are used together to estimate the FCB at the receiver and satellite ends. Equations (10) and (11) are used to summarize the conventional ambiguity fixing process and to explain the impact of abnormal floating-point ambiguity on FCB estimation and ambiguity reference selection. After fixing the ambiguity, a carrier phase distance consistency check is further performed. The fixed integer ambiguity is substituted back into the carrier phase observation to obtain the ambiguity-free carrier phase distance. In this check, the prior simplified dynamic trajectory is fixed, and only the receiver clock error and carrier phase residual are re-estimated. The test statistic is: ， In the formula, and These are the carrier phase distance residual vector and its weight matrix, respectively. For reference variance, For the corresponding degrees of freedom, at the significance level Under the following conditions: ， Accept the ambiguity fixation result for this epoch. If the test fails, remove the satellite that contributes the most to the weighted sum of squared residuals and re-perform the ambiguity fixation and consistency test. Only one satellite is removed in each round. Stop the iteration after the test passes. If there are fewer than 4 satellites remaining, the ambiguity will not be fixed for that epoch.

4. The method for enhancing orbit determination stability based on system deviation detection and receiver clock bias constraint as described in claim 3, characterized in that, In step 3, for adjacent epochs and The receiver clock bias increment is used to construct a pseudo-observation with a mean of zero: ， In the formula, For the corresponding arc segment The receiver clock error of an epoch, measured in seconds; For the corresponding arc segment The receiver clock error of an epoch, measured in seconds; Let be the random clock bias increment expressed in terms of distance; under the condition of continuous clock bias, The expected value is 0, and its variance varies with the epoch interval. change; The variance of the nominal clock bias increment is calculated using the following formula according to the original procedure: ， ， In the formula, Pick empirical coefficient Take 3, based on the estimated clock bias increment. Perform clock jump detection when: ， At that time, the epoch was determined. When a clock jump occurs, the detection threshold is directly taken as the nominal standard deviation given by equation (15), without setting a sliding time window or local scale estimation. After detecting a clock jump, only the standard deviation is reduced. and Weights of clock bias constraints between: ， Normal epochs use nominally constrained variance; The clock jump epoch uses a larger constraint variance so that the receiver clock error parameter can absorb discontinuities and prevent clock jump errors from being transmitted to the orbit coordinates. The original GNSS observations of the clock jump epoch still participate in the solution, while the nominal clock error constraint continues to be used for the remaining epochs.

Images