BDS cycle slip detection method based on GM model

Abstract

The invention discloses a BDS cycle slip detection method based on a GM model, and relates to the technical field of Beidou navigation positioning data processing. Firstly, a GNSS observation equation is constructed to obtain a double-frequency carrier phase; carrying out rough detection by utilizing a Melbourne-Wubbena wide lane combination, and positioning a suspicious epoch; calculating a total electron content change rate (TECR) sequence, and establishing a GM (1, 1) model in a moving window to predict the TECR after level ratio test and translation transformation; and finally, outputting a cycle slip mark by fusing rough detection and fine detection results. Actually measured BDS2 / 3-GEO, IGSO and MEO1s sampling data show that the method can accurately detect equal-cycle slip, frequency-ratio cycle slip and one-cycle small cycle slip, does not need an external sensor or large-scale training, and has the advantages of real-time performance, high precision and low cost.

Description

Technical Field

[0001] The present invention belongs to the technical field of Beidou navigation and positioning data processing, and particularly relates to a BDS cycle-slip detection method based on a GM model. Background Art

[0002] With the completion of the global networking of the Beidou satellite navigation system (BDS) and the official provision of services to users, the global satellite navigation system has welcomed a new important member. The carrier phase observation value is a key observable for BDS precise positioning, navigation, and timing (PNT). High-quality and continuous carrier phase directly determines the performance of high-precision applications. However, in actual observations, affected by factors such as signal loss of lock, tree or building occlusion, and receiver failure, the carrier phase often cannot be continuously collected, and cycle-slips will occur in the observation sequence. Cycle-slips will interrupt the continuity of the integer ambiguity, thereby reducing the accuracy and reliability of GNSS positioning and timing. Therefore, cycle-slip detection and repair have always been the core link in the GNSS data preprocessing chain.

[0003] Traditional BDS cycle-slip detection methods are mainly based on models such as the Melbourne-Wübbena (MW) wide-lane combination, geometry-free (GF), ionosphere-free (IF), and total electron content ratio (TECR) of the free ionosphere combination. Some researchers have proposed improvements on this basis. For example, by studying the characteristics of the MW combination of different constellations through static observations; combining TECR with the MW combination to uniquely determine the cycle-slips at the L1 and L2 frequency points; using the phase geometry-free combination in the second-order time difference model to jointly combine with the MW combination to improve the detection sensitivity. These methods have improved the detection rate to a certain extent, but most still take the MW combination as the core. When the same amplitude cycle-slips (equal cycle-slips) occur at two frequency points, the MW combination cannot judge the cycle-slip amplitude or determine the frequency point where it occurs. At this time, it is often necessary to further discriminate by means of the TECR method. If stable detection is to be maintained at a higher sampling rate or in a multi-constellation environment, reliable prediction of the TECR sequence is also required.

[0004] Traditional BDS cycle slip detection methods mainly rely on the application of models such as the MW combination, the Geometry-Free (GF) combination, the Ionosphere-Free (IF) combination, and the TECR method. Many scholars have improved traditional methods. For example, the MW combination characteristics of different types of BDS satellites have been studied through static observations; the TECR method and the Melbourne-Wübbena wide lane (MWWL) linear combination have been jointly used to uniquely determine cycle slips at the L1 and L2 frequencies; a new combination of two test parameters, the phase off-sphere free combination and the MW combination in the second-order time difference model, has been used, but this method depends on the stability of the ionosphere; some scholars have also proposed algorithmic innovations, such as the circular arc segmentation method using the GF and MW combinations to mark cycle slip values by dividing the observed circular arc; a new cycle slip detection and repair method based on polynomial fitting of the IF combination TECR method, especially the TECR method, which proves that the new algorithm can detect and repair cycle slips with a very high success rate; an improved MW combination cycle slip detection method based on complete ensemble empirical mode decomposition, permutation entropy, and wavelet denoising. To improve the detection accuracy, some scholars have proposed first comparing the MW and PIR (Phase Ionospheric Residual) combinations with the moving window average value, using these two combinations to detect cycle slips. Then, the epoch-differenced wide lane combination is used to estimate the changes in the smartphone position and clock bias. Finally, the epoch-differenced wide lane and IF combinations are used to estimate the cycle slip values at each frequency, and the least squares ambiguity decorrelation adjustment method is further used to obtain integer solutions. Based on the multi-frequency observation combination theory, some scholars have derived a three-frequency MW and geometry-free combination model and verified its effectiveness at a 1-second sampling interval. However, as the sampling interval increases, this method will be affected to a certain extent.

[0005] In recent years, several related patent solutions have been disclosed for the characteristics of BDS multi-frequency observations: The invention patent with the patent number and publication number CN115932922A discloses a cycle slip detection method based on BDS four-frequency data, which proposes to construct ultra-wide lane and narrow lane detection quantities through four-frequency pseudorange-phase combinations, and jointly perform four-frequency cycle slip detection with a geometry-free combination; this solution makes up for the detection blind spots of a single method and improves the detection accuracy, but still relies on the optimization of combination coefficients and has limited adaptability to rapidly changing ionospheric scenarios. The invention patent with the patent number and publication number CN115826003A discloses a BDS three-frequency cycle slip detection method assisted by Doppler integration, which uses three-frequency Doppler integration to calculate the change in pseudorange, and then jointly detects it with phase difference and the second difference of ionospheric residuals; this method can weaken the influence of Doppler integration error at low sampling rates, but still has insufficient ability to predict carrier phase with high frequency updates. The invention patent with the patent number and publication number CN107505642A / B discloses an INS-assisted real-time BDS single-frequency cycle slip detection method, which uses an inertial navigation system (INS) to extrapolate the double-difference phase prediction value and compare it with the measured double-difference to identify cycle slips; although the real-time performance is improved, it relies on an external INS sensor and is only applicable to single-frequency scenarios. The invention patent with the patent number and publication number CN112346093A discloses a method for repairing BDS cycle slips, which uses an improved BP neural network and NAR network to repair cycle slips, and selects the better one for output after comparison; this solution focuses on cycle slip repair rather than detection, and the neural network relies on a large number of training samples and is difficult to handle short-arc observations or data loss.

[0006] GM (Grey Model) is one of the most commonly used prediction models in grey system theory. Its characteristic is that it can establish an approximate differential equation to extrapolate the sequence trend with only a very small amount of partial known data. Its basic idea is to perform cumulative smoothing on the original data to amplify the overall change trend and weaken random fluctuations, and then fit the cumulative sequence with a first-order linear equation and back-calculate the predicted value. Due to its extremely low requirements for sample size, distribution type, and prior information, the GM model is especially suitable for rapid prediction in environments with small samples, short time periods, and incomplete information.

[0007] In summary, existing technologies have made progress in multi-frequency BDS cycle slip detection, but they still face the following common problems: Limitations of the MW-dominated method: When equal cycle slips occur at two frequencies, the MW combination cannot accurately locate the cycle slip and frequency attribution, and other combinations (such as TECR) are required to assist. TECR sequence prediction gap: Although the TECR method can make up for the blind spot of the MW combination, its effectiveness depends on the accurate prediction of the rate of change of ionospheric content; existing schemes mostly use polynomials or empirical models, which makes it difficult to obtain stable predictions under conditions of small data volume and unknown models. Real-time and adaptive requirements: With the popularization of high-frequency, short sampling interval observations, there is a need for a cycle slip detection method that does not require a large number of training samples, can be updated over time, and can work in real time in a multi-constellation and multi-frequency environment.

[0008] Therefore, how to construct a method that combines the MW combination and the TECR method and uses the grey GM(1,1) prediction model to adaptively predict the TECR sequence within a moving window, and then achieve high-precision real-time detection of equal-cycle cycle slips, frequency ratio cycle slips and even small cycle slips, has become a key issue that needs to be urgently solved in BDS cycle slip processing technology. Summary of the invention

[0009] In view of the above limitations and needs, the present invention proposes a BDS cycle slip detection method based on the GM model, which adopts the technical scheme of MW wide-lane combination coarse detection, TECR combination fine detection and moving window grey prediction adaptive compensation, and overcomes the problems of MW combination missed detection when equal cycle slips occur at two frequency points, difficulty in real-time prediction of TECR sequence and model instability under small sample conditions; the applicability of the grey model is ensured by level ratio test and translation transformation, and the accurate positioning of cycle slips at all frequency points is achieved by residual judgment and set fusion, so as to complete the high-precision real-time detection of equal cycle slips, frequency ratio cycle slips and small cycle slips of dual-frequency observation values ​​of BDS GEO, IGSO and MEO satellites at a sampling interval of 1s.

[0010] The present invention provides a BDS cycle slip detection method based on the GM model, which specifically comprises the following steps:

[0011] S1: Original observation value acquisition: construct the GNSS observation equation, collect the carrier phase observation values ​​of the target BDS satellite at the first carrier frequency point and the second carrier frequency point, and obtain the dual-frequency carrier phase observation sequence;

[0012] S2: Wide-lane combination rough inspection: Based on the dual-frequency carrier phase observation sequence described in step S1, a Melbourne-Wübbena combination observation quantity is constructed, and the first round of cycle slip rough judgment is performed on each epoch according to the wide-lane distance mutation threshold to obtain a first suspicious epoch set;

[0013] S3: TECR sequence formation: Based on the first set of suspicious epochs described in step S2, calculate the TECR sequence of the free ionosphere combination using the dual-frequency carrier phase observation sequence described in step S1;

[0014] S4: TECR sequence preprocessing: Perform a cumulative generation operation on the TECR sequence described in step S3 to obtain a first-generation sequence, and perform a ratio test on the first-generation sequence; if the ratio test fails, perform a translation transformation on the TECR sequence until the ratio test passes to obtain a corrected TECR sequence and its first-generation sequence;

[0015] S5: GM(1,1) prediction: In a moving window containing the most recent K epochs, establish a grey GM(1,1) model with the first-generation sequence described in step S4 as the input to predict the TECR value of the next epoch;

[0016] S6: Residual determination: Calculate the residual ε between the TECR predicted value described in step S5 and the corresponding measured TECR value; when the residual is greater than the preset threshold, record the corresponding epoch in the second set of suspicious epochs;

[0017] S7: Result fusion and output: Fusion the first set of suspicious epochs described in step S2 and the second set of suspicious epochs described in step S6 to generate and output the final cycle slip determination result.

[0018] As a preferred technical solution of the present invention, constructing the GNSS observation equation in step S1 specifically includes constructing a carrier phase observation equation and a pseudorange observation equation, and the calculation formulas for constructing the equations are respectively:

[0019] Construction of the carrier phase observation equation:

[0020]

[0021] where λ i is the wavelength of the carrier phase observation value , i represents the number of frequencies; ρ is the distance between the satellite and the ground; c is the speed of light in a vacuum; dt r , dt s are the clock errors of the receiver and the satellite respectively; ρ ion,i is the ionospheric delay; ρ trop is the tropospheric delay; N i is the integer ambiguity; is the carrier phase observation noise.

[0022] Construction of the pseudorange observation equation:

[0023]

[0024] where, is the pseudorange observation noise.

[0025] As a preferred technical solution of the present invention, the calculation of the first set of suspicious epochs in step S2 includes the following steps:

[0026] S2-1: Wide-lane phase combination construction: In the dual-frequency carrier phase observation sequence, according to formula (1), a wide-lane phase combination is established through dual-frequency GNSS observations. The specific calculation formula is:

[0027]

[0028] where, assuming the carrier phase observations of frequency points f1 and f2 are respectively and their wavelengths are λ1 and λ2 respectively;

[0029] S2-2: Narrow-lane pseudorange combination construction: In the corresponding pseudorange sequence P, according to formula (2), a narrow-lane pseudorange combination is established. The specific calculation formula is:

[0030]

[0031] where, assuming the pseudorange observations of frequency points f1 and f2 are P1 and P2 respectively;

[0032] S2-3: MW combination construction: The MW combination observation is generated through formula (3) and formula (4). The specific calculation formula is:

[0033]

[0034] S2-4: Wide-lane normalization and cycle slip detection calculation: The cycle slip detection is obtained by normalizing formula (5) to obtain the cycle slip detection quantity of the combination. The specific calculation formula is:

[0035]

[0036] where, the wide-lane combination wavelength is λ MW = c / (f1 - f2);

[0037] S2-5: Wide-lane normalization and cycle slip detection calculation: Differencing the N MW between adjacent epochs to weaken the multipath effect and receiver noise. The specific calculation formula is:

[0038] ΔN MW = N MW (t + 1) - N MW (t) (7)

[0039] S2-6: Threshold determination and generation of suspicious epochs: The specific calculation formula is:

[0040] ΔN MW= ΔN1 - ΔN2 (8)

[0041] Wherein, ΔN1 and ΔN2 represent cycle slip values at two frequencies; when ΔN MW ≥ γMW, the corresponding epoch is recorded in the first set of suspicious epochs M, where γMW is a preset threshold value.

[0042] As a preferred technical solution of the present invention, the formation of the TECR sequence in step S3 includes the following steps:

[0043] S3-1: Calculation of the total ionospheric electron content: Using two sets of carrier phase observation values obtained by a dual-frequency receiver, calculate the total ionospheric electron content at a certain epoch k, and the specific calculation formula is:

[0044]

[0045] Wherein, b1 and b2 are the hardware time delays at the receiver end and the satellite end, f1 and f2 are the frequency points in the dual-frequency carrier phase observation sequence, and are the observed quantities, λ1 and λ2 are the wavelengths,

[0046] S3-2: Calculation of the ionospheric change rate: Obtain the ionospheric change rate according to the difference between epochs, and the specific calculation formula is:

[0047]

[0048] Wherein, Δt is the sampling interval;

[0049] S3-3: Formation of cycle slip detection quantity: At epoch k, deduce the cycle slip detection quantity according to formula (9) and formula (10), and the specific calculation formula is:

[0050]

[0051] S3-4: Auxiliary cycle slip determination, and the specific calculation formula is:

[0052]

[0053] Wherein, when there is no cycle slip or the cycle slip has been repaired at epoch (k–1), use the ΔN shown in formula (12) TECR and the relationship between the frequency point and the wavelength to deduce the suspicious cycle slip value at epoch k, and participate in the subsequent residual determination process as a part of the free ionospheric combination TECR sequence.

[0054] As a preferred technical solution of the present invention, the formation of the TECR sequence in step S3 further includes a prediction process of the ionospheric total electron content change rate, which specifically includes the following steps:

[0055] At epoch k, first calculate the instantaneous change rate using the measured TECR values at epochs (k - 1) and (k - 2). The specific calculation formula is:

[0056]

[0057] Extrapolate the change rate obtained from formula (13) to epoch k to obtain the TECR prediction value. The specific calculation formula is:

[0058]

[0059] Substitute the TECR prediction value prediction into formulas (8) and (12), and solve the equations simultaneously to obtain ΔN1 and ΔN2. Then round ΔN1 and ΔN2 to obtain the cycle slips at the two frequencies.

[0060] As a preferred technical solution of the present invention, the specific calculation formula for the ratio test in step S4 is:

[0061]

[0062] If ω(t) is within the interval then it is determined that the TECR sequence can directly adopt the GM(1,1) model;

[0063] If the TECR value sequence is not within the interval corresponding to ω(t) after inspection, this TECR value sequence cannot adopt the GM(1,1) model; perform a translation transformation on this TECR value, specifically add an arbitrary constant c to each TECR value, solve using the GM(1,1) model and then subtract c, and repeat the ratio test operation until the TECR value sequence is within the interval corresponding to ω(t) after inspection.

[0064] As a preferred technical solution of the present invention, step S5 further includes performing grey GM(1,1) modeling and feasibility test on the free ionosphere combined TECR sequence, specifically:

[0065] x (0) (1) is the calculated TECR value of the first epoch, and the original TECR sequence arranged in epoch order is: X (0) =(x (0) (1), x (0) (2), …, x (0) (n)). Perform an accumulative generating operation on X (0) to obtain the first - order generated sequence:

[0066]

[0067] Assume that the first - order generated sequence satisfies the first - order grey differential equation:

[0068]

[0069] Discretize formula (17) at the sampling interval Δt=(t + 1)-t = 1, and rewrite the in formula (17) as Find the best function match by minimizing the sum of the squares of the errors and The discretized equation is obtained as:

[0070] x (0) (t)+ax (1) (t)=u (18)

[0071] where, Δx (1) =x (1) (t)-x (1) (t - 1)=x (0) (t), then formula (18) is rewritten as:

[0072] x (0) (t)=-ax (1) (t)+u (19)

[0073] Perform adjacent mean generation on the accumulated generating sequence. Adjacent mean generation is for an equidistant time series, and new data is constructed using the average of adjacent data. Its calculation formula is:

[0074] z (1) (t)=(x (1) (t)+x (1) (t - 1)) / 2, t = 2, 3, …, n (20)

[0075] Then formula (19) is rewritten as:

[0076] x (0) (t)=-az (1) (t)+u (21)

[0077] Using the least squares method to minimize the sum of the squares of the differences between the values of the fitting function and the known data, then formula (21) is written in matrix form:

[0078] Y = BU (22) where,

[0079]

[0080] Then the calculation formula for U is:

[0081] U = [a u] T =(B T B) -1 (B T Y) (23)

[0082] The parameters a and u have been obtained. Substituting them into Equation (17), we get:

[0083]

[0084] Predict the TECR value of the next epoch according to Equation (24)

[0085] Perform a residual test on the model. The residual test formula is:

[0086]

[0087] where the residual ε(t) < 0.2; if the residual is greater than or equal to 0.2, it is determined that the grey model's prediction of TECR for this epoch is unreliable, and return to step S4 to re-adjust the TECR sequence and then re-build the model.

[0088] As a preferred technical solution of the present invention, the moving window GM(1,1) prediction in step S5 includes:

[0089] At epoch k, select the first-generated sequence data of the nearest m epochs as the modeling sample;

[0090] As the epoch progresses, every time a new epoch k + 1 arrives, move the window forward by one epoch as a whole: eliminate the data of the oldest epoch (k - m + 1), introduce the data of the latest epoch (k + 1), and keep the sample size unchanged at m;

[0091] After each window update, re-build the grey GM(1,1) model based on the m samples, obtain the TECR prediction value of the new epoch, and use it for the residual determination in step S6.

[0092] Compared with the related prior art, the beneficial effects of the present invention are:

[0093] Through the three-layer structure of "MW wide-lane combination rough detection + TECR combination fine detection + GM(1,1) moving window prediction", the present invention can simultaneously detect equidistant cycle slips, frequency ratio cycle slips, and 1-week small cycle slips, making up for the undetected blind area of the traditional MW combination in the equidistant cycle slip scenario.

[0094] The GM(1,1) model of the present invention uses the accumulated generating operator (AGO) to transform the original TECR sequence into an approximate exponential (linear) sequence, then solves the grey differential equation parameters by the least square method, and finally obtains the prediction value by the inverse accumulated generating operator (IAGO); this process only needs a small amount of adjacent data to converge, and is especially suitable for the scenarios of 1-second sampling, short-arc observation, or data missing.

[0095] Before GM modeling, the present invention performs a ratio test on the TECR sequence; if the test is not satisfied, the sequence is re - fallen into the modelable interval through an overall translation transformation and then modeled to ensure the reliability of grey prediction.

[0096] The grey model of the present invention is repeatedly reconstructed within a sliding window containing the most recent K epochs, which can adaptively capture the local changes of ionospheric delay, improve the real - time response ability to sudden cycle slips, and at the same time suppress the cumulative error of prediction caused by long - term drift.

[0097] The present invention can be realized without relying on Doppler observations, inertial sensors or large - sample training of deep learning, and only using conventional dual - frequency carrier and pseudorange observations; it is easy to be embedded in the existing GNSS pre - processing chain, providing a simple and efficient new cycle slip detection solution for high - precision positioning, navigation and timing of BDS. BRIEF DESCRIPTION OF THE DRAWINGS

[0098] Figure 1 is a flowchart of a BDS cycle slip detection method based on the GM model provided by the present invention;

[0099] Figure 2 is a TECR value diagram of L2I and L7I / L7Z / L7D combined observations of the embodiment provided by the present invention;

[0100] Figure 3 is a TECR cumulative value diagram of L2I and L7I / L7Z / L7D combined observations of the embodiment provided by the present invention;

[0101] Figure 4 is a ratio test result diagram of the original TECR values of L2I and L7I / L7Z / L7D of the embodiment provided by the present invention;

[0102] Figure 5 is a horizontal ratio test result diagram of the TECR values after translation transformation of L2I and L7I / L7Z / L7D of the embodiment provided by the present invention;

[0103] Figure 6 is a cycle slip detection result diagram of L2I and L7I of BDS2 of the embodiment provided by the present invention;

[0104] Figure 7 is a cycle slip detection result diagram of L2I and L7Z / L7D of BDS3 of the embodiment provided by the present invention;

[0105] Figure 8 is a cycle slip detection result diagram of some satellites of BDS2 of the embodiment provided by the present invention;

[0106] Figure 9 is a cycle slip detection result diagram of some satellites of BDS3 of the embodiment provided by the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0107] The present invention will be further described below in conjunction with 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; on the contrary, these embodiments provide implementations that meet the applicable legal requirements for those skilled in the art.

[0108] Embodiment 1: As Figure 1 shown, the present invention specifically includes the following steps:

[0109] S1: Acquisition of original observations: Construct a GNSS observation equation, collect the carrier phase observations of the target BDS satellite at the first carrier frequency and the second carrier frequency, and obtain a dual-frequency carrier phase observation sequence; the construction of the GNSS observation equation specifically includes the construction of the carrier phase observation equation and the construction of the pseudorange observation equation, and the calculation formulas for the equation construction are respectively:

[0110] Construction of the carrier phase observation equation:

[0111]

[0112] In the formula, λ i is the wavelength of the carrier phase observation value , i represents the number of frequencies; ρ is the distance between the satellite and the ground; c is the speed of light in vacuum; dt r , dt s are respectively the clock errors of the receiver and the satellite; ρ ion,i is the ionospheric delay; ρ trop is the tropospheric delay; N i is the integer ambiguity; is the carrier phase observation noise.

[0113] Construction of the pseudorange observation equation:

[0114]

[0115] Among them, is the pseudorange observation noise.

[0116] S2: Coarse detection of the wide-lane combination: Based on the dual-frequency carrier phase observation sequence described in step S1, construct a Melbourne-Wübbena combined observable, and perform the first-round cycle slip coarse determination on each epoch according to the wide-lane distance mutation threshold to obtain the first set of suspicious epochs; the calculation of the first set of suspicious epochs includes the following steps:

[0117] S2-1: Construction of the wide-lane phase combination: In the dual-frequency carrier phase observation sequence, according to formula (1), establish a wide-lane phase combination through dual-frequency GNSS observations, and the specific calculation formula is:

[0118]

[0119] Among them, let the carrier phase observables of frequency points f1 and f2 be Their wavelengths are λ1 and λ2 respectively;

[0120] S2-2: Narrow-lane pseudorange combination construction: In the corresponding pseudorange sequence P, according to formula (2), establish the narrow-lane pseudorange combination, and the specific calculation formula is:

[0121]

[0122] Among them, let the pseudorange observables of frequency points f1 and f2 be P1 and P2 respectively;

[0123] S2-3: MW combination construction: Generate the MW combination observable through formula (3) and formula (4), and the specific calculation formula is:

[0124]

[0125] S2-4: Wide-lane normalization and cycle slip detection calculation: Normalize formula (5) to obtain the cycle slip detection quantity and obtain the cycle slip detection quantity of the combination. The specific calculation formula is:

[0126]

[0127] Among them, the wide-lane combination wavelength is λ MW = c / (f1 - f2);

[0128] S2-5: Wide-lane normalization and cycle slip detection calculation: Differentiate N MW between adjacent epochs to weaken the multipath effect and receiver noise. The specific calculation formula is:

[0129] ΔN MW = N MW (t + 1) - N MW (t) (7)

[0130] S2-6: Threshold determination and suspicious epoch generation: The specific calculation formula is:

[0131] ΔN MW = ΔN1 - ΔN2 (8)

[0132] Among them, ΔN1 and ΔN2 represent the cycle slip values at two frequencies; when ΔN MW ≥ γMW, record the corresponding epoch into the first suspicious epoch set M, where γMW is a preset threshold.

[0133] S3: Formation of TECR sequence: Based on the first set of suspicious epochs described in step S2, use the dual-frequency carrier phase observation sequence described in step S1 to calculate the free ionosphere combined TECR sequence; the formation of the TECR sequence includes the following steps:

[0134] S3-1: Calculation of total ionospheric electron content: Through two sets of carrier phase observations obtained by a dual-frequency receiver, calculate the total ionospheric electron content at a certain epoch k. The specific calculation formula is:

[0135]

[0136] where b1 and b2 are the hardware time delays at the receiver end and the satellite end, f1 and f2 are the frequency points in the dual-frequency carrier phase observation sequence, and are the observables, λ1 and λ2 are the wavelengths,

[0137] S3-2: Calculation of ionospheric change rate: Obtain the ionospheric change rate according to the difference between epochs. The specific calculation formula is:

[0138]

[0139] where Δt is the sampling interval;

[0140] S3-3: Formation of cycle slip detection quantity: At epoch k, derive the cycle slip detection quantity according to formulas (9) and (10). The specific calculation formula is:

[0141]

[0142] S3-4: Auxiliary cycle slip determination, the specific calculation formula is:

[0143]

[0144] where when there is no cycle slip or the cycle slip has been repaired at epoch (k–1), use the ΔN shown in formula (12) TECR and the relationship between the frequency point and the wavelength to derive the suspicious cycle slip value at epoch k, and participate in the subsequent residual determination process as part of the free ionosphere combined TECR sequence.

[0145] The formation of the TECR sequence also includes the prediction process of the ionospheric total electron content change rate, which specifically includes the following steps:

[0146] At epoch k, first calculate the instantaneous change rate using the measured TECR values of epochs (k-1) and (k-2). The specific calculation formula is:

[0147]

[0148] Extrapolate the rate of change obtained from Equation (13) to epoch k to obtain the predicted value of TECR. The specific calculation formula is as follows:

[0149]

[0150] Substitute the predicted value of TECR into Equations (8) and (12) and solve them simultaneously to obtain ΔN1 and ΔN2. Then, round ΔN1 and ΔN2 to obtain the cycle slips at the two frequencies.

[0151] S4: Preprocessing of the TECR sequence: Perform an accumulative generation operation on the TECR sequence to obtain a first-generation sequence, and perform a ratio test on the first-generation sequence. If the ratio test fails, perform a translation transformation on the TECR sequence until the ratio test passes to obtain a corrected TECR sequence and its first-generation sequence. The specific calculation formula for the ratio test described in step S4 is as follows:

[0152]

[0153] If ω(t) is within the interval , it is determined that the GM(1,1) model can be directly applied to the TECR sequence;

[0154] If the TECR value sequence fails the test and is not within the interval corresponding to ω(t), the GM(1,1) model cannot be applied to this TECR value sequence. Perform a translation transformation on the TECR values, specifically, add an arbitrary constant c to each TECR value, solve using the GM(1,1) model, and then subtract c. Repeat the ratio test operation until the TECR value sequence passes the test and is within the interval corresponding to ω(t).

[0155] S5: GM(1,1) prediction: In a moving window containing the most recent K epochs, establish a grey GM(1,1) model with the first-generation sequence as the input to predict the TECR value of the next epoch. Step S5 further includes grey GM(1,1) modeling and feasibility testing for the free ionosphere combined TECR sequence, specifically:

[0156] x (0) (1) is the calculated value of TECR for the first epoch. The original TECR sequence arranged in epoch order is: X (0) =(x (0) (1), x (0) (2), …, x (0) (n)). Perform an accumulative generation operation on X (0) to obtain the first-generation sequence:

[0157]

[0158] Suppose the sequence generated at one time satisfies the first-order grey differential equation as follows:

[0159]

[0160] Discretize formula (17) at the sampling interval Δt = (t + 1) - t = 1, and rewrite the in the formula as Find the best function matching by minimizing the sum of the squares of the errors and where Δx (1) = x (1) (t) - x (1) (t - 1) = x (0) (t), and the discretized equation is obtained as:

[0161] x (0) (t) + ax (1) (t) = u (18)

[0162] Then formula (18) is rewritten as:

[0163] x (0) (t) = -ax (1) (t) + u (19)

[0164] Perform adjacent mean generation on the accumulated generating sequence. Adjacent mean generation is for an equidistant time series, and new data is constructed using the average of adjacent data. Its calculation formula is:

[0165] z (1) (t) = (x (1) (t) + x (1) (t - 1)) / 2, t = 2, 3, …, n (20)

[0166] Then formula (19) is rewritten as:

[0167] x (0) (t) = -az (1) (t) + u (21)

[0168] Using the least squares method to minimize the sum of the squares of the differences between the values obtained from the fitting function and the known data, then formula (21) is written in matrix form:

[0169] Y = BU (22)

[0170] where

[0171]

[0172] Then the calculation formula for U is:

[0173] U = [a u] T=(B T B) -1 (B T Y) (23)

[0174] The parameters a and u have been obtained and can be substituted into formula (17) to obtain:

[0175]

[0176] According to formula (24), the TECR value of the next epoch is predicted

[0177] The residual test is performed on the model, and the residual test formula is:

[0178]

[0179] in, Residual ε(t)<0.2. If the residual is greater than or equal to 0.2, it is determined that the grey model is unreliable in predicting the TECR of this epoch, and the model is built again after the TECR sequence is readjusted in step S4.

[0180] The moving window GM(1,1) prediction in step S5 includes:

[0181] At epoch k, select the generated sequence data of the most recent m epochs as the modeling sample;

[0182] As the epoch advances, when the new epoch k+1 comes, the window is moved forward by one epoch: the data of the oldest epoch (k-m+1) is eliminated, and the data of the latest epoch (k+1) is introduced, keeping the sample size m unchanged;

[0183] After each window update, the grey GM (1,1) model is rebuilt based on the m samples to obtain the TECR prediction value of the new epoch, which is used for residual determination in step S6.

[0184] S6: residual error determination: calculating the residual error ε between the predicted TECR value and the corresponding measured TECR value, and when the residual error is greater than a preset threshold, recording the corresponding epoch into the second suspicious epoch set;

[0185] S7: Result fusion and output: The first suspicious epoch set and the second suspicious epoch set are fused to generate and output a final cycle slip determination result.

[0186] Considering that there are three types of satellites in BDS2 and BDS3: Geostationary Earth Orbit Satellites (GEO), Inclined Geosynchronous Orbit Satellites (IGSO), and Medium Earth Orbit Satellites (MEO), we respectively select the dual-frequency observation data of C02 (BDS2-GEO), C11 (BDS2-MEO), C13 (BDS2-IGSO), C38 (BDS3-IGSO), C43 (BDS3-MEO), and C60 (BDS3-GEO) with a sampling interval of 1 s at the URUM station from 0:00:00 to 0:14:59 on January 10, 2024. The dual-frequency observation data of 900 epochs for each satellite are used as experimental samples.

[0187] In this embodiment, a BDS cycle slip detection method based on the GM model is detailed. The GM model is used to predict the TECR value. Before prediction, the regularity of the TECR value and the feasibility of GM prediction are studied. The specific implementation process is as follows:

[0188] Calculate the combined observations of L2I and L7I / L7Z / L7D for each epoch of each satellite, that is, the TECR value of the experimental data. The results are as Figure 2 shown.

[0189] It can be seen from the figure that there is no pattern among the values after TECR combination, and the functional relationship between the data is also unknown. Although there is an internal relationship between ionospheric delays in adjacent epochs, the amount of data after TECR combination is not large and cannot be predicted by neural networks. Based on this, a grey prediction model is introduced.

[0190] Assume that x (0) (1) represents the TECR calculated value of the first epoch, then the original sequence can be obtained as: X (0) =(x (0) (1), x (0) (2), …, x (0) (n)). Then the formula for the accumulated generating sequence is:

[0191]

[0192] Thus, the accumulated TECR calculated values are as Figure 3 shown:

[0193] From Figure 3 it can be seen that the TECR accumulated values are basically linearly distributed, that is, a special exponential curve e α , that is, a new sequence of TECR accumulated values can be approximated by an expression of an exponential curve or even a straight line, and the functional expression of the fitting curve can be solved by constructing a first-order ordinary differential equation.

[0194] Solve the functional expression of the fitting curve by constructing a first-order ordinary differential equation. Let x (1)Satisfy:

[0195]

[0196] If \(a\) and \(u\) in formula (17) are known, the cumulative value of TECR at the next epoch can be predicted by solving this differential equation. However, the previous cumulative values of TECR are discrete rather than continuous. Therefore, the in the above formula is rewritten as Then, by minimizing the sum of the squares of the errors, the and \(\Delta t=(t + 1)-t = 1\), always being 1, while \(\Delta x\) (1) \(=x\) (1) (t)-x (1) (t - 1)=x (0) (t). Therefore, formula (17) can be transformed into:

[0197] x (0) (t)+ax (1) (t)=u (18)

[0198] That is:

[0199] x (0) (t)=-ax (1) (t)+u (19)

[0200] Considering that there is in the original equation, so it is more reasonable to change \(x\) (1) (t) in formula (19) to the mean value of the previous and the next moments. Next, adjacent mean generation is performed on the cumulative generation sequence. Adjacent mean generation is for an equidistant time series, and new data is constructed using the average value of adjacent data. Its calculation formula is:

[0201] z (1) (t)=(x (1) (t)+x (1) (t - 1)) / 2, t = 2, 3, …, n (20)

[0202] Then formula (19) is rewritten as:

[0203] x (0) (t)=-az (1) (t)+u (21)

[0204] Thus, the least squares method can be used to minimize the sum of the squares of the differences between the values obtained from the fitting function and the known data. Formula (21) written in matrix form is:

[0205] Y = BU (22)

[0206] Where

[0207]

[0208] Thus

[0209] U = [a u] T = (B T B) -1 (B T Y) (23)

[0210] The parameters a and u have been obtained. Substituting them into formula (17), we can get:

[0211]

[0212] According to formula (24), the TECR value of the next epoch can be predicted

[0213] Next, the residual test of the model needs to be carried out to test the rationality of the model. The residual test formula is:

[0214]

[0215] Among them,

[0216] In this paper, the residual ε(t) < 0.2 is taken. If the residual is greater than or equal to 0.2, it means that it is not very reasonable to use the grey model to predict the TECR value. This requires the ratio test of the data before starting the modeling.

[0217] Example 2: As Figure 4 , Figure 5 and Table 1 show, this example details a BDS cycle slip detection method based on the GM model. Before establishing the model, the ratio test is carried out. Only by passing the ratio test can the TECR value sequence adopt the GM(1,1) model. The specific implementation process is as follows:

[0218] Calculate the ratio parameter:

[0219]

[0220] If ω(t) is in the interval then it means that the GM(1,1) model can be used.

[0221] The ratio test of the TECR value obtained from the experimental data is as Figure 4 shown:

[0222] According to the 900 epoch data in this time period, ω(t) should be in the interval (0.9978, 1.0022), and from Figure 4 it can be concluded that this TECR value sequence cannot adopt the GM(1,1) model.

[0223] To solve this problem, we tried to perform a translation transformation on the TECR value, that is, add an arbitrary constant c to each TECR value, check if it is within the interval (0.9978, 1.0022), and then subtract c after solving using the GM(1,1) model.

[0224] Table 1: Translation constants for L2I and L7I / L7Z / L7D

[0225] Type Translation Constant Type Translation Constant C022I-7I 66 C382I-7Z 25 C112I-7I 30 C432I-7Z 27 C132I-7I 34 C602I-7D 32

[0226] Based on the data in Table 1, perform a translation transformation by adding the minimum value in Table 1 to the TECR value of each epoch, thus obtaining the ratio inspection as shown in Figure 5 shown.

[0227] From Figure 5 it can be seen that after the TECR value is translated according to Table 1, the ratio inspection result is within the interval (0.9978, 1.0022), and the GM(1,1) model can be used to predict the TECR value of the next epoch.

[0228] Example 3: As shown in Figure 6 、 Figure 7 and Table 2, this example details a grey prediction GM(1,1) model based on a moving window in a BDS cycle slip detection method based on the GM model. The specific implementation process is as follows:

[0229] The TECR value (or the translated TECR value) after the ratio inspection can adopt the GM(1,1) model. Different numbers of modeling data can be used when constructing the model, and the TECR value closer to the epoch being studied has a better prediction effect. Thus, the data for constructing the model moves over time. Add 1 cycle slip to the 80th epoch of L2I and L7I in C02, C11, and C13 respectively, add 1 cycle slip to the 80th epoch of L2I and L7Z in C38 and C43 respectively, and add 1 cycle slip to the 80th epoch of L2I and L7D in C60 respectively. The accuracy of different models for the modeling data also varies, as shown in Figure 6 shown.

[0230] Table 2: Detection values for L2I and L7I / L7Z / L7D

[0231]

[0232]

[0233] From Figure 6 and Figure 7It can be seen that for different satellites, as the amount of modeling data increases, the cycle slip detection results are either close to or deviate from the true values. However, as the amount of modeling data further increases, the cycle slip detection results remain unchanged. In this example, when GM(1,1) models are constructed using different data, the highest cycle slip detection accuracy values on L2I and L7I / 7Z / 7D are both rounded to 1 cycle, which is the same as the artificially added cycle slip detection results. The specific situation is shown in Table 2.

[0234] Example 4: As Figure 8 、 Figure 9 、shown in Table 3 and Table 4, this example details the cycle slip detection results of different BDS constellations in a BDS cycle slip detection method based on the GM model. The specific implementation process is as follows:

[0235] To verify the practicality of the GM model in cycle slip detection, different cycle slips were artificially added. Based on the experimental data, for each satellite, 9 epoch pairs with 9 groups of cycle slips were randomly selected and added at epochs 10, 70, 100, 130, 280, 310, 320, 350, and 390, adding (1,1), (-1,-1), (-1,1), (1,-1), (1,0), (0,1), (-1,0), (0,-1), and (f i / f j ,1) cycle slip pairs. For example, the cycle slip pair (1,1) means adding 1 cycle of cycle slip on L1 and 1 cycle of cycle slip on L2 respectively. These 9 groups of cycle slip pairs fully consider the cycle slips such as those in dual-frequency carrier phase observations, that is, ΔN i =ΔN j (i and j refer to different carrier phases), the frequency ratio cycle slip, that is, ΔN i / ΔN j =f i / f j , as well as 1-cycle small cycle slips. The cycle slip detection results for different constellations are as follows.

[0236] Table 3: Cycle Slip Detection Results of BDS2 Partial Satellites (L2I, L7I) and (L6I, L7I) Observations

[0237]

[0238]

[0239] Due to space limitations, Figure 8Only the cycle slip detection results of 9 groups of cycle slips of L6I and L7I for C02, C11, and C13 are shown. For clearer presentation, Table 3 lists the cycle slip detection results. If the artificially added cycle slip value is equal to the t detection value, the cycle slip detection is considered successful; otherwise, it is considered a failure. Through statistical analysis, equal-cycle cycle slips, frequency-ratio cycle slips, and small cycle slips artificially added to the dual-frequency carrier phase observations (L2I, L7I) and (L6I, L7I) of C02, C11, and C13 satellites can all be successfully detected. If the method adopted can detect small cycle slips of 1 cycle, there is no doubt that large cycle slips can also be detected. This shows that the method adopted can successfully detect equal-cycle cycle slips, frequency-ratio cycle slips, and small cycle slips of the dual-frequency observations of GEO, IGSO, and MEO satellites of BDS2 with a 1-second sampling interval.

[0240] Table 4: Cycle slip detection results of the observations of (L1X, L5X), (L1X, L6I), (L1X, L8X), (L2I, L5X), (L2I, L7Z), (L5X, L6I), (L5X, L8X), (L6I, L7Z), (L6I, L8X),

[0241] (L7Z, L8X) of C38 / C43 of BDS3 and the observations of (L2I, L7D), (L6I, L7D) of C60

[0242]

[0243] Figure 9 The cycle slip detection results of 9 groups of cycle slips of L6I and L7Z / L7D for C38, C43, and C60 are shown. For clearer presentation, Table 4 lists the percentage of the cycle slip detection results, and the data processing strategy is the same as that of BDS2. Through statistics, equal-cycle cycle slips, frequency-ratio cycle slips, and small cycle slips artificially added to the dual-frequency carrier phase observations (L1X, L5X), (L1X, L6I), (L1X, L8X), (L2I, L5X), (L2I, L7Z), (L5X, L6I), (L5X, L8X), (L6I, L7Z), (L6I, L8X), and (L7Z, L8X) of C38 and C43 satellites can all be successfully detected, and equal-cycle cycle slips, frequency-ratio cycle slips, and small cycle slips artificially added to the dual-frequency carrier phase observations (L2I, L7D) and (L6I, L7D) of C60 satellite can also be successfully detected. This shows that the method adopted can successfully detect equal-cycle cycle slips, frequency-ratio cycle slips, and small cycle slips of the dual-frequency observations of GEO, IGSO, and MEO satellites of BDS3 with a 1-second sampling interval.

[0244] Through the system verification of Embodiment 1 to Embodiment 4, it can be seen that the cycle slip detection method of "MW wide-lane combination rough detection + TECR combination fine detection + moving window GM(1,1) prediction" proposed by the present invention shows excellent performance in the dual-frequency carrier phase observation data of each constellation (GEO, IGSO, MEO) of BDS2 and BDS3. The specific conclusions are as follows:

[0245] Predictability of TECR sequence: After performing the cumulative generating operation on 900 epochs of measured data, the TECR cumulative sequence shows a linear (exponential) distribution; after ratio test and necessary translation transformation, all sequences meet the GM(1,1) modeling conditions, proving that the grey model has stable feasibility in small sample and high noise scenarios.

[0246] Cycle slip detection accuracy: In the control experiments of injecting 1 cycle slip per epoch and injecting 9 groups of equal cycle / frequency ratio / small cycle slips at random epochs, 100% of the 54 cycle slip events of six satellites were successfully detected, and the false alarm rate was 0; the GM prediction residual was always lower than the 0.2 threshold, verifying the effectiveness of the residual judgment criterion.

[0247] Applicability and robustness: The method maintains a 100% success rate for different frequency point combinations (L2I / L7I, L2I / L7Z, L2I / L7D, L1X / L5X, etc.); as the length of the moving window increases, the detection results converge rapidly and remain stable, indicating that the algorithm is not sensitive to window parameters and has good robustness and real-time performance.

[0248] Comprehensive advantages: Compared with the existing technologies that only rely on MW, GF, IF or Doppler integration, INS assistance, and neural network prediction, the present method does not require additional sensors or large-scale training samples, and can detect various types of cycle slips in real time and reliably at a 1s sampling rate, providing an efficient and easily engineering-implemented new solution for high-precision positioning, navigation and timing of BDS.

[0249] The above embodiments only represent the implementation manners of the present invention, and the description thereof is relatively specific and detailed, but it should not be construed as a limitation to the scope of the invention. It should be noted that for those of ordinary skill in the art, without departing from the concept of the present invention, several modifications and improvements can still be made, and these all belong to the protection scope of the present invention.

Claims

1. A BDS cycle slip detection method based on the GM model, characterized in that: It includes the following steps: S1: Acquisition of original observations: Construct a GNSS observation equation, collect the carrier phase observations of the target BDS satellite at the first carrier frequency and the second carrier frequency, and obtain a dual-frequency carrier phase observation sequence; S2: Coarse detection of wide-lane combination: Based on the dual-frequency carrier phase observation sequence in step S1, construct Melbourne-Wübbena combined observables, and perform the first-round cycle slip coarse determination for each epoch according to the wide-lane distance mutation threshold to obtain the first set of suspicious epochs; S3: Formation of TECR sequence: Based on the first set of suspicious epochs obtained in step S2, use the dual-frequency carrier phase observation sequence obtained in step S1 to calculate the free ionosphere combined TECR sequence; S4: Preprocessing of TECR sequence: Perform a cumulative generation operation on the TECR sequence in step S3 to obtain a first-generation sequence, and perform a ratio test on the first-generation sequence; if the ratio test is not passed, perform a translation transformation on the TECR sequence until the ratio test is passed to obtain the corrected TECR sequence and its first-generation sequence; S5: GM(1,1) prediction: In a moving window containing the most recent K epochs, establish a grey GM(1,1) model with the first-generation sequence in step S4 as the input to predict the TECR value of the next epoch; S6: Residual determination: Calculate the residual ε between the predicted TECR value and the corresponding measured TECR value; when the residual ε is greater than the preset threshold, record the corresponding epoch in the second set of suspicious epochs; S7: Result fusion and output: Fusion of the first set of suspicious epochs in step S2 and the second set of suspicious epochs in step S6 to generate and output the final cycle slip determination result.

2. The BDS cycle slip detection method based on the GM model according to claim 1, wherein: The construction of the GNSS observation equation in step S1 specifically includes the construction of the carrier phase observation equation and the construction of the pseudorange observation equation, and the calculation formulas are as follows: Construction of carrier phase observation equation: where λ i is the carrier phase observation value wavelength, i represents the number of frequencies; ρ is the distance between the satellite and the ground; c is the speed of light in vacuum; dt r and dt s are the clock errors of the receiver and the satellite respectively; ρ ion,i is the ionospheric delay; ρ trop is the tropospheric delay; N i is the integer ambiguity; is the carrier phase observation noise; Construction of pseudorange observation equation: Among them, is the pseudorange observation noise.

3. A BDS cycle slip detection method based on the GM model according to claim 2, characterized in that: The calculation of the first set of suspicious epochs in step S2 includes the following steps: S2-1: Construction of wide-lane phase combination: In the dual-frequency carrier phase observation sequence in step S1, according to formula (1), establish a wide-lane phase combination through dual-frequency GNSS observations, and the specific calculation formula is: Among them, let the carrier phase observables of frequency points f1 and f2 be respectively and their wavelengths be λ1 and λ2 respectively; S2-2: Construction of narrow-lane pseudorange combination: In the corresponding pseudorange sequence P, according to formula (2), establish a narrow-lane pseudorange combination, and the specific calculation formula is: where the pseudorange observables of frequency points f1 and f2 are P1 and P2 respectively; S2-3: Construction of MW combination: Generate MW combined observables through formula (3) and formula (4), and the specific calculation formula is: S2-4: Wide-lane normalization and cycle slip detection calculation: Normalize formula (5) to obtain the cycle slip detection quantity and obtain the combined cycle slip detection quantity, and the specific calculation formula is: wherein, the wavelength of the wide-lane combination is λ MW = c / (f1 - f2); S2-5: Wide-lane normalization and cycle slip detection calculation: Differencing the N of adjacent epochs MW to weaken multipath effects and receiver noise. The specific calculation formula is as follows: ΔN MW = N MW (t + 1)-N MW (t)(7) S2-6: Threshold determination and generation of suspicious epochs: The specific calculation formula is: ΔN MW = ΔN1 - ΔN2 (8) where ΔN1 and ΔN2 represent cycle slip values at two frequencies; when ΔN MW ≥ γMW, the corresponding epoch is recorded in the first set of suspicious epochs M; where γMW is a preset threshold.

4. A BDS cycle slip detection method based on the GM model according to claim 1, characterized in that: The formation of the TECR sequence in step S3 includes the following steps: S3-1: Calculation of total ionospheric electron content: Through two sets of carrier phase observations obtained by a dual-frequency receiver, calculate the total ionospheric electron content at a certain epoch k, and the specific calculation formula is: where b1 and b2 are the hardware delays at the receiver and satellite ends, f1 and f2 are the frequency points in the dual-frequency carrier phase observation sequence, and are the observables, λ1 and λ2 are the wavelengths, S3-2: Calculation of ionospheric change rate: The change rate of the total electron content in the ionosphere is obtained based on the difference between epochs. The specific calculation formula is: Where, Δt is the sampling interval; S3-3: Formation of cycle slip detection value: At epoch k, the cycle slip detection value is derived according to formula (9) and formula (10). The specific calculation formula is: S3-4: Auxiliary cycle slip determination, the specific calculation formula is: Among them, when no cycle slips are confirmed or cycle slips have been repaired in epoch (k–1), the suspected cycle slip value of epoch k is derived using ΔN shown in formula (12), TECR TECR and the relationship between the frequency point wavelength is used as part of the TECR sequence of the free ionosphere combination to participate in the subsequent residual determination process.

5. A BDS cycle slip detection method based on the GM model according to claim 4, characterized in that: The formation of the TECR sequence also includes the prediction process of the rate of change of the total electron content of the ionosphere, which specifically includes the following steps: At epoch k, the instantaneous rate of change is calculated using the measured TECR values ​​of epochs (k-1) and (k-2). The specific calculation formula is: The change rate obtained by formula (13) is extrapolated to epoch k to obtain the TECR prediction value. The specific calculation formula is: Substitute the predicted value of the TECR prediction into equations (8) and (12), solve the equations simultaneously to obtain ΔN1 and ΔN2, and then round ΔN1 and ΔN2 to get the cycle slip values at the two frequencies.

6. The BDS cycle slip detection method based on the GM model according to claim 1, wherein: The specific calculation formula of the level ratio test in step S4 is: If ω(t) is within the interval , it is determined that the TECR sequence can directly adopt the GM(1,1) model; If the TECR value sequence is not in the interval corresponding to ω(t) after verification, the GM(1,1) model cannot be used for this TECR value sequence; a translation transformation is performed on this TECR value, specifically, an arbitrary constant c is added to each TECR value, and then c is subtracted after the solution is obtained using the GM(1,1) model, and the level ratio test operation is repeated until the TECR value sequence is in the interval corresponding to ω(t) after verification.

7. A BDS cycle slip detection method based on the GM model according to claim 1, characterized in that: Step S5 further includes performing grey GM (1,1) modeling and feasibility test on the free ionosphere combined TECR sequence described in step S3, specifically: x (0) (1) is the calculated value of TECR for the first epoch, and the original sequence of TECR arranged in epoch order is: X (0) =(x (0) (1), x (0) (2), …, x (0) (n)), and performing an accumulative generation operation on X (0) results in a first-generation sequence: Assume that the generated sequence satisfies the first-order grey differential equation: Discretize Equation (17) at the sampling interval Δt = (t + 1) - t = 1, and rewrite the in Equation (17) as Find the best function match by minimizing the sum of the squares of the errors of and The discretized equation is obtained as follows: x (0) (t) + ax (1) (t) = u (18) where, Δx (1) = x (1) (t) - x (1) (t - 1) = x (0) (t), then formula (18) is rewritten as: x (0) x(t) = -ax (1) x(t) + u (19) Generate the neighboring mean of the cumulatively generated sequence, and construct new data using the average value of adjacent data. The calculation formula is: z (1) (t) = (x (1) (t) + x (1) (t - 1)) / 2, t = 2, 3, …, n (20) Then formula (19) can be rewritten as: x (0) z(t) = -az (1) z(t) + u (21) The least squares method is used to minimize the square difference between the value of the fitting function and the known data, and formula (21) is written in matrix form: Y=BU (22)where, Then the calculation formula of U is: U = [a u] T = (B T B) -1 (B T Y) (23) Calculate the parameters a and U and substitute them into formula (17) to obtain: Predict the TECR value of the next epoch according to formula (24) The residual test is performed on the model, and the residual test formula is: Among them, the residual ε(t) < 0.2; if the residual is greater than or equal to 0.2, it is determined that the gray model's prediction of TECR for this epoch is unreliable, and the process returns to step S4 to readjust the TECR sequence and then perform modeling again.

8. A BDS cycle slip detection method based on the GM model according to claim 1, characterized in that: The moving window GM(1,1) prediction in step S5 includes: At epoch k, select the generated sequence data of the most recent m epochs as the modeling sample; As the epoch advances, when the new epoch k+1 comes, the window is moved forward by one epoch: the data of the oldest epoch (k-m+1) is eliminated, and the data of the latest epoch (k+1) is introduced, keeping the sample size m unchanged; After each window update, the grey GM (1,1) model is rebuilt based on the m samples to obtain the TECR prediction value of the new epoch, which is used for residual determination in step S6.

Images