A fast encryption method for navigation satellite clock errors based on M estimation

Abstract

The application discloses a kind of navigation satellite clock difference fast encryption methods based on M estimation, it is related to satellite navigation positioning field, comprising the following steps: S1: consistency check is carried out to data, reads navigation satellite orbit and clock difference data, station coordinates and receiver clock difference data, GNSS auxiliary data, carries out data initialization and data preprocessing;S2: epoch difference observation equation is constructed;S3: relative clock difference parameter robust estimation is carried out;S4: high sampling relative clock difference is obtained;S5: clock difference is synthesized according to low sampling absolute clock difference and high sampling relative clock difference, and high sampling absolute clock difference is output;It can effectively improve the calculation efficiency and robustness of clock difference encryption, while ensuring high precision, significantly reduce the calculation overhead, meet the actual application needs such as low-orbit satellite precise orbit determination and high-precision dynamic positioning.

Description

TECHNICAL FIELD

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

[0002] High-precision estimation of navigation satellite clock error is one of the key problems in global navigation satellite system applications. The existing satellite clock error estimation method is usually based on the traditional least squares method, and its typical sampling interval is 300 seconds. However, low-orbit satellite precise orbit determination and real-time dynamic positioning require clock error with a sampling rate of 30 seconds or even higher to ensure data processing accuracy. The existing 300-second sampling interval cannot meet the actual application requirements.

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

[0004] However, when there are gross errors in the observation 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 observation values to eliminate abnormal data as much as possible. This purification process is very time-consuming and will significantly reduce the 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, making it difficult to adapt to the demand for high-sampling-rate, high-timeliness and high-precision clock error products in modern navigation and positioning.

[0006] Therefore, a navigation satellite clock error fast encryption method based on M estimation is provided to solve the above problems. SUMMARY

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

[0008] To achieve the above purpose, the present application provides a navigation satellite clock error fast encryption method based on M estimation, comprising the following steps:

[0009] S1: Perform consistency check on data, read 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 interepoch difference observation equation;

[0011] S3: performing relative clock bias parameter estimation;

[0012] S4: obtaining high-sampling relative clock bias;

[0013] S5: performing clock bias synthesis according to low-sampling absolute clock bias and high-sampling relative clock bias, and outputting high-sampling clock bias.

[0014] Preferably, in step S1, the GNSS auxiliary data comprises earth rotation parameters and troposphere parameters, and the data preprocessing process comprises cycle slip detection and gross error elimination on GNSS observation values, and the GNSS observation values comprise pseudo-range observation values and carrier phase observation values.

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

[0016] S21: calculating pseudo-range observation values pseudo-range observation values Specifically represented as:

[0017]

[0018] Wherein, r represents a station, s represents a satellite, represents a geometric distance, c represents the speed of light, δt r (t) represents a receiver clock bias, δt s (t) represents a satellite clock bias, represents ionospheric delay, represents tropospheric delay, ε P represents pseudo-range observation noise;

[0019] S22: performing difference on pseudo-range observation values between adjacent epochs t and t+1, and the pseudo-range observation difference result Specifically represented as:

[0020]

[0021] Wherein, represents a change amount of geometric distance caused by relative motion of the satellite and the receiver, Δδt r represents a change amount of receiver clock bias, Δδt s represents a change amount of satellite clock bias, represents a change amount of ionospheric delay, represents a change amount of tropospheric delay, Δε P represents a change amount of pseudo-range observation noise;

[0022] S23: calculating carrier phase observation values carrier phase observation values Specifically represented as:

[0023]

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

[0025] S24: between adjacent epochs t and t+1, the carrier phase observation value is differenced, and the carrier phase observation difference result is specifically represented as:

[0026]

[0027] where Δε Φ represents the change amount of carrier phase observation noise.

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

[0029] S31: based on the posteriori residual r i M-estimation weighting is performed;

[0030] S32: relative clock difference parameter estimation is performed through M-estimation and iterative weighted least squares method.

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

[0032] Step 1: calculate the posteriori residual r i , the posteriori residual r i is specifically represented as:

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

[0034] where O i represents an observation value, X represents an estimated parameter, and F(X) represents an estimated value;

[0035] Step 2: weighting adjustment is performed through M-estimation objective function, and the M-estimation objective function is specifically represented as:

[0036]

[0037] where ρ(x) represents an M-estimation loss function, and σ represents an observation noise standard deviation;

[0038] Step 3: the robust weight w i is calculated through an M-estimation weight function, and the robust weight w i is specifically represented as:

[0039]

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

[0041] Step one: construct a linearization model of the observation equation, specifically represented as:

[0042] y = Ax + ε

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

[0044] Wherein, 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 the zero mean normal distribution;

[0045] Step two: solve the parameters by minimizing the residual sum of squares through the least square method, and the parameter solving method is specifically represented as:

[0046]

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

[0048] Wherein, represents the estimated parameter vector, v i represents the residual, specifically the deviation between the observation value and the fitting value, y i represents the observation value;

[0049] Step three: introduce the loss function ρ(v i ) through M estimation, so that the loss function ρ(v i ) replaces the square loss Get a new parameter solving method, and the new parameter solving method is specifically represented as:

[0050]

[0051] Wherein, the loss function ρ(v i ) is set as the loss function ρ W (v) of Welsch M estimation, and the Welsch M estimation loss function ρ W (v) is specifically represented as:

[0052]

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

[0054] Step four: through the weight function w W(v) Assign weights to the observations, the weight function w of the Welsch M-estimator W (v) is specifically expressed as:

[0055]

[0056] Step five: Solve the parameters based on M-estimation by the iterative weighted least squares method, and the solving process is specifically expressed as:

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

[0058]

[0059] Wherein, N k represents the observation noise covariance matrix at the current time, N c,k represents the carrier phase observation noise covariance matrix, N P,k represents the pseudo-range observation noise covariance matrix, A c represents the design matrix of the carrier phase observation, A P represents the design matrix of the pseudo-range 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 pseudo-range observation residual weight matrix, represents the normalized carrier phase residual of the i-th satellite, represents the normalized carrier phase residual of the i-th satellite, represents the estimated value of the pseudo-range 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 as the deviation between the observation value and the fitting value, Δy P represents the residual vector of the pseudo-range observation, represents the estimated parameter vector at the next time, v c,k+1 represents the carrier phase observation correction value at the next time, v P,k+1 represents the pseudo-range observation correction value at the next time, represents the weighting function based on the normalized residual, represents the weighting function of the pseudo-range observation, represents the trace of the noise covariance matrix.

[0060] Preferably, in step five, the iteration termination condition of the iterative weighted 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 pseudo-range observation value The M estimation weights, Indicates the carrier phase observation value at the kth iteration 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 Satellite clock error, δ(t i ) represents the time t i 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 low sampling rate data, δ(t n ) represents the time t n 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) The application adopts M estimation method in the relative clock difference estimation process between epochs, estimates parameters through M estimation, effectively suppresses the influence of abnormal values, simultaneously adjusts the observation data weight dynamically by using iterative weighted least square method, realizes the rapid convergence of the estimated parameters, greatly improves the calculation efficiency, compared with the traditional LS method, does not need multiple iteration residual editing, greatly reduces the calculation time consumption, and meets the high efficiency requirement;

[0075] (2) In the clock difference synthesis process, the M estimation method is used to estimate the dynamic weighting of observation data, realize low sampling absolute clock difference and high sampling relative clock difference, and improve the robustness of the algorithm;

[0076] (3) In the clock difference synthesis process, the high sampling relative clock difference and the low sampling absolute clock difference are synthesized through time series, realize the conversion of low sampling clock difference to high sampling clock difference, and ensure the continuity and precision of the time series;

[0077] (4) The application fuses the M estimation robustness technology and the clock difference encryption strategy, ensures high precision, significantly improves the calculation efficiency and robustness, and provides an efficient and reliable technical solution for high-frequency satellite navigation application.

[0078] The method scheme of the application will be further described in detail below by means of the accompanying drawings and examples. BRIEF DESCRIPTION OF DRAWINGS

[0079] Figure 1 The flow chart of the navigation satellite clock difference fast encryption method based on M estimation of the application;

[0080] Figure 2 The schematic diagram of residual difference comparison result of the navigation satellite clock difference fast encryption method based on M estimation and LS method in 30 second clock difference sampling task of the application. DETAILED DESCRIPTION

[0081] The method scheme of the application will be further described in detail below by means of the accompanying drawings and examples.

[0082] Unless otherwise defined, the method terms or scientific terms used in the application should be the usual meanings understood by those skilled in the art in the field of the application.

[0083] The terms such as "comprising" or "including" or similar words in the present invention mean that the elements before the word encompass the elements listed after the word, and do not exclude the possibility of also encompassing other elements. The terms "in", "out", "up", "down", etc. indicate the orientation or positional relationship based on the orientation or positional relationship shown in the drawings, and are only for the convenience of describing the present invention and simplifying the description, and do not indicate or imply that the indicated device or element must have a particular orientation, be constructed and operated in a particular orientation, and therefore cannot be understood as a limitation on the present invention. When the absolute position of the described object changes, the relative positional relationship may also change accordingly. In the present invention, unless otherwise explicitly specified and limited, the term "attached" and other terms should be understood broadly, for example, it can be fixedly connected, or it can be detachably connected, or it can be integrated; it can be directly connected, or it can be indirectly connected through an intermediate medium; it can be the internal connection of two elements or the interaction relationship between two elements. For those skilled in the art, the specific meaning of the above terms in the present invention can be understood according to the specific circumstances.

[0084] Embodiment

[0085] As Figure 1 shown, the present application provides a navigation satellite clock error fast encryption method based on M estimation, comprising the following steps:

[0086] S1: consistency check on data, reading navigation satellite orbit and clock error data, station coordinates and receiver clock error data, GNSS auxiliary data, data initialization and data preprocessing;

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

[0088] S2: constructing epoch difference observation equation, epoch difference ionosphere-free combined observation;

[0089] Step S2 specifically includes the following steps:

[0090] S21: calculating pseudorange observation values pseudorange observation values Specifically represented as:

[0091]

[0092] Wherein, r represents the station, 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 satellite clock error, represents ionospheric delay, represents tropospheric delay, ε P represents pseudo-range observation noise.

[0093] S22: between adjacent epochs t and t+1, the pseudo-range observation value is differenced, and the pseudo-range observation value difference result is specifically represented as:

[0094]

[0095] wherein, represents the change amount of geometric distance caused by relative motion of the satellite and the receiver, Δδt r represents the change amount of receiver clock error, since the receiver clock error is linearly changed between adjacent epochs, Δδt r can be approximately eliminated, Δδt s represents the change amount of satellite clock error, represents the change amount of ionospheric delay, represents the change amount of tropospheric delay, Δε P represents the change amount of pseudo-range observation noise.

[0096] S23: the carrier phase observation value is calculated, and the carrier phase observation value is specifically represented as:

[0097]

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

[0099] S24: between adjacent epochs t and t+1, the carrier phase observation value is differenced, and the carrier phase observation value difference result is specifically represented as:

[0100]

[0101] wherein, Δε Φ represents the change amount of carrier phase observation noise, since the integer ambiguity is not changed in time, thus the integer ambiguity can be eliminated, which makes the carrier phase observation more stable than the pseudo-range observation, and becomes the main data source for high-precision clock error estimation, and in the embodiment, the carrier phase epoch difference observation is also used for clock error encryption.

[0102] S3: Perform robust estimation of relative clock difference 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 posterior residual r i refers to the residual obtained by comparing the calculated value of the observation value using the estimated parameters with the original observation value after the parameter estimation is completed, in M-estimation, the posterior residual r i is used for weighted adjustment to reduce the influence of abnormal observation values, the posterior residual r i is specifically expressed as:

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

[0108] Wherein, O i represents the observation value, X represents the estimated parameter, and F(X) represents the estimated value.

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

[0110]

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

[0112] Step 3: Calculate the robust weight w i through the M-estimation weight function, common M-estimation weight functions include Huber loss function and Welsch loss function, the robust weight w i is specifically expressed as:

[0113]

[0114] S32: Perform relative clock difference parameter estimation through M-estimation and iterative weighted least squares method, dynamically adjust the weight of observation data to ensure that the estimated clock difference result is more robust.

[0115] Step S32 specifically includes the following steps:

[0116] Step one: M estimation is to weaken the influence of outliers on the estimation results by assigning different weights to the observed data. In order to solve the parameter estimation problem, the embodiment combines M estimation and iterative weighted least squares method to estimate clock error, and constructs a linearized model of observation equation, which is specifically represented as:

[0117] y = Ax + ε

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

[0119] Wherein, 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 zero mean normal distribution.

[0120] Step two: the least squares method is used to minimize the residual sum of squares to solve the parameters, and the parameter solving method is specifically represented as:

[0121]

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

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

[0124] Step three: when there are outliers in the observed data, the LS method is sensitive to them, so M estimation is introduced to reduce the influence of outliers. By introducing the loss function ρ(v i ) through M estimation, the loss function ρ(v i ) is replaced by the square loss A new parameter solving method is obtained, and the new parameter solving method is specifically represented as:

[0125]

[0126] Wherein, the loss function ρ(v i ) is set as the loss function ρ W (v) of Welsch M estimation, and the loss function ρ W (v) of Welsch M estimation is specifically represented as:

[0127]

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

[0129] The objective of M-estimation is to minimize the sum of the loss function above, thus reducing the influence of outliers while maintaining the stability of the estimation. The core idea of M-estimation is to assign different weights to different observations, thus weakening the influence of outliers.

[0130] Step four: weight function w of Welsch M-estimation W (v) Assign weights to observations, weight function w of Welsch M-estimation W (v) Specifically represented as:

[0131]

[0132] For observations with large outliers, the weight tends to zero, reducing the influence on parameter estimation.

[0133] Step five: solve parameters based on M-estimation by iterative weighted least squares method, the solving process is specifically represented as:

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

[0135]

[0136] Where, N k represents the observation noise covariance matrix at the current time, N c,k represents the carrier phase observation noise covariance matrix, N P,k represents the pseudo-range observation noise covariance matrix, A c represents the design matrix of carrier phase observation, A P represents the design matrix of pseudo-range phase observation, W c,k represents the carrier phase residual weight matrix, defined as the diagonal matrix after weighting the residual, W P,k represents the pseudo-range observation residual weight matrix, represents the normalized carrier phase residual of the i-th satellite, represents the normalized carrier phase residual of the i-th satellite, represents the estimated value of pseudo-range observation noise variance, represents the estimated value of carrier phase observation noise variance, Δy c represents the residual vector of carrier phase observation, specifically set as the deviation between the observation value and the fitting value, Δy P represents the residual vector of pseudo-range observation, represents the estimated parameter vector at the next time, v c,k+1 represents the carrier phase observation correction value at the next time, v P,k+1 represents the pseudo-range observation correction value at the next time, denotes a weighting function based on normalized residuals, commonly used in robust estimation, denotes a weighting function of pseudorange observations, denotes the trace of the noise covariance matrix, used to calculate the normalization factor of the weighted residuals.

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

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

[0139]

[0140] When the relative change of the variance estimate value of the two consecutive iterations is less than 1%, it means that the noise variance has stabilized and will not be updated further.

[0141] In step five, the convergence condition of the M-estimation weight is specifically expressed as:

[0142]

[0143] wherein, denotes the M-estimation weight for the pseudorange observation value at the kth iteration, denotes the M-estimation weight for the carrier phase observation value at the kth iteration.

[0144] When the change of the M-estimation weight of all observations is less than 0.01, it means that the weight has stabilized and will not be updated further.

[0145] S4: Obtain high-sampling relative clock difference;

[0146] S5: Integrate the low-sampling absolute clock difference and the high-sampling relative clock difference, and output the high-sampling clock difference. By integrating the low-sampling clock difference and the high-sampling clock difference, the low-sampling clock difference is encrypted to the high-sampling clock difference, for example, the 300-second interval low-sampling clock difference is encrypted to the 30-second or 5-second high-sampling clock difference. In the clock encryption strategy, the low-sampling absolute clock difference is provided based on the carrier phase observation value, and the clock change value, i.e. the high-sampling relative clock difference, is provided based on the clock estimation process. Combining the absolute clock difference information and the relative clock difference information, a high-precision clock difference sequence is obtained. The low-sampling clock difference and the high-sampling clock difference are fused to ensure the continuity and accuracy of the time sequence, to meet the application requirements of real-time kinematic positioning and low-orbit satellite precise orbit determination, and to provide more stable and efficient clock difference information support for high-precision navigation applications.

[0147] Step S5 specifically comprises the following steps:

[0148] S51: Obtain a clock error change equation through phase observation, and the clock error change equation is specifically expressed as:

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

[0150] Wherein, Δδ(t i+1 ,t i ) represents a clock error change amount, δ(t i+1 ) represents a satellite clock error at time t i+1 , and δ(t i ) represents a satellite clock error at time t i .

[0151] S52: Establish a linear observation equation set, and fuse the low-sampling absolute clock error δ fix (t) and the high-sampling relative clock error δ(t), and the linear observation equation set is specifically expressed as:

[0152]

[0153] Wherein, δ fix (t n ) represents an absolute clock error obtained by low-sampling rate data at time t n , and δ(t n ) represents a satellite clock error at time t n .

[0154] S53: Solve the linear observation equation set through the M estimation and the iterative weighted least square method.

[0155] Based on the time sequence characteristic, the clock error parameter dimension is optimized, the calculation efficiency is improved, the observation data weight is dynamically adjusted by combining the robust estimation method, and therefore the precision and reliability of the clock error information fusion are further improved.

[0156] The embodiment aims at the key technical problem of fast solving of navigation satellite clock error, carries out a comparative study on the algorithm of high-sampling clock error estimation, and carries out clock error solving on four satellite navigation systems of GPS, Galileo, GLONASS and BDS based on observation data of 100 global monitoring stations, sets the sampling navigation satellite to 113, sets the sampling rate to 30 seconds, uses a high-performance computing server with an Intel Xeon Gold 6326 processor as the hardware platform, sets the main frequency of the Intel Xeon Gold 6326 processor to 2.90 GHz, and sets the computing core to 28.

[0157] Under the same experimental conditions, the traditional least squares method needs to go through two stages of parameter estimation and residual editing, and the time consumption is 2558 seconds and 276 seconds respectively, and the total calculation time is 2836 seconds. However, the navigation satellite clock error fast encryption method based on M estimation proposed in the embodiment can dynamically adjust the observation weight by introducing M estimation, does not need to edit the residual multiple times, and realizes fast convergence combined with the iterative weighted least squares method, so that the calculation time consumption is greatly reduced, and the calculation time is successfully shortened to 276 seconds.

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

[0159] As shown in Figure 2 The precision of the embodiment method is better than that of the traditional LS method, and the calculation efficiency is improved by 90%, and the operation speed is improved by about 10 times, which realizes a magnitude breakthrough. This breakthrough is mainly due to the dynamic adjustment of the observation weight by M estimation, which effectively suppresses the influence of abnormal values and greatly improves the calculation efficiency. Compared with the traditional LS method, the residual editing does not need to be iterated multiple times, and the iterative weighted least squares method is combined to realize fast convergence, so that the calculation time consumption is greatly reduced, and the real-time requirement is met.

[0160] The experimental results show that the navigation satellite clock error fast encryption method based on M estimation proposed in the embodiment not only meets the high-precision GNSS data processing requirement, but also has excellent calculation efficiency and robustness characteristics, and is especially suitable for real-time precise point positioning, star-based enhancement system and other navigation application scenarios with strict time requirement, thereby providing reliable technical support for efficient data processing of the new generation of navigation system.

[0161] Therefore, the navigation satellite clock error fast encryption method based on M estimation can effectively process abnormal values and noise in the observation data, improve the calculation efficiency, has strong robustness, can reduce the iteration times while maintaining high precision, thereby greatly improving the timeliness, and meeting the demand of real-time navigation and positioning application.

[0162] Finally, it should be noted that: the above embodiments are only used to illustrate the method scheme of the present application, but not to limit it. Although the present application has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the method scheme of the present application can still be modified or equivalently replaced, and these modifications or equivalent replacements cannot make the modified method scheme deviate from the spirit and scope of the method scheme of the present application.

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 data consistency, read the navigation satellite orbit and clock error data, station coordinates and receiver clock error data, GNSS assistance data, and perform data initialization and data preprocessing; S2: Construct the inter-epoch differential observation equation; S3: Estimation of relative clock error parameters; Step S3 specifically includes the following steps: S31: Based on the posterior residual Perform M estimation weighting; S32: Estimation of relative clock error parameters using M estimation and iterative weighted least squares method; 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; 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: ； in, represents the change in clock error, Indicates time The satellite clock error, Indicates time Satellite clock error; S52: Establish a linear observation equation system to convert the low-sampled absolute clock error and high sampling relative clock difference For fusion, the linear observation equations are specifically expressed as: ； in, Indicates at time The absolute clock error obtained from low sampling rate data, Indicates time Satellite clock error; S53: Solve the linear observation equations by M estimation and iterative weighted least squares method.

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 assistance data includes Earth rotation parameters and tropospheric parameters. The data preprocessing process includes cycle slip detection and gross error removal of GNSS observation values. 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 observation value Specifically expressed as: ； in, Indicates the site, Indicates satellite, represents the geometric distance, represents the speed of light, represents the receiver clock error, represents the satellite clock error, represents the ionospheric delay, represents the tropospheric delay, represents the pseudorange observation noise; S22: In adjacent epochs and Between the pseudo-range observations Perform differential analysis, pseudorange observation value differential results Specifically expressed as: ； ； ； ； ； in, It represents the change in geometric distance caused by the relative motion between the satellite and the receiver. Indicates the change in receiver clock error, represents the change in satellite clock error, represents the change in ionospheric delay, represents the change in tropospheric delay, represents the variation of pseudorange observation noise; S23: Calculate carrier phase observation value , carrier phase observation value Specifically expressed as: ； in, represents the carrier wavelength, represents the integer ambiguity, represents the carrier phase observation noise; S24: In adjacent epochs and Between, the carrier phase observation value Perform differential analysis, carrier phase observation value differential result Specifically expressed as: ； ； in, Indicates the variation of 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 S31 specifically includes the following steps: Step 1: Calculate the post-test residuals , post-test residual Specifically expressed as: ； in, represents the observed value, represents the estimated parameters, Indicates an estimated value; Step 2: Perform weighted adjustment through the M-estimation objective function. The M-estimation objective function is specifically expressed as: ； in, represents the M estimation loss function, represents the standard deviation of observation noise; Step 3: Calculate the robustness weight through the M estimation weight function , robustness weight Specifically expressed as: 。 5. The method for fast encryption of navigation satellite clock errors based on M estimation according to claim 1, 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: ； ； in, represents the observation vector, represents the design matrix, represents the parameter vector to be estimated, represents the observation noise vector, 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: ； ； in, represents the estimated parameter vector, Represents the residual, specifically the deviation between the observed value and the fitted value, represents the observed value; Step 3: Introduce loss function through M estimation , so that the loss function Instead of square loss , a new parameter solving method is obtained, which is specifically expressed as: ； Among them, the loss function Set to the loss function estimated by WelschM , WelschM estimates the loss function Specifically expressed as: ； in, represents the residual, represents the scale parameter; Step 4: Weight function estimated by WelschM Assign weights to observations, the weight function estimated by WelschM Specifically expressed as: ； Step 5: Solve the parameters based on M estimation through iterative weighted least squares method. The solution process is specifically expressed as: ； in, represents the observation noise covariance matrix at the current moment, represents the carrier phase observation noise covariance matrix, represents the pseudorange observation noise covariance matrix, represents the design matrix of carrier phase observation, represents the design matrix of pseudorange phase observations, Represents the carrier phase residual weight matrix, which is defined as the diagonal matrix after weighting the residual. represents the pseudorange observation residual weight matrix, Represents the normalized The carrier phase residual of the satellites, Represents the normalized The pseudorange observation residuals of the satellites, represents the estimated value of the pseudorange observation noise variance, represents the estimated value of the carrier phase observation noise variance, represents the residual vector of the carrier phase observation, which is specifically set to the deviation between the observed value and the fitted value, represents the residual vector of pseudorange observations, represents the parameter vector estimated at the next moment, Indicates the carrier phase observation correction value at the next moment, represents the pseudorange observation correction value at the next moment, represents the weighting function based on the normalized carrier phase residual, represents the weighting function of the pseudorange observation residual, represents the trace of the carrier phase noise covariance matrix.

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

Images