Satellite-borne GPS cycle slip detection method for low-orbit satellite

Abstract

The invention is suitable for the field of satellite navigation signal processing, and provides a low-orbit satellite satellite-borne GPS cycle slip detection method, which comprises the following steps: acquiring observation value data of a low-orbit satellite satellite-borne GPS, and calculating an MW combination wide lane ambiguity value of the satellite-borne GPS; an IAO algorithm is adopted to optimize the set average number of times and the noise standard deviation of the ICEEMDAN algorithm; decomposing the wide lane ambiguity value into a plurality of intrinsic mode function (IMF) components and a residual component by adopting an optimized ICEEMDAN algorithm; a permutation entropy algorithm is adopted to calculate permutation entropy values of all the IMF components, the IMF components needing to be denoised are screened out, and the IMF components needing to be denoised are denoised in combination with a wavelet threshold denoising method; and performing signal reconstruction on the IMF component which is denoised by the wavelet threshold denoising method and the IMF component which does not need to be denoised to obtain denoised wide lane ambiguity, calculating the difference value of the denoised wide lane ambiguity between the current epoch and the adjacent epoch, and judging whether cycle slip exists or not according to a preset threshold. According to the method, the cycle slip detection accuracy and reliability of the low-orbit satellite in a high-dynamic and strong-noise environment are remarkably improved, and the misjudgment rate and the missed judgment rate are reduced.

Description

Technical Field

[0001] This application belongs to the field of satellite navigation signal processing, and in particular relates to a method for detecting cycle slips on low-Earth orbit satellite-borne GPS. Background Technology

[0002] Low Earth Orbit (LEO) satellite-borne GPS technology has become an important means of achieving high-precision orbit determination and real-time navigation and positioning. Compared with GPS observations acquired from ground stations, LEO satellites, due to their high-speed motion and complex space environment, often experience cycle slips in their carrier phase observations, which significantly affect phase ambiguity resolution and orbit accuracy. A cycle slip refers to a sudden change in integer phase during carrier phase observation caused by factors such as signal interruption, multipath effects, receiver noise, or high dynamic conditions. To maintain the continuity of phase observations, cycle slips must be accurately detected and corrected. Therefore, research on satellite-borne GPS cycle slip detection methods applicable to the dynamic characteristics of LEO satellites is of great significance.

[0003] With the increasing application of Global Navigation Satellite System (GNSS) in LEO precise orbit determination, researchers both domestically and internationally have proposed various cycle slip detection strategies for dual-frequency and multi-frequency GPS signals. Currently common methods include the higher-order difference method, pseudorange-phase combination method, geometry-independent combination method, ionospheric residual method, wavelet transform method, and the Melbourne-Wunnema (MW) combination method. The MW combination method utilizes pseudorange and carrier phase observations of dual-frequency signals for cycle slip detection. This method detects cycle slips by constructing a wide-lane combination of carrier phase values ​​and a narrow-lane combination of pseudorange values. The MW combination method effectively eliminates the influence of geometric distance and ionospheric delay, and is simple in structure and widely applicable. However, because it relies on pseudorange observations, its results are highly sensitive to pseudorange observation noise. If the observation noise is large or the signal quality is poor, small cycle slips of 1-2 cycles may be masked by noise, leading to misjudgments and missed detections, making it difficult to meet the requirements of real-time precise orbit determination for LEO satellites. Summary of the Invention

[0004] The purpose of this application is to provide a method for detecting cycle slips on low-Earth orbit (LEO) satellite-borne GPS. This method aims to address the problem that the MW combination method relies on pseudorange observations, and its results are highly sensitive to pseudorange observation noise. If the observation noise is large or the signal quality is poor, small cycle slips of 1-2 cycles may be masked by noise, leading to misjudgment and missed detection, which makes it difficult to meet the requirements of real-time and precise orbit determination for LEO satellites.

[0005] This application provides a method for detecting cycle slips on a low-Earth orbit satellite-borne GPS system, including the following steps: S101. Obtain observation data of the low-orbit satellite-borne GPS. The observation data includes carrier phase observations and pseudorange observations corresponding to two different carrier frequencies. Calculate the MW combined wide-lane ambiguity value of the satellite-borne GPS based on the two different carrier frequencies and the corresponding carrier phase observations and pseudorange observations. S102. The IAO algorithm is used to optimize the ensemble mean number and noise standard deviation of the ICEEMDAN algorithm in order to minimize the fitness function. S103. The optimized ICEEMDAN algorithm is used to decompose the wide aisle ambiguity value into multiple intrinsic mode function (IMF) components and a residual component. S104. The permutation entropy algorithm is used to calculate the permutation entropy value of all IMF components, the IMF components that need to be denoised are selected, and the wavelet threshold denoising method is used to denoise the IMF components that need to be denoised. S105. The IMF component after denoising by wavelet threshold denoising method and the IMF component that does not need to be denoised are reconstructed to obtain the denoised wide-lane ambiguity. The difference between the denoised wide-lane ambiguity between the current epoch and the adjacent epoch is calculated, and the existence of cycle slip is determined according to the preset threshold.

[0006] In this application, the limitations of traditional empirical parameter selection are overcome by employing the IAO algorithm to adaptively optimize the key parameters of ICEEMDAN (ensemble mean degree and noise standard deviation). The optimized ICEEMDAN adaptively decomposes the noisy wide-lane ambiguity value into multiple IMF components and residual components. High-noise components are identified using the permutation entropy criterion, and then multi-level fine denoising is performed using a wavelet thresholding method, effectively suppressing the noise amplification effect of the MW combination. Finally, high-precision cycle slip detection is achieved through signal reconstruction and adjacent epoch difference determination. This application significantly improves the accuracy and reliability of cycle slip detection in high-dynamic, high-noise environments for low-Earth orbit satellites, and reduces the false positive and false negative rates. Attached Figure Description

[0007] Figure 1 This is a flowchart of a low-orbit satellite-borne GPS cycle slip detection method provided in one embodiment of this application.

[0008] Figure 2 This is a diagram showing the results of cycle slip detection using a traditional MW combination.

[0009] Figure 3 This is a fitness curve diagram of the IAO algorithm optimizing the ICEEMDAN parameters in one embodiment of this application.

[0010] Figure 4 This is a diagram showing the result of decomposing the MW combined wide aisle ambiguity using the ICEEMDAN algorithm in one embodiment of this application.

[0011] Figure 5 This is a schematic diagram of the entropy values ​​of each IMF component in one embodiment of this application.

[0012] Figure 6 This is a comparison diagram of the effects of wavelet threshold denoising on each IMF component in one embodiment of this application.

[0013] Figure 7 It is the epoch difference of the MW detection value after denoising in one embodiment of this application. Detailed Implementation

[0014] To make the objectives, technical solutions, and beneficial effects of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0015] To illustrate the technical solution described in this application, specific embodiments are provided below.

[0016] Please see Figure 1 The low-Earth orbit satellite-borne GPS cycle slip detection method provided in one embodiment of this application includes the following steps: S101. Obtain observation data of the low-orbit satellite-borne GPS. The observation data includes carrier phase observations and pseudorange observations corresponding to two different carrier frequencies. Calculate the MW combined wide-lane ambiguity value of the satellite-borne GPS based on the two different carrier frequencies and the corresponding carrier phase observations and pseudorange observations.

[0017] In one embodiment of this application, S101 may further include the following steps: preprocessing the observed data to remove outliers.

[0018] In one embodiment of this application, the calculation formula for calculating the MW combined wide-lane ambiguity value of the satellite-borne GPS based on the two different carrier frequencies and the corresponding carrier phase observations and pseudorange observations is as follows: in, The MW combined wide-lane ambiguity value, These are two different carrier frequencies. Carrier frequency The corresponding carrier phase observation value, Carrier frequency The corresponding pseudorange observations, For MW combined wavelength, y is the speed at which electromagnetic waves propagate in a vacuum.

[0019] S102. The IAO algorithm is used to optimize the ensemble mean number and noise standard deviation of the ICEEMDAN algorithm in order to minimize the fitness function.

[0020] In one embodiment of this application, S102 may specifically include the following steps: The IAO algorithm is used to adaptively optimize the key parameters of the ICEEMDAN algorithm, and the iteration steps, population size and parameter search range of the IAO algorithm are determined. The ensemble mean number NR and noise standard deviation Nstd of the ICEEMDAN algorithm are used as optimization variables, and the fitness function is the minimization of permutation entropy. The original signal is decomposed into multiple Intrinsic Mode Function (IMF) components and residual terms using the preset parameters of the ICEEMDAN algorithm. The fitness of each IMF component is evaluated one by one to obtain the fitness value of each IMF component. The IAO algorithm simulates four hunting strategies of eagles: high-altitude soaring expansion search, low-altitude soaring contraction search, low-flying grasping, and walking capture. It iteratively updates the set mean number of iterations (NR) and noise standard deviation (Nstd) through these strategies to minimize the fitness function and retain the best individuals based on an elite selection strategy. Finally, it determines whether the algorithm has converged. If it has not converged, it continues to iterate. When the fitness value of the IMF component changes less than a preset threshold or the maximum number of iterations is reached, the algorithm is considered to have converged, and the optimized parameters of the ICEEMDAN algorithm are output as the global optimal solution.

[0021] In one embodiment of this application, the step of decomposing the original signal into multiple Intrinsic Mode Function (IMF) components and residual terms using preset ICEEMDAN algorithm parameters, and evaluating the fitness of each IMF component to obtain the fitness value of each IMF component may specifically include the following steps: Input the MW combined wide-lane ambiguity value as the signal value, initialize the ICEEMDAN algorithm parameters and IAO algorithm parameters, and set the fitness function to minimize the permutation entropy; The signal value is decomposed using the ensemble mean number NR and noise standard deviation Nstd of the current ICEEMDAN algorithm to obtain multiple IMF components and residual terms; Calculate the fitness value for each IMF component based on the fitness function.

[0022] The Aquila Optimizer (AO) is a swarm intelligence optimization algorithm based on the hunting behavior of eagles, characterized by strong global search capabilities and fast convergence speed. The IAO algorithm is an improved version of the Aquila Optimizer. Based on the AO algorithm, the IAO algorithm mainly makes two optimizations: (1) using the Tent chaotic mapping method to complete population initialization; and (2) adaptively updating the global optimum with weights. These improvements aim to enhance the algorithm's diversity and convergence stability, avoiding getting trapped in local optima.

[0023] The ICEEMDAN (Intrinsic Complete Ensemble Empirical Mode Decomposition with Adaptive Noise) algorithm is an improved signal decomposition algorithm that combines the advantages of adaptive noise and complete ensemble empirical mode decomposition for denoising and feature extraction of nonlinear and non-stationary signals.

[0024] Ensemble averaging (NR) refers to the number of times the same signal (or data) is averaged after repeated measurements / sampling. Noise standard deviation (Nstd) is a statistic describing the dispersion of noise data (i.e., the standard deviation of noise), representing the amplitude of noise fluctuations around the mean.

[0025] In one embodiment of this application, the population initialization using the Tent chaotic mapping method specifically includes the following steps: The initial population is generated using the Tent chaotic mapping. For the j-th dimension variable of the i-th individual, the iterative formula for the Tent mapping is: in, is a chaotic variable, k=2 is the chaos control parameter, at which point the system is in a completely chaotic state; D is the problem dimension. Mapping chaotic variables to the actual search space is as follows: in, Let j be the position of the i-th individual. and are the lower and upper bounds of the j-th dimension, respectively, and N is the population size.

[0026] In one embodiment of this application, to enhance the exploration and development capabilities of dynamic balancing algorithms, an adaptive weight function is proposed. The adaptive weight update of the global optimal solution specifically includes the following steps: The adaptive weighting function is: Where t is the current iteration number and T is the maximum iteration number. For adaptive weighting coefficients, It increases non-linearly with the number of iterations; Calculate the fitness of each individual ; Adaptive adjustment of the globally optimal position yields the globally optimal solution. and optimal fitness value : if ,but , in, Let be the globally optimal position obtained before the t-th iteration. This is the position after adaptive weight adjustment. The objective function; For the exploration phase , The first strategy is to expand the search by soaring at high altitudes. This involves a vertical descent from a high altitude, where the eagle identifies the prey area and selects the best hunting ground by soaring vertically at high altitudes. in, Let be the globally optimal position obtained before the t-th iteration. This represents the average position of the current solution being connected at the t-th iteration. A random number in the interval [0, 1] This is a time decay factor that decreases with each iteration to avoid premature convergence. For the first i The position of each individual in the t-th iteration; The second strategy is to conduct low-altitude soaring and narrowing searches, namely contour flight and Levy flight. The Skyhawk explores the divergent search space by using contour flight with short gliding attacks, and combines the Levy flight mechanism to enhance the global search capability. in, Let Levy's flight step size vector be in D-dimensional space; For randomly selected individual locations, , This is the parameter vector for the spiral search trajectory.

[0027] The Levy flight stride length is calculated as follows: in, =1.5 is the stability index. Let g be the Gamma function, and u and v be random vectors that follow a normal distribution. express Follows a mean of 0 and a variance of The normal distribution; express It follows a standard normal distribution with a mean of 0 and a variance of 1.

[0028] The parameters for the spiral search trajectory are as follows: in, Let the initial radius be , The spiral growth coefficient is... Angular velocity, For the initial phase, and Indicates a spiral shape during the search. For the first The polar radius of each point For the first The polar angle of each point.

[0029] For the development phase , The third strategy: low-altitude capture. After successfully locking onto the target's specific location and completing the necessary landing and attack deployment, the Skyhawk will execute a vertical descent strategy. Its attack characteristic is that it first uses a probing strike method, using this initial contact to assess and analyze the prey's stress response; this is known as a low-altitude slow descent attack. in, The mean position of the population at the t-th iteration; To develop parameters, control the extent to which it approaches the global optimum; To develop parameters, the range of random walks is limited; UB and LB are the upper and lower bounds of the search space, respectively. It is a random number in the interval [0, 1].

[0030] Fourth strategy: Foot capture. The eagle pounces on its prey on foot and seizes it, using a mass function for the final attack. The formula is as follows: Where QF is the quality function used to balance the search strategy in the t-th iteration. For flight speed parameters, These are parameters representing the direction of motion.

[0031] After each new solution is generated, a greedy selection is performed based on the fitness value: in, To pass the strategy or The generated new solution, The objective function value is used; this mechanism ensures that the population fitness improves monotonically.

[0032] S103. The optimized ICEEMDAN algorithm is used to decompose the wide lane ambiguity value into multiple intrinsic mode function (IMF) components and a residual component.

[0033] In one embodiment of this application, S103 may specifically include the following steps: The optimized ICEEMDAN algorithm is used to decompose the wide-lane ambiguity values. By adaptively adding Gaussian white noise to the wide-lane ambiguity values ​​during each filtering process, non-stationary wide-lane ambiguity values ​​are decomposed. It is decomposed into multiple IMF components and a residual component; the multiple IMF components are arranged in descending order of frequency to extract the local features of the signal and the oscillation components at different frequencies, while the residual component represents the overall trend of the signal or the low-frequency components.

[0034] The non-stationary wide-lane ambiguity value This refers to the ambiguity of the wide alley. It is not a fixed integer, but a time-varying parameter that changes with time, space, or environmental factors.

[0035] In one embodiment of this application, the adaptive addition of Gaussian white noise to the wide-lane ambiguity value reduces the non-stationary wide-lane ambiguity value. The decomposition into multiple IMF components and one residual component specifically includes the following steps: The formula for adaptively adding Gaussian white noise to the wide-lane ambiguity value is as follows: in, The original signal, i.e., the wide-lane ambiguity value. To implement the j-th Gaussian white noise, The first IMF for white noise, Where is the noise figure and N is the number of integration iterations; Calculate the local mean and obtain the first residual component. The formula is as follows: in, This indicates a local mean operation; Extract the first IMF component The formula is as follows: Add adaptive Gaussian white noise in the q-th stage, as shown in the following formula: in, Let be the noise figure for the (q-1)th stage. The q-th IMF component of the white noise; Calculate the q-th residual component The formula is as follows: Recursive usage formula Extracting the q-th IMF component ; Continue extracting until the residual component is reached. Once the termination condition is met, the original signal is decomposed into q IMF components and one residual component. The final decomposition result is shown in the following formula: Where y is the original wide-lane ambiguity value, without artificially added noise. For the q-th IMF component, K represents the final residual components, and K is the total number of IMF components.

[0036] S104. The permutation entropy algorithm is used to calculate the permutation entropy value of all IMF components, the IMF components that need to be denoised are selected, and the wavelet threshold denoising method is used to denoise the IMF components that need to be denoised.

[0037] In one embodiment of this application, S104 may specifically include the following steps: The permutation entropy algorithm is used to calculate the permutation entropy value of all IMF components. If the permutation entropy value is greater than the preset denoising threshold, it is determined that the IMF component needs to be denoised. For IMF components that require denoising, wavelet thresholding is used to eliminate noise.

[0038] Permutation entropy is a time series complexity measurement method based on symbolic dynamics. The core idea of ​​this method is to transform a continuous time series into a discrete symbolic sequence, and then quantify the complexity of the system by analyzing the distribution of permutation patterns in the symbolic sequence.

[0039] In one embodiment of this application, the step of calculating the permutation entropy value of all IMF components using the permutation entropy algorithm, and determining that the IMF component needs to be denoised if the permutation entropy value is greater than a preset denoising threshold, may specifically include the following steps: For a signal sequence X of length N, i.e., q IMF components obtained by decomposing the wide-lane ambiguity value using the optimized ICEEMDAN algorithm, phase space reconstruction is performed using Takens' delay embedding theorem, based on the preset embedding dimension m and delay time. With embedding dimension = 3 and delay time = 1, the resulting matrix Y is: in, Let K be the phase space reconstruction matrix, where each row of the matrix represents an m-dimensional reconstruction vector, and there are a total of K m-dimensional reconstruction vectors. The number of reconstruction vectors is 1. ; The elements within each reconstructed vector are arranged in ascending order as follows: in, A column representing the elements within the reconstructed vector; After arranging the elements within each reconstructed vector in ascending order, a symbol sequence is constructed using the indices from their positions before sorting. This symbol sequence can represent each reconstructed component as follows: The probability of each permutation pattern is calculated using the frequency estimation method: Based on the probability distribution of the permutation pattern, and according to Shannon's definition of information entropy, the permutation entropy value is calculated as follows: To facilitate comparison and standardization under different embedding dimensions, a normalized permutation entropy value is defined: Based on the preset noise reduction threshold =0.7, then determine The IMF components need to be denoised.

[0040] In one embodiment of this application, the wavelet threshold denoising method performs threshold processing on wavelet coefficients to achieve signal denoising. This application uses db4 as the wavelet basis, which performs excellently in noise removal while preserving abrupt changes and step transitions in the signal, ensuring that the main information and features of the signal are retained. The wavelet decomposition has three layers, which can ensure both denoising effect and computational efficiency and signal integrity.

[0041] This application uses VisuShrink as the threshold. VisuShrink thresholding achieves denoising by removing all wavelet coefficients smaller than the threshold value, thus applying a uniform global threshold. This method is suitable for handling various signal and noise levels and can effectively reduce the risk of overfitting. The VisuShrink thresholding formula is as follows: ,in Let N be the standard deviation of the noise, and N be the signal length.

[0042] This application employs soft thresholding as a denoising strategy. By shrinking the wavelet coefficients according to a preset threshold and reconstructing the signal using the updated coefficients, the denoising result is obtained. Compared with other methods, the signal processed by soft thresholding exhibits better smoothness and continuity. The denoising formula is as follows: in, Let T represent the k-th wavelet coefficient of the j-th layer, where T is the VisuShrink threshold. When the wavelet coefficients are considered to contain a useful signal, they are retained but shrunken T units towards zero; when At that time, it was assumed that the wavelet coefficients were mainly noise, and the wavelet coefficients were set to zero.

[0043] S105. The IMF component after denoising by wavelet threshold denoising method and the IMF component that does not need to be denoised are reconstructed to obtain the denoised wide-lane ambiguity. The difference between the denoised wide-lane ambiguity between the current epoch and the adjacent epoch is calculated, and the existence of cycle slip is determined according to the preset threshold.

[0044] In one embodiment of this application, determining whether a cycle slip exists based on a preset threshold specifically involves: in, The standard deviation of the noise is given. Based on research, the difference between adjacent epochs without cycle slips is around 0.3. Therefore, the detection threshold is set as follows: Set to 0.3, the current epoch i and the previous epoch Denoising of the wide aisle blur When the difference exceeds a preset threshold, if the next epoch... With the current epoch Denoising the width of the wide aisle If the difference is less than 0.1, it is judged that a cycle slip has occurred; otherwise, it is judged as a gross error.

[0045] The following describes a low-orbit satellite-borne GPS cycle slip detection method provided in an embodiment of this application through specific examples.

[0046] The following analysis examines the cycle slip detection performance at a sampling rate of 1:10s.

[0047] To verify the effectiveness of the low-Earth orbit satellite-borne GPS cycle slip detection method provided in one embodiment of this application, a simulation experiment was conducted using satellite-borne GPS data from the GRACE-FO satellite on the 152nd day of 2025. The original observation data consisted of 223 epochs of cycle slip-free satellite-borne GPS data with a sampling rate of 10 seconds. Cycle slips were artificially added to the original observation data, and then the low-Earth orbit satellite-borne GPS cycle slip detection method provided in one embodiment of this application was used to detect cycle slips and verify the performance of the cycle slip detection method.

[0048] To verify the accuracy of the low-orbit satellite-borne GPS cycle slip detection provided in one embodiment of this application, cycle slips of different sizes were artificially added at some epochs of the original carrier phase when the sampling rate was 10 seconds, as shown in Table 1.

[0049] Table 1: Cases of manually added cycle jumps of different sizes like Figure 2 As shown, the traditional MW combination detection capability is significantly insufficient for satellite-borne GPS data. In the introduced cycle slips across eight epochs, the traditional MW combination method can only detect one cycle slip, resulting in obvious missed detections.

[0050] like Figure 3 As shown, the IAO algorithm is used to optimize the ensemble mean number of iterations and the noise standard deviation of the ICEEMDAN algorithm. The algorithm reaches stable convergence after 15 iterations, and the fitness value decreases from 4.842 to approximately 4.812. The convergence process is smooth and monotonically decreasing. The optimal parameters obtained are NR=75 and Nstd=0.5. The IAO algorithm achieves significant performance improvement and stable global optimization within a finite number of iterations.

[0051] like Figure 4 As shown, the MW combined wide-lane ambiguity value is decomposed by ICEEMDAN into 5 IMF components (IMF1~IMF5) and 1 residual component Res. It can be clearly seen from the figure that at epochs 30, 60, 90, 120, 150, and 180, the IMF1~IMF3 components all show obvious jump characteristics, corresponding to artificially added cycle slip positions.

[0052] like Figure 5 As shown, the permutation entropy value is calculated for each IMF component. The PE threshold is set to 0.7, and IMF1 and IMF2 with PE > 0.7 are identified as noise-dominant components.

[0053] like Figure 6As shown, wavelet threshold denoising is used to further reduce noise in the IMF1 and IMF2 components. Wavelet threshold denoising significantly reduces the noise level of the MW detection value while preserving the signal jump that occurs during cycle slip.

[0054] like Figure 7 As shown, the upper and lower horizontal lines represent the set upper and lower threshold values, respectively, and the waveform in the middle represents the difference between denoised MW detection values ​​at each epoch. The cycle slip detection method, which combines improved MW combination ICEEMDAN decomposition with wavelet threshold denoising, compensates for the shortcomings of the MW combination. The set ±0.3 threshold successfully identified all six artificial cycle slips, proving the reliability of the method. Specifically, cycle slips at epochs 30 and 60 exhibit positive jumps of approximately +0.7, while those at epochs 90 and 120 exhibit negative jumps of approximately -0.6. The largest cycle slip amplitudes are at epochs 150 and 180, reaching +1.6 and -1.6 respectively, showing a clear alternation of positive and negative values. MW detection values ​​during normal observation periods remain stable within the threshold range with very small fluctuations, indicating that the denoising process effectively suppresses the interference of observation noise. Experimental results show that the improved MW combination method has good detection capability for cycle slips of different sizes, providing an effective technical means for cycle slip detection and repair in actual GNSS data processing.

[0055] In this application, the limitations of traditional empirical parameter selection are overcome by employing the IAO algorithm to adaptively optimize the key parameters of ICEEMDAN (ensemble mean degree and noise standard deviation). The optimized ICEEMDAN adaptively decomposes the noisy wide-lane ambiguity value into multiple IMF components and residual components. High-noise components are identified using the permutation entropy criterion, and then multi-level fine denoising is performed using a wavelet thresholding method, effectively suppressing the noise amplification effect of the MW combination. Finally, high-precision cycle slip detection is achieved through signal reconstruction and adjacent epoch difference determination. This application significantly improves the accuracy and reliability of cycle slip detection in high-dynamic, high-noise environments for low-Earth orbit satellites, and reduces the false positive and false negative rates.

[0056] It should be understood that the steps in the various embodiments of this application are not necessarily executed sequentially according to the order indicated by the step numbers. Unless explicitly stated herein, there is no strict order restriction on the execution of these steps, and they can be executed in other orders. Moreover, at least some steps in each embodiment may include multiple sub-steps or multiple stages. These sub-steps or stages are not necessarily completed at the same time, but can be executed at different times. The execution order of these sub-steps or stages is not necessarily sequential, but can be performed alternately or in turn with other steps or at least a portion of the sub-steps or stages of other steps.

[0057] Those skilled in the art will understand that all or part of the processes in the above embodiments can be implemented by a computer program instructing related hardware. The program can be stored in a non-volatile computer-readable storage medium, and when executed, it can include the processes of the embodiments described above. Any references to memory, storage, databases, or other media used in the embodiments provided in this application can include non-volatile and / or volatile memory. Non-volatile memory can include read-only memory (ROM), programmable ROM (PROM), electrically programmable ROM (EPROM), electrically erasable programmable ROM (EEPROM), or flash memory. Volatile memory can include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM is available in various forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), dual data rate SDRAM (DDRSDRAM), enhanced SDRAM (ESDRAM), synchronous link DRAM (SLDRAM), RAMbus direct RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and RAMbus dynamic RAM (RDRAM), etc.

[0058] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0059] The above embodiments merely illustrate several implementation methods of the present invention, and their descriptions are relatively specific and detailed, but they should not be construed as limiting the scope of the invention patent. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these all fall within the protection scope of the present invention. Therefore, the protection scope of this invention patent should be determined by the appended claims.

Claims

1. A method for detecting cycle slips on a low-Earth orbit satellite-borne GPS system, characterized in that, Includes the following steps: S101. Obtain observation data of the low-orbit satellite-borne GPS. The observation data includes carrier phase observations and pseudorange observations corresponding to two different carrier frequencies. Calculate the MW combined wide-lane ambiguity value of the satellite-borne GPS based on the two different carrier frequencies and the corresponding carrier phase observations and pseudorange observations. S102. The IAO algorithm is used to optimize the ensemble mean number and noise standard deviation of the ICEEMDAN algorithm in order to minimize the fitness function. S103. The optimized ICEEMDAN algorithm is used to decompose the wide aisle ambiguity value into multiple intrinsic mode function (IMF) components and a residual component. S104. The permutation entropy algorithm is used to calculate the permutation entropy value of all IMF components, the IMF components that need to be denoised are selected, and the wavelet threshold denoising method is used to denoise the IMF components that need to be denoised. S105. The IMF component after denoising by wavelet threshold denoising method and the IMF component that does not need to be denoised are reconstructed to obtain the denoised wide-lane ambiguity. The difference between the denoised wide-lane ambiguity between the current epoch and the adjacent epoch is calculated, and the existence of cycle slip is determined according to the preset threshold.

2. The low-orbit satellite-borne GPS cycle slip detection method as described in claim 1, characterized in that, S101 further includes the following steps: preprocessing the observed data to remove outliers; The formula for calculating the MW combined wide-lane ambiguity value of the satellite-borne GPS based on the two different carrier frequencies and the corresponding carrier phase observations and pseudorange observations is as follows: in, The MW combined wide-lane ambiguity value, These are two different carrier frequencies. Carrier frequency The corresponding carrier phase observation value, Carrier frequency The corresponding pseudorange observations, For MW combined wavelength, y is the speed at which electromagnetic waves propagate in a vacuum.

3. The low-orbit satellite-borne GPS cycle slip detection method as described in claim 1, characterized in that, S102 specifically includes the following steps: The IAO algorithm is used to adaptively optimize the key parameters of the ICEEMDAN algorithm, and the iteration steps, population size and parameter search range of the IAO algorithm are determined. The ensemble mean number NR and noise standard deviation Nstd of the ICEEMDAN algorithm are used as optimization variables, and the fitness function is the minimization of permutation entropy. The original signal is decomposed into multiple Intrinsic Mode Function (IMF) components and residual terms using the preset parameters of the ICEEMDAN algorithm. The fitness of each IMF component is evaluated one by one to obtain the fitness value of each IMF component. The IAO algorithm simulates four hunting strategies of eagles: high-altitude soaring expansion search, low-altitude soaring contraction search, low-flying grasping, and walking capture. It iteratively updates the set mean number of iterations (NR) and noise standard deviation (Nstd) through these strategies to minimize the fitness function and retain the best individuals based on an elite selection strategy. Finally, it determines whether the algorithm has converged. If it has not converged, it continues to iterate. When the fitness value of the IMF component changes less than a preset threshold or the maximum number of iterations is reached, the algorithm is considered to have converged, and the optimized parameters of the ICEEMDAN algorithm are output as the global optimal solution.

4. The low-orbit satellite-borne GPS cycle slip detection method as described in claim 3, characterized in that, The process of decomposing the original signal into multiple Intrinsic Mode Function (IMF) components and residual terms using preset ICEEMDAN algorithm parameters, and evaluating the fitness of each IMF component to obtain its fitness value, specifically includes the following steps: Input the MW combined wide-lane ambiguity value as the signal value, initialize the ICEEMDAN algorithm parameters and IAO algorithm parameters, and set the fitness function to minimize the permutation entropy; The signal value is decomposed using the ensemble mean number NR and noise standard deviation Nstd of the current ICEEMDAN algorithm to obtain multiple IMF components and residual terms; Calculate the fitness value for each IMF component based on the fitness function.

5. The low-orbit satellite-borne GPS cycle slip detection method as described in claim 1, characterized in that, The IAO algorithm is an improved version of the Skyhawk optimization algorithm. Based on the Skyhawk optimization algorithm, the IAO algorithm mainly makes two optimizations: it uses the Tent chaotic mapping method to complete the population initialization. Adaptive weight update of the global optimal solution; The population initialization using the Tent chaotic mapping method specifically includes the following steps: The initial population is generated using the Tent chaotic mapping. For the j-th dimension variable of the i-th individual, the iterative formula for the Tent mapping is: in, is a chaotic variable, k=2 is the chaos control parameter, at which point the system is in a completely chaotic state; D is the problem dimension. Mapping chaotic variables to the actual search space is as follows: in, Let j be the position of the i-th individual. and These are the lower and upper bounds of the j-th dimension, respectively, and N is the population size; The adaptive weight update of the global optimal solution specifically includes the following steps: The adaptive weighting function is: Where t is the current iteration number and T is the maximum iteration number. For adaptive weighting coefficients, It increases non-linearly with the number of iterations; Calculate the fitness of each individual ; Adaptive adjustment of the globally optimal position yields the globally optimal solution. and optimal fitness value : if ,but , in, Let be the globally optimal position obtained before the t-th iteration. This is the position after adaptive weight adjustment. The objective function; For the exploration phase , First strategy: Soaring high to expand the search, that is, vertical descent from high altitude, the eagle identifies the prey area, and selects the best hunting area by soaring vertically at high altitude. in, Let be the globally optimal position obtained before the t-th iteration. This represents the average position of the current solution being connected at the t-th iteration. A random number in the interval [0, 1] This is a time decay factor that decreases with each iteration to avoid premature convergence. For the first i The position of each individual in the t-th iteration; The second strategy is to soar and shrink the search, namely contour flight and Levy flight. The Sky Eagle explores the divergent search space by using contour flight with short gliding attacks, and combines the Levy flight mechanism to enhance the global search capability. in, Let Levy's flight step size vector be in D-dimensional space; For randomly selected individual locations, , This is the parameter vector for the spiral search trajectory; The Levy flight stride length is calculated as follows: in, =1.5 is the stability index. Let g be the Gamma function, and u and v be random vectors that follow a normal distribution. express Follows a mean of 0 and a variance of The normal distribution; express It follows a standard normal distribution with a mean of 0 and a variance of 1; The parameters for the spiral search trajectory are as follows: in, Let the initial radius be , The spiral growth coefficient is... Angular velocity, For the initial phase, and Indicates a spiral shape during the search. For the first The polar radius of each point For the first The polar angle of each point; For the development phase , The third strategy: low-flying capture, that is, after successfully locking onto the specific location of the target and completing the necessary landing and attack deployment work, the Skyhawk will execute a vertical descent strategy; in, The mean position of the population at the t-th iteration; To develop parameters, control the extent to which it approaches the global optimum; To develop parameters, the range of random walks is limited; UB and LB are the upper and lower bounds of the search space, respectively. A random number in the interval [0, 1]; Fourth strategy: Foot capture. The eagle pounces on its prey on foot and seizes it, using a mass function for the final attack; the formula is as follows: Where QF is the quality function used to balance the search strategy in the t-th iteration. For flight speed parameters, These are parameters related to the direction of motion. After each new solution is generated, a greedy selection is performed based on the fitness value: in, To pass the strategy or The generated new solution, The objective function value is used; this mechanism ensures that the population fitness improves monotonically.

6. The low-orbit satellite-borne GPS cycle slip detection method as described in claim 1, characterized in that, S103 specifically includes the following steps: The optimized ICEEMDAN algorithm is used to decompose the wide-lane ambiguity values. By adaptively adding Gaussian white noise to the wide-lane ambiguity values ​​during each filtering process, non-stationary wide-lane ambiguity values ​​are decomposed. It is decomposed into multiple IMF components and a residual component; the multiple IMF components are arranged in descending order of frequency to extract the local features of the signal and the oscillation components at different frequencies, while the residual component represents the overall trend of the signal or the low-frequency components.

7. The low-orbit satellite-borne GPS cycle slip detection method as described in claim 6, characterized in that, The method of adaptively adding Gaussian white noise to the wide-lane ambiguity value reduces the non-stationary wide-lane ambiguity value. The decomposition into multiple IMF components and one residual component specifically includes the following steps: The formula for adaptively adding Gaussian white noise to the wide-lane ambiguity value is as follows: in, The original signal, i.e., the wide-lane ambiguity value. To implement the j-th Gaussian white noise, The first IMF for white noise, Where is the noise figure and N is the number of integration iterations; Calculate the local mean and obtain the first residual component. The formula is as follows: in, This indicates a local mean operation; Extract the first IMF component The formula is as follows: Add adaptive Gaussian white noise in the q-th stage, as shown in the following formula: in, Let be the noise figure for the (q-1)th stage. The q-th IMF component of the white noise; Calculate the q-th residual component The formula is as follows: Recursive usage formula Extracting the q-th IMF component ; Continue extracting until the residual component is reached. Once the termination condition is met, the original signal is decomposed into q IMF components and one residual component. The final decomposition result is shown in the following formula: Where y is the original wide-lane ambiguity value, without artificially added noise. For the q-th IMF component, K represents the final residual components, and K is the total number of IMF components.

8. The low-orbit satellite-borne GPS cycle slip detection method as described in claim 1, characterized in that, S104 specifically includes the following steps: The permutation entropy algorithm is used to calculate the permutation entropy value of all IMF components. If the permutation entropy value is greater than the preset denoising threshold, it is determined that the IMF component needs to be denoised. For IMF components that require denoising, wavelet thresholding is used to eliminate noise. The step of calculating the permutation entropy value of all IMF components using the permutation entropy algorithm, and determining that the IMF component needs to be denoised if the permutation entropy value is greater than a preset denoising threshold, specifically includes the following steps: For a signal sequence X of length N, i.e., q IMF components obtained by decomposing the wide-lane ambiguity value using the optimized ICEEMDAN algorithm, phase space reconstruction is performed using Takens' delay embedding theorem, based on the preset embedding dimension m and delay time. With embedding dimension = 3 and delay time = 1, the resulting matrix Y is: in, Let K be the phase space reconstruction matrix, where each row of the matrix represents an m-dimensional reconstruction vector, and there are a total of K m-dimensional reconstruction vectors. The number of reconstruction vectors is 1. ; The elements within each reconstructed vector are arranged in ascending order as follows: in, A column representing the elements within the reconstructed vector; After arranging the elements within each reconstructed vector in ascending order, a symbol sequence is constructed using the indices from their positions before sorting. This symbol sequence can represent each reconstructed component as follows: The probability of each permutation pattern is calculated using the frequency estimation method: Based on the probability distribution of the permutation pattern, and according to Shannon's definition of information entropy, the permutation entropy value is calculated as follows: To facilitate comparison and standardization under different embedding dimensions, a normalized permutation entropy value is defined: Based on the preset noise reduction threshold =0.7, then determine The IMF components need to be denoised.

9. The low-orbit satellite-borne GPS cycle slip detection method as described in claim 8, characterized in that, The wavelet threshold denoising method is to perform threshold processing on wavelet coefficients to achieve signal denoising. It uses db4 as the wavelet basis and has 3 wavelet decomposition layers. VisuShrink is used as the threshold, and soft thresholding is used as the denoising strategy. The denoising result is obtained by shrinking the wavelet coefficients according to the preset threshold and reconstructing the signal using the updated coefficients. The noise reduction formula is as follows: in, Let T represent the k-th wavelet coefficient of the j-th layer, where T is the VisuShrink threshold. When the wavelet coefficients are considered to contain a useful signal, they are retained but shrunken T units towards zero; when At that time, it was assumed that the wavelet coefficients were mainly noise, and the wavelet coefficients were set to zero.

10. The low-orbit satellite-borne GPS cycle slip detection method as described in claim 1, characterized in that, The specific steps for determining whether a cycle slip exists based on a preset threshold are as follows: in, Let $i$ be the standard deviation of the noise, when the current epoch $i$ is the same as the previous epoch $i$. Denoising of the wide aisle blur When the difference exceeds a preset threshold, if the next epoch... With the current epoch Denoising the width of the wide aisle If the difference is less than 0.1, it is judged that a cycle slip has occurred; otherwise, it is judged as a gross error.

Images