Real-time cycle slip detection method based on ionosphere change trend constraint

Abstract

The invention belongs to the field of global navigation satellite system (GNSS) data processing, and particularly relates to a real-time cycle slip detection method based on ionosphere variation trend constraint, which comprises the following steps of: firstly, fitting an ionosphere variation trend by using all observation satellite data at the current moment, and selecting a non-geometric combination observation value according to a sliding window; then ionosphere variation trend constraint is applied on the basis of a polynomial function, geomagnetic latitude is used as an independent variable to fit a polynomial coefficient, the ionosphere variation of each satellite observed value at the current moment is forecasted, and then the ionosphere variation and the non-geometric combination are differentiated to obtain ambiguity variation; and finally, a detection threshold is determined through the adaptive expansion factor changing along with latitude, and the ambiguity variation is detected to judge whether cycle slip occurs or not. The method is convenient to calculate, fully considers the change trend of the ionosphere, is suitable for real-time cycle slip detection application, and is more suitable for cycle slip detection application in the ionosphere disturbance period.

Description

Technical Field

[0001] This invention relates to the field of GNSS data processing for satellite positioning and navigation, and in particular to a real-time cycle slip detection method based on ionospheric change trend constraints for applications such as navigation and positioning. Background Technology

[0002] With the arrival of the 25th solar cycle, increasingly intense ionospheric disturbances significantly impact GNSS positioning performance, especially in low-latitude regions. In GNSS navigation and positioning, ionospheric delay is one of the main sources of error in GNSS signals. As a crucial indicator in cycle slip detection, geometry-free combination (GFS) observations, which include ionospheric delay, will fluctuate dramatically during ionospheric disturbances, causing cycle slip detection errors and over-initialization of ambiguity parameters, ultimately leading to a decrease in GNSS positioning accuracy.

[0003] Currently, commonly used cycle slip detection methods include the TurboEdit algorithm, STPIR method, FBMWA method, and empirical threshold method. The TurboEdit algorithm is widely used for cycle slip detection because it is unaffected by clock instability, selectivity availability, receiver-satellite kinematics, and tropospheric conditions. However, the TurboEdit algorithm uses pseudorange observations in both Melbourne–Wübbena (MW) and Geometry-Free (GF) combinations, and the large pseudorange noise limits its cycle slip sensitivity. Based on the TurboEdit algorithm, FBMWA uses forward and reverse sliding windows to smooth the MW combination, reducing MW combination noise and improving cycle slip sensitivity. However, when two frequencies produce cycle slips of the same magnitude, the MW combination cannot detect the cycle slip. To address this, the STPIR algorithm uses the second-order time-domain difference of the carrier phase GF, which not only avoids the influence of pseudorange noise but also weakens the effect of ionospheric delay. Although the above methods can achieve efficient cycle slip detection, because they use information after the detection time, they are only suitable for post-event cycle slip detection and cannot be applied to real-time scenarios. The empirical threshold method determines the cycle slip detection threshold function as a function of satellite elevation angle based on measured data. It is easy to calculate, highly efficient, and can be used for both real-time and post-detection purposes. However, this algorithm frequently produces false alarms during periods of ionospheric calm, leading to decreased positioning accuracy. Therefore, developing a real-time cycle slip detection method applicable to various periods of ionospheric activity is crucial for GNSS navigation and positioning. Summary of the Invention

[0004] A real-time cycle slip detection method based on ionospheric variation trend constraints improves the accuracy of real-time cycle slip detection, solves the problem of frequent false detections of cycle slips during ionospheric disturbances, and can handle real-time cycle slip detection under various levels of ionospheric activity. This invention has certain application value for real-time cycle slip detection.

[0005] To achieve the above objectives, this invention provides a real-time cycle slip detection method based on ionospheric variation trend constraints, comprising: Step S1: Obtain all satellite observations at the current moment and construct a geometrically independent combination; Step S2: Obtain the information from the previous time step and calculate the first and second differences without geometric combinations; Step S3: Fit the ionospheric variation trend using a geometrically non-geometrically combined quadratic difference; Step S4: For each satellite, select the observation values ​​within the sliding window, calculate the geometric combination and its first difference, apply the ionospheric change trend constraint on the basis of the polynomial function, fit the polynomial coefficients and predict the ionospheric change in the current epoch. Step S5: Calculate the difference between the predicted ionospheric change and the first difference without geometric combination at the current time, and take the absolute value to obtain the ambiguity change. Step S6: Calculate the unit weighted mean error of the fitting residual, calculate the expansion factor based on the geomagnetic latitude of the target, and then calculate the detection threshold; Step S7: Compare the change in ambiguity with the detection threshold to determine whether the change in ambiguity is greater than the detection threshold; If so, determine that a cycle slip has occurred in the current satellite observation value at the current moment; If not, it is determined that no cycle slip has occurred in the current satellite observation value at the current moment; This exploration is complete.

[0006] This invention utilizes all satellite observations to fit the ionospheric variation trend. By applying constraints on the ionospheric variation trend to a polynomial function based on geomagnetic latitude, it improves the ionospheric forecasting performance, accurately determines the magnitude of ambiguity changes, and realizes real-time cycle slip detection under various ionospheric activity levels (calm, active, and disturbed). This enhances the real-time cycle slip detection performance and broadens the application scenarios of real-time cycle slip detection in various space weather environments.

[0007] Beneficial effects Compared with the prior art, the present invention has the following advantages: Existing cycle slip detection methods are generally affected by ionospheric activity or are not specifically adapted to different ionospheric states. During ionospheric disturbances, cycle slip detection performance is low, leading to a decrease in GNSS navigation and positioning accuracy. In contrast, this invention is designed to address the characteristic that ionospheric activity varies with geomagnetic latitude, and can adapt to various ionospheric states, including calm, active, and disturbed states, making it applicable to a wider range of scenarios.

[0008] Existing cycle slip detection methods have significant limitations and struggle to meet the real-time dynamic requirements of GNSS. Methods such as TurboEdit, STPIR, and FBMWA rely on retrospective verification using observational information after the detection time, essentially representing post-processing. They cannot perform cycle slip judgment simultaneously with observational data generation, making them unsuitable for GNSS dynamic scenarios with extremely high real-time requirements, such as vehicle navigation and UAV positioning. While empirical threshold methods can support real-time detection, they suffer from frequent false detections during periods of ionospheric activity. This invention focuses on utilizing only geometrically uncombined observations from all satellites at the current moment. By selecting data through a sliding window, fitting ionospheric change trends, and predicting ionospheric changes, the calculation process requires no subsequent data and can complete cycle slip detection immediately after observational data acquisition. This satisfies the high-precision requirements of post-processing data while adapting to real-time dynamic positioning scenarios, filling the gap in existing technologies for "efficient real-time detection" and better meeting the practical application needs of GNSS navigation and positioning.

[0009] Existing cycle slip detection methods are susceptible to external interference. For example, TurboEdit, relying on pseudorange observations, is weakened by significant pseudorange noise, reducing its cycle slip sensitivity. While the empirical threshold method is simple to calculate, it only sets a fixed threshold based on the satellite elevation angle, neglecting ionospheric variations. During periods of high ionospheric activity, the threshold often deviates significantly from the actual ionospheric delay, frequently leading to misjudgments of cycles even when none exist. This invention improves accuracy and stability through multiple technological designs: First, it optimizes observation processing by constraining ionospheric variation trends. Based on polynomial fitting, it combines the correlation between geomagnetic latitude and ionospheric activity, making the prediction of ionospheric variations more closely match the actual atmospheric environment and reducing detection errors caused by ionospheric model biases. Second, it uses an adaptive expansion factor that varies with latitude to set the threshold. This factor dynamically adjusts the threshold range according to the ionospheric activity characteristics at different geomagnetic latitudes. This invention can accurately identify cycle slips while avoiding invalid corrections, effectively ensuring the accuracy and stability of GNSS positioning calculations. Attached Figure Description

[0010] Figure 1 This is a schematic diagram of the real-time cycle slip detection method according to an embodiment of the present invention; Figure 2 This is a schematic diagram of the specific process of step S1 in an embodiment of the present invention; Figure 3 This is a schematic diagram of the specific process of step S3 in an embodiment of the present invention; Figure 4 This is a schematic diagram illustrating the specific process of step S4 in an embodiment of the present invention; Figure 5 This is a schematic diagram of the specific process of step S6 in an embodiment of the present invention; Figure 6This is a distribution chart of the daily average weekly false alarms in 2022, based on an embodiment of the present invention. Detailed Implementation

[0011] The specific embodiments of the present invention will now be described in more detail with reference to the accompanying drawings. The advantages and features of the present invention will become clearer from the following description and claims. It should be noted that the drawings are all in a very simplified form and use non-precise proportions, and are only used to facilitate and clarify the illustration of the embodiments of the present invention.

[0012] To ensure accurate real-time cycle slip detection, this invention proposes a real-time cycle slip detection method based on ionospheric variation trend constraints. The basic idea is as follows: First, using all observed satellite data at the current moment, the ionospheric variation trend is fitted using a weighted average method based on satellite elevation angles. Observations without geometric combinations are selected according to a sliding window. Then, an observation equation based on a polynomial function is constructed, and ionospheric variation trend constraints are applied. Polynomial coefficients are fitted with geomagnetic latitude as the independent variable to predict the ionospheric variation of each satellite's observation at the current moment. Subsequently, the ionospheric variation is differiated from the geometrically unconstrained variation to obtain the ambiguity variation. Finally, a detection threshold is determined using an adaptive dilation factor that varies with latitude, and the ambiguity variation is detected to determine whether a cycle slip has occurred. If the ambiguity variation is less than the detection threshold, it is determined that the carrier phase observation of that satellite at the current moment has not experienced a cycle slip; otherwise, it is determined that the carrier phase observation of that satellite at the current moment has experienced a cycle slip.

[0013] refer to Figure 1 In a preferred embodiment of the present invention, a real-time cycle slip detection method based on ionospheric change trend constraints includes the following steps: Step S1: Obtain all satellite observations at the current moment and construct a geometrically independent combination; For details, please refer to Figure 2 Step S1 includes: S11: Perform data preprocessing on satellite observation data; The satellite observation data includes carrier phase observations; Data preprocessing includes, but is not limited to, setting the satellite cutoff elevation angle and iterative calculation of satellite signal transmission time; this invention does not impose any limitations on these aspects.

[0014] Preferably, the satellite signal transmission time iteration formula used is: Among them, subscript r , s These represent the receiver and the satellite, respectively. and When the clock faces of the receiver and satellite are displayed respectively,c Represents the speed of light. P Indicates the distance from its satellite to the ground. It is the time when the satellite signal is transmitted. Indicates satellite clock bias.

[0015] S12: Construct observations without geometric combinations; Geometric-free observations are constructed using carrier phase observations, as shown in the formula. in, For satellite s The geometric combination without and Satellites s Carrier phase observations at the first and second frequencies; Step S2: Obtain the information on the geometric combination without geometric combination in the previous time step, and calculate the first and second differences without geometric combination; Step S3: Fit the ionospheric variation trend using a geometrically non-geometrically combined quadratic difference; For details, please refer to Figure 3 Step S3 includes: S31: Calculate weights using satellite elevation angles; The weight calculation formula is as follows: in, For satellite s Weights without geometric combination quadratic difference For satellite s The corresponding elevation angle; S32: Fitting the trend of ionospheric changes using the weighted average method; The fitting formula for the ionosphere variation trend is: in, This represents the trend of ionospheric changes. The number of satellites observed at the current moment. For satellite s Weights without geometric combination quadratic difference For satellite s The quadratic difference without geometric combinations; Step S4: For each satellite, select the observation values ​​within the sliding window, calculate the geometric combination and its first difference, apply the ionospheric change trend constraint on the basis of the polynomial function, fit the polynomial coefficients and predict the ionospheric change in the current epoch. The sliding window contains data from a period preceding the ionospheric change trend observation points in steps S1-S3, and its length is [missing information]. mThis invention utilizes geometrically unobserved values ​​within a sliding window for fitting and forecasting to obtain the predicted value of ionospheric changes at the current moment.

[0016] For details, please refer to Figure 4 Step S4 includes: S41: Calculate the geometrically undefined combination and its first difference within the sliding window; The formula for calculating the geometric combination-free combination is: in, For satellite s The geometric combination without and Satellites s Carrier phase observations at the first and second frequencies; The formula for calculating the first difference without geometric combination is: in, For satellite s The first difference of geometric combinations, For satellite s At the present moment Observations without geometric combination For satellite s In the previous moment Observations without geometric combinations; S42: Construct the observation equations for the polynomial function; The formula for calculating the polynomial function is: in, It is a polynomial function with a first-order difference and no geometric combination. As the independent variable, this invention uses geomagnetic latitude as the independent variable. Let be the order of the polynomial, ( , , , , ) represents the polynomial coefficients; For satellites s In the sliding window m The observation equation for the first difference without geometric combinations within the range is: in, The observation matrix, To design the matrix, Let be the parameter matrix to be estimated. Here is the error matrix. For satellite s In the sliding windowm The Middle j The difference between the first and second ionosphere combinations at a given time. For sliding windows m The Middle j Moment x of n Power of 1 For satellite s In the sliding window m The Middle j Observational noise at any given moment; S43: Apply constraints on the trend of ionospheric change; The constraint formula for the ionospheric change trend is: in, It is a polynomial function with no geometric combination of first-order difference forecast values. For satellite s At any moment The first difference of observations without geometric combination. This represents the trend of ionospheric changes; Combining the polynomial function observation equation (6), the observation equation with constraints on ionospheric variation trends is as follows: in, The observation matrix, To design the matrix, Let be the parameter matrix to be estimated. The error matrix; This is the constraint matrix. The design matrix for the constraints. The error matrix represents the constraints. The above equation can be written in matrix form: in, For observational equations constrained by ionospheric variation trends, This is a design matrix constrained by ionospheric variation trends. Let be the parameter matrix to be estimated. This is the error matrix constrained by the trend of ionospheric change; S44: Parameter fitting and prediction of ionospheric changes; The parameter fitting method is the least squares method, and the calculation formula is: in, The least squares solution. For observational equations constrained by ionospheric variation trends, This is a design matrix constrained by ionospheric variation trends. The error matrix is ​​constrained by the trend of ionospheric change. For the power formation; Power Array Specifically, it is expressed as follows: in, For the power formation, Given a diagonal matrix function, construct a square matrix whose non-diagonal elements are all 0; For sliding windows m The weights of the first difference of geometric combinations are not included. ,in For the first in the sliding window k Elevation angle of epoch; The weights of the constraints are specifically represented as follows: in, The variance of the first difference without geometric combinations. The variance of the ionospheric variation trend. The number of satellites participating in the trend of ionospheric changes. For satellite s Weights without geometric combination quadratic difference; Using the least-squares solution of the polynomial coefficients, the formula for predicting the ionospheric change is as follows: in, For the predicted changes in the ionosphere, The least squares solution for the polynomial coefficients. This represents the geomagnetic latitude at the next moment (the current moment) within the sliding window. The order of the polynomial; Step S5: Calculate the difference between the predicted ionospheric change and the first difference without geometric combination at the current time, and take the absolute value to obtain the ambiguity change. The formula for calculating the change in ambiguity is: in, For the predicted changes in the ionosphere, Let be the first difference of the geometric combination at the current time (i.e., the next time step of the sliding window). It is an absolute value function; Step S6: Calculate the unit weighted mean error of the fitting residual, calculate the expansion factor based on the geomagnetic latitude of the target, and then calculate the cycle slip detection threshold; For details, please refer to Figure 5 Step S6 includes: S61: Calculate the unit weight error; The formula for calculating the unit weight mean square error is: in, The error is the unit weight. For the residual vector, For the observation weight matrix, m The length of the sliding window. n The order of the polynomial fitting is given. S62: Calculate the expansion factor based on geomagnetic latitude; The formula for calculating the expansion factor is: in, It is the expansion factor. For the absolute value function, It is a cosine function. This represents the current geomagnetic latitude. S63: Cycle slip detection threshold determination; The formula for calculating the cycle slip detection threshold is: in, The cycle slip detection threshold, It is the expansion factor. This represents the unit weight error.

[0017] Step S7: Compare the change in ambiguity with the detection threshold to determine whether the change in ambiguity is greater than the detection threshold; If so, determine that a cycle slip has occurred in the current satellite observation value at the current moment; If not, it is determined that no cycle slip has occurred in the current satellite observation value at the current moment; This exploration is complete.

[0018] In summary, this invention utilizes all satellite observations to fit the ionospheric variation trend. By applying constraints on the ionospheric variation trend to a polynomial function based on geomagnetic latitude, it improves the ionospheric forecasting performance, accurately determines the magnitude of ambiguity changes, and achieves real-time cycle slip detection under various ionospheric activity levels (calm, active, and disturbed). This enhances the real-time cycle slip detection performance and broadens the application scenarios of real-time cycle slip detection in various space weather environments.

[0019] Specifically, compared with the prior art, the present invention has the following advantages: Compared with existing polynomial fitting forecasts, this invention can accurately predict the amount of ionospheric change under conditions of drastic ionospheric variations, making the calculation of ambiguity changes more accurate. Specific data comparisons are shown in Table 1: Table 1 Comparison of prediction errors for ionospheric changes During periods of ionospheric calm and activity, the prediction error of ionospheric changes in this invention is superior to that of existing polynomial fitting prediction methods. During periods of ionospheric calm, the prediction error of ionospheric changes in this invention is 0.0069m, which is 9% lower than that of polynomial fitting prediction methods. During periods of ionospheric activity, the prediction error of ionospheric changes in this invention is less than 0.1m, which is 28% lower than that of polynomial fitting prediction methods.

[0020] Compared with existing polynomial fitting predictions, this invention is designed specifically to address the variation of ionospheric activity with geomagnetic latitude, effectively improving the false detection rate of cycle slip predictions in ionospherically active regions (low latitudes and the Earth's poles). Specific results are as follows... Figure 6 As shown.

[0021] As shown in the figure, the empirical threshold method frequently misdetects cycle slips in low-latitude regions and at the Earth's poles, while the polynomial fitting method has fewer cycle slip misdetections at the Earth's poles, but still exhibits frequent cycle slip misdetections in low-latitude regions. However, the present invention achieves a cycle slip misdetection count close to 0 in both ionospherically active regions (low-latitude regions and the Earth's poles) and ionospherically calm regions (mid- to high-latitude regions).

[0022] The above description is merely a description of preferred embodiments of this application and is not intended to limit the scope of this application in any way. Any changes or modifications made by those skilled in the art based on the above-disclosed technical content should be considered as equivalent and valid embodiments and fall within the scope of protection of the technical solution of this application.

Claims

1. A real-time cycle slip detection method based on ionospheric variation trend constraints, characterized in that, Includes the following steps: Step S1: Obtain all satellite observations at the current moment and construct a geometrically independent combination; Step S2: Obtain the information from the previous time step and calculate the first and second differences without geometric combinations; Step S3: Fit the ionospheric variation trend using a geometrically non-geometrically combined quadratic difference; Step S4: For each satellite, select the observation values ​​within the sliding window, calculate the geometric combination and its first difference, apply the ionospheric change trend constraint on the basis of the polynomial function, fit the polynomial coefficients and predict the ionospheric change in the current epoch. Step S5: Calculate the difference between the predicted ionospheric change and the first difference without geometric combination at the current time, and take the absolute value to obtain the ambiguity change. Step S6: Calculate the unit weighted mean error of the fitting residual, calculate the expansion factor based on the geomagnetic latitude of the target, and then calculate the detection threshold; Step S7: Compare the change in ambiguity with the detection threshold to determine whether the change in ambiguity is greater than the detection threshold; If so, determine that a cycle slip has occurred in the current satellite observation value at the current moment; If not, it is determined that no cycle slip has occurred in the current satellite observation value at the current moment; This exploration is complete.

2. The real-time cycle slip detection method as described in claim 1, characterized in that, Step S1 includes: S11: Perform data preprocessing on satellite observation data; The satellite observation data includes carrier phase observations; Data preprocessing includes: setting the satellite cutoff elevation angle and iteratively calculating the satellite signal transmission time; S12: Construct observations without geometric combinations.

3. The real-time cycle slip detection method as described in claim 1, characterized in that, Step S3 includes: S31: Calculate weights using satellite elevation angles; The weight calculation formula is as follows: in, For satellite s Weights without geometric combination quadratic difference For satellite s The corresponding elevation angle; S32: Fitting the trend of ionospheric changes using the weighted average method; The fitting formula for the ionosphere variation trend is: in, This represents the trend of ionospheric changes. The number of satellites observed at the current moment. For satellite s Weights without geometric combination quadratic difference For satellite s The quadratic difference without geometric combination.

4. The real-time cycle slip detection method as described in claim 1, characterized in that, Step S4 includes: S41: Calculate the geometrically undefined combination and its first difference within the sliding window; The formula for calculating the geometric combination-free combination is: in, For satellite s The geometric combination without and Satellites s Carrier phase observations at the first and second frequencies; The formula for calculating the first difference without geometric combination is: in, For satellite s The first difference of geometric combinations, For satellite s At the present moment Observations without geometric combination For satellite s In the previous moment Observations without geometric combinations; S42: Construct the observation equations for the polynomial function; The formula for calculating the polynomial function is: in, It is a polynomial function with a first-order difference and no geometric combination. Using geomagnetic latitude as the independent variable, Let be the order of the polynomial, ( , , , , ) represents the polynomial coefficients; For satellites s In the sliding window m The observation equation for the first difference without geometric combinations within the range is: in, The observation matrix, To design the matrix, Let be the parameter matrix to be estimated. Here is the error matrix. For satellite s In the sliding window m The Middle j The difference between the first and second ionosphere combinations at a given time. For sliding windows m The Middle j Moment x of n Power of 1 For satellite s In the sliding window m The Middle j Observational noise at any given moment; S43: Apply constraints on the trend of ionospheric change; The constraint formula for the ionospheric change trend is: in, It is a polynomial function with no geometric combination of first-order difference forecast values. For satellite s At any moment The first difference of observations without geometric combination. This represents the trend of ionospheric changes; Combining the polynomial function observation equation (6), the observation equation with constraints on ionospheric variation trends is as follows: in, The observation matrix, To design the matrix, Let be the parameter matrix to be estimated. The error matrix; This is the constraint matrix. The design matrix for the constraints. The error matrix represents the constraints. The above equation can be written in matrix form: in, For observational equations constrained by ionospheric variation trends, This is a design matrix constrained by ionospheric variation trends. Let be the parameter matrix to be estimated. This is the error matrix constrained by the trend of ionospheric change; S44: Parameter fitting and prediction of ionospheric changes; The parameter fitting method is the least squares method, and the calculation formula is: in, The least squares solution. For observational equations constrained by ionospheric variation trends, This is a design matrix constrained by ionospheric variation trends. The error matrix is ​​constrained by the trend of ionospheric change. For the power formation; Power Array Specifically, it is expressed as follows: in, For the power formation, Given a diagonal matrix function, construct a square matrix whose non-diagonal elements are all 0; For sliding windows m The weights of the first difference of geometric combinations are not included. ,in For the first in the sliding window k Elevation angle of epoch; The weights of the constraints are specifically represented as follows: in, The variance of the first difference without geometric combinations. The variance of the ionospheric variation trend. The number of satellites participating in the trend of ionospheric changes. For satellite s Weights without geometric combination quadratic difference; Using the least-squares solution of the polynomial coefficients, the formula for predicting the ionospheric change is as follows: in, For the predicted changes in the ionosphere, The least squares solution for the polynomial coefficients. This represents the geomagnetic latitude at the next moment of the sliding window. Let be the order of the polynomial.

5. The real-time cycle slip detection method as described in claim 1, characterized in that, In step S5, the formula for calculating the change in ambiguity is: in, For the predicted changes in the ionosphere, For the first difference of the geometric combination at the current moment, This is the function for taking the absolute value.

6. The real-time cycle slip detection method as described in claim 1, characterized in that, Step S6 includes the following steps: S61: Calculate the unit weight error; The formula for calculating the unit weight mean square error is: in, The error is the unit weight. For the residual vector, For the observation weight matrix, m The length of the sliding window. n The order of the polynomial fitting is given. S62: Calculate the expansion factor based on geomagnetic latitude; The formula for calculating the expansion factor is: in, It is the expansion factor. For the absolute value function, It is a cosine function. This represents the current geomagnetic latitude. S63: Cycle slip detection threshold determination; The formula for calculating the cycle slip detection threshold is: in, The cycle slip detection threshold, It is the expansion factor. This represents the unit weight error.

Images