Navigation satellite clock error rapid encryption method based on M estimation

Abstract

The invention discloses a navigation satellite clock error rapid encryption method based on M estimation, and relates to the field of satellite navigation positioning, and the method comprises the following steps: S1, carrying out the consistency check of data, reading navigation satellite orbit and clock error data, observation station coordinates and receiver clock error data, and GNSS auxiliary data, and carrying out the data initialization and data preprocessing; s2, constructing an inter-epoch difference observation equation; s3, carrying out relative clock error parameter robust estimation; s4, obtaining a high-sampling relative clock difference; s5, performing clock correction synthesis according to the low-sampling absolute clock correction and the high-sampling relative clock correction, and outputting the high-sampling absolute clock correction; the calculation efficiency and robustness of clock correction encryption can be effectively improved, the calculation overhead is remarkably reduced while high precision is guaranteed, and the actual application requirements of low-orbit satellite precision orbit determination, high-precision dynamic positioning and the like are met.

Description

Technical Field

[0001] The invention relates to the field of satellite navigation and positioning, and in particular to a navigation satellite clock error fast encryption method based on M estimation. Background Art

[0002] High-precision estimation of navigation satellite clock errors is one of the key issues in the application of global navigation satellite systems. Existing satellite clock error estimation methods are usually based on the traditional least squares method with a typical sampling interval of 300 seconds. However, applications such as precise orbit determination of low-orbit satellites and real-time dynamic positioning require clock errors of 30 seconds or even higher sampling rates to ensure data processing accuracy. The existing 300-second sampling interval can no longer meet the needs of actual applications.

[0003] Therefore, there is an urgent need for clock error data with a higher sampling rate. The traditional satellite clock error estimation method is based on least squares estimation, the core idea of ​​which is to estimate parameters by minimizing the sum of squares of the residuals between the observed values ​​and the model predicted values.

[0004] However, when there are gross errors in the observed data, the estimation accuracy of this method will be significantly reduced. In order to overcome this limitation, multiple rounds of residual editing and outlier detection must be performed on the observed values ​​to eliminate abnormal data as much as possible. This purification process is very time-consuming and will significantly reduce computational efficiency.

[0005] In addition, when the clock error sampling interval is reduced from 300 seconds to 30 seconds, the amount of data also increases by an order of magnitude, which further reduces the computational efficiency of the traditional method and makes it difficult to adapt to the modern navigation and positioning requirements for high sampling rate, high timeliness and high-precision clock error products.

[0006] Therefore, a method for fast encryption of navigation satellite clock errors based on M estimation is provided to solve the above problems. Summary of the invention

[0007] The purpose of the present invention is to provide a method for fast encryption of navigation satellite clock errors based on M estimation, which can effectively deal with outliers and noise in observation data, improve computational efficiency, have strong robustness, and can reduce the number of iterations while maintaining high precision, thereby greatly improving timeliness and meeting the needs of near real-time navigation and positioning applications.

[0008] To achieve the above object, the present invention provides a method for fast encryption of navigation satellite clock errors based on M estimation, comprising the following steps:

[0009] S1: Check the consistency of data, read the navigation satellite orbit and clock error data, station coordinates and receiver clock error data, GNSS auxiliary data, and perform data initialization and data preprocessing;

[0010] S2: construct the inter-epoch differential observation equation;

[0011] S3: Estimation of relative clock error parameters;

[0012] S4: obtain high sampling relative clock difference;

[0013] S5: Perform clock error synthesis based on the low sampling absolute clock error and the high sampling relative clock error, and output the high sampling clock error.

[0014] Preferably, in step S1, the GNSS auxiliary data includes earth rotation parameters and tropospheric parameters, the data preprocessing process includes cycle slip detection and gross error elimination of GNSS observation values, and the GNSS observation values ​​include pseudorange observation values ​​and carrier phase observation values.

[0015] Preferably, step S2 specifically includes the following steps:

[0016] S21: Calculate pseudorange observations Pseudorange observations Specifically expressed as:

[0017]

[0018] Among them, r represents the site, s represents the satellite, represents the geometric distance, c represents the speed of light, δt r (t) represents the receiver clock error, δt s (t) represents the satellite clock error, represents the ionospheric delay, represents the tropospheric delay, ε P represents the pseudorange observation noise;

[0019] S22: Between adjacent epochs t and t+1, the pseudorange observations Perform differential analysis, pseudorange observation value differential results Specifically expressed as:

[0020]

[0021] in, Indicates the change in geometric distance caused by the relative motion between the satellite and the receiver, Δδt r Indicates the change in receiver clock error, Δδt s represents the change in satellite clock error, represents the change in ionospheric delay, represents the change in tropospheric delay, Δε P Represents the change in pseudorange observation noise;

[0022] S23: Calculate carrier phase observation value Carrier phase observations Specifically expressed as:

[0023]

[0024] Where λ represents the carrier wavelength, represents the integer ambiguity, ε Φ represents the carrier phase observation noise;

[0025] S24: Between adjacent epochs t and t+1, the carrier phase observation value Perform differential analysis, and the carrier phase observation value differential result Specifically expressed as:

[0026]

[0027] Among them, Δε Φ Represents the change in carrier phase observation noise.

[0028] Preferably, step S3 specifically includes the following steps:

[0029] S31: Based on the posterior residual r i Perform M estimation weighting;

[0030] S32: Relative clock error parameter estimation is performed using M estimation and iterative weighted least squares method.

[0031] Preferably, step S31 specifically includes the following steps:

[0032] Step 1: Calculate the posterior residual r i , the post-test residual r i Specifically expressed as:

[0033] r i =O i -F(X)

[0034] Among them, O i represents the observed value, X represents the estimated parameter, and F(X) represents the estimated value;

[0035] Step 2: Perform weighted adjustment through the M-estimation objective function. The M-estimation objective function is specifically expressed as:

[0036]

[0037] Where ρ(x) represents the M estimation loss function, and σ represents the standard deviation of the observation noise;

[0038] Step 3: Calculate the robustness weight w through M estimation weight function i , robustness weight w i Specifically expressed as:

[0039]

[0040] Preferably, step S32 specifically includes the following steps:

[0041] Step 1: Construct a linear model of the observation equation, which is specifically expressed as:

[0042] y=Ax+ε

[0043] ε~N(0,σ 2 I)

[0044] Where y represents the observation value vector, A represents the design matrix, x represents the parameter vector to be estimated, ε represents the observation noise vector, N(0,σ 2 I) represents a zero-mean normal distribution;

[0045] Step 2: Solve the parameters by minimizing the residual sum of squares using the least squares method. The parameter solution method is specifically expressed as:

[0046]

[0047] v i =y i -(Ax) i

[0048] in, represents the estimated parameter vector, v i Represents the residual, specifically the deviation between the observed value and the fitted value, y i represents the observed value;

[0049] Step 3: Introduce the loss function ρ(v i ), so that the loss function ρ(v i ) instead of square loss A new parameter solving method is obtained, which is specifically expressed as:

[0050]

[0051] Among them, the loss function ρ(v i ) is set to the loss function ρ estimated by WelschM W (v), WelschM estimates the loss function ρ W (v) is specifically expressed as:

[0052]

[0053] Among them, v represents the residual, and C represents the scale parameter;

[0054] Step 4: The weight function w estimated by WelschM W(v) Assign weights to the observations, and use the weight function w estimated by WelschM W (v) is specifically expressed as:

[0055]

[0056] Step 5: Solve the parameters based on M estimation through iterative weighted least squares method. The specific solution process is as follows:

[0057] N k =N c,k +N P,k

[0058]

[0059] Among them, N k Represents the observation noise covariance matrix at the current moment, N c,k represents the carrier phase observation noise covariance matrix, N P,k represents the pseudorange observation noise covariance matrix, A c represents the design matrix of carrier phase observation, A P represents the design matrix of pseudorange phase observation, W c,k Represents the carrier phase residual weight matrix, which is defined as the diagonal matrix after weighting the residual, W P,k represents the pseudorange observation residual weight matrix, represents the normalized carrier phase residual of the ith satellite, represents the normalized carrier phase residual of the ith satellite, represents the estimated value of the pseudorange observation noise variance, represents the estimated value of the carrier phase observation noise variance, Δy c Represents the residual vector of the carrier phase observation, which is specifically set to the deviation between the observed value and the fitted value, Δy P represents the residual vector of pseudorange observation, represents the parameter vector estimated at the next moment, v c,k+1 Indicates the carrier phase observation correction value at the next moment, v P,k+1 represents the pseudorange observation correction value at the next moment, represents the weighting function based on the normalized residual, represents the weighting function of pseudorange observations, represents the trace of the noise covariance matrix.

[0060] Preferably, in step 5, the iterative termination condition of the iterative weighting is specifically expressed as:

[0061]

[0062] Preferably, in step 5, the convergence condition of the M estimated weight is specifically expressed as:

[0063]

[0064] in, Indicates that in the kth iteration, for the pseudorange observation value The M estimation weights, Indicates that in the kth iteration, for the carrier phase observation value The M estimation weights.

[0065] Preferably, step S5 specifically includes the following steps:

[0066] S51: The clock error variation equation is obtained through satellite observation data. The clock error variation equation is specifically expressed as:

[0067] Δδ(t i+1 ,t i )=δ(t i+1 )-δ(t i )

[0068] Among them, Δδ(ti +1 ,ti) represents the change in clock error, δ(t i+1 ) represents the time ti +1 The satellite clock error, δ(t i ) represents the time t i The satellite clock error;

[0069] S52: Establish a linear observation equation system, and transform the low-sampled absolute clock error δ fix (t) and the high-sampling relative clock error δ(t) are fused, and the linear observation equations are specifically expressed as:

[0070]

[0071] Among them, δ fix (t n ) indicates that at time t n The absolute clock error obtained from the low sampling rate data, δ(t n ) represents the time t n The satellite clock error;

[0072] S53: Solve the linear observation equations by M estimation and iterative weighted least squares method.

[0073] Therefore, the present invention adopts the above-mentioned navigation satellite clock error fast encryption method based on M estimation, which has the following beneficial effects:

[0074] (1) In the process of estimating the relative clock error between epochs, the present invention adopts the M estimation method to perform parameter estimation through M estimation, effectively suppressing the influence of outliers. At the same time, the iterative weighted least squares method is used to dynamically adjust the weights of the observation data to achieve rapid convergence of the estimated parameters, greatly improving the calculation efficiency. Compared with the traditional LS method, it does not require multiple iterative residual editing, greatly reducing the calculation time, and meeting the high effectiveness requirements;

[0075] (2) In the clock error synthesis process, the present invention adopts the M estimation method to estimate the dynamic weighting of the observation data, so as to achieve low sampling absolute clock error and high sampling relative clock error, thereby improving the algorithm's robustness;

[0076] (3) In the clock error synthesis process, the present invention synthesizes the high-sampling relative clock error and the low-sampling absolute clock error through the time series to achieve the conversion of the low-sampling clock error to the high-sampling clock error, thereby ensuring the continuity and accuracy of the time series;

[0077] (4) The present invention integrates the M-estimation anti-error technology and the clock error encryption strategy, which significantly improves the computational efficiency and robustness while ensuring high accuracy, and provides an efficient and reliable technical solution for high-frequency satellite navigation applications.

[0078] The method scheme of the present invention is further described in detail below through the drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0079] Figure 1 This is a flow chart of a method for fast encryption of navigation satellite clock errors based on M estimation of the present invention;

[0080] Figure 2 It is a schematic diagram of the residual comparison results of a navigation satellite clock error fast encryption method based on M estimation and the LS method in a 30-second difference sampling task. DETAILED DESCRIPTION

[0081] The method scheme of the present invention is further described below through drawings and embodiments.

[0082] Unless otherwise defined, method terms or scientific terms used in the present invention shall have the common meanings understood by one of ordinary skill in the art to which the present invention belongs.

[0083] The words "include" or "comprises" and the like used in the present invention mean that the elements before the word include the elements listed after the word, and do not exclude the possibility of also including other elements. The orientation or position relationship indicated by the terms "inside", "outside", "upper", "lower", etc. is based on the orientation or position relationship shown in the drawings, which is only for the convenience of describing the present invention and simplifying the description, and does not indicate or imply that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation. Therefore, it cannot be understood as a limitation of the present invention. When the absolute position of the described object changes, the relative position relationship may also change accordingly. In the present invention, unless otherwise clearly specified and limited, the terms "attachment" and the like should be understood in a broad sense, for example, it can be a fixed connection, a detachable connection, or an integral body; it can be directly connected, or indirectly connected through an intermediate medium, and it can be the internal connection of two elements or the interaction relationship between two elements. For ordinary technicians in this field, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0084] Example

[0085] like Figure 1 As shown, the present invention provides a method for fast encryption of navigation satellite clock errors based on M estimation, comprising the following steps:

[0086] S1: Check the consistency of data, read the navigation satellite orbit and clock error data, station coordinates and receiver clock error data, GNSS auxiliary data, and perform data initialization and data preprocessing;

[0087] In step S1, the GNSS auxiliary data includes the earth's rotation parameters and tropospheric parameters. The data preprocessing process includes cycle slip detection and gross error elimination of GNSS observations to provide relatively clean observations and ensure data consistency and availability. GNSS observations include pseudorange observations and carrier phase observations.

[0088] S2: construct the inter-epoch differential observation equation and the inter-epoch differential ionospheric-free combined observation quantity;

[0089] Step S2 specifically includes the following steps:

[0090] S21: Calculate pseudorange observations Pseudorange observations Specifically expressed as:

[0091]

[0092] Among them, r represents the site, s represents the satellite, represents the geometric distance, c represents the speed of light, δt r (t) represents the receiver clock error, δt s(t) represents the satellite clock error, represents the ionospheric delay, represents the tropospheric delay, ε P represents the pseudorange observation noise.

[0093] S22: Between adjacent epochs t and t+1, the pseudorange observations Perform differential analysis, pseudorange observation value differential results Specifically expressed as:

[0094]

[0095] in, Indicates the change in geometric distance caused by the relative motion between the satellite and the receiver, Δδt r It represents the change of the receiver clock error. Since the receiver clock error changes linearly between adjacent epochs, Δδt r Can be approximately eliminated, Δδt s represents the change in satellite clock error, represents the change in ionospheric delay, represents the change in tropospheric delay, Δε P Represents the change in pseudorange observation noise.

[0096] S23: Calculate carrier phase observation value Carrier phase observations Specifically expressed as:

[0097]

[0098] Where λ represents the carrier wavelength, represents the integer ambiguity, ε Φ represents the carrier phase observation noise.

[0099] S24: Between adjacent epochs t and t+1, the carrier phase observation value Perform differential analysis, and the carrier phase observation value differential result Specifically expressed as:

[0100]

[0101] Among them, Δε Φ represents the change in carrier phase observation noise due to integer ambiguity is constant in time, so the integer ambiguity can be eliminated, which makes the carrier phase observation more stable than the pseudorange observation and becomes the main data source for high-precision clock error estimation. In this embodiment, the carrier phase epoch-to-epoch differential observation is also used for clock error encryption.

[0102] S3: perform robust estimation of relative clock error parameters;

[0103] Step S3 specifically includes the following steps:

[0104] S31: Based on the posterior residual r i Perform M estimation weighting;

[0105] Step S31 specifically includes the following steps:

[0106] Step 1: Calculate the posterior residual r i , the post-test residual r i It refers to the residual obtained by calculating the observed value using the estimated parameters and comparing it with the original observed value after the parameter estimation is completed. In M estimation, the posterior residual r i Used for weighted adjustment to reduce the impact of abnormal observations, the post-test residual r i Specifically expressed as:

[0107] r i =O i -F(X)

[0108] Among them, O i represents the observed value, X represents the estimated parameter, and F(X) represents the estimated value.

[0109] Step 2: Perform weighted adjustment through the M-estimation objective function. The M-estimation objective function is specifically expressed as:

[0110]

[0111] Where ρ(x) represents the M estimation loss function and σ represents the standard deviation of the observation noise.

[0112] Step 3: Calculate the robustness weight w through M estimation weight function i Common M estimation weight functions include Huber loss function and Welsch loss function, and robust weight w i Specifically expressed as:

[0113]

[0114] S32: Relative clock error parameters are estimated through M estimation and iterative weighted least squares method, and the weights of observation data are dynamically adjusted to ensure that the estimated clock error results are more robust.

[0115] Step S32 specifically includes the following steps:

[0116] Step 1: M estimation assigns different weights to the observation data to weaken the influence of outliers on the estimation results. In order to solve the parameter estimation problem, this embodiment combines M estimation and iterative weighted least squares method to estimate the clock error and construct a linear model of the observation equation, which is specifically expressed as follows:

[0117] y=Ax+ε

[0118] ε~N(0,σ 2 I)

[0119] Where y represents the observation value vector, A represents the design matrix, x represents the parameter vector to be estimated, ε represents the observation noise vector, N(0,σ 2 I) represents a normal distribution with zero mean.

[0120] Step 2: Solve the parameters by minimizing the residual sum of squares using the least squares method. The parameter solution method is specifically expressed as:

[0121]

[0122] v i =y i -(Ax) i

[0123] in, represents the estimated parameter vector, v i Represents the residual, that is, the deviation between the observed value and the fitted value, y i Represents the observed value.

[0124] Step 3: When there are outliers in the observed data, the LS method is sensitive to them, so the M estimate is introduced to reduce the impact of outliers. The loss function ρ(v i ), so that the loss function ρ(v i ) instead of square loss A new parameter solving method is obtained, which is specifically expressed as:

[0125]

[0126] Among them, the loss function ρ(v i ) is set to the loss function ρ estimated by WelschM W (v), loss function ρ estimated by WelschM W (v) is specifically expressed as:

[0127]

[0128] Among them, v represents the residual; C represents the scale parameter, which controls the sensitivity of the weight to the error.

[0129] The goal of M estimation is to minimize the sum of the above loss functions, thereby reducing the impact of outliers while maintaining the stability of the estimate. The core idea of ​​M estimation is to assign weights to different observations to weaken the impact of outliers.

[0130] Step 4: The weight function w estimated by WelschM W (v) Assign weights to the observations, and use the weight function w estimated by WelschM W (v) is specifically expressed as:

[0131]

[0132] For observations with large outliers, the weights approach zero, reducing their impact on parameter estimates.

[0133] Step 5: Solve the parameters based on M estimation through iterative weighted least squares method. The specific solution process is as follows:

[0134] N k =N c,k +N P,k

[0135]

[0136] Among them, N k Represents the observation noise covariance matrix at the current moment, N c,k represents the carrier phase observation noise covariance matrix, N P,k represents the pseudorange observation noise covariance matrix, A c represents the design matrix of carrier phase observation, A P represents the design matrix of pseudorange phase observation, W c,k Represents the carrier phase residual weight matrix, which is defined as the diagonal matrix after weighting the residual, W P,k represents the pseudorange observation residual weight matrix, represents the normalized carrier phase residual of the ith satellite, represents the normalized carrier phase residual of the ith satellite, represents the estimated value of the pseudorange observation noise variance, represents the estimated value of the carrier phase observation noise variance, Δy c Represents the residual vector of the carrier phase observation, which is specifically set to the deviation between the observed value and the fitted value, Δy P represents the residual vector of pseudorange observation, represents the parameter vector estimated at the next moment, v c,k+1 Indicates the carrier phase observation correction value at the next moment, v P,k+1 represents the pseudorange observation correction value at the next moment, represents a weighting function based on the normalized residual, usually used for robust estimation, represents the weighting function of pseudorange observations, Represents the trace of the noise covariance matrix, which is used to compute the normalization factor for the weighted residual.

[0137] In the iterative weighted least squares method of M estimation, it is necessary to set a suitable convergence criterion to ensure that the iteration converges within a reasonable error range. Therefore, an iterative convergence criterion is expressed by mathematical conditions. The core idea is that the relative change of the variance estimate should be small enough, and the change of the weight function should be small enough.

[0138] In step 5, the iterative termination condition of the iterative weighting is specifically expressed as:

[0139]

[0140] When the relative change of the variance estimate between two consecutive iterations is less than 1%, it indicates that the noise variance has become stable and no further updates are performed.

[0141] In step 5, the convergence condition of the M estimated weight is specifically expressed as:

[0142]

[0143] in, Indicates that in the kth iteration, for the pseudorange observation value The M estimation weights, Indicates that in the kth iteration, for the carrier phase observation value The M estimation weights.

[0144] When the change in the M-estimated weights of all observations is less than 0.01, it means that the weights have stabilized and no further updates are required.

[0145] S4: obtain high sampling relative clock difference;

[0146] S5: Perform clock error synthesis based on low-sampling absolute clock error and high-sampling relative clock error, and output high-sampling clock error. Encrypt low-sampling clock error to high-sampling clock error by synthesizing high-sampling relative clock error with low-sampling absolute clock error. For example, encrypt low-sampling clock error of 300 seconds interval to high-sampling clock error of 30 seconds or 5 seconds. In the clock error encryption strategy, low-sampling absolute clock error is provided based on carrier phase observation value, and clock error variation, i.e. high-sampling relative clock error, is provided based on clock error estimation process. Absolute clock error information and relative clock error information are combined to obtain high-precision clock error sequence. The low-sampling clock error is merged with the high-sampling clock error to ensure the continuity and accuracy of the time series, so as to meet the application requirements of real-time dynamic positioning and precise orbit determination of low-orbit satellites, and provide more stable and efficient clock error information support for high-precision navigation applications.

[0147] Step S5 specifically includes the following steps:

[0148] S51: The clock error variation equation is obtained through phase observation. The clock error variation equation is specifically expressed as:

[0149] Δδ(t i+1 ,t i )=δ(t i+1 )-δ(t i )

[0150] Among them, Δδ(t i+1 ,t i ) represents the change in clock error, δ(t i+1 ) represents the time t i+1 The satellite clock error, δ(t i ) represents the time t i The satellite clock error.

[0151] S52: Establish a linear observation equation system, and transform the low-sampled absolute clock error δ fix (t) and the high-sampling relative clock error δ(t) are fused, and the linear observation equations are specifically expressed as:

[0152]

[0153] Among them, δ fix (t n ) indicates that at time t n The absolute clock error obtained from the low sampling rate data, δ(t n ) represents the time t n The satellite clock error.

[0154] S53: Solve the linear observation equations by M estimation and iterative weighted least squares method.

[0155] The clock error parameter dimensions are optimized based on time series characteristics to improve computational efficiency. At the same time, the weights of observation data are dynamically adjusted in combination with robust estimation methods to further improve the accuracy and reliability of clock error information fusion.

[0156] This embodiment aims at the key technical problem of fast solution of navigation satellite clock error, and carries out comparative study on algorithms of high sampling clock error estimation. The experiment is based on observation data of 100 monitoring stations around the world, and solves clock errors of four major satellite navigation systems, namely GPS, Galileo, GLONASS and BDS. The sampling navigation satellites are set to 113, the sampling rate is set to 30 seconds, and the hardware platform adopts a high-performance computing server equipped with Intel Xeon Gold 6326 processor. The main frequency of Intel Xeon Gold 6326 processor is set to 2.90 GHz, and the computing core is set to 28.

[0157] Under the same experimental conditions, the traditional least squares method needs to go through two stages: parameter estimation and residual editing, which take 2558 seconds and 276 seconds respectively, and the total calculation time is 2836 seconds. The fast encryption method of navigation satellite clock error based on M estimation proposed in this embodiment dynamically adjusts the observation weights by introducing M estimation, does not require multiple iterative residual editing, and combines the iterative weighted least squares method to achieve rapid convergence, greatly reducing the calculation time, and successfully shortening the calculation time to 276 seconds.

[0158] The residuals of the method of this embodiment and the LS method in the 30-second difference sampling task are compared. The residuals are obtained by subtracting the method of this embodiment from the official IGS clock error, and the residuals are also obtained by subtracting the LS method from the official IGS clock error. The residual distributions of the two are further compared, and a bar chart is drawn to analyze the differences in clock error estimation accuracy between different methods.

[0159] like Figure 2 As shown in the figure, the accuracy of this implementation method is better than the traditional LS method. At the same time, the calculation efficiency is significantly improved by 90%, and the operation speed is increased by about 10 times, achieving an order of magnitude breakthrough. This breakthrough is mainly due to the dynamic adjustment of observation weights through M estimation, which effectively suppresses the influence of outliers and greatly improves the calculation efficiency. Compared with the traditional LS method, there is no need for multiple iterative residual editing. Combined with the iterative weighted least squares method, rapid convergence can be achieved, and the calculation time is greatly reduced to meet the real-time requirements.

[0160] Experimental results show that the fast encryption method of navigation satellite clock error based on M estimation proposed in this embodiment not only meets the requirements of high-precision GNSS data processing, but also has excellent computing performance and robustness characteristics. It is particularly suitable for real-time precise single-point positioning in high dynamic environments, satellite-based augmentation systems and other navigation application scenarios with strict timeliness requirements, providing reliable technical support for efficient data processing of the new generation of navigation systems.

[0161] Therefore, the present invention adopts the above-mentioned M-estimation-based navigation satellite clock error fast encryption method, which can effectively process outliers and noise in the observation data, improve calculation efficiency, have strong robustness, and can reduce the number of iterations while maintaining high precision, thereby greatly improving timeliness and meeting the needs of real-time navigation and positioning applications.

[0162] Finally, it should be noted that the above embodiments are only used to illustrate the method scheme of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, ordinary method personnel in the field should understand that they can still modify or replace the method scheme of the present invention with equivalents, and these modifications or equivalent replacements cannot cause the modified method scheme to deviate from the spirit and scope of the method scheme of the present invention.

Claims

1. A method for fast encryption of navigation satellite clock errors based on M estimation, characterized in that: The following steps are involved: S1: Check the consistency of data, read the navigation satellite orbit and clock error data, station coordinates and receiver clock error data, GNSS auxiliary data, and perform data initialization and data preprocessing; S2: construct the inter-epoch differential observation equation; S3: Estimation of relative clock error parameters; S4: obtain high sampling relative clock difference; S5: Perform clock error synthesis based on the low sampling absolute clock error and the high sampling relative clock error, and output the high sampling clock error.

2. The method for fast encryption of navigation satellite clock errors based on M estimation according to claim 1, characterized in that: In step S1, the GNSS auxiliary data includes earth rotation parameters and tropospheric parameters, and the data preprocessing process includes cycle slip detection and gross error elimination of GNSS observation values, and the GNSS observation values ​​include pseudorange observation values ​​and carrier phase observation values.

3. The method for fast encryption of navigation satellite clock errors based on M estimation according to claim 1, characterized in that: Step S2 specifically includes the following steps: S21: Calculate pseudorange observations Pseudorange observations Specifically expressed as: Among them, r represents the site, s represents the satellite, represents the geometric distance, c represents the speed of light, δt r (t) represents the receiver clock error, δt s (t) represents the satellite clock error, represents the ionospheric delay, represents the tropospheric delay, ε P represents the pseudorange observation noise; S22: Between adjacent epochs t and t+1, the pseudorange observations Perform differential analysis, pseudorange observation value differential results Specifically expressed as: in, Indicates the change in geometric distance caused by the relative motion between the satellite and the receiver, Δδt r Indicates the change in receiver clock error, Δδt s represents the change in satellite clock error, represents the change in ionospheric delay, represents the change in tropospheric delay, Δε P Represents the change in pseudorange observation noise; S23: Calculate carrier phase observation value Carrier phase observations Specifically expressed as: Where λ represents the carrier wavelength, represents the integer ambiguity, ε Φ represents the carrier phase observation noise; S24: Between adjacent epochs t and t+1, the carrier phase observation value Perform differential analysis, and the carrier phase observation value differential result Specifically expressed as: Among them, Δε Φ Represents the change in carrier phase observation noise.

4. The method for fast encryption of navigation satellite clock errors based on M estimation according to claim 1, characterized in that: Step S3 specifically includes the following steps: S31: Based on the posterior residual r i Perform M estimation weighting; S32: Relative clock error parameter estimation is performed using M estimation and iterative weighted least squares method.

5. The method for fast encryption of navigation satellite clock errors based on M estimation according to claim 4, characterized in that: Step S31 specifically includes the following steps: Step 1: Calculate the posterior residual r i , the post-test residual r i Specifically expressed as: r i =O i -F(X) Among them, O i represents the observed value, X represents the estimated parameter, and F(X) represents the estimated value; Step 2: Perform weighted adjustment through the M-estimation objective function. The M-estimation objective function is specifically expressed as: Where ρ(x) represents the M estimation loss function, and σ represents the standard deviation of the observation noise; Step 3: Calculate the robustness weight w through M estimation weight function i , robustness weight w i Specifically expressed as:

6. The method for fast encryption of navigation satellite clock errors based on M estimation according to claim 4, characterized in that: Step S32 specifically includes the following steps: Step 1: Construct a linear model of the observation equation, which is specifically expressed as: y=Ax+ε ε~N(0,σ 2 I) Where y represents the observation value vector, A represents the design matrix, x represents the parameter vector to be estimated, ε represents the observation noise vector, N(0,σ 2 I) represents a normal distribution with zero mean; Step 2: Solve the parameters by minimizing the residual sum of squares using the least squares method. The parameter solution method is specifically expressed as: v i =y i -(Ax) i in, represents the estimated parameter vector, v i Represents the residual, specifically the deviation between the observed value and the fitted value, y i represents the observed value; Step 3: Introduce the loss function ρ(v i ), so that the loss function ρ(v i ) instead of square loss A new parameter solving method is obtained, which is specifically expressed as: Among them, the loss function ρ(v i ) is set to the loss function ρ estimated by WelschM W (v), WelschM estimates the loss function ρ W (v) is specifically expressed as: Among them, v represents the residual, C represents the scale parameter; Step 4: The weight function w estimated by WelschM W (v) Assign weights to the observations, and use the weight function w estimated by WelschM W (v) is specifically expressed as: Step 5: Solve the parameters based on M estimation through iterative weighted least squares method. The specific solution process is as follows: N k =N c,k +N P,k Among them, N k Represents the observation noise covariance matrix at the current moment, N c,k represents the carrier phase observation noise covariance matrix, N P,k represents the pseudorange observation noise covariance matrix, A c represents the design matrix of carrier phase observation, A P represents the design matrix of pseudorange phase observation, W c,k Represents the carrier phase residual weight matrix, which is defined as the diagonal matrix after weighting the residual, W P,k represents the pseudorange observation residual weight matrix, represents the normalized carrier phase residual of the ith satellite, represents the normalized carrier phase residual of the ith satellite, represents the estimated value of the pseudorange observation noise variance, represents the estimated value of the carrier phase observation noise variance, Δy c Represents the residual vector of the carrier phase observation, which is specifically set to the deviation between the observed value and the fitted value, Δy P represents the residual vector of pseudorange observation, represents the parameter vector estimated at the next moment, v c,k+1 Indicates the carrier phase observation correction value at the next moment, v P,k+1 represents the pseudorange observation correction value at the next moment, represents the weighting function based on the normalized residual, represents the weighting function of pseudorange observations, represents the trace of the noise covariance matrix.

7. The method for fast encryption of navigation satellite clock errors based on M estimation according to claim 6, characterized in that: In step 5, the iterative termination condition of the iterative weighting is specifically expressed as:

8. The method for fast encryption of navigation satellite clock errors based on M estimation according to claim 6, characterized in that: In step 5, the convergence condition of the M estimated weight is specifically expressed as: in, Indicates that in the kth iteration, for the pseudorange observation value The M estimation weights, Indicates that in the kth iteration, for the carrier phase observation value The M estimation weights.

9. The method for fast encryption of navigation satellite clock errors based on M estimation according to claim 1, characterized in that: Step S5 specifically includes the following steps: S51: The clock error variation equation is obtained through satellite observation data. The clock error variation equation is specifically expressed as: Δδ(t i+1 ,t i )=δ(t i+1 )-δ(t i ) Among them, Δδ(t i+1 ,t i ) represents the change in clock error, δ(t i+1 ) represents the time t i+1 The satellite clock error, δ(t i ) represents the time t i The satellite clock error; S52: Establish a linear observation equation system, and transform the low-sampled absolute clock error δ fix (t) and the high-sampling relative clock error δ(t) are fused, and the linear observation equations are specifically expressed as: Among them, δ fix (t n ) indicates that at time t n The absolute clock error obtained from the low sampling rate data, δ(t n ) represents the time t n The satellite clock error; S53: Solve the linear observation equations by M estimation and iterative weighted least squares method.

Images