Single-frequency BDS cycle slip detection and repair method based on wavelet transform and model switching

Abstract

The invention relates to the technical field of Beidou navigation positioning data processing, and particularly discloses a single-frequency BDS cycle slip detection and restoration method based on wavelet transform and model switching, and the method comprises the steps: firstly reading single-frequency carrier phase observation data, and constructing a cycle slip detection amount; performing three-layer decomposition on the double-difference observed quantity by adopting a db4 wavelet to obtain coefficients D1, D2, D3 and A3; for the time-frequency characteristic of each coefficient, designing adaptive multi-threshold abnormal interval detection, and respectively matching CNN / ELM, BP / CNN, ELM / SVR and ARIMA prediction models for D1, D2, D3 and A3; through abnormal interval identification and recursive sliding window prediction, high-sensitivity detection and centimeter-level repair of cycle slip are realized. Experiments show that compared with an optimal single CNN model, the RMSE is improved by 17.57%, the RMSE is stable at the sampling rate of 1-10 s and on satellites of all orbits, and the method is suitable for high-precision positioning application of a low-cost single-frequency receiver.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of Beidou navigation and positioning data processing, and particularly relates to a single-frequency BDS cycle slip detection and repair method based on wavelet transform and model switching. BACKGROUND

[0002] Cycle-slip detection and repair of carrier phase observations is a key prerequisite for achieving centimeter-level or even millimeter-level GNSS positioning accuracy. For dual-frequency or multi-frequency receivers, cycle slips can be effectively distinguished and repaired by constructing geometrically independent, wide-lane or ionospheric independent combinations; however, in single-frequency receiver applications with limited cost, power consumption and volume, due to the lack of multi-frequency observations, ionospheric delay cannot be eliminated, wide-lane combinations cannot be established, and cycle slip detection is significantly more difficult. Although single-frequency receivers still occupy a considerable market share in mass navigation and Internet of Things scenarios, efficient cycle slip processing algorithms for this scenario are very scarce.

[0003] Traditional single-frequency cycle slip processing ideas can be roughly summarized as follows: polynomial fitting residual method, carrier phase sequence is modeled by a low-order polynomial, and the residual mutation is regarded as cycle slip. However, carrier phase usually exhibits strong non-stationary characteristics, especially when satellite motion is fast or signal is blocked, the fitting accuracy of the polynomial decreases sharply, and small cycle slips of less than 1 cycle are almost undetectable. High-order phase difference method, through epoch differencing to eliminate low-frequency trends step by step, but the differencing process amplifies observation noise synchronously, and high-order differencing will significantly deteriorate the signal-to-noise ratio, small amplitude cycle slips are easily covered by amplified noise. Kalman filter method, by means of state prediction-update framework, detects the residual between predicted value and observation value; this method relies on accurate dynamic model and noise statistical characteristics, but these priors are often difficult to obtain, resulting in uncontrollable false detection rate and missed detection rate. Pseudorange-phase combination method, which takes advantage of the complementarity of carrier and pseudorange observations, but the pseudorange noise itself is large, limiting the sensitivity to small cycle slips. Wavelet analysis method, from time-frequency domain to capture the mutation point, can locate the cycle slip occurrence epoch, but it is difficult to directly give the specific cycle slip size; if combined with machine learning, it can improve the detection rate, but it is still insufficient under the conditions of BDS single frequency, high speed dynamic or low sampling.

[0004] The Chinese invention patent with patent publication number CN109799521A discloses a BDS / GPS triple-difference combination cycle slip detection and repair method, which eliminates clock error, ambiguity and system bias by triple-difference combination of multi-system and multi-frequency observations, but it cannot be directly migrated to a single frequency scene because it relies on multi-frequency observations; the Chinese invention patent with patent publication number CN113189628A discloses a BDS multi-frequency observation cycle slip detection and repair method, which constructs an optimal linear combination for BDS multi-frequency observations, improves the cycle slip detection success rate and repair accuracy, and also takes multi-frequency signals as the premise; the Chinese invention patent with patent publication number CN115826003A discloses a BDS triple-frequency cycle slip detection method based on Doppler integral assistance, which enhances the detection reliability by combining ultra-wide lane combination and ionospheric residual difference, but it is not suitable for single frequency receivers; the Chinese invention patent with patent publication number CN107505642B discloses an INS-assisted real-time BDS single-frequency cycle slip detection method, which realizes real-time BDS single-frequency double-difference cycle slip detection using INS prediction information, solves the problem of insufficient prior model, but this scheme is sensitive to the quality of inertial navigation and coupling accuracy, and does not consider the fine repair of small cycle slips.

[0005] For cost-sensitive BDS single-frequency receivers that can only receive a single carrier frequency, lack of combined observations, and are significantly affected by noise and ionospheric delay, the existing cycle slip detection and repair techniques have the following shortcomings: low detection rate for small and short cycle slips, easily affected by noise amplification and model mismatch; insufficient cycle slip size estimation accuracy, difficult to meet the centimeter-level positioning requirements; poor adaptability to dynamic environments or obstructed working conditions, lack of adaptive adjustment mechanism for parameters and thresholds; excessive dependence on multi-frequency, inertial navigation and other external information, unable to be efficiently implemented on a pure single-frequency hardware platform. Therefore, there is an urgent need for a new wavelet-intelligent fusion cycle slip processing method that does not require multi-frequency observations, can comprehensively utilize signal characteristics in time-frequency domain, has adaptive threshold and model switching capability, and can detect with high sensitivity and repair with high accuracy, to solve the technical bottleneck of reliable discrimination and accurate compensation of single-frequency BDS carrier phase cycle slips. SUMMARY

[0006] The purpose of the present application is to address the defects of existing single-frequency BDS carrier phase cycle slips that are difficult to detect with high sensitivity and accurately repair, and to propose a single-frequency BDS cycle slip detection and repair method based on wavelet transform and model switching. By db4 three-layer wavelet decomposition of cycle slip detection in time-frequency domain, constructing an adaptive multi-threshold multi-scale anomaly interval detection model, and combining dynamic optimal model switching and recursive prediction strategies for D1, D2, D3 and A3 coefficients, high-reliable detection and centimeter-level repair accuracy of small and short cycle slips are achieved.

[0007] The application provides a single-frequency BDS cycle slip detection and repair method based on wavelet transform and model switching, and comprises the following steps:

[0008] S1: reading carrier data, wherein the carrier data comprises carrier phase observation values of a Beidou satellite system;

[0009] S2: constructing a cycle slip detection quantity, wherein the cycle slip detection quantity is formed by one or more of undifference, single difference, double difference and triple difference;

[0010] S3: selecting db4 wavelet as a wavelet base, and performing three-layer discrete wavelet decomposition on the cycle slip detection quantity in step S2 to obtain high-frequency detail coefficients D1, D2 and D3 and a low-frequency approximation coefficient A3;

[0011] S4: performing adaptive multi-threshold abnormal interval detection based on the response characteristics of multi-scale coefficients; wherein for the coefficients of wavelet decomposition at different levels in step S3, the calculation formula of the abnormal interval judgment criterion is:

[0012]

[0013] wherein, CS i represents a cycle slip detection result of the i-th coefficient, D1, D2 and D3 represent the high-frequency coefficients of the first, second and third levels of wavelet respectively, A3 represents the low-frequency coefficient of wavelet, μ i and σ i respectively represent the mean value and standard deviation of the absolute value of the i-th coefficient after robust statistical processing; a i represents the fluctuation amplitude in the abnormal interval Ω i , and is defined as the difference between the maximum value and the minimum value; represents the first-order difference of the approximation coefficient A3, and is defined as the difference between the A3 value at the current time t and the A3 value at the previous time t-1; f(D1) and g(D1) are adaptive functions based on the maximum value of the D1 coefficient, and are used to adjust the detection sensitivity of the A3 coefficient;

[0014] S5: when the cycle slip is detected, performing multi-prediction model switching repair on the time-frequency characteristics of D1, D2, D3 and A3;

[0015] S6: performing inverse wavelet reconstruction on the repaired D1, D2, D3 and A3 to obtain the repaired carrier phase observation value as the final output.

[0016] As a preferred technical scheme of the application, after reading the carrier data, step S1 further comprises pre-processing of the carrier phase observation value, and the pre-processing at least comprises any one or more of the following:

[0017] Epoch time alignment: synchronously correcting the observation epochs of different satellite channels to ensure consistent time labels;

[0018] Weighting processing: weights are assigned to observations according to satellite elevation angle and signal-to-noise ratio, and the weight coefficient is positively correlated with the sine value of the elevation angle and the signal-to-noise ratio;

[0019] Outlier rejection: the gross error in the observation sequence is identified and rejected by using the median absolute deviation or 3σ criterion.

[0020] As a preferred technical solution of the present application, the abnormal interval determination in step S4 further comprises: in the cycle slip detection result, 1 represents that cycle slip is detected, and 0 represents that cycle slip is not detected; and the robust statistical processing is specifically deleting the first 5% of abnormal values.

[0021] As a preferred technical solution of the present application, the adaptive multi-threshold multi-scale wavelet coefficient abnormal interval detection in step S4 comprises the following steps:

[0022] S4-1: data points with absolute values exceeding the statistical threshold T i are marked as an abnormal point set The specific calculation formula is: Wherein, C i (t) represents the i-th coefficient value at time t, and T i is the corresponding statistical threshold.

[0023] S4-2: the abnormal point set is connected and processed by using morphological closing operation Close to obtain Wherein, the structure element B selects a third-order structure element for D1 and A3 coefficients, and selects a second-order structure element for D2 and D3 coefficients.

[0024] S4-3: all connected regions in are identified, and the largest area is selected as a candidate abnormal interval; if the coefficients are A3 and D1 and there is a specified epoch point, the connected region closest to the epoch point is selected.

[0025] S4-4: the edge points are respectively and The interval boundary is extended outward by 3 epochs for D1 coefficients, and the interval boundary is extended outward by 2 epochs for D3 coefficients.

[0026] S4-5: the final abnormal interval is Ω i ={t:t start ≤t≤t end}; wherein the number of data points contained in the abnormal interval is k=t end -t start +1.

[0027] In the adaptive multi-threshold multi-scale wavelet coefficient anomaly interval detection in step S4, the low-frequency approximation coefficient A3 is first subjected to first-order differencing to enhance its sensitivity to abrupt changes. Based on the positive correlation between the high-frequency detail coefficient D1 and the cycle slip amplitude, the detection threshold and minimum fluctuation amplitude of A3 are adjusted in real time. When the anomaly interval corresponding to A3 deviates significantly from the anomaly intervals of the high-frequency coefficients D1, D2, and D3, the A3 coefficient is not considered. The processing using the low-frequency approximation coefficient A3 is divided into three cases:

[0028] When the maximum modulus of D1 is greater than 9:

[0029]

[0030] When the maximum modulus of D1 is less than or equal to 9 and greater than 3:

[0031]

[0032] When the maximum modulus of D1 is less than or equal to 3 and greater than 0:

[0033]

[0034] Where T represents the threshold, Δ m in represents the minimum fluctuation range, μ a3 and σ a3 These represent the mean and standard deviation of the A3 coefficient, respectively.

[0035] Based on the characteristics of different wavelet coefficients, the final formula for calculating cycle slip is:

[0036]

[0037] The specific conditions for determining the existence of cycle slips are: at least two coefficients detect cycle slips, or only the D1 coefficient detects a cycle slip.

[0038] As a preferred embodiment of the present invention, step S5, multi-prediction model switching repair, includes the following steps:

[0039] S5-1: Establish a candidate model set that includes at least convolutional neural networks, backpropagation neural networks, radial basis function neural networks, extreme learning machines, support vector regression, and autoregressive integral moving average;

[0040] S5-2: Divide the historical samples into training subsets and test subsets, and train the coefficient predictors of each model;

[0041] S5-4: Select the optimal model based on the evaluation index for recursive prediction and repair of the current coefficients within the abnormal interval;

[0042] S5-5: Adopting sliding window point-by-point extrapolation until the missing segment of length k is repaired.

[0043] The evaluation indexes of the above model at least include root mean square error (RMSE), mean absolute error (MAE), standard deviation (SD) of prediction error and variance (Var), and the calculation formula of each evaluation index is:

[0044]

[0045] Wherein, n represents the sample quantity, i.e. the total number of observation values in the data set; y i represents the actual observation value of the i-th sample, represents the prediction value of the i-th sample by the model, is the mean value of the error.

[0046] Compared with the related prior art, the beneficial effects of the present application are:

[0047] High sensitivity and high reliability of cycle slip detection, through three-layer decomposition of db4 wavelet, small amplitude and short-time carrier phase mutation can be captured on different scales at the same time; adaptive multi-threshold and multi-scale anomaly interval detection model, combined with statistical threshold and fluctuation amplitude criterion, significantly reduces the missed detection and false detection probability under noise interference;

[0048] Centimeter-level repair accuracy, in the cycle slip occurrence interval, combined with multiple prediction models such as convolutional neural network, support vector regression and ARIMA, the optimal model is dynamically selected through prediction error indicators; recursive sliding window prediction method is adopted to extrapolate and repair the missing segment point by point, realizing high-precision compensation of carrier phase observation value and meeting the centimeter-level positioning requirement;

[0049] High calculation efficiency, which can be deployed on resource-limited devices; in the model evaluation stage, only the limited starting point and its subsequent k samples in the test set are recursively predicted and error calculated, greatly reducing the calculation overhead of offline training and online switching; the wavelet decomposition and model prediction module can run in real time, which is suitable for single-frequency and low-power BDS receivers;

[0050] Strong adaptability and robustness, the adaptive threshold function adjusts the detection sensitivity of low-frequency approximation coefficients in real time according to the amplitude of high-frequency coefficients, and has good adaptability to signal characteristics under different ionospheres, shielding or dynamic working conditions; the multi-model switching mechanism avoids over-reliance on a single prediction model, improving the robustness of the overall algorithm. BRIEF DESCRIPTION OF DRAWINGS

[0051] Figure 1 The wavelet transform and model switching-based single-frequency BDS cycle slip detection and repair method flowchart provided by the present application;

[0052] Figure 2 A D1 high frequency coefficient contrast chart provided by an exemplary embodiment of the present disclosure is provided at a sampling rate of 1s;

[0053] Figure 3 A single-difference-double-difference low frequency coefficient chart provided by an exemplary embodiment of the present disclosure is provided;

[0054] Figure 4 A db4 wavelet reconstruction cycle slip detection result chart provided by an exemplary embodiment of the present disclosure is provided;

[0055] Figure 5 A db4 wavelet reconstruction 1 cycle slip detection result chart provided by an exemplary embodiment of the present disclosure is provided;

[0056] Figure 6 A cycle slip influence in wavelet transform chart provided by an exemplary embodiment of the present disclosure is provided;

[0057] Figure 7 A multi-prediction model switching flow chart provided by an exemplary embodiment of the present disclosure is provided;

[0058] Figure 8 A D1 wavelet coefficient prediction result contrast chart provided by an exemplary embodiment of the present disclosure is provided;

[0059] Figure 9 A D1 prediction error contrast chart provided by an exemplary embodiment of the present disclosure is provided;

[0060] Figure 10 A D2 wavelet coefficient prediction result contrast chart provided by an exemplary embodiment of the present disclosure is provided;

[0061] Figure 11 A D2 prediction error contrast chart provided by an exemplary embodiment of the present disclosure is provided;

[0062] Figure 12 A D3 wavelet coefficient prediction result contrast chart provided by an exemplary embodiment of the present disclosure is provided;

[0063] Figure 13 A D3 prediction error contrast chart provided by an exemplary embodiment of the present disclosure is provided;

[0064] Figure 14 A A3 wavelet coefficient prediction result contrast chart provided by an exemplary embodiment of the present disclosure is provided;

[0065] Figure 15 A A3 prediction error contrast chart provided by an exemplary embodiment of the present disclosure is provided;

[0066] Figure 16 A different size cycle slip detection chart provided by an exemplary embodiment of the present disclosure is provided. DETAILED DESCRIPTION

[0067] The present invention will be further described below with reference to the accompanying drawings. However, the present invention can be implemented in many different ways and should not be construed as limited to the embodiments shown; rather, these embodiments provide those skilled in the art with implementation methods that meet applicable legal requirements.

[0068] Example 1: As Figure 1 As shown, the single-frequency BDS cycle slip detection and repair method based on wavelet transform and model switching specifically includes the following steps: S1: Read carrier data, which includes carrier phase observation values ​​of the Beidou satellite system;

[0069] S2: Structural circumference probe measurement;

[0070] S3: Select the db4 wavelet as the wavelet basis, and decompose the cycle slip probe through three layers of wavelet decomposition to decompose the complex signal into wavelet coefficients D1, D2, D3 and A3 with different frequency characteristics;

[0071] S4: Based on the different oscillation interval characteristics of wavelet coefficients at different scales for satellite signal cycle slips, an adaptive multi-threshold multi-scale wavelet coefficient abnormal interval detection method is used to determine whether cycle slips occur.

[0072] S5: When cycle slips occur, a multi-prediction model switching mechanism is designed based on the time-frequency characteristics of each coefficient. The optimal prediction models of CNN / ELM, BP / CNN, ELM / SVR and ARIMA are matched for coefficients D1, D2, D3 and A3 respectively. Through multi-threshold abnormal interval detection and recursive prediction strategies, accurate repair of cycle slips is achieved.

[0073] Example 2: As Figures 2-3 As shown in the figure, this embodiment details the implementation process of the single-frequency BDS cycle slip detection and repair method based on wavelet transform and model switching, according to the constructed cycle slip detection measurement:

[0074] First, a non-differential probe measurement is constructed using L2I carrier phase observations from the C01 satellite at the WUH2 observation station. Then, a single-differential probe measurement is constructed by subtracting the L2I carrier phase observations from the C02 satellite. Using the inter-epoch single-differential probe measurement from the C01 satellite and the single-differential probe measurement constructed by subtracting the previous epoch from the next epoch of the C02 satellite, a double-differential probe measurement is then constructed between satellites C01 and C02. Finally, a triple-differential probe measurement is obtained by subtracting the previous probe measurement from the subsequent double-differential probe measurement. The constructed non-differential, single-differential, double-differential, and triple-differential probe measurements are then reconstructed using wavelet decomposition. Figure 2 This is the first layer of high-frequency signal after decomposition.

[0075] High-frequency coefficients directly reflect abrupt changes in the signal, such as... Figure 2As shown, the single-difference and double-difference measurements are obviously better than the non-difference and triple-difference measurements in the decomposed high-frequency coefficients D1. Through the analysis of the first layer high-frequency coefficients of the non-difference, single-difference, double-difference and triple-difference coefficients, it can be found that the single-difference and double-difference exhibit better signal characteristics and noise suppression capability, making the signal structure clearer. Although the single-difference has the smallest fluctuation in the high-frequency coefficients, from the perspective of the prediction error, the double-difference is better than the single-difference. Although the single-difference has the smallest fluctuation in the high-frequency coefficients, from the perspective of the prediction error, the double-difference is better than the single-difference. Therefore, the double-difference is used for cycle slip detection and repair in the present application. Figure 3 As can be seen, the single-difference base is large, and the prediction difficulty of its low-frequency coefficient (A3) is obviously greater than that of the double-difference. In the prediction process, it will cause the prediction result to deviate from the actual value more than the double-difference. The double-difference keeps good signal characteristics while having smaller prediction error, which is beneficial to improve the accuracy and reliability of cycle slip detection. Therefore, the double-difference measurement is used for cycle slip detection and repair in the present application.

[0076] Embodiment 3: as Figures 4-5 As shown in Table 1 and Table 1, this embodiment details the single-frequency BDS cycle slip detection and repair method based on wavelet transform and model switching. According to the selection of db4 wavelet as the wavelet base, the cycle slip detection quantity is decomposed into wavelet coefficients D1, D2, D3 and A3 with different frequency characteristics, and the specific implementation process is as follows:

[0077] Based on the comparison of various wavelet characteristics shown in Table 1, db wavelet is selected as the decomposition wavelet in this study, mainly based on the following reasons: db wavelet has both orthogonal and biorthogonal characteristics, which ensures the accuracy of the signal decomposition and reconstruction process. Compared with the single support length of haar wavelet, the support length of db wavelet is 2N-1, which provides better time-frequency positioning capability and is more suitable for capturing multi-scale features in GNSS carrier phase signals. Although db wavelet does not have symmetry, its N-order vanishing moment feature can effectively eliminate the polynomial trend in the signal, which has a significant advantage in dealing with complex noise in GNSS carrier phase data. Considering the characteristics of orthogonality, support length and vanishing moment number, db wavelet provides better signal decomposition performance while maintaining computational efficiency, so it is selected as the optimal wavelet base function in this study.

[0078] Table 1: Common index diagram of wavelet base

[0079] haar db N sym N coif N orthogonality yes yes yes yes biorthogonality yes yes yes yes compact support yes yes yes yes support length 1 2N-1 2N-1 6N-1 symmetry symmetric asymmetric approximately symmetric approximately symmetric regularity no no no no vanishing moments 1 N N 2N

[0080] The non-cycle-slip detection quantity and the 800 epoch added 1 cycle-slip detection quantity are decomposed and reconstructed into four wavelet coefficients using wavelet transform, and db4 wavelet base is used for three-layer wavelet decomposition. In order to effectively deal with the edge effect that may occur in the wavelet transform process, the signal extension mode is specially set to "symmetric padding" (spd) mode. This filling strategy effectively reduces the distortion phenomenon at the boundary while maintaining the continuity of the signal, improving the accuracy of the decomposition. The complete results of the reconstruction are shown in Figure 4 and Figure 5 .

[0081] By analyzing Figure 4 and Figure 5 In-depth analysis can find that the low-frequency coefficient A3 accurately captures and reflects the macroscopic change trend and overall characteristics of the cycle slip detection, and presents the main skeleton of the data; and the high-frequency coefficients D1, D2 and D3 clearly show the subtle fluctuations and transient characteristics of the cycle slip detection in each frequency band in different scales. The advantage of this multi-scale decomposition is that, compared with the original data, the component data not only has better stationarity and smoothness, but also makes the details information contained in the data more prominent, concentrated, and presents certain periodic characteristics. This structured data representation form greatly reduces the complexity of subsequent analysis, provides a more ideal data basis for cycle slip detection and repair, and effectively improves the overall processing efficiency and accuracy.

[0082] Embodiment 4: as Figure 6 shown, this embodiment details a single-frequency BDS cycle slip detection and repair method based on wavelet transform and model switching, which adopts an adaptive multi-threshold multi-scale wavelet coefficient abnormal interval detection method according to the different oscillation interval characteristics of different scale wavelet coefficients in response to satellite signal cycle slip, and the specific implementation process is as follows:

[0083] Figure 6 The comparative analysis of different scale wavelet coefficients in response to satellite signal cycle slip is shown. Through wavelet transform, the coefficients show obvious different characteristics under the conditions of cycle slip and non-cycle slip: the coefficient curve (solid line) under the condition of non-cycle slip shows small amplitude steady fluctuation, and the coefficient curve (red dotted line) under the condition of cycle slip produces significant disturbance near the cycle slip point. It can be observed from the decomposition result that the high-frequency detail coefficient D1 shows the highest sensitivity to cycle slip, generating a "Mexican hat" structure oscillation with significant amplitude (about ±0.5); with the increase of decomposition scale, the response amplitude of D2 and D3 coefficients to cycle slip decreases gradually, but still maintains similar oscillation mode; and the approximation coefficient A3 presents typical non-oscillating step characteristics. It can be found from Figure 7 that compared with D1, the influence epoch of D2 and D3 is also increased layer by layer. The different response characteristics of frequency domain components to cycle slip in multi-scale analysis provide an effective means for the detection of cycle slip in satellite observation data.

[0084] In view of Figure 6 the different oscillation interval characteristics of different coefficients, the application adopts an adaptive multi-threshold multi-scale wavelet coefficient abnormal interval detection method. For each level wavelet coefficient, the abnormal interval judgment criterion can be formally expressed as:

[0085]

[0086] Among them,

[0087] CS i indicates the cycle slip detection result of the i-th coefficient (1 indicates detecting cycle slip, 0 indicates not detecting cycle slip)

[0088] D1, D2, D3 indicate wavelet high-frequency coefficients of different decomposition scales

[0089] A3 indicates a wavelet low-frequency coefficient

[0090] μ i and σ i respectively indicate the mean and standard deviation of the absolute value of the i-th coefficient after robust statistical processing (deleting the top 5% outliers)

[0091] a i indicates the fluctuation amplitude in the abnormal interval Ω i , defined as the difference between the maximum and minimum values

[0092] indicates the first-order difference of the approximate coefficient A3, defined as the difference between the A3 values at the current time t and the previous time t-1

[0093] f(D1) and g(D1) are adaptive functions based on the maximum value of the D1 coefficient, used to adjust the detection sensitivity of the A3 coefficient.

[0094] The k data points in the abnormal interval Ω i are determined by the following algorithm:

[0095] Step 1: Mark the data points with absolute values exceeding the statistical threshold T i as the abnormal point set , i.e. where C i (t) indicates the i-th coefficient value at time t, and T i is the corresponding statistical threshold.

[0096] Step 2: Apply morphological structure closing operation Close to the abnormal point set to connect adjacent abnormal points into patches: where B is a morphological structure element. Based on the varying characteristics of different coefficients, a third-order structure element (morphological closing operator) is used for D1 and A3 coefficients to enhance the connection ability; a second-order structure element is used for D2 and D3 coefficients. This differentiated processing can effectively adapt to the varying characteristics of different wavelet coefficients.

[0097] Step 3: Identify all connected regions in , and select the largest region as the candidate abnormal interval. For A3 and D1 coefficients, when there is a specified epoch point, preferentially select the abnormal interval closest to the epoch point rather than the largest interval.

[0098] Step 4: First, by and Determine edge points, 3 points for D1 coefficient, 2 points for D3 coefficient, to capture cycle slip characteristics more comprehensively.

[0099] Step 5: Final abnormal interval Ω i = {t: t start ≤ t ≤ t end} contains k = t end -t start +1 data points.

[0100] For A3 coefficients, a differential processing method is used to enhance the mutation detection capability. According to the characteristics of D1 mutation coefficients and cycle slip values, the detection threshold and minimum fluctuation amplitude of A3 are dynamically adjusted. However, for large phase changes in a short time, cycle slip detection may not be able to judge for small cycle slip low frequency parts. In order to avoid this situation, when the abnormal interval of A3 coefficient and the abnormal interval of high frequency coefficient (d1, d2, d3) appear large deviation, then A3 coefficient is not considered. According to experience, the following three kinds of judgment conditions are divided:

[0101] When the modulus maximum value of D1 is greater than 9:

[0102]

[0103] When the modulus maximum value of D1 is less than or equal to 9 and greater than 3:

[0104]

[0105] When the modulus maximum value of D1 is less than or equal to 3 and greater than 0:

[0106]

[0107] Equations (1)-(4) represent three threshold and interval fluctuation settings, where T represents the threshold, Δ m in represents the minimum fluctuation amplitude, μ a3 and σ a3 represent the mean and standard deviation of A3 coefficient, respectively.

[0108] Based on the characteristics of different wavelet coefficients, the final condition for judging the existence of cycle slip is:

[0109]

[0110] i.e. at least two coefficients detect cycle slip, or only D1 coefficient detects cycle slip. This comprehensive judgment strategy makes full use of the characteristic response of different wavelet coefficients to cycle slip, especially the high sensitivity of high-frequency detail coefficient D1 to cycle slip, and improves the reliability of cycle slip detection, laying a foundation for subsequent cycle slip repair.

[0111] Embodiment 5: as Figures 7-15 As shown in Tables 2-7, this embodiment details the multi-prediction model switching mechanism designed according to the time-frequency characteristics of each coefficient when cycle slip occurs in the single-frequency BDS cycle slip detection and repair method based on wavelet transform and model switching, which matches the CNN / ELM, BP / CNN, ELM / SVR and ARIMA optimal prediction models for D1, D2, D3 and A3 coefficients respectively, and the specific implementation process is as follows:

[0112] Multi-prediction model switching mechanism:

[0113] Step 1: Divide the training set into a training subset and a test subset, use the training subset to establish cycle slip repair models for the initial model set, and input the test subset into the initial model set for prediction to obtain the prediction data of each model.

[0114] Step 2: Evaluate each model using model prediction evaluation parameters, find the model with good model evaluation parameters and high prediction accuracy for each wavelet coefficient, and form a candidate model set for subsequent prediction.

[0115] Step 3: Randomly select 6 points from the test set, use the candidate model set to predict the next k data using the data before these 6 points, and calculate the average root mean square of the 6 points of different models according to the root mean square of the 6 points.

[0116] Step 4: Select the model with the smallest average root mean square error, predict the next point before the cycle slip occurs, remove the earliest data point each time and add the newly predicted point until the k data prediction is completed.

[0117] Model prediction evaluation parameters: select root mean square error (RMSE), mean absolute error (MAE), standard deviation (SD) and variance (Variance) as model evaluation indexes, where RMSE and MAE range from [0, +∞), the smaller the better. The present application does not use the coefficient of determination (R 2 ) and mean absolute percentage error (MAPE) as evaluation indexes for wavelet coefficient prediction, because high-frequency wavelet coefficients usually exhibit small positive and negative fluctuations around zero. In this case, R2 may produce misleading results. Since the denominator (total sum of squares) is close to zero, even if the absolute error of prediction is small, R 2 may be abnormally low or unstable, or even negative. Similarly, MAPE also has serious problems when evaluating high-frequency wavelet coefficients, because its calculation formula contains an operation of dividing by the actual value (|actual value-predicted value| / |actual value|). When the wavelet coefficient is close to zero or crosses zero, the MAPE will appear a numerical maximum or even infinity, completely losing the evaluation significance. Especially when the data contains positive and negative alternating fluctuations, this instability is more obvious.

[0118] The calculation formulas of each index are as follows:

[0119]

[0120] Formulas (5)-(8) are the calculation formulas of the evaluation indexes, wherein n represents the sample number, i.e. the total number of observation values in the data set, y i represents the actual observation value of the i th sample, represents the predicted value of the i th sample by the model, is the average value of the error.

[0121] The back propagation (BP) neural network, radial basis function (RBF) neural network, extreme learning machine (ELM), support vector regression (SVR), convolutional neural network (CNN), and auto-regressive integrated moving average (ARIMA) are selected as the initial model set. The normal carrier phase observation values before the cycle slip are taken as the data set. Assuming that the cycle slip occurs at 1001 epochs, the first 1000 epochs are taken as the data set, wherein 700 epochs are taken as the training set and 300 epochs are taken as the test set. Since the generation of cycle slip will only cause the wavelet coefficient to fluctuate at the cycle slip, the wavelet coefficient interval affected by the cycle slip is k, so the first 700-k data of the training set are taken as the training subset and the remaining k data are taken as the test subset. The sequence length of the input and output is set to 15 and 1 respectively, and the prediction result is obtained by the recursive prediction method. The training subset is used to establish the cycle slip repair model of the initial model set, and the parameter settings of each initial model set are shown in Table 2.

[0122] Table 2: Parameter settings of the initial model set

[0123]

[0124] The recursive prediction is first performed on the test subset, which has a significant advantage in anomaly interval repair. Through wavelet decomposition, the original signal is separated into wavelet signals of different frequency bands, so that the prediction error of high-frequency wavelet coefficients significantly reduces the influence on the reconstructed overall signal. Although the high-frequency wavelet coefficients have small positive and negative fluctuation characteristics, which may lead to directional error and error accumulation problems in the recursive prediction process, since these high-frequency coefficients contribute relatively small in the signal reconstruction process, their influence will be effectively suppressed. Therefore, even if the recursive prediction has certain limitations in the processing of high-frequency coefficients, in the actual application of cycle slip repair, due to the multi-resolution characteristics of wavelet decomposition effectively isolating error propagation in different frequency bands, the overall repair effect is still quite reliable. This method is particularly suitable for processing cycle slip anomalies, because it can not only retain the overall trend and main features of the signal, but also effectively repair local abnormal points, achieving accurate restoration of the original signal.

[0125] As shown in Figure 8 , Figure 9 and Table 3, through comprehensive analysis of the prediction results and errors of the D1 high-frequency coefficient test subset, this study selects CNN and ELM as the candidate model set. As shown in Table 3, the CNN model performs the most outstanding, with an RMSE of 0.0142, an MAE of 0.0097, an SD of 0.015, and a Var of only 0.000225, all of which are the best; the ELM model follows closely, with an RMSE of 0.0158, an MAE of 0.0119, an SD of 0.0166, and a Var of 0.000277, which is significantly better than other traditional models. As shown in Figure 8 , CNN and ELM have stronger ability to capture the high-frequency detail features of D1 wavelet coefficients, and can more accurately reflect the fluctuation trend of the coefficients. As shown in Figure 9 , these two models have relatively stable and small errors at most epoch points, especially in handling the positive and negative fluctuation of high-frequency wavelet coefficients, and the error suppression effect is obvious. Considering the prediction accuracy, stability and adaptability to the characteristics of wavelet coefficients, CNN and ELM become the candidate model set for predicting D1 wavelet coefficients.

[0126] Table 3: Prediction effect of the initial model set on the D1 high-frequency coefficient test subset

[0127]

[0128] As shown in Figure 10 , Figure 11As shown in Table 4, by comprehensive analysis of the prediction results and errors of D2 wavelet coefficients, it can be observed that the overall prediction effect of D2 is indeed lower than that of D1 due to the increase in the number of prediction steps. The data shows that the BP model performs best in D2 prediction, with an RMSE of 0.011, followed by the CNN model, with an RMSE of 0.0112, and the MAE value of BP is slightly better than that of CNN. Figure 10 As shown in Table 4, by comprehensive analysis of the prediction results and errors of D2 wavelet coefficients, it can be observed that the overall prediction effect of D2 is indeed lower than that of D1 due to the increase in the number of prediction steps. The data shows that the BP model performs best in D2 prediction, with an RMSE of 0.011, followed by the CNN model, with an RMSE of 0.0112, and the MAE value of BP is slightly better than that of CNN. Figure 11 As shown in Table 4, by comprehensive analysis of the prediction results and errors of D2 wavelet coefficients, it can be observed that the overall prediction effect of D2 is indeed lower than that of D1 due to the increase in the number of prediction steps. The data shows that the BP model performs best in D2 prediction, with an RMSE of 0.011, followed by the CNN model, with an RMSE of 0.0112, and the MAE value of BP is slightly better than that of CNN.

[0129] Table 4: Prediction effect of the initial model set on the test subset of D2 high-frequency coefficients

[0130]

[0131] As Figure 12 , Figure 13 As shown in Table 5, by comprehensive analysis of the prediction results and errors of D3 wavelet coefficients, it can be observed that each model has a prediction effect on D3 wavelet coefficients. As shown in Table 5, in D3 prediction, the SVR model performs best, with an RMSE of 0.004, followed by the ELM model, with an RMSE of 0.0048, and the MAE value of SVR is also better than that of ELM. As Figure 12 shown, with the increase in the number of prediction steps, the error of recursive prediction gradually accumulates. Figure 13 Further, it is shown that although the error fluctuation trend of each model is similar, the prediction deviation of SVR and ELM models is relatively small compared to other models in subsequent epochs. Compared to the prediction results of D2, the fluctuation amplitude of D3 compared to D2 is reduced by about half. Considering the prediction accuracy, stability, and variance, the SVR model exhibits the best prediction performance, followed by the ELM model, so ELM and SVR are selected as the candidate model set for predicting D3 wavelet coefficients.

[0132] Table 5: Prediction effect of the initial model set on the test subset of D3 high-frequency coefficients

[0133]

[0134] As Figure 14 、 Figure 15 and Table 6, through comprehensive analysis of the prediction results and errors of A3 low frequency coefficients, it can be clearly observed that the ARIMA model performs well in predicting A3 low frequency coefficients, far exceeding other models. The data shows that the RMSE of the ARIMA model is only 0.0034, and the MAE is 0.0029, which is significantly better than the RBF model. From Figure 14 It can be seen that the prediction curve of the ARIMA model is almost completely coincident with the true value, while the ELM model has large fluctuations, and the prediction curve deviates from the actual value by more than 0.5. Figure 15 Further confirmed that the error of the ARIMA model always remains at a level close to zero with minimal fluctuations, so the ARIMA model is selected as the prediction model for predicting A3 wavelet coefficients.

[0135] Table 6: Prediction effect of initial model set on A3 high frequency coefficient test subset

[0136]

[0137] According to the analysis of the prediction of the candidate model set on different wavelet coefficients, the selected candidate model set is shown in Table 7.

[0138] Table 7: Determination of prediction model for different wavelet coefficients

[0139]

[0140]

[0141] Example 6: As Figure 16 and Tables 8-12, this embodiment details the specific implementation process according to the cycle slip repair results in the single-frequency BDS cycle slip detection and repair method based on wavelet transform and model switching:

[0142] In order to verify the effectiveness of the method of the application, the application uses the GNSS observation data of November 25, 2024 published by the CDDIS website for research. Satellite information collected at WUH2 station and URUM station is selected, data integrity is verified before the experiment, Turbo-Edit algorithm is used to detect and repair cycle slip by using GNSSer software, and the correctness of the experimental data is ensured.

[0143] (1) Model comparison results

[0144] In order to verify the effectiveness of the method, four prediction models without data decomposition reconstruction, BP / CNN / ELM / SVR, and four prediction models with data decomposition reconstruction, WT-BP / WT-CNN / WT-ELM / WT-SVR, are selected for comparison between prediction models. The experimental data is 900 epochs collected from 0:00 to 0:15 of WUH2 station with a sampling rate of 1s. Recursive prediction is used for multi-step prediction, and the prediction results of different methods are compared. Without data decomposition reconstruction means directly using the model to predict the double difference detection. With data decomposition reconstruction means decomposing the data into different wavelet coefficients, predicting different wavelet coefficients, and finally reconstructing different wavelet coefficients as the final prediction result. The comparison results of different methods are shown in Table 8.

[0145] Table 8: Comparison of prediction results of different prediction models

[0146] prediction model RMSE MAE SD Var BP 0.0349 0.0325 0.0291 0.000845 CNN 0.0239 0.0214 0.0231 0.000533 ELM 0.0355 0.0307 0.03 0.000898 SVR 0.0293 0.0248 0.0303 0.000918 WT-BP 0.0512 0.0402 0.0475 0.002254 WT-CNN 0.0316 0.0273 0.0316 0.000998 WT-ELM 0.0521 0.0409 0.0324 0.00105 WT-SVR 0.0298 0.0254 0.0297 0.000882 the method 0.0197 0.0164 0.0194 0.000376

[0147] From Table 8, under the RMSE index, taking the BP model as an example, the RMSE index of the model without data decomposition reconstruction is 0.0349, and the RMSE index of the model with data decomposition reconstruction is 0.0512, and the prediction accuracy is reduced by 46.7%. This is because the low frequency coefficients and high frequency coefficients of the double difference detection have different regular characteristics, and the same model cannot accurately predict the data with different characteristics, so different models are needed to predict different coefficients.

[0148] In addition, according to the comparison of prediction results without data decomposition reconstruction and with data decomposition reconstruction, the prediction results of CNN, SVR and WT-SVR are better, and the best prediction result is CNN with an RMSE of 0.0239. The prediction results of SVR model and WT-SVR are similar, and compared with other models, the prediction accuracy of SVR model is improved by 1.68% than WT-SVR model. This further illustrates that different wavelet coefficients need to select the most suitable prediction model.

[0149] Compared with the CNN model, the RMSE index of the method is 0.0197, the accuracy is improved by 17.57%; the MAE index is 0.0164, the accuracy is improved by 23.36%; the SD index is 0.0194, the accuracy is improved by 16.02%; and the Var index is 0.000376, the accuracy is improved by 29.46%. From the perspective of prediction performance index, the effectiveness of the method is fully verified.

[0150] (2) Different sampling rate detection results

[0151] In order to verify the effectiveness of the method of the application under different sampling rates, the carrier phase observation values of WUH2 station from 0:00 to 0:30 are selected, the double-difference probing quantity is constructed by using the carrier phase observation values of C01 and C02 satellites, and the detection repair analysis is carried out by adding cycle slip of C01 satellite. The 1s sampling rate is 1800 epochs, the 3s sampling rate is 600 epochs, the 5s sampling rate is 360 epochs, the 10s sampling rate is 180 epochs, and the 15s sampling rate is 120 epochs.

[0152] Table 9: Cycle slip detection results under different sampling rates

[0153]

[0154] The experimental results of Table 9 show that the method of the application has good cycle slip detection ability in the range of 1-10s sampling rate. Through the analysis of the detection results under different sampling rates, it can be seen that under the condition of 1s sampling rate, the detection value range is 0.96-1.03, which is closest to the theoretical value 1.0, the average detection value is 0.988, and the relative error is only 1.2%; under the conditions of 3s and 5s sampling rates, the detection results still maintain high accuracy, and the relative errors are 5.8% and 5.3% respectively; under the condition of 10s sampling rate, although the detection value is slightly higher than the theoretical value, the relative error is 5.5%, and the coefficient of variation is only 2.6%, which shows that the results have good consistency. The coefficient of variation analysis shows that the detection results under the conditions of 1s to 10s sampling rate have good stability, and the coefficients of variation are all controlled within 5% (1s is 3.0%, 3s is 4.9%, 5s is 4.2%, and 10s is 2.6%). This shows that in the application of cycle slip detection, the results with a coefficient of variation less than 5% can be considered as having good reliability and stability. Only when the sampling rate is reduced to 5s, the detection performance decreases significantly, the average detection value decreases to 0.805, the relative error reaches 19.5%, and the coefficient of variation rises to 8.5%, which shows that under such low sampling rate, the extraction ability of wavelet transform to signal characteristics is seriously affected. This performance attenuation is mainly due to the loss of high-frequency information caused by too large sampling interval, which makes the key D1 wavelet coefficient obviously decrease in the capture ability of cycle slip characteristics. Comprehensive analysis shows that the method of the application has good sampling rate adaptability, and can maintain stable detection effect in the commonly used range of 1-10s sampling rate, and even under the condition of 15s sampling rate, it can also detect 1 cycle slip through rounding, which provides strong support for the optimization of data acquisition strategy in different application scenarios, so that users can flexibly select appropriate sampling rate according to actual precision requirements and data storage conditions.

[0155] (3) Detection results of different size cycle slips

[0156] To verify the detection ability of the method of the application under different sizes of cycle slips, 1800 carrier phase observations of 1s sampling rate from 0:00 to 0:30 of WUH2 station were selected, double-difference detection quantities were constructed using carrier phase observations of C01 and C02 satellites, and different sizes of cycle slips were added to C01 satellite for detection repair analysis. A 0.5-cycle cycle slip was added at the 200th epoch, a 1-cycle cycle slip was added at the 400th epoch, a 5-cycle cycle slip was added at the 600th epoch, a 10-cycle cycle slip was added at the 800th epoch, a 50-cycle cycle slip was added at the 1000th epoch, and a 100-cycle cycle slip was added at the 1500th epoch. The detection results are shown in Figure 16 .

[0157] Figure 16 The three left-hand side illustrations show the detection details of smaller cycle slips. The first illustration shows that the peak value of the 0.5-cycle cycle slip at the 200th epoch is 0.5164 cycles, and the relative error is 3.28%; the second illustration shows that the peak value of the 1-cycle cycle slip at the 400th epoch is 0.991 cycles, and the relative error is only 0.9%; and the third illustration shows that the peak value of the 5-cycle cycle slip detection is 4.9765 cycles, and the relative error is 0.47%. The right-hand side illustrations show the detection results of the 10-cycle cycle slip at the 800th epoch, the peak value is 9.9848 cycles, and the relative error is 0.152%. The peak value at the 1000th epoch is 50.0125 cycles, and the relative error is 0.025%. The peak value at the 1500th epoch is 99.867 cycles, and the relative error is 0.133%.

[0158] It can be seen from Figure 16 that the method shows excellent detection ability for cycle slips from 0.5 cycles to 100 cycles, especially for the effective detection of half-cycle slips, with a relative error of 3.28%, which is considered a difficult point in cycle slip detection. The relative errors of various sizes of cycle slips are kept at a low level, from 3.28% of 0.5-cycle slip to 0.025% of 50-cycle slip, indicating that the method has good adaptability and stability. The experiment proves the strong ability of wavelet transform and multi-prediction model switching strategy in feature extraction and anomaly detection, enabling the algorithm to effectively distinguish cycle slips of different sizes and positions. This comprehensive detection ability for various cycle slips provides a reliable solution for the BDS cycle slip detection problem, significantly improving the usability and positioning accuracy of carrier phase observations.

[0159] (4) Model generalization verification

[0160] To verify the generalization performance of the method, this section systematically verifies it from two dimensions of different BDS satellite combinations and different stations. This experiment selects 900 carrier phase observations of 1 s sampling rate of WUH2 station, URUM station and ULAB station from 0:00 to 0:15. First, at WUH2 station, multiple groups of different BDS satellites are selected, and double-difference measurements are constructed by combining A satellite with a cycle slip and B satellite to verify the adaptability of the algorithm to different satellite signal characteristics. Second, stations with different geographical locations and observation conditions (WUH2, URUM and ULAB) are selected, and double-difference measurements are constructed using carrier phase observations of C01 satellite and C02 satellite, and a cycle slip is added to C01 satellite to verify the stability of the algorithm under different observation environments. All experiments use 1 s sampling rate data, and a 1-cycle cycle slip is artificially added at 400 epochs for detection and repair.

[0161] Table 10: Detection and repair performance under different BDS satellite combinations

[0162]

[0163] The satellite combination verification experiment selects five groups of satellites with different characteristics from the BDS constellation: C01-C02, C03-C08, C06-C09, C48-C59 and C59-C62. Double-difference measurements are constructed using carrier phase observations of two satellites in each group. These combinations cover satellites with different orbit types, different elevation angles and different signal strength characteristics. Table 10 shows the detection and repair performance under different satellite combinations. From the results in Table 10, it can be seen that the performance of different satellite combinations differs significantly. The C59-C62 combination performs best, with a relative error of only 0.2%, RMSE and MAE of 0.0162 and 0.0127, respectively, demonstrating excellent detection and repair capabilities; the C01-C02 combination is second, with a relative error of 0.4% and stable repair accuracy; while the C03-C08 combination has relatively poor performance, with a relative error of up to 10.5% and a repair RMSE of 0.1363, which is about 3-8 times higher than other combinations, but after data correction, it does not affect the final repair result.

[0164] Table 11: Cycle slip detection and repair performance under different stations

[0165]

[0166] In the station verification experiment, three stations with wide geographical distribution and different observation conditions were selected: WUH2 station located in central China, URUM station located in the arid region of Xinjiang, and ULAB station located in the Mongolian plateau. These three stations have significant differences in climate conditions, terrain environment and ionospheric characteristics, providing a good basis for the environmental adaptability verification of the algorithm. Table 11 shows the cycle slip detection and repair results at different stations. The results in Table 11 show that the method of the present application performs well in different station environments. WUH2 station performs best with a relative error of only 0.43%; ULAB station is second with a relative error of 1.93%; and URUM station has a relative error of 2.68%. The repair accuracy indicators of the three stations are also very close, with RMSE and MAE values remaining at a low level. This result shows that although there are differences in geographical environment and observation conditions at different stations, the method of the present application can still maintain stable detection and repair performance.

[0167] Combining the verification results of the two dimensions, the method of the present application shows good generalization ability in different satellite combinations and different geographical environments. The generalization verification results prove that the method of the present application has wide practical value and can be stably operated in various satellite combinations and different geographical environments, providing a reliable solution for cycle slip processing of BDS carrier phase observation data. This generalization ability is mainly due to the robust feature extraction of wavelet transform and the adaptive learning ability of multi-model switching strategy for different data characteristics.

[0168] (5) Detection results of different satellite types

[0169] In order to verify the detection ability of the method of the present application for different satellite types of Beidou satellites, this section experiment uses the 1s sampling rate carrier phase observation value of WUH2 station from 0:00 to 0:15, and carries out experiments on carrier phase observation values of GEO (Geostationary Earth Orbit), IGSO (Inclined Geosynchronous Orbit) and MEO (Medium Earth Orbit) three different types of Beidou satellites. A and B satellites of the same type are combined, and cycle slip is added at 500 epochs of satellite A to build double difference detection for cycle slip detection. Due to the difference in motion characteristics of satellites of different types, their observation value characteristics and noise levels may be significantly different, which poses a challenge to the robustness of the method.

[0170] Table 12: Comparison of cycle slip repair accuracy of different types of satellites

[0171]

[0172] From the experimental results of Table 12, it can be seen that the method of the present application has good cycle slip detection and repair capability on Beidou satellites of different orbit types. For 0.5 cycle slip, the detection values of satellites of the three orbit types are in the range of 0.476-0.481, which is very close to the theoretical value of 0.5 cycle, and the relative error is controlled within 4.8%; in 1 cycle slip detection, the detection values of GEO, IGSO and MEO satellites are 1.011 cycles, 0.950 cycles and 0.918 cycles respectively, which can accurately reflect the true cycle slip size; the 5 cycle slip detection results show that the detection values of satellites of each orbit type are in the range of 4.997-5.052 cycles, and the relative error is less than 1.1%; in 10 cycle slip detection, the detection values of the three satellite types are 9.932 cycles, 10.189 cycles and 10.125 cycles respectively, and the relative error is controlled within 1.9%, which can be accurately repaired to the theoretical value after integer processing. Due to the differences in motion characteristics, the observation values of satellites of different orbit types have significant differences in signal strength, multipath effect and other aspects, but the experimental results show that the method of the present application can effectively adapt to the observation characteristics of satellites of various orbit types, which is mainly due to the multi-resolution analysis capability of wavelet transform and the adaptive characteristics of the multi-prediction model switching strategy, which verifies the applicability of the method to accurately detect and repair the satellites of all orbit types of Beidou satellite system in the range of 0.5-10 cycle slips.

[0173] The above examples only express the embodiments of the present application, and the description is more specific and detailed, but it should not be understood as limiting the scope of the application. It should be noted that for ordinary skilled in the art, without departing from the concept of the present application, a number of modifications and improvements can be made, which are within the scope of protection of the present application.

Claims

1. A single-frequency BDS cycle slip detection and repair method based on wavelet transform and model switching, characterized in that: The method comprises the following steps: S1: reading carrier data, the carrier data comprising carrier phase observations of a Beidou satellite system; S2: constructing cycle slip detection quantities, the cycle slip detection quantities being formed by one or more of undifferenced, single-difference, double-difference and triple-difference; S3: selecting db4 wavelet as a wavelet base, performing three-layer discrete wavelet decomposition on the cycle slip detection quantities in step S2 to obtain high-frequency detail coefficients D1, D2 and D3 and low-frequency approximation coefficient A3; S4: performing adaptive multi-threshold abnormal interval detection based on the response characteristics of multi-scale coefficients; wherein, for the coefficients of wavelet decomposition at different levels in step S3, the calculation formula of the abnormal interval determination criterion is: where CS i denotes the cycle slip detection result of the ith coefficient, D1, D2, and D3 denote the first, second, and third level wavelet high-frequency coefficients, respectively, A3 denotes the wavelet low-frequency coefficient, μ i and σ i denote the mean and standard deviation of the absolute value of the ith coefficient after robust statistical processing, respectively; a i denotes the fluctuation amplitude within the abnormal interval Ω i , defined as the difference between the maximum and minimum values; denotes the first-order difference of the approximate coefficient A3, defined as the difference between the A3 values at the current time t and the previous time t-1; f(D1) and g(D1) are adaptive functions based on the maximum value of the D1 coefficient, used to adjust the detection sensitivity of the A3 coefficient; S5: when cycle slip is detected, performing multi-prediction model switching repair on the time-frequency characteristics of D1, D2, D3 and A3; S6: performing inverse wavelet reconstruction on the repaired D1, D2, D3 and A3 to obtain the repaired carrier phase observations as the final output. 2.The single-frequency BDS cycle slip detection and repair method based on wavelet transform and model switching according to claim 1, characterized in that: After reading the carrier data in step S1, the method further comprises preprocessing the carrier phase observations, which at least comprises any one or more of the following: Epoch time alignment: synchronously correcting the observation epochs of different satellite channels to ensure consistent time labels; Weighted processing: assigning weights to the observations according to the satellite elevation angle and signal-to-noise ratio, the weight coefficients being positively correlated with the sine value of the elevation angle and the signal-to-noise ratio; Outlier rejection: identifying and rejecting gross errors in the observation sequence using the median absolute deviation or 3σ criterion. 3.The single-frequency BDS cycle slip detection and repair method based on wavelet transform and model switching according to claim 1, characterized in that: The abnormal interval determination in step S4 further comprises: in the cycle slip detection result, 1 represents detecting cycle slip, and 0 represents not detecting cycle slip; the robust statistical processing specifically comprises deleting the first 5% of abnormal values.

4. The single-frequency BDS cycle slip detection and repair method based on wavelet transform and model switching according to claim 1, characterized in that: The adaptive multi-threshold multi-scale wavelet coefficient abnormal interval detection in step S4 comprises the following steps: S4-1: Mark data points with absolute values exceeding a statistical threshold T as an anomaly point set i The specific calculation formula is: Wherein, C i (t) represents the i-th coefficient value at time t, T i is the corresponding statistical threshold;​ S4-2: to the abnormal point set Adopt morphological close operation Close to handle in one piece, obtain Among them, the structure element B selects the third order structure element to D1 and A3 coefficient, selects the second order structure element to D2 and D3 coefficient; S4-3: In all connected regions are identified, and the largest region is selected as the candidate abnormal interval If the coefficients are A3 and D1 and there is a specified epoch point, the connected region closest to the epoch point is selected S4-4: Let the edge points be and For D1 coefficients, expand the interval boundaries outwards by 3 epochs each, and for D3 coefficients, expand the interval boundaries outwards by 2 epochs each; S4-5: The final abnormal interval is Ω i = {t: t start ≤ t ≤ t end}; where the number of data points contained in the abnormal interval is k = t end - t start + 1.

5. The single-frequency BDS cycle slip detection and repair method based on wavelet transform and model switching according to claim 4, characterized in that: In the adaptive multi-threshold multi-scale wavelet coefficient abnormal interval detection of step S4, the low-frequency approximation coefficient A3 is first subjected to first-order difference processing; And according to the characteristic that the high-frequency detail coefficient D1 is positively correlated with the cycle slip amplitude, the detection threshold and the minimum fluctuation amplitude of A3 are adjusted in real time; when there is a significant deviation between the abnormal interval corresponding to A3 and the abnormal intervals of high-frequency coefficients D1, D2 and D3, A3 is not considered.

6. The single-frequency BDS cycle slip detection and repair method based on wavelet transform and model switching according to claim 5, characterized in that: There are three cases for processing using the low-frequency approximation coefficient A3: When the modulus maximum value of D1 is greater than 9: When the modulus maximum value of D1 is less than or equal to 9 and greater than 3: When the modulus maximum value of D1 is less than or equal to 3 and greater than 0: where T denotes a threshold, Δ m in denotes a minimum fluctuation amplitude, μ a3 and σ a3 denote the mean and standard deviation of the A3 coefficients, respectively. Based on the characteristics of different wavelet coefficients, the calculation formula for finally judging cycle slip is: The specific condition for judging the existence of cycle slip is that at least two coefficients detect cycle slip, or only the D1 coefficient detects cycle slip.

7. The single-frequency BDS cycle slip detection and repair method based on wavelet transform and model switching according to claim 1, characterized in that: The multi-prediction model switching repair of step S5 comprises the following steps: S5-1: establishing a candidate model set comprising at least convolutional neural network, back propagation neural network, radial basis neural network, extreme learning machine, support vector regression and autoregressive integrated moving average; S5-2: dividing the historical samples into a training subset and a test subset to train the coefficient predictors of each model; S5-3: randomly extracting reference points from the test set to calculate the evaluation indexes of each model at the points; S5-4: Select the optimal model for the recursive prediction repair of the current coefficient in the abnormal interval according to the evaluation index; S5-5: Use a sliding window to extrapolate point by point until the missing segment repair of length k is completed.

8. The single-frequency BDS cycle slip detection and repair method based on wavelet transform and model switching according to claim 7, characterized in that: The evaluation index at least includes root mean square error RMSE, mean absolute error MAE, standard deviation SD of prediction error, and variance Var, and the calculation formula of each evaluation index is: where n represents the number of samples, i.e., the total number of observations in the dataset; y i represents the actual observation value of the i-th sample, represents the predicted value of the i-th sample by the model; is the mean of the errors.

Images