BDS cycle slip repairing method based on wavelet and generalized extrapolation

Abstract

The invention provides a BDS cycle slip repair method based on wavelet and generalized extrapolation, which comprises the following steps: reading carrier data and constructing a single-difference detection quantity, and obtaining a high-frequency coefficient by using wavelet transform decomposition to judge cycle slip; and when cycle slip occurs, extrapolation is carried out by using a historical correct carrier phase observation value, and an extrapolation value is compared with an actual value and interpolation restoration is carried out. Experiments combining a generalized continuation extrapolation method and a polynomial fitting extrapolation method show that the generalized continuation extrapolation method and the polynomial fitting extrapolation method are better in precision and stability, and the single-frequency BDS cycle slip detection and repair capability can be remarkably improved.

Description

Technical Field

[0001] The invention belongs to the technical field of Beidou navigation and positioning data processing, and in particular relates to a BDS cycle slip repair method based on wavelet and generalized extrapolation. Background Art

[0002] The BeiDou satellite system has been widely used in many fields such as navigation and positioning. At present, the system has realized the function of broadcasting multi-frequency satellite signals. Through multi-frequency receivers, it can receive signals at multiple frequencies, thereby performing more precise data processing. However, due to many factors such as receiver cost, power consumption and size, single-frequency receivers still have important advantages in the mass market. When using single-frequency receivers for precise single-point positioning, the continuity of carrier phase observations is crucial, otherwise, once a cycle slip occurs, the positioning accuracy will be significantly reduced. Single-frequency receivers are not like multi-frequency receivers that can use multi-frequency signals to form effective combined observations. When performing cycle slip preprocessing on single-frequency observations, only the traditional single-frequency cycle slip detection and repair method can be used. Therefore, the study of cycle slips of single-frequency observations should be further deepened.

[0003] In view of the impact of cycle slips on positioning accuracy, domestic and foreign scholars have studied many methods to avoid errors caused by cycle slips on positioning, including the ionospheric phase delay method, the application of Kalman filtering technology in cycle slip research, and the polynomial fitting method. The ionospheric phase delay method is only applicable to specific observation environments where the ionospheric activity is calm and the residuals between adjacent epochs change slowly; the Kalman filtering technology is applied to cycle slip research. This method compares the carrier phase measurement value with the estimated value obtained by the dynamic model and Kalman filtering. If the difference is obvious, it indicates that a cycle slip has occurred. However, its performance depends to a large extent on the selection of Kalman filtering statistical parameters, and the determination of these parameters often requires rich engineering experience as support; the polynomial fitting method determines the occurrence of cycle slips by analyzing the residuals between the observation value sequence and the fitting polynomial. This method performs well in dealing with large cycle slips and lays an important foundation for subsequent research, but its ability to identify small cycle slips is still insufficient.

[0004] In recent years, many scholars have improved the traditional methods. These methods mainly include single-frequency cycle slip detection and repair methods based on Doppler residual and inertial assistance, ambiguity resolution and cycle slip processing methods using inertial navigation assistance, methods based on unmodeled error constraints and geometric model-free (GF) combined with geometric model (GB), and a new single-frequency cycle slip detection method that uses a dual-board design method based on a shared antenna and uses two circuit boards to provide two independent phase measurement values ​​for each satellite at the same frequency.

[0005] A single-frequency cycle slip detection and repair method based on Doppler residual and inertial assistance uses a micro-electromechanical system inertial measurement unit (MEMS-IMU) for tight coupling integration to accurately predict the Doppler measurement value of the receiver for cycle slip detection and repair. If the equipment itself has noise and drift errors, it will affect the accuracy of Doppler prediction; the ambiguity resolution and cycle slip processing method assisted by inertial navigation is used to perform ambiguity resolution and satellite selection, and the short-term high-precision position predicted by INS is used as the constraint equation for cycle slip detection. The unit weighted mean square error of the cycle slip detection equation is used to determine the epoch of the cycle slip, and the residual calculated by robust estimation is used for cycle slip research. However, the system performance of this method depends largely on the accuracy level of INS; based on unmodeled error constraints A method without geometric model (GF) combined with geometric model (GB) is used to accurately predict the short-term GF and GB unmodeled errors, and the predicted GF and GB unmodeled errors are used to correct the GF and GB models that may have deviations. Finally, the cycle slips are detected and repaired by combining the optimized GF and GB models, but this method does not highlight the role of the latest data in the prediction process; a new single-frequency cycle slip detection method is adopted with a dual-board design based on a shared antenna, using two circuit boards to provide two independent phase measurement values ​​for each satellite at the same frequency. Based on the consistency criterion of the measurement values, triple differences between boards, epochs and satellites are constructed as statistics to achieve cycle slip detection and repair, but it requires the support of special hardware equipment, which will increase the cost of the hardware equipment.

[0006] RINEX is a well-known universal data exchange format for satellite navigation receivers, proposed and maintained by the International GNSS Service (IGS), and is compatible with satellite navigation receiver data produced by various manufacturers. By parsing the RINEX observation file, complete satellite observation information including carrier phase, pseudorange, Doppler frequency shift, etc. can be obtained, providing important data source and format support for subsequent data processing, analysis and precision calculation. However, the carrier observation values ​​in the RINEX data may have cycle slips in high-precision positioning applications, especially when the signal is blocked, interfered, or the receiver switches satellites. Traditional cycle slip detection and repair methods include signal processing methods such as carrier phase combination method, carrier phase and pseudorange combination method, and empirical mode decomposition (EMD). However, in the case of complex noise environment, high sampling rate and huge amount of observation data, traditional methods are often interfered by noise when extracting mutation information, resulting in detection efficiency and accuracy that are difficult to meet the requirements.

[0007] In the scenario of cycle slip detection, one of the commonly used methods is "single difference detection between epochs", that is, the difference between the carrier phase observation values ​​of the same satellite in two adjacent observation epochs. Ideally, since the satellite motion and user motion between adjacent epochs change little in a short period of time, and the clock error change is effectively weakened or completely eliminated in the single difference, the single difference between epochs should show a relatively stable change characteristic. If the observation value has a sudden change of integer wavelength (cycle slip), it will cause the single difference between epochs to produce obvious abnormal deviation, so that the abnormal position can be quickly identified. Compared with higher-order differences (such as double differences, triple differences, etc.), the construction of single difference detection between epochs is simpler and the calculation time is lower. It is suitable for cycle slip detection in real-time or quasi-real-time environments with high sampling rates. In conventional GNSS data processing, single difference detection between epochs is usually combined with threshold methods or other statistical test methods to judge whether a cycle slip has occurred by comparing indicators such as the change amplitude of the single difference and the residual distribution.

[0008] As a multi-resolution analysis tool, wavelet analysis has unique advantages in mutation point detection, singularity analysis and signal denoising due to its good time-frequency localization characteristics. The selection of wavelet functions is crucial to the analysis of different types of signals. Common and well-known types of wavelets include: Haar wavelet (Haar): the earliest discrete wavelet proposed, with the advantages of simplicity and low computational complexity, suitable for scenes that are most sensitive to mutation responses, and can well detect jump points in signals; Daubechies wavelet (Db): proposed by Ingrid Daubechies proposed that wavelets of different orders with characteristics such as compact support and orthogonality can take into account both smoothness and localization capabilities, and are suitable for detailed analysis of various signals; Symlets wavelet (Sym): is an improved version of the Daubechies wavelet, with approximate symmetry characteristics, reduced phase distortion, and is often used in situations where signal symmetry or low phase distortion is required; Coiflets wavelet (Coif): is another compactly supported wavelet with both high vanishing moment and good orthogonality, and shows good accuracy in singularity detection and feature extraction; Biorthogonal wavelet (Bior): has a biorthogonal basis, and different wavelets can be used in the decomposition and reconstruction stages respectively, allowing more flexible signal processing while having both linear phase characteristics and symmetry.

[0009] Although the existing methods have achieved certain results, they still face problems such as hardware dependence or high cost in practical applications. The cycle slip problem has always been a difficult problem in the field of satellite navigation when processing GNSS data. Wavelet transform can reflect the localized signal characteristics of the signal, and the cycle slip itself is a signal mutation. In the field of cycle slip detection, wavelet transform itself has its own unique advantages. Since the high-frequency coefficients after wavelet transform decomposition can only reflect the location where the satellite carrier signal has a mutation, it cannot reflect the specific size of the cycle slip value. Therefore, most scholars combine wavelet transform with other methods to detect cycle slips. Using a limited number of error-free data to extrapolate the cycle slip is a common method to solve the cycle slip problem, and the generalized extension extrapolation method has obvious advantages over the traditional polynomial fitting method in terms of accuracy and stability, which has opened up new research ideas for effectively improving the ability to detect cycle slips. Summary of the invention

[0010] In view of the problems of the prior art in terms of cycle slip detection accuracy and repair stability, the present invention proposes a BDS cycle slip repair method based on wavelet and generalized extrapolation, which uses wavelet transform to identify the mutation point of the signal and combines the generalized extension extrapolation method to accurately repair the cycle slip value, thereby effectively improving the accuracy of cycle slip detection and the stability of repair, and solving the shortcomings of traditional methods in high noise environments. The method includes the following steps:

[0011] S1: Read carrier data, where the carrier data includes carrier phase observation values ​​of the Beidou satellite system;

[0012] S2: construct a single difference detection quantity according to the carrier phase observation value obtained in step S1;

[0013] S3: Decomposing the single difference detection amount by wavelet transform, reconstructing and obtaining the first layer high frequency coefficient, specifically comprising the following steps:

[0014] The signal f(t) is subjected to translation and scaling operations using the wavelet function family generated by the mother wavelet ψ(t) to perform continuous wavelet transform. The specific calculation formula is:

[0015]

[0016] Where a represents the scaling parameter, b represents the translation parameter, (a, b∈R and a≠0); ψ*(t) is the complex conjugate of ψ(t);

[0017] The specific calculation formula for signal reconstruction is:

[0018]

[0019] Among them, the constant C ψ The following conditions are met:

[0020]

[0021] Among them, Ψ(ω) is the Fourier transform of ψ(t), and ω is the signal frequency;

[0022] S4: Use the first-layer high-frequency coefficients obtained in step S3 to determine whether cycle slips occur in the carrier data through the local singularity detection characteristics of wavelet transform, which specifically includes the following steps:

[0023] At scale s 0 , if there exists a point (s 0 , τ 0 ) that satisfies the local limit condition, then the point (s 0 , τ 0 ) is a local extreme point, and the local limit condition is where W f (s 0 , τ 0 ) is the wavelet transform result of the signal at scale s 0 and time τ 0 ;

[0024] If there exists a zero-crossing point at τ 0 , and for any point τ in a certain neighborhood of τ 0 , there is a maximum condition satisfied, then the point (s 0 , τ 0 ) is a modulus maximum point of the wavelet transform, and the maximum condition is |W f (s 0 , τ)| ≤ |W f (s 0 , τ 0 )|;

[0025] Use the Lipschitz exponent to extract the singularity characteristics of the signal. For the signal function f(x), let n be a non-integer (n ≤ a < n + 1). If there exist constants A and h 0 (h 0 > 0) and an nth-degree polynomial P n (h), such that for any h ≤ h 0 the following calculation conditions are satisfied:

[0026] |f(x 0 + h) - P n (h)| ≤ A|h| a (5)

[0027] then f(x) has a Lipschitz of α at point x 0 . Specifically, if the above conditions are satisfied for any x 0 in the interval (a, b) 0+h∈(a,b), then f(x) is uniformly Lipschitzα on the interval (a,b);

[0028] By performing multi-scale analysis on the signal, when the signal has a mutation, its wavelet transformed coefficient has the modulus maximum feature, and by detecting the modulus maximum point, it is determined whether the carrier data has a cycle slip.

[0029] S5: When step S4 determines that there is a cycle slip in the carrier data, the following method is used to repair the cycle slip:

[0030] S5-1: Select the Haar wavelet with the most significant response to mutations, decompose and reconstruct the satellite data according to different sampling rates, and select the correct carrier phase observation values ​​in several epochs before the cycle slip occurs;

[0031] S5-2: Use the generalized extension extrapolation method to perform extrapolation repair, including:

[0032] Assume that the known data sequence is {(t 1 , x 1 ), (t 2 , x 2 ), (t 3 , x 3 )…(t n , x n )}, need to extrapolate t n+1 The extrapolated value at time i∈[1, n], when the data points are dense, the average value of the latest point m is selected as the constraint condition to establish a generalized extension model:

[0033]

[0034] Where: (a 0 , a 1 , a 2 , a 3 ) is the constant coefficient to be solved, min I(a 0 , a 1 , a 2 , a 3 ) is the optimal objective function to be minimized, and the difference between the extrapolated result and the actual carrier phase observation value is used to repair the cycle slip, thereby achieving accurate detection and repair of cycle slips in single-frequency BDS carrier data.

[0035] As a technical preferred solution of the present invention, the reading of carrier data in step S1 includes parsing a Rinex observation file to obtain the carrier phase observation value.

[0036] As a technical preferred solution of the present invention, the construction of the single-difference detection amount described in step S2 adopts the single-difference detection amount between epochs.

[0037] As a technical preferred solution of the present invention, step S3 specifically includes: using five different wavelet functions, namely Haar, Db, Sym, Coif and Bior, to perform wavelet transform decomposition and reconstruction on the single difference detection quantity of the carrier data, obtain the first layer of high-frequency coefficients, and evaluate the response characteristics of the above five wavelets to the mutation signal.

[0038] As a technical preferred solution of the present invention, the wavelet transform decomposition described in step S3 is accelerated by using discrete wavelet transform, which specifically includes the following steps:

[0039] Select the scaling parameter and translation parameter in binary form, and the calculation formula is:

[0040] a=2 m (7)

[0041] b=n2 m (8)

[0042] where m and n are integers; for a particular choice of ψ(t), there exists a corresponding discrete wavelet ψ that satisfies the good time-frequency localization condition n，n , the calculation formula for good time-frequency localization conditions is:

[0043] ψ m，n (t) = 2 -m / 2 ψ(2 -m tn) (9)

[0044] Discrete wavelet is the structure of L 2 (R) is an orthogonal basis. Using this orthogonal basis, L 2 The calculation formula for any f(t) in (R) space is:

[0045]

[0046] Among them, the discrete wavelet coefficient a m，n The calculation formula is:

[0047]

[0048] Among them, the signal can be decomposed into an approximate part and a detail part, and the detail calculation formula of the M layer is:

[0049]

[0050] Among them, U is a set of positive integers, and the approximation of the M layer is defined as the sum of all details before this layer, and its calculation formula is:

[0051]

[0052] Among them, the calculation formula of the signal f(t) is:

[0053]

[0054] Formula (14) further yields:

[0055] A M-1 (t) = A M (t)+D M (t) (15)

[0056] Through the above signal reconstruction, detail components are obtained at different scales to detect and analyze local feature changes in the signal and achieve accurate detection of cycle slips.

[0057] As a technical preferred solution of the present invention, step S4 determines whether a cycle slip occurs in the carrier data by detecting a modulus maximum point, which specifically includes the following steps:

[0058] The multi-scale analysis of the signal is performed by using wavelet transform. By detecting the local singularity of the signal, when the signal changes suddenly, the coefficients after wavelet transform have the characteristics of modulus maximum value.

[0059] Perform modulus maximum point detection on the signal coefficients after wavelet transformation, and judge whether there is abnormality or mutation in the signal by identifying the modulus maximum point;

[0060] Combining the concept of Lipschitz index, when the signal is x 0 When the Lipschitz index α of a point is less than 1, the point is considered to be a singular point of the signal, and further confirms whether the signal has a cycle slip;

[0061] According to the local characteristics of the coefficients after wavelet transformation, the mutation points in the signal are detected, so as to effectively identify whether the cycle slip phenomenon occurs in the carrier data.

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

[0063] Accurate and reliable cycle slip detection, using the sensitivity of wavelet transform to the high-frequency components of the signal, can quickly locate the mutation points in the carrier phase at multiple scales, thereby accurately identifying the epochs where cycle slips may occur and reducing the misjudgment rate caused by noise or other interference.

[0064] High-precision cycle slip repair uses a generalized extension extrapolation method to predict and repair the cycle slip value based on the correct carrier phase observations within several epochs before the cycle slip occurs. Compared with the traditional polynomial fitting method, this extrapolation method has more advantages in prediction accuracy and stability, and effectively reduces the accumulation of repair errors.

[0065] No additional hardware support is required. The method of the present invention only relies on software algorithm implementation and has no additional requirements for receiver hardware. It can be directly applied on the basis of existing GNSS receiving equipment and can greatly reduce system cost and deployment difficulty.

[0066] It is adaptable to complex environments and has high reliability. Wavelet analysis has good performance in time-frequency localization and can effectively cope with high noise and dynamic environments. It improves the ability to capture weak mutation signals through multi-scale decomposition. Combined with generalized extension and extrapolation, it can enhance the stability of the repair results and achieve robust detection and repair of cycle slips in different scenarios. BRIEF DESCRIPTION OF THE DRAWINGS

[0067] Figure 1 It is a flow chart of a BDS cycle slip repair method based on wavelet and generalized extrapolation of the present invention;

[0068] Figure 2 It is an experimental data diagram of an embodiment provided by the present invention;

[0069] Figure 3 is a wavelet decomposition diagram of a single difference detection quantity according to an embodiment of the present invention, wherein a3 is an approximation coefficient, and d1, d2, and d3 are detail coefficients;

[0070] Figure 4 is a wavelet decomposition diagram of the original data of the embodiment provided by the present invention, wherein a3 is an approximation coefficient, and d1, d2, and d3 are detail coefficients;

[0071] Figure 5 It is a detection effect diagram of different wavelets provided in the embodiment of the present invention;

[0072] Figure 6 It is a first-layer reconstructed signal diagram of single-difference detection quantities with different sampling rates according to an embodiment of the present invention;

[0073] Figure 7 It is a difference diagram between the true value and the extrapolated value of the carrier phase according to the embodiment provided by the present invention. DETAILED DESCRIPTION

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

[0075] Example 1: According to Figure 1 As shown, this embodiment provides a specific implementation process of a BDS cycle slip repair method based on wavelet and generalized extrapolation. The steps are as follows:

[0076] S1: Read carrier data, where the carrier data includes carrier phase observation values ​​of the Beidou satellite system; wherein reading carrier data includes parsing a Rinex observation file to obtain the carrier phase observation values.

[0077] S2: construct a single-difference detection value according to the carrier phase observation value obtained in step S1; wherein, the single-difference detection value constructed uses the single-difference detection value between epochs.

[0078] S3: Decomposing the single difference detection amount by wavelet transform, reconstructing and obtaining the first layer high frequency coefficient, specifically comprising the following steps:

[0079] The signal f(t) is subjected to translation and scaling operations using the wavelet function family generated by the mother wavelet ψ(t) to perform continuous wavelet transform. The specific calculation formula is:

[0080]

[0081] Where a represents the scaling parameter, b represents the translation parameter, (a, b∈R and a≠0); ψ*(t) is the complex conjugate of ψ(t);

[0082] The specific calculation formula for signal reconstruction is:

[0083]

[0084] Among them, the constant C ψ The following conditions are met:

[0085]

[0086] in, is the Fourier transform of ψ(t), ω is the signal frequency;

[0087] Five different wavelet functions, Haar, Db, Sym, Coif and Bior, are used to decompose and reconstruct the single difference detection of the carrier data to obtain the first layer of high frequency coefficients, and the response characteristics of the above five wavelets to the mutation signal are evaluated.

[0088] The wavelet transform decomposition is accelerated by using discrete wavelet transform, which specifically includes the following steps:

[0089] Select the scaling parameter and translation parameter in binary form, and the calculation formula is:

[0090] a=2 m (7)

[0091] b=n2 m (8)

[0092] where m and n are integers; for a particular choice of ψ(t), there exists a corresponding discrete wavelet ψ that satisfies the good time-frequency localization condition m，n , the calculation formula for good time-frequency localization conditions is:

[0093] ψ m，n (t) = 2 -m / 2 ψ(2 -mt -n) (9)

[0094] Discrete wavelet is the structure of L 2 (R) is an orthogonal basis. Using this orthogonal basis, L 2 The calculation formula for any f(t) in (R) space is:

[0095]

[0096] Among them, the discrete wavelet coefficient a m，n The calculation formula is:

[0097]

[0098] Among them, the signal can be decomposed into an approximate part and a detail part, and the detail calculation formula of the M layer is:

[0099]

[0100] Among them, U is a set of positive integers, and the approximation of the M layer is defined as the sum of all details before this layer, and its calculation formula is:

[0101]

[0102] Among them, the calculation formula of the signal f(t) is:

[0103]

[0104] Formula (14) further yields:

[0105] A M-1 (t) = A M (t)+D M (t) (15)

[0106] Through the above signal reconstruction, detail components are obtained at different scales to detect and analyze local feature changes in the signal and achieve accurate detection of cycle slips.

[0107] S4: using the first layer high frequency coefficients obtained in step S3, judging whether the carrier data has cycle slips by using the local singularity detection characteristics of wavelet transform, specifically comprising the following steps:

[0108] In scale s 0 If there is a point (s0 If τ 0 satisfies the local limit condition, then the point (s 0 , τ 0 ) is a local extreme point, and the local limit condition is where W f (s 0 , τ 0 ) is the result of the wavelet transform of the signal at scale s 0 and time τ 0 ;

[0109] If at τ 0 , there is a zero-crossing point, and for any point τ in a certain neighborhood of τ 0 , the maximum value condition is satisfied, then the point (s 0 , τ 0 ) is a modulus maximum point of the wavelet transform, and the maximum value condition is |W f (s 0 , τ)| ≤ |W f (s 0 , τ 0 )|;

[0110] The Lipschitz exponent is used to extract the singularity characteristics of the signal. For the signal function f(x), let n be a non-integer (n ≤ a < n + 1). If there exist constants A and h 0 (h 0 > 0) and an nth-degree polynomial P n (h), such that for any h ≤ h 0 the following calculation conditions are satisfied:

[0111] |f(x 0 + h) - P n (h)| ≤ A|h| a (5)

[0112] then f(x) is Lipschitz α at the point x 0 . Specifically, if the above conditions are satisfied for any x 0 in the interval (a, b) such that x 0 + h ∈ (a, b), then f(x) is uniformly Lipschitz α in the interval (a, b);

[0113] By performing multi-scale analysis on the signal, when the signal has a mutation, the coefficients after its wavelet transform have modulus maximum characteristics. By detecting the modulus maximum points, it is judged whether the carrier data has a cycle slip; among them, the detection steps of the cycle slip include:

[0114] The multi-scale analysis of the signal is performed by using wavelet transform. By detecting the local singularity of the signal, when the signal changes suddenly, the coefficients after wavelet transform have the characteristics of modulus maximum value.

[0115] Perform modulus maximum point detection on the signal coefficients after wavelet transformation, and judge whether there is abnormality or mutation in the signal by identifying the modulus maximum point;

[0116] Combining the concept of Lipschitz index, when the signal is x 0 When the Lipschitz index α of a point is less than 1, the point is considered to be a singular point of the signal, and further confirms whether the signal has a cycle slip;

[0117] According to the local characteristics of the coefficients after wavelet transformation, the mutation points in the signal are detected, so as to effectively identify whether the cycle slip phenomenon occurs in the carrier data.

[0118] S5: When step S4 determines that there is a cycle slip in the carrier data, the following method is used to repair the cycle slip:

[0119] S5-1: Select the Haar wavelet with the most significant response to mutations, decompose and reconstruct the satellite data according to different sampling rates, and select the correct carrier phase observation values ​​in several epochs before the cycle slip occurs;

[0120] S5-2: Use the generalized extension extrapolation method to perform extrapolation repair, including:

[0121] Assume that the known data sequence is {(t 1 , x 1 ), (t 2 , x 2 ), (t 3 , x 3 )…(t n , x n )}, need to extrapolate t n+1 The extrapolated value at time i∈[1, n], when the data points are dense, the average value of the latest point m is selected as the constraint condition to establish a generalized extension model:

[0122]

[0123] Where: (a 0 , a 1 , a 2 , a 3 ) is the constant coefficient to be solved, min I(a 0 , a 1 , a 2 , a 3) is the optimal objective function to be minimized, and the difference between the extrapolated result and the actual carrier phase observation value is used to repair the cycle slip, thereby achieving accurate detection and repair of cycle slips in single-frequency BDS carrier data.

[0124] Example 2: In order to verify the beneficial effects of the method of the present invention, this experiment is based on the GNSS observation data published on the CDDIS website. The carrier phase observation values ​​of the L2 frequency point of the C01 satellite collected by the URUM station are selected as the research object, and the observation data of 600 epochs are continuously extracted as experimental samples.

[0125] like Figures 2 to 4 As shown, this embodiment shows in detail a single-frequency BDS cycle slip detection and repair method based on wavelet transform combined with generalized extension extrapolation method, in which a single difference detection amount is constructed according to the carrier phase observation value, and the specific implementation process is:

[0126] The observation data with a sampling rate of 1s is used to extract the carrier phase observation values ​​of 600 consecutive epochs. In order to verify the detection ability of the algorithm, a cycle slip of one week is artificially added at the 300th epoch. The experiment uses a variety of wavelet functions to perform three-layer wavelet decomposition on the single-difference observation signal between epochs, and reconstructs the obtained low-frequency and high-frequency components. The focus is on analyzing the cycle slip detection effect after the reconstruction of the first layer of high-frequency coefficients. Figure 2 The comparison between the original observation data and the constructed single-difference detection between epochs is shown. Taking the db4 wavelet as an example, the two types of data are decomposed by three layers of wavelet respectively. The results are as follows Figure 3 , Figure 4 As shown in the figure, it can be seen from the comparison that due to the large value of the original data of the carrier phase, significant edge effects will be generated at both ends during the signal decomposition process, which will interfere with the precise positioning of the cycle slip location. By constructing the single difference detection quantity between epochs, the edge effects at both ends of the signal can be effectively eliminated, and the reliability of cycle slip detection can be improved.

[0127] This experiment uses five different wavelet functions, bior3.5, coif4, db4, haar and sym4, to perform wavelet analysis experiments on BDS carrier phase observation data with a sampling rate of 1s. Since signal mutation characteristics such as cycle slips and noise are mainly reflected in high-frequency components, the analysis of the high-frequency part is of great significance for discovering the location of cycle slips. Figure 5 Table 1 and Table 2 are the comparison results of different wavelets. Through comparative analysis, it is found that when different wavelet functions are used to decompose and reconstruct the same data sequence, significant "burr" mutations appear in the high-frequency components. These mutations are mainly caused by various error sources such as ionospheric influence, tropospheric delay, and multipath effect. Figure 5It can be clearly observed from the distribution of different wavelet coefficients that at the moment when the signal changes suddenly, the wavelet transform coefficients show significant modulus maximum characteristics. This feature provides a reliable basis for accurately locating the specific epoch when the cycle slip occurs. By setting a suitable threshold, these modulus maximum points can be effectively identified, thereby realizing automatic detection of cycle slips.

[0128] Table 1: First layer reconstruction coefficients of different wavelets

[0129]

[0130] Analysis of the experimental results in Table 1 shows that different wavelet functions exhibit significantly different response characteristics in cycle slip detection. By analyzing the experimental results, it was found that the modulus maximum appeared at 300 epochs. Although all types of wavelet functions can detect this feature, they show obvious differences in response patterns. In particular, in terms of detection sensitivity, the three wavelets Bior3.5, DB4 and Haar performed well. It was observed that Bior3.5 and DB4 produced significant coefficient mutations at 299 and 300 epochs, and quickly decayed to zero at 301 epochs. This feature makes them very suitable for accurately locating the occurrence of cycle slips.

[0131] The haar wavelet presents a unique characteristic. It not only produces the strongest mutation response when a cycle slip occurs, but also maintains strong oscillation in the next epoch. This reflection feature is particularly helpful for the initial identification of cycle slips. In comparison, the responses of coif4 and sym4 are relatively mild. Although they are not as sensitive to mutation signals as the first three wavelets, they can still effectively identify a cycle slip of 1 week. The characteristics of these two types of wavelets are good smoothness and relatively gentle fluctuations.

[0132] In actual cycle slip detection applications, it is found that wavelet functions with stronger mutation response capabilities often achieve better results. Specifically, Bior3.5, DB4, and Haar perform best in capturing signal mutation characteristics. Although Coif4 and Sym4 perform smoothly overall, their relatively weak fluctuation amplitudes may affect the reliability of detection, which is particularly evident when dealing with small cycle slips. Therefore, when choosing a wavelet function, it is necessary to decide based on the specific application scenario. If you pay special attention to the significance of the amplitude response, then Bior3.5, DB4, or Haar would be a better choice.

[0133] Example 3: Figure 6 As shown in Table 2, this embodiment shows in detail how to implement the influence of different sampling rates on the cycle slip detection accuracy in a single-frequency BDS cycle slip detection and repair method based on wavelet transform combined with generalized extension extrapolation method. The specific implementation process is as follows:

[0134] The gfzrnx software is used to extract 600 epochs of data with a sampling rate of 1s into data sets of 200 epochs with a sampling rate of 3s and 120 epochs with a sampling rate of 5s. The Haar wavelet transform is used to decompose and reconstruct the single-difference carrier phase observations between epochs with different sampling rates. The experimental results show that Figure 6 It can be seen from Table 2 that the sampling rate has a significant impact on the accuracy of cycle slip detection. Through systematic experiments, it is found that the cycle slip detection performance shows an obvious positive correlation with the sampling rate: the higher the sampling rate, the better the cycle slip detection effect, and as the sampling rate decreases, the detection accuracy gradually decreases. Specifically, under the condition of 1s sampling rate, the system can accurately locate and identify cycle slips greater than or equal to 1 cycle; when the sampling rate drops to 3s and 5s, the system's detection capability drops to 2 cycles and 3 cycles, respectively. This phenomenon is mainly due to the relationship between the sampling rate and the phase change between adjacent epochs: under high sampling rate conditions, the phase change between adjacent epochs is small, and the discontinuity characteristics caused by cycle slips are more prominent, which is conducive to the identification and positioning of cycle slips; on the contrary, under low sampling rate conditions, there may be large phase changes between adjacent epochs, making it difficult for the system to effectively distinguish between normal phase changes and cycle slips.

[0135] Table 2: Cycle slip detection based on wavelet transform combined with generalized extension extrapolation at different sampling rates

[0136]

[0137] like Figure 7 As shown in Table 3 and Table 4, this embodiment shows in detail how to verify the accuracy of repairing cycle slips by the preferred generalized extension method in a single-frequency BDS cycle slip detection and repair method based on wavelet transform combined with generalized extension extrapolation. The specific implementation process is:

[0138] Although the existence of cycle slips can be identified by singular values ​​at the signal mutation points, the size of the cycle slips cannot be detected, which makes it difficult to directly repair the cycle slips. Therefore, wavelet transform needs to be combined with other methods to determine the size of the cycle slips using the generalized extension extrapolation method. This method is an extension of the generalized extension interpolation theory in the application of extrapolation and has significant advantages. Compared with the general interpolation model, the generalized extension extrapolation method can use more prior data, which makes it have better predictions than the general interpolation model. This feature enables the model to better reflect the law of signal changes. It also highlights the importance of the latest data. This method can obtain more accurate estimation results when determining the size of the cycle slip, laying the foundation for accurate detection and reliable repair of the cycle slip.

[0139] By comparing the performance of the polynomial fitting extrapolation method using 20 prior data and the generalized extension extrapolation method, the detection accuracy of the cycle slip size is evaluated by rounding the extrapolation results to one decimal place. The results in Table 3 are rounded to one decimal place to obtain Table 4. The experimental results are shown in Figure 7 , Table 3 and Table 4 show that the two extrapolation methods are stably compared. Figure 7 It can be seen that the black curve represents the gap between the generalized extension extrapolation value and the true value, and the red curve represents the gap between the polynomial fitting extrapolation value and the true value. When the sampling rate is 1 second, the standard deviation of the generalized extension extrapolation method is 0.18, while the standard deviation of the polynomial fitting extrapolation method is 0.258. Therefore, the accuracy of the generalized extension extrapolation method is improved by 30.23% compared with the polynomial fitting extrapolation method. The smaller standard deviation indicates that the data fluctuation of the generalized extension extrapolation method is more stable, and the prediction trend is more stable.

[0140] Table 3: Comparison of cycle slip repair accuracy using different extrapolation methods

[0141]

[0142] Table 4: Comparison of cycle slip repair results by different extrapolation methods

[0143]

[0144] The results in Tables 3 and 4 show that the generalized extension extrapolation method exhibits better performance characteristics than the polynomial fitting extrapolation method. Under the sampling rate of 1s, the repair range of the generalized extension extrapolation method for 1 cycle jump is (0.5908, 1.1693), and all 1 cycle jumps are successfully detected. The repair range of the polynomial fitting extrapolation method is (0.3168, 1.4688), and a misjudgment occurs at the 250th epoch, failing to effectively repair the 1-cycle small cycle jump. When the sampling rate is reduced to 3s and 5s, although the error values ​​of both methods gradually increase, the generalized extension extrapolation method exhibits more stable numerical characteristics, and its prediction results fluctuate less. The study found that at a sampling rate of 1s, this method predicts the cycle jump epoch best. Even under the condition of reduced sampling rate, its prediction results are more accurate than polynomial fitting.

[0145] Through comparative analysis, the method combining wavelet transform and generalized extension extrapolation has obvious advantages in cycle slip detection. Experiments show that this method can successfully detect and repair a small cycle slip of 1 cycle under the condition of 1s sampling rate; even under the condition of 5s sampling rate, it still maintains the position of detecting 3 cycles and above cycle slips, and can repair most cycle slips. Compared with the traditional combination of wavelet transform and polynomial fitting extrapolation, the generalized extension extrapolation method significantly improves the accuracy of cycle slip detection and shows stronger robustness.

[0146] The above embodiments only express several implementation methods of the present invention, and the description is relatively specific and detailed, but it cannot be understood as limiting the scope of the invention. It should be pointed out that for ordinary technicians in this field, several modifications and improvements can be made without departing from the concept of the present invention, which all belong to the protection scope of the present invention.

Claims

1. A BDS cycle slip repair method based on wavelet and generalized extrapolation, characterized by: The method comprises the following steps: S1: Read carrier data, where the carrier data includes carrier phase observation values ​​of the Beidou satellite system; S2: construct a single difference detection quantity according to the carrier phase observation value obtained in step S1; S3: Decomposing the single difference detection amount by wavelet transform, reconstructing and obtaining the first layer high frequency coefficient, specifically comprising the following steps: The signal f(t) is subjected to translation and scaling operations using the wavelet function family generated by the mother wavelet ψ(t) to perform continuous wavelet transform. The specific calculation formula is: Among them, W f represents the transformed wavelet, a represents the scaling parameter, b represents the translation parameter, (a, b∈R and a≠0); ψ*(t) is the complex conjugate of ψ(t); The specific calculation formula for signal reconstruction is: Among them, ψ a，b represents the wavelet basis, the constant C ψ The following conditions are met: C ψ ＞0 (3) Where Ψ(ω) is the Fourier transform of ψ(t), and ω is the signal frequency; S4: using the first layer high frequency coefficients obtained in step S3, judging whether the carrier data has cycle slips by using the local singularity detection characteristics of wavelet transform, specifically comprising the following steps: At scale s0, if there exists a point (s0, τ0) that satisfies the local limit condition, then the point (s0, τ0) is a local extreme point, and the local limit condition is Where W f (s0, τ0) is the wavelet transform result of the signal at scale s0 and time τ0; If there is a zero crossing point at τ0, and for any point τ in a certain neighborhood of τ0, the maximum condition is satisfied, then the point (s0, τ0) is the modulus maximum point of the wavelet transform, and the maximum condition is |W f (s0, τ)|≤|W f (s0, τ0)|; The Lipschitz index is used to extract the singularity characteristics of the signal. For the signal function f(x), if n is a non-integer (n≤a<n+1), if there are constants A and h0 (h0>0) and an n-order polynomial P n (h), so that for any h≤h0, the following calculation conditions are satisfied: |f(x0+h)-P n (h)|≤A|h| a (5) Then f(x) is Lipschitz α at point x0. Specifically, if the above condition is satisfied for any x0 on the interval (a, b), x0+h∈(a, b), then f(x) is uniformly Lipschitz α on the interval (a, b); By performing multi-scale analysis on the signal, when the signal has a mutation, its wavelet transformed coefficient has the modulus maximum feature, and by detecting the modulus maximum point, it is determined whether the carrier data has a cycle slip. S5: When step S4 determines that there is a cycle slip in the carrier data, the following method is used to repair the cycle slip: S5-1: Select the Haar wavelet with the most significant response to mutations, decompose and reconstruct the satellite data according to different sampling rates, and select the correct carrier phase observation values ​​in several epochs before the cycle slip occurs; S5-2: Use the generalized extension extrapolation method to perform extrapolation repair, including: Assume that the known data sequence is {(t1, x1), (t2, x2), (t3, x3)...(t n , x n )}, need to extrapolate t n+1 The extrapolated value at time When the data points are dense, the average value of the latest point m is selected as the constraint condition to establish a generalized extension model: Where: (a0, a1, a2, a3) are the constant coefficients to be solved, minI(a0, a1, a2, a3) is the minimization optimal objective function, and the difference between the extrapolated result and the actual carrier phase observation value is used to repair the cycle slip, thereby realizing the accurate detection and repair of the cycle slip in the single-frequency BDS carrier data.

2. A BDS cycle slip repair method based on wavelet and generalized extrapolation according to claim 1, characterized in that: The reading of carrier data in step S1 includes parsing the Rinex observation file to obtain the carrier phase observation value.

3. The BDS cycle slip repair method based on wavelet and generalized extrapolation according to claim 1, characterized in that: The single-difference detection amount constructed in step S2 adopts the single-difference detection amount between epochs.

4. The BDS cycle slip repair method based on wavelet and generalized extrapolation according to claim 1, characterized in that: Step S3 specifically includes: using five different wavelet functions, namely Haar, Db, Sym, Coif and Bior, to perform wavelet transform decomposition and reconstruction on the single difference detection quantity of the carrier data, obtain the first layer of high frequency coefficients, and evaluate the response characteristics of the above five wavelets to the mutation signal.

5. The BDS cycle slip repair method based on wavelet and generalized extrapolation according to claim 1, characterized in that: The wavelet transform decomposition described in step S3 is accelerated by using discrete wavelet transform, which specifically includes the following steps: Select the scaling parameter and translation parameter in binary form, and the calculation formula is: a＝2 m (7) b=n2 m (8) where m and n are integers; for a particular choice of ψ(t), there exists a corresponding discrete wavelet ψ that satisfies the good time-frequency localization condition m，n , the calculation formula for good time-frequency localization conditions is: ψ m，n (t)=2 -m / 2 ψ(2 -m (9) Discrete wavelet is the structure of L 2 (R) is an orthogonal basis. Using this orthogonal basis, L 2 The calculation formula for any f(t) in (R) space is: Among them, the discrete wavelet coefficient a m，n The calculation formula is: a m，n =∫ R f(t)ψ* m，n (t)dt (11) Among them, ψ* m，n represents the complex conjugate of ψ(t) under m and n. The signal can be decomposed into an approximate part and a detail part. The detail calculation formula of the M layer is: Among them, U is a set of positive integers, and the approximation of the M layer is defined as the sum of all details before this layer, and its calculation formula is: Among them, the calculation formula of the signal f(t) is: Formula (14) further yields: A M-1 (t)=A M (t)+D M (t) (15) Through the above signal reconstruction, detail components are obtained at different scales to detect and analyze local feature changes in the signal and achieve accurate detection of cycle slips.

6. A BDS cycle slip repair method based on wavelet and generalized extrapolation according to claim 1, characterized in that: Step S4 determines whether a cycle slip occurs in the carrier data by detecting the modulus maximum point, and specifically includes the following steps: The multi-scale analysis of the signal is performed by using wavelet transform. By detecting the local singularity of the signal, when the signal changes suddenly, the coefficients after wavelet transform have the characteristics of modulus maximum value. Perform modulus maximum point detection on the signal coefficients after wavelet transformation, and judge whether there is abnormality or mutation in the signal by identifying the modulus maximum point; Combined with the concept of Lipschitz index, when the Lipschitz index α of the signal at point x0 is less than 1, this point is considered to be a singular point of the signal, further confirming whether the signal has a cycle slip; According to the local characteristics of the coefficients after wavelet transformation, the mutation points in the signal are detected, so as to effectively identify whether the cycle slip phenomenon occurs in the carrier data.

Images