Method, device and medium for removing abnormal values of GNSS precise clock error product

Abstract

The application provides a GNSS precise clock difference product outlier elimination method, equipment and medium, and belongs to the technical field of satellite navigation. Each analysis center is paired two by two, and the clock difference data of each pair of analysis centers is used to generate a transient difference sequence by using a clock difference per epoch alignment method; a clock difference residual sequence is generated by using a linear transformation residual calculation method; a second difference sequence is generated by using a clock difference second difference method; outliers are eliminated based on the transient difference sequence, the clock difference residual sequence and the second difference sequence by using a weighted fusion algorithm, consistency verification is performed on the final clock difference product, and a precise clock difference product after the outliers are eliminated is output. The method combines linear transformation between analysis centers and the clock difference second difference method on the basis of the clock difference per epoch estimation method, comprehensively utilizes advantages and disadvantages of each method, and completes clock difference product outlier elimination through joint verification of multiple analysis centers.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of satellite navigation, and particularly relates to a GNSS precise clock difference product outlier elimination method, equipment and medium. BACKGROUND

[0002] Currently, there are mainly two methods for data preprocessing before clock difference product synthesis: clock difference quadratic difference method and analysis center inter-similarity transformation. Among them, the full-system satellite linear transformation method involves all satellite data of the constellation, and if the clock difference calculation of a certain satellite at a specific epoch is wrong, it will cause a large jump in the clock difference difference of all satellites of the system on that day, thereby affecting the clock difference outlier detection result. Although the per-satellite linear transformation method can avoid global error propagation, it cannot effectively detect these abnormal situations when facing weak fluctuations caused by signal interruption. In contrast, the clock difference quadratic difference method eliminates common errors through inter-satellite single difference, which can capture these small fluctuations, but its effectiveness is highly dependent on the integrity and stability of the clock difference data of the selected reference satellite. If the clock difference of the reference satellite at a certain epoch is missing or abnormal, it will cause a chain jump in the quadratic difference of all satellites of the system at that epoch, thereby affecting the accuracy of the overall analysis result. In view of the above limitations, we propose a clock difference per-epoch estimation method, which removes the diurnal bias of the clock difference of each satellite of the two analysis centers, and then removes the epoch bias of the system, and calculates the clock difference residual after removing the bias. Its sensitivity is significantly better than that of the traditional method. Due to its single-epoch isolated calculation characteristics, when a satellite has a large anomaly in a single epoch, it will affect the outlier detection of all satellites in this epoch.

[0003] The prior art with publication number CN119471737A discloses a GNSS precise satellite clock difference correction number calculation method based on state space representation, including the following steps: step one: select a data segment of 10 epochs as the fitting object, thereby effectively reducing the calculation amount, improving the speed of real-time solution, and ensuring the efficiency and accuracy of clock difference correction number fitting; step two: for the clock difference data of each epoch, deduct the average value of the clock difference of each system at the current epoch; step three: use the SSR data before and after detection to restore the satellite clock difference after a 5-second delay in broadcast time, and compare the clock difference data after detection of gross errors with the original clock difference data; step four: generate SSR information to restore real-time precise clock difference; step five: after saving the real-time data stream data obtained by BNC, perform pseudo-real-time clock difference estimation. The calculation steps are complex and the precision is low.

[0004] The prior art disclosed in CN108107455A describes a real-time satellite clock bias prediction method based on phase jumps. The method involves acquiring the original satellite clock bias data sequence, preprocessing it, performing polynomial fitting on the processed satellite clock bias phase sequence, removing periodic terms from the phase sequence, and calculating the fitting accuracy. If the fitting accuracy is less than a threshold, the satellite clock bias is predicted, and the prediction ends. If it is greater than the threshold, frequency anomalies are detected from the end to the beginning, obtaining several data points before and after the anomalies. If there are no frequency anomalies and the data is stable, the satellite clock bias is predicted, and the prediction ends. If there are frequency anomalies and / or data jumps, the predicted satellite clock bias is used, and the prediction ends. This method requires significant network resources, processes large amounts of data, and has a slow response time. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this paper provides a method, equipment, and medium for outlier removal in GNSS precision clock error products. By progressively identifying and removing outliers, clock error data among multiple analysis centers is optimized, resulting in more accurate and reliable final clock error products, and enabling precise satellite positioning and monitoring.

[0006] To address the aforementioned technical objectives, this invention discloses a multi-center jointly verified method for eliminating outliers in GNSS precision clock error products, comprising the following steps:

[0007] S1. Preprocess the GNSS clock bias historical data provided by each analysis center, and pair each analysis center with the grouped preprocessed clock bias historical data.

[0008] S2. Take three times the standard deviation of the original clock difference of each pair of analysis centers as the threshold a. By performing joint detection on the paired analysis centers, outliers in the clock difference data of the analysis centers that exceed the threshold a are eliminated.

[0009] S3. The clock difference epoch-by-epoch alignment method is used to align each pair of central clock differences in the analysis epoch-by-epoch to obtain the instantaneous differences of clock differences. Outliers are removed with a threshold of 3 times the standard deviation as b, and an instantaneous difference sequence is generated.

[0010] S4. Use the linear transformation method between analysis centers to obtain the linear transformation parameters between the clock errors of each satellite in each pair of analysis centers, calculate the clock error residuals after linear transformation, remove outliers with a threshold of 3 times the standard deviation c, and generate a clock error residual sequence.

[0011] S5. Using the clock difference quadratic difference method, calculate the quadratic difference for each pair of clock difference data from the analysis center, and remove outliers with a threshold of 3 times the standard deviation d; generate a quadratic difference sequence.

[0012] S6. Repeat S2-S5 until the change rate of the number of abnormal epochs removed in two consecutive iterations is less than 5%, or the preset maximum number of iterations is reached, and output the final clock difference product.

[0013] S7, consistency verification is performed on all final clock error products, and root mean square errors of clock errors between analysis centers are calculated;

[0014] S8, when the root mean square errors of clock errors between analysis centers all meet the requirements, output the precise clock error product after removing abnormal values.

[0015] Further, the process of preprocessing is: packing and numbering the GNSS clock error historical data provided by all analysis centers, and then traversing and combining each GNSS clock error historical data packet with all the other GNSS clock error historical data packets.

[0016] Further, the mathematical expression of GNSS clock error historical data is constructed: the original pseudo-range and carrier phase observation value included in the GNSS clock error historical data are expressed as:

[0017] ,

[0018] In the formula, s, T, r and j respectively represent satellite, satellite system, receiver and carrier frequency; and respectively represent pseudo-range and carrier phase observation; represents the coefficient matrix; represents the three-dimensional coordinate increment of the receiver; c represents the speed of light in vacuum; and respectively represent receiver information and satellite clock error information; and respectively represent wet mapping function and zenith wet delay; represents the slant ionospheric delay on the signal frequency; represents the frequency-dependent multiplier factor; and respectively represent the frequency-dependent receiver uncalibrated code delay and the frequency-dependent satellite uncalibrated code delay relative to satellite s; represents the carrier wavelength; N is the integer phase ambiguity; and respectively represent the frequency-dependent receiver uncalibrated phase delay and the frequency-dependent satellite uncalibrated phase delay relative to satellite s; and respectively represent the sum of measurement noise and multipath error of pseudo-range and carrier phase observation.

[0019] Further, the clock error per epoch estimation method is used to estimate each pair of analysis center clock error per epoch, and the difference between the clock error data of each pair of analysis center after pairing is calculated as follows:

[0020] In the clock difference data of the pairwise analysis centers, the clock difference data of one analysis center is defined as the clock difference data of the analysis center ac, and the clock difference data of the other analysis center is defined as the clock difference data of the reference center ref;

[0021] The original difference of the clock difference of the satellite s of the analysis center ac and the reference center ref at the k ephemeris is calculated by using the following formula:

[0022] ,

[0023] In the formula, Dac(k) represents the clock correction number of the satellite s of the analysis center ac at the k ephemeris, Dref(k) represents the clock correction number of the satellite s of the reference center ref at the k ephemeris; by subtracting the clock difference between the analysis center ac and the reference center ref, the clock difference change characteristics of the satellite clock itself are eliminated; The clock difference difference of the satellite s after removing the daily bias at the k ephemeris is calculated by using the following formula:

[0024]

[0025] , In the formula, Dac(k) represents the clock correction number of the satellite s of the analysis center ac at the k ephemeris,

[0026] Dref(k) represents the clock correction number of the satellite s of the reference center ref at the k ephemeris; by subtracting the clock difference between the analysis center ac and the reference center ref, the clock difference change characteristics of the satellite clock itself are eliminated; The clock difference difference of the satellite s at the k ephemeris is calculated by using the following formula:

[0027]

[0028] ,

[0029] In the formula, Dac(k) represents the clock correction number of the satellite s of the analysis center ac at the k ephemeris, Dref(k) represents the clock correction number of the satellite s of the reference center ref at the k ephemeris; by subtracting the clock difference between the analysis center ac and the reference center ref, the clock difference change characteristics of the satellite clock itself are eliminated;

[0030] Further, the linear transformation parameters between the clock differences of each satellite of each pair of analysis centers are obtained by using the linear transformation method between the analysis centers, and the clock difference residual after linear transformation is calculated:

[0031] The clock differences of the analysis center ac and the reference center ref are aligned by using the linear transformation, and a set of linear transformation parameters is estimated for each satellite system T, and the conversion formula of any satellite s of any satellite system T at the ephemeris k is represented as:

[0032] ,

[0033] In the formula, Dac(k) represents the clock correction number of the satellite s of the analysis center ac at the k ephemeris, Dref(k) represents the clock correction number of the satellite s of the reference center ref at the k ephemeris; by subtracting the clock difference between the analysis center ac and the reference center ref, the clock difference change characteristics of the satellite clock itself are eliminated; , , respectively represent linear transformation conversion parameters,​​​ GPS time representing the clock error at epoch k, is the reference time of ephemeris; let the satellite system T contain n satellites, each satellite has m epochs of data;

[0034] The linear transformation conversion parameters are solved by using the following formula 、 、 :

[0035] ,

[0036] ,

[0037] ,

[0038] ,

[0039] ,

[0040] The clock error residual between the analysis center ac linear transformation and the reference center ref is solved by using the following formula :

[0041] ,

[0042] When the clock error is integrated, a set of linear transformation parameters is estimated for each satellite clock error, eliminating the influence of the middle linear part: the transformation formula of the analysis center ac clock error relative to the reference center ref clock error is rewritten as

[0043] ,

[0044] Wherein is a set of linear transformation parameters calculated for each satellite of the system .

[0045] Further, the second difference method of clock error is used to calculate the second difference of each pair of analysis center clock error products:

[0046] Let the satellite clock clock error correction number of GNSS data processing is represented as:

[0047] ,

[0048] In the formula, indicates the reference clock used by the analysis center in data processing; indicates the deviation related to the analysis center and the satellite clock, indicates the determined part in the satellite clock itself characteristics; indicates the random part in the other observation noise and the clock itself characteristics;

[0049] The first difference of satellite clock error is obtained by eliminating the clock error variation characteristics of satellite clock itself by subtracting the analysis center ac clock correction number from the reference center ref using the following formula:

[0050] ,

[0051] In the formula, and are the clock correction numbers of the reference center and the analysis center respectively; represents the difference between the reference clock of the analysis center ac and the reference clock of the reference center ref; represents the difference between the reference clock of the analysis center ac and the reference clock of the reference center ref; represents the difference between the reference clock of the analysis center ac and the reference clock of the reference center ref; represents the difference between the reference clock of the analysis center ac and the reference clock of the reference center ref; represents the difference between the reference clock of the analysis center ac and the reference clock of the reference center ref; represents the difference between the reference clock of the analysis center ac and the reference clock of the reference center ref; represents the difference between the reference clock of the analysis center ac and the reference clock of the reference center ref; represents the difference between the reference clock of the analysis center ac and the reference clock of the reference center ref;

[0052] A reference satellite is selected for each satellite system, and the first difference of satellite clock error is subtracted between the reference satellite si and other satellites sj of the satellite system to obtain the second difference of satellite clock error, eliminate the reference clock related to the analysis center, and only leave the observation noise term and the term related to the data processing strategy of the analysis center, i and j are satellite numbers, using the following formula:

[0053] ,

[0054] In the formula, respectively represent the difference between the clock of satellite sj and the clock of satellite si; is the clock difference second difference result; represents the difference between the reference clock of the analysis center ac and the reference clock of the reference center ref; respectively represent the clock difference first difference between satellite sj and satellite si. Further, the epoch data to be eliminated needs to be marked in the data of each analysis center for more than twice; when the consistency of the final clock difference product is verified, the root mean square error of the clock difference between each analysis center is calculated, and when the root mean square error of the satellite clock difference of each system is less than 0.3 ns, the verification is passed, otherwise the S2 loop is returned.

[0055] A computer device comprising a processor and a memory, the processor and the memory being electrically connected, the memory being used to store instructions and data, and the processor being used to execute the GNSS precise clock difference product outlier elimination method.

[0056]

[0057] ​​A computer-readable storage medium storing a computer program adapted to be loaded by a processor and executed as a method for removing outliers from GNSS precision clock products.

[0058] Beneficial Effects: This invention effectively improves the accuracy and consistency of the final clock difference product by integrating three clock difference consistency analysis methods and combining clock difference data from multiple analysis centers to jointly identify and eliminate outliers in the clock difference products from each analysis center. Its implementation process covers the entire workflow from data preprocessing to outlier removal, consistency verification, and final product output. Through this series of optimized steps, the reliability and accuracy of the final clock difference data can be ensured in multiple fields such as precise positioning, satellite orbit determination, and crustal monitoring, demonstrating significant application value. Attached Figure Description

[0059] Figure 1 This is a flowchart illustrating the outlier removal method for GNSS precision clock error products according to the present invention. Detailed Implementation

[0060] The embodiments of the present invention will be further described below with reference to the accompanying drawings:

[0061] like Figure 1 As shown, this invention discloses a method for removing outliers in GNSS precision clock error products, the steps of which are as follows:

[0062] S1, for multiple analysis centers AC1, AC2 The GNSS clock bias historical data provided by ACn is preprocessed, grouped by satellite system (GPS, GLONASS, Galileo, BDS), and then matched with each other in pairs.

[0063] S2. For each pair of analysis centers and each satellite, calculate the standard deviation of the clock error sequence difference, using 3 times the standard deviation. The threshold 'a' is used as the basis for determining the clock difference. If the abnormal fluctuation of a clock difference in a certain epoch exceeds the threshold 'a', then that epoch is marked in the corresponding analysis center pair. When a single epoch is marked more than or equal to 2 times, the first outlier removal is performed. This step mainly eliminates outliers that significantly deviate from the population consistency.

[0064] S3. Using the epoch-wise clock bias estimation method, after removing the daily bias, epoch-wise alignment estimation is performed on each pair of analysis center clock biases to calculate the instantaneous difference in clock bias. Considering the characteristics of the satellite system, a value of 3 standard deviations is used. Generate a differential threshold b. Epochs exceeding threshold b are marked. If a single epoch is marked two or more times, a second rejection is performed. This step focuses on anomalies such as short-term jumps or high-frequency noise.

[0065] S4. Using the linear transformation method between analysis centers, calculate the linear transformation parameters between the clock errors of each satellite for each pair of analysis centers, and calculate the clock error residuals after the linear transformation. (Using 3 times the standard deviation...) A threshold 'c' is generated for the linearly transformed clock error residuals. For the transformed clock error residuals, a threshold of three times the standard deviation is recalculated according to the satellite system. Epochs exceeding threshold 'c' are marked. If an epoch in the analysis center is marked two or more times, a third rejection is performed. This step can suppress systematic anomalies caused by baseline differences or equipment drift.

[0066] S5. Using the clock difference quadratic difference method, calculate the quadratic difference for each pair of clock difference products from the analysis center. The difference is expressed as three times the standard deviation. A threshold 'd' is generated for the epoch-by-epoch clock error difference. When the fluctuation of the clock error residual exceeds this threshold 'd', the clock error epoch corresponding to the analysis center is marked. Epochs exceeding the threshold 'd' are marked, and a fourth removal is performed when a single epoch is marked two or more times. This step effectively suppresses the interference of single-star clock error estimation errors on the overall results.

[0067] S6. For the data after removing anomalies in the first four steps, recalculate the thresholds a, b, c, and d for clock anomalies in each analysis center, and iterate from S1 to S6 until one of the following conditions is met: 1) The change rate of the number of anomaly epochs removed in two consecutive iterations is <5%; 2) The preset maximum number of iterations (5 times) is reached.

[0068] S7. Perform consistency verification on the final clock error product, calculate the root mean square error (RMS) of clock errors between each analysis center, and require the system RMS < 0.3n. Satellite clock errors that fail verification need to be returned for review.

[0069] S8 outputs a precise clock error product after outlier removal, simultaneously providing an anomaly marking log and a consistency assessment report. This product can be directly used in scenarios such as multi-center clock error product integration, precise single-point positioning, low-Earth orbit satellite orbit determination, and crustal deformation monitoring.

[0070] Constructing a mathematical representation of GNSS clock bias history data: The raw pseudoranges included in the GNSS clock bias history data and carrier phase observations , is represented as:

[0071] ,

[0072] In the formula, s, T, r, and j represent the satellite, satellite system, receiver, and carrier frequency, respectively; and These represent pseudorange and carrier phase observable measurements, respectively. Represents the coefficient matrix; The increment represents the receiver's three-dimensional coordinates; c represents the speed of light in a vacuum. and denote receiver information and satellite clock bias information, respectively; and denote wet mapping function and zenith wet delay, respectively; denotes slant ionospheric delay at signal frequency; denotes a frequency-dependent multiplier factor; and denote frequency-dependent receiver uncalibrated code delay and frequency-dependent satellite uncalibrated code delay with respect to satellite s, respectively; denotes carrier wavelength; N is an integer phase ambiguity; and denote frequency-dependent receiver uncalibrated phase delay and frequency-dependent satellite uncalibrated phase delay with respect to satellite s, respectively; and denote the sum of measurement noise and multipath error for pseudorange and carrier phase observations, respectively.

[0073] Calculation of raw clock difference between analysis centers: In the process of clock anomaly value monitoring, the difference of clock data of each pair of analysis centers needs to be calculated. In the clock data of the pairwise analysis centers, the clock data of one analysis center is defined as the clock data of analysis center ac, and the clock data of the other analysis center is defined as the clock data of reference center ref.

[0074] The raw difference of clock of analysis center and reference center k epochs satellite s is calculated by the following formula:

[0075] ,

[0076] In the formula, is the satellite s clock correction number of analysis center ac at k epochs, is the satellite s clock correction number of reference center ref at k epochs; by subtracting the clock difference between the analysis center and the reference center, the clock difference of the satellite clock itself is eliminated.

[0077] Abnormal value detection method of clock difference per epoch estimation: after calculating the clock difference of two analysis centers, the clock difference of satellite s after removing the daily bias at k epochs is calculated by the following formula:

[0078] ,

[0079] In the formula, is the average value of the raw difference of satellite s daily clock.

[0080] The per-epoch estimated clock difference of satellite s at k epochs is calculated by the following formula: ​

[0081] ,

[0082] where, is the average of the differences of the clock corrections of all satellites after removing the daily bias of the clock corrections.

[0083] Detecting outliers by the method of quadratic difference of clock corrections: Let the clock corrections of the satellites processed by GNSS data be be expressed as:

[0084] ,

[0085] where, denotes the reference clock used by the analysis center in data processing; denotes the bias related to the analysis center and the satellite clock, denotes the determined part in the characteristics of the satellite clock itself; denotes the random part in the other observation noise and the characteristics of the clock itself;

[0086] The first difference of the satellite clock corrections is obtained by eliminating the clock variation characteristics of the satellite clock itself by differencing the clock correction of the analysis center ac with that of the reference center ref using the following formula:

[0087] ,

[0088] where, and are the clock correction of the reference center and the analysis center, respectively; denotes the difference between the reference clock of the analysis center ac and the reference clock of the reference center ref ; denotes the difference between the bias related to the analysis center and the satellite clock and the bias related to the reference center and the satellite clock ; denotes the difference between the determined bias part in the characteristics of the satellite clock itself and the reference center ;

[0089] A reference satellite is selected for each system, and the first difference of the satellite clock corrections is differenced between the reference satellite si and the other satellites sj of the system using the following formula to obtain the second difference of the satellite clock corrections, eliminating the influence of the reference clock related to the analysis center and leaving only the observation noise term and the term related to the data processing strategy of the analysis center:

[0090] ,

[0091] where, denote the difference between the clock correction of the satellite sj and the clock correction of the satellite si, respectively; is the second difference of clock error; represents the reference clock error value of the analysis center ac and the reference center ref; respectively represent the first difference of clock error of satellite sj and satellite si.

[0092] The linear transformation of clock error is calculated to detect outliers, including full-satellite linear transformation detection of outliers and single-satellite linear transformation detection of outliers:

[0093] The full-satellite linear transformation detection of outliers is calculated: the analysis center ac and the reference center ref are aligned by linear transformation, and a set of linear transformation parameters is estimated for each system. The conversion formula of any satellite s of a certain system at epoch k can be written as:

[0094] ,

[0095] In the formula, and are the clock errors of satellite s at epoch k in the reference center and the analysis center, , , represent the linear transformation conversion parameters, represents the GPS time of the clock error at epoch k, is the ephemeris reference time. Considering that a satellite system contains n satellites and each satellite has m epochs of data, an error equation is established to analyze and correct the clock error. The following formula is used to solve , , :

[0096] ,

[0097] ,

[0098] ,

[0099] ,

[0100] .

[0101] Single-satellite linear transformation detection of outliers: after estimating a set of linear transformation parameters for all satellite clock error data of a single system, it is found that there is still a relatively obvious linear deviation between satellites, which is closely related to the systematic deviation related to the analysis center and the satellite clock. By estimating a set of linear transformation parameters for each satellite clock error, the influence of the linear deviation part of the systematic deviation related to the analysis center and the satellite clock can be effectively eliminated. Further, the conversion formula of the analysis center clock error relative to the reference center clock error is re-expressed as:

[0102] ,

[0103] After linear transformation, the linear deviation between satellites is no longer obvious.

[0104] The above is only the preferred embodiment of the present application, and is not intended to limit the present application. Any slight modification, equivalent replacement and improvement made according to the technical essence of the present application to the above embodiment shall be included in the protection scope of the technical scheme of the present application.

Claims

1. A method for removing outliers of a GNSS precise clock product, characterized in that, The steps are as follows: S1, preprocessing the GNSS clock error historical data provided by each analysis center, pairing each analysis center with another according to the grouped preprocessed clock error historical data; S2, taking 3 times the standard deviation of the original clock error difference of each pair of analysis centers as a threshold value a, and eliminating abnormal values of the clock error data of the analysis centers exceeding the threshold value a through joint detection of the two-by-two paired analysis centers; S3, using a clock error per epoch alignment method to align the clock error of each pair of analysis centers per epoch, obtaining the instantaneous difference of the clock error, and eliminating abnormal values by taking 3 times the standard deviation as a threshold value b to generate an instantaneous difference sequence; S4, using an inter-analysis center linear transformation method to obtain the linear transformation parameters between the clock errors of each satellite of each pair of analysis centers, calculating the clock error residual after linear transformation, and eliminating abnormal values by taking 3 times the standard deviation as a threshold value c to generate a clock error residual sequence; S5, using a clock error second difference method, calculating the second difference of the clock error data of each pair of analysis centers, and eliminating abnormal values by taking 3 times the standard deviation as a threshold value d; to generate a second difference sequence; S6, repeating S2-S5 until the change rate of the number of abnormal epochs eliminated in the last two iterations is less than 5% or the maximum preset iteration number is reached, and outputting the final clock error product; S7, performing consistency verification on all final clock error products, and calculating the root mean square error of the clock error between each analysis center; S8, when the root mean square error of the clock error between each analysis center meets the requirements, outputting the precise clock error product after eliminating abnormal values; Constructing the mathematical expression of GNSS clock error historical data: the original pseudorange and carrier phase observation values included in the GNSS clock error historical data are expressed as: are expressed as: ， where s, T, r and j represent satellite, satellite system, receiver and carrier frequency, respectively; and represent pseudo-range and carrier phase observations, respectively; represents the coefficient matrix; represents the receiver three-dimensional coordinate increment; c represents the speed of light in vacuum; and represent receiver information and satellite clock error information, respectively; and represent wet mapping function and zenith wet delay, respectively; represents the slant ionospheric delay at the signal frequency; represents the frequency-dependent multiplicative factor; and represent the frequency-dependent receiver uncalibrated code delay and the frequency-dependent satellite uncalibrated code delay with respect to satellite s, respectively; represents the carrier wavelength; N is the integer phase ambiguity; and represent the frequency-dependent receiver uncalibrated phase delay and the frequency-dependent satellite uncalibrated phase delay with respect to satellite s, respectively; and represent the sum of the measurement noise and multipath error for the pseudo-range and carrier phase observations, respectively; The steps for calculating the difference of the clock error data of each pair of analysis centers are as follows: In the clock error data of the two-by-two paired analysis centers, the clock error data of one analysis center is defined as the clock error data of the analysis center ac, and the clock error data of the other analysis center is defined as the clock error data of the reference center ref. The original difference of the clock error of satellite s at k epochs between the analysis center ac and the reference center ref is calculated by the following formula: ， wherein denotes the clock correction number of satellite s at epoch k for the analysis center ac, is the clock correction number of satellite s at epoch k for the reference center ref; by differencing the clock between the analysis center ac and the reference center ref, the clock variation characteristics of the satellite clock itself are eliminated; The clock error difference of satellite s at k epochs after removing the daily bias is calculated by the following formula: ， In the formula, is the average of the original difference of the satellite s daily clock difference; The satellite s clock error difference is calculated at each ephemeris k using the following formula : ， wherein is the average of the clock error differences of all satellites after removal of the daily bias of the clock error at the epoch; The linear transformation parameters between the clock errors of each satellite of each pair of analysis centers are obtained by using an inter-analysis center linear transformation method, and the clock error residual after linear transformation is calculated: The linear transformation is used to align the clock error of the analysis center ac with the clock error of the reference center ref, and a set of linear transformation parameters is estimated for each satellite system T, and the conversion formula of any satellite s of any satellite system T at epoch k is represented as: ， wherein and are the clock errors of the satellite s at the epoch k in the reference center ref and in the analysis center ac, respectively, , , denote the linear transformation conversion parameters, denotes the GPS time at the epoch k of the clock error, is the ephemeris reference time; let the satellite system T comprise n satellites, each having data for m epochs; The linear transformation conversion parameters are solved using the following equation , , : ， ， ， ， The clock error residual between the analysis center ac linearly transformed and the reference center ref is calculated using the following equation : ， A set of linear transformation parameters is estimated for each satellite clock error in the clock error synthesis, and the influence of the linear part is eliminated: the transformation formula of the clock error of the analysis center ac relative to the clock error of the reference center ref is rewritten as: ， wherein is a set of linear transformation parameters calculated for each satellite of the system of the system.

2. The GNSS precise clock product outlier rejection method according to claim 1, characterized in that, The preprocessing process is as follows: all GNSS clock error historical data provided by all analysis centers are packaged and numbered, and then each GNSS clock error historical data package is combined with all other GNSS clock error historical data packages two by two.

3. The GNSS precise clock product outlier rejection method according to claim 1, wherein, Using a clock error second difference method, the second difference of each pair of analysis center clock error products is calculated: Satellite clock correction number for processing GNSS data is represented as: ， wherein denotes the reference clock used by the analysis center in data processing; denotes the bias associated with the analysis center and satellite clocks, denotes the part of the satellite clock's own characteristics determined; denotes the random part of the other observation noise and clock's own characteristics; The first difference of satellite clock error is obtained by eliminating the clock change characteristics of satellite clock itself by subtracting the analysis center ac clock correction number from the reference center ref using the following formula: ， wherein and are the reference center and analysis center clock corrections, respectively; denotes the difference between the reference clock of the analysis center ac and the reference clock of the reference center ref ; denotes the difference between the analysis center and satellite clock related biases and the reference center and satellite clock related biases ; denotes the difference between the satellite clock self-characteristic determined bias part and the reference center ; A reference satellite is selected for each satellite system, and the first difference of satellite clock error is subtracted between the reference satellite si and the satellite sj to obtain the second difference of satellite clock error using the following formula, i and j being satellite numbers, thereby eliminating the influence of the reference clock related to the analysis center, leaving only the observation noise term and the term related to the data processing strategy of the analysis center: ， wherein respectively represent the clock difference quadratic difference result of the satellite sj; the difference between the satellite sj and the satellite si; is the clock difference quadratic difference result; respectively represent the clock difference quadratic difference result of the satellite sj; respectively represent the clock difference quadratic difference result of the satellite sj.

4. The GNSS precise clock product outlier rejection method of claim 1, wherein, The epoch data to be eliminated needs to be marked in the data of each analysis center for more than or equal to 2 times; when the consistency of the final clock difference product is verified, the root mean square error of the clock difference between each system is calculated, and when the root mean square error of the satellite clock difference of each system is less than 0.3 ns, the verification is passed, otherwise the S2 loop is returned.

5. A computer device, comprising: The GNSS precise clock difference product outlier elimination method of any one of claims 1-4 is executed by the processor.

6. A computer-readable storage medium, characterized in that, The computer program stored in the computer readable storage medium is adapted to be loaded and executed by the processor to execute the GNSS precise clock difference product outlier elimination method of any one of claims 1-4.

Images