DI-based adaptive GF cycle slip detection threshold model method and system

Abstract

The invention belongs to the technical field of GNSS precise point positioning cycle slip detection, and discloses a DI-based adaptive GF cycle slip detection threshold model method, and the method comprises the steps: obtaining IGS tracking station data during an extra-large magnetic storm; preprocessing the data, and calculating an inter-epoch difference GF value and a DI value; obtaining 99.7 percentile in each interval section according to a 3 sigma criterion by using distribution of a time domain difference GF value and a DI value; performing quadratic polynomial fitting on the percentile to obtain a final threshold model; the GF adaptive threshold model based on the DI avoids the limitation that the ROTI is influenced by cycle slip, and effectively solves the problem of false detection.

Description

Technical Field

[0001] This invention belongs to the field of GNSS precise single-point positioning cycle slip detection technology, and particularly relates to an adaptive GF cycle slip detection threshold model method and system based on DI. Background Technology

[0002] Precise Point Positioning (PPP) is an advanced Global Navigation Satellite System (GNSS) positioning technology that enables high-precision location-based services using a single receiver. Due to its numerous advantages, PPP has been widely applied in geodesy and geophysics, ionospheric monitoring, meteorology, and other space-related fields over the past 20 years. However, PPP still faces many challenges in practical applications, such as service degradation or interruption due to ionospheric activity. Under complex ionospheric conditions, the total electron content (TEC) of the ionosphere varies drastically, significantly affecting raw GNSS observations and thus reducing the positioning performance of PPP. In low-latitude regions where TEC changes significantly, the positioning error of a single epoch PPP can reach several meters.

[0003] Geomagnetic storms (GSM) are global, violent disturbances to the Earth's magnetic field caused by high-speed solar wind plasma clouds or enhanced solar wind flows generated by solar activity (such as coronal mass ejections and solar flares). The physical process involves the interaction of the interplanetary magnetic field carried by the solar wind with the Earth's magnetosphere magnetic field, resulting in magnetic reconnection at the magnetopause. This drives energy coupling in the magnetosphere-ionosphere system, leading to enhanced loop currents and distortion of the magnetosphere's magnetic field structure. Consequently, the intensity and direction of the Earth's magnetic field on the surface change significantly within minutes to hours. The space environment disturbances caused by GSM can interfere with satellite orbit operations, high-frequency communications, power transmission networks, and other technological systems. They also induce auroras in high-latitude regions, making them one of the core physical processes in space weather research. As one of the most representative complex ionospheric events, GSM has a significant adverse impact on PPP (Potentially Parallel Sphere Analysis) results.

[0004] During geomagnetic storms, due to severe ionospheric disturbances, traditional geometry-free (GF) cycle slip detection models suffer from numerous false detections due to their sensitivity to the ionosphere. Currently, there are GF dynamic threshold models based on the Rate of Change of Total Electron Content (ROTI) index. ROTI is calculated from observed TEC sequences. High-precision TECs are generally derived from GNSS dual-frequency carrier phase observations; however, GNSS carrier phase observations are affected by cycle slips, leading to discontinuous measurements. Therefore, cycle slip detection and correction are required before calculating TECs, which contradicts the introduction of ROTI to establish a cycle slip detection model. In fact, another commonly used type of GNSS observation is Doppler measurement, which has the advantage of being unaffected by cycle slips. Therefore, a GF adaptive threshold model constructed based on the Doppler Index (DI) calculated from Doppler observations can effectively avoid the problem of the conflict between ROTI and cycle slip detection logic.

[0005] Based on the above analysis, the problems and shortcomings of the existing technology are as follows:

[0006] (1) Traditional GF cycle slip detection models are prone to false detections due to their sensitivity to the ionosphere. Currently, there are GF dynamic threshold models based on ROTI.

[0007] (2) ROTI is calculated from the observed TEC sequence. High-precision TEC generally comes from GNSS dual-frequency carrier phase observations. However, GNSS carrier phase observations are affected by cycle slips, resulting in discontinuous measurement values. Summary of the Invention

[0008] To address the problems existing in the prior art, this invention provides an adaptive GF cycle slip detection threshold model method based on DI.

[0009] This invention is implemented as follows: A DI-based adaptive GF cycle slip detection threshold model method includes:

[0010] Step S1: Obtain IGS tracking station data during the severe geomagnetic storm;

[0011] Step S2: Preprocess the data and calculate the inter-epoch difference GF value and DI value;

[0012] Step S3: Based on the distribution of the time-domain difference GF value and DI value, obtain the boundary point between the piecewise model constant threshold and the function threshold interval, and obtain the 99.7 percentile of the time-domain difference GF corresponding to DI in each step of the function threshold interval according to the 3σ criterion;

[0013] Step S4: Perform quadratic polynomial fitting on the percentiles and take a reasonable constant threshold to obtain the final piecewise threshold model;

[0014] Wherein, DI is the Doppler index, GF is the geometric cycle slip detection model, and IGS is the International GNSS Service Organization.

[0015] Furthermore, step S1, the process of acquiring IGS tracking station data during severe geomagnetic storms, includes:

[0016] Raw observation data from IGS tracking stations with dual-frequency 30-second sampling rate throughout the day were selected during typical severe geomagnetic storms to support Doppler observations.

[0017] Furthermore, step S2, the process of preprocessing the data and calculating the inter-epoch difference GF value and DI value, includes:

[0018] The expression for dual-frequency GF combination (Geometry-Free) phase combination is:

[0019]

[0020] in, and They are frequencies and Phase measurement value in meters. and It is the phase wavelength at the corresponding frequency. and It is phase ambiguity; For ionospheric delay term, This is GF combined phase noise;

[0021] For consecutive epochs Perform the difference to obtain the time-domain difference. ,Right now The expression is:

[0022]

[0023] in, and These represent the GF phase combination values ​​for the next epoch and the previous epoch, respectively;

[0024] By calculating consecutive epochs The difference, without detecting cycle slips, eliminates the ambiguity term in the equation; therefore, the time-domain difference... The size, that is Dependence on time-domain differential residual ionospheric delay ,in It is a time-domain difference operator;

[0025] Under normal circumstances, The values ​​are very small, and the time intervals are short, such as 1 to 30 seconds. Therefore, when using GF combination to detect cycle slips, an empirical threshold of 0.05 m is usually set. This method performs well during the ionospheric calm period. However, studies have shown that the empirical threshold cannot effectively handle ionospheric anomalies, leading to false cycle slip detections and generating a large number of unnecessary ambiguity initializations, thereby reducing positioning accuracy. Therefore, it is necessary to construct a GF cycle slip detection model that adaptively sets the threshold according to ionospheric changes.

[0026] The threshold model is fitted using the GNSS Doppler index (DI), which is defined as follows:

[0027]

[0028]

[0029]

[0030] in, Indicates frequency, These are Doppler measurements; This represents the rate of change of the geometric distance between the satellite and the receiver. At the speed of light, and These are the clock bias rates of the satellite and receiver, respectively. The ionospheric delay rate, The rate of change of tropospheric delay. For noise, This indicates the average calculation; the unit of DI is TECU / min, consistent with ROTI; in this study, a 5-minute sliding window was used to calculate DI, just like ROTI.

[0031] Further, step S3, which involves obtaining the boundary points of the piecewise model constant threshold and function threshold intervals based on the distribution of the time-domain difference GF and DI values, and obtaining the 99.7 percentile of the time-domain difference GF corresponding to DI within each step of the function threshold interval according to the 3σ criterion, includes:

[0032] Obtained from raw GNSS data of IGS stations selected during typical severe geomagnetic storms. Based on the distribution of DI, the boundary points between the constant threshold and the function threshold intervals should be reasonably set. For the function threshold model interval, a step size of 0.05 TECU / min should be used, and within each step size... The 99.7th percentile (3σ).

[0033] Furthermore, in step S4, a quadratic polynomial fitting is performed on the percentiles, and a reasonable constant threshold is taken to obtain the final piecewise threshold model.

[0034] For the function threshold interval, for DI within each step... A quadratic polynomial was fitted using the 99.7 percentile; for the constant threshold interval, a reasonable constant threshold was selected. This yielded the final piecewise threshold model.

[0035] Another objective of this invention is to provide a DI-based adaptive GF cycle slip detection threshold model system comprising:

[0036] The acquisition module is used to acquire data from IGS tracking stations during severe geomagnetic storms.

[0037] The preprocessing module is used to preprocess the data and calculate the interepoch difference GF value and DI value;

[0038] The distribution module is used to obtain the 99.7 percentile within each interval segment based on the 3σ criterion by using the distribution of time-domain difference GF and DI values.

[0039] The fitting module is used to perform quadratic polynomial fitting on percentiles to obtain the final threshold model.

[0040] Another object of the present invention is to provide a computer device including a memory and a processor, the memory storing a computer program, which, when executed by the processor, causes the processor to perform the steps of the DI-based adaptive GF cycle slip detection threshold model method.

[0041] Another object of the present invention is to provide a computer-readable storage medium storing a computer program, which, when executed by a processor, causes the processor to perform the steps of the DI-based adaptive GF cycle slip detection threshold model method.

[0042] Another objective of this invention is to provide an information data processing terminal for implementing the DI-based adaptive GF cycle slip detection threshold model system.

[0043] Based on the above technical solutions and the technical problems solved, the advantages and positive effects of the technical solution to be protected by this invention are as follows:

[0044] In an active ionosphere, the source of false detections in the GF fixed threshold method is not simply "increased noise," but rather the significant thick tails and skewness of the ΔGF (inter-epoch GF) distribution under scintillation-driven conditions. The thresholds traditionally estimated using constants or global variance no longer correspond to stable upper bounds for Type I errors. Many adaptive schemes attempt to characterize perturbation strength using ROTI (standard deviation of STEC rate of change), but ROTI itself depends on the TEC sequence constructed from phase / code. Once a full-cycle or half-cycle slip occurs, the STEC time series is "contaminated" by discrete jumps, immediately amplifying ROTI and feeding error information back into the threshold update, forming a closed loop of "cycle slip → increased ROTI → threshold distortion → increased false detections." This invention selects the Doppler exponent DI (based on carrier frequency shift / phase rate and PLL dynamics), decoupled from the integer term of the cycle slip, as the perturbation strength characterization measure. This shifts the adaptive basis for cycle slip judgment to the "mechanical quantity of the tracking loop being perturbed," avoiding the inherent coupling limitations of ROTI.

[0045] At the methodological level, we empirically validated the joint statistical relationship between <ΔGF, DI> using samples from extreme geomagnetic storms: We divided the DI axis into steps, taking the 99.7 percentile of ΔGF within each step according to the 3σ criterion, to obtain the data-driven profile of "disturbance intensity → safety boundary." Then, we constructed a function threshold T(DI) using a quadratic polynomial fit, and introduced a physically reasonable upper limit constant as a guardrail at the sparse tail of high DI. This piecewise fitting structure ensures that the threshold can scale finely with DI in the low-to-medium disturbance region, maintaining stable false positive control; on the other hand, it avoids overfitting and suppresses missed detections in the sample-poor region of extreme disturbances. The key innovations are: replacing ROTI with DI as the adaptive independent variable, isolating the cycle slip integer term from the observation construction level; establishing the threshold using "conditional percentiles" instead of global variance, refining statistical control to each disturbance intensity level; and achieving a real-time computable, cross-constellation / frequency-band transferable engineering model through piecewise continuous splicing and low-order fitting.

[0046] The engineering benefits are reflected in the continuity and robustness of the solution link: when a geomagnetic storm arrives or scintillation suddenly increases, T(DI) expands in real time with DI, significantly reducing unnecessary cycle slip marking and ambiguity reconstruction, and avoiding frequent reconvergence and coordinate drift in PPP / PPP-AR; when the disturbance subsides, the threshold automatically tightens, suppressing abnormal phase leakage into the filter, thereby maintaining the consistency of residual statistics. In our validation based on IGS multi-station multi-frequency samples, we observed that the false detection rate during peak periods was significantly lower than that of the fixed threshold method, and the reconvergence time and coordinate repeatability improved accordingly; at the same time, because a constant guardrail was set in the high DI segment, the above improvements did not come at the cost of a surge in the false negative rate.

[0047] This invention explicitly incorporates the "external intensity of ionospheric space weather" into the cycle slip detection criterion via DI (Discriminator Indicator), achieving a discriminator that can reproduce the upper limit of experimental statistics even in extreme environments. This adaptive-robust balance is the source of its innovative technical effect. By introducing DI calculated from Doppler measurements and constructing an adaptive GF (Geothermal Field) threshold model for detecting cycle slips during geomagnetic storms, this invention overcomes the limitations of traditional fixed-threshold GF models during periods of ionospheric activity. Furthermore, due to the fact that Doppler measurements are unaffected by cycle slips, it effectively distinguishes itself from threshold models based on ROTI (Reactive Origin of Optical Time Indicators), achieving better cycle slip detection and significantly improving the positioning accuracy of PPP (Potential Point of Detection) during geomagnetic storms. Attached Figure Description

[0048] Figure 1 This is a flowchart of the adaptive GF cycle slip detection threshold model method based on DI provided in the embodiments of the present invention.

[0049] Figure 2 This is a system structure block diagram of the adaptive GF cycle slip detection threshold model method based on DI provided in the embodiments of the present invention.

[0050] Figure 3 This is a graph of Dst, IMF-Bz, and Kp indices from May 10 to 12, 2024, provided in an embodiment of the present invention.

[0051] Figure 4 This is a correlation analysis diagram between GNSS data from 110 IGS tracking stations and DI, provided in an embodiment of the present invention.

[0052] Figure 5 This is a comparison chart of dynamic PPP results based on Scheme 1 (blue) and Scheme 2 (red) using GPS data from FRDN, PRDS and HOFN stations on May 11, 2024, provided by an embodiment of the present invention.

[0053] Figure 6 This is an embodiment of the present invention that compares the dynamic precise point positioning results based on scheme 1 (blue) and scheme 2 (red) using GPS data from FRDN, PRDS and HOFN stations on October 11, 2024. Detailed Implementation

[0054] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0055] The engineering challenge addressed by this invention stems from the instability of high-precision GNSS applications under severe space weather events. Ionospheric irregularities and phase scintillation triggered by severe geomagnetic storms significantly amplify instantaneous fluctuations in carrier observations. This causes the traditional geometrically independent (GF) cycle slip discrimination threshold, set with a fixed constant, to become "too small" under strong disturbances, leading to frequent false alarms; conversely, the same threshold becomes "too large" during calm periods, causing missed detections. For real-time PPP / RTK, deformation monitoring, and timing systems, false alarms trigger unnecessary ambiguity reconstruction and convergence restarts, while missed detections introduce abnormal phases into the solution, causing coordinate / time-frequency solution drift. Therefore, the industry urgently needs a threshold model that can dynamically adapt to ionospheric and signal conditions, maintaining stable and reproducible cycle slip detection performance even under extreme geomagnetic storm environments.

[0056] The mechanism of this method is based on the joint characterization of GF combination and Doppler index (DI). GF combination eliminates geometric and clock errors through cross-frequency linearization, highlighting the ionospheric term and integer ambiguity in the carrier phase. Epoch-level differencing of GF can essentially suppress the integer term and manifest rapid ionospheric changes, receiver noise, and cycle slip abrupt changes as time-domain pulses. DI, as a dimensionless intensity index constructed based on carrier Doppler / phase rate, can characterize the "stress level" of the tracking loop under disturbance. When ionospheric activity intensifies, DI rises and the thick tail of the ΔGF distribution becomes obvious. The separability of cycle slips and noise depends on whether the threshold can expand and contract synchronously with DI. This invention replaces the static threshold with a "DI→threshold" mapping, enabling the detector to have a flexible response to space weather.

[0057] At the data level, the method first selects multi-constellation, multi-frequency data from IGS tracking stations during severe geomagnetic storms to ensure that the training samples cover extreme scenarios and have regional diversity. The preprocessing stage completes observation synchronization, intra-weekly bias correction, and outlier removal, and calculates the inter-epoch difference GF and corresponding DI to form a binary statistical sample. Subsequently, the boundary between the function threshold interval and the constant threshold interval is defined on the DI axis based on the actual distribution density: in the main density interval with low to moderate DI fluctuations, a function-type threshold is used to finely track perturbations; in the extremely high DI sparse region, a constant upper limit is used to control variance. To obtain a robust empirical interface for the "DI-ΔGF amplitude," the 99.7 percentile of ΔGF within each DI step size is calculated according to the 3σ criterion as an approximate "safety boundary."

[0058] During the model building phase, the percentile curves mentioned above are fitted with a quadratic polynomial within the function threshold range. The quadratic form captures the nonlinear expansion relationship between DI and ΔGF while avoiding fluctuations introduced by high-order overfitting. Continuity and (if necessary) monotonicity constraints are applied at the piecewise junctions to ensure a smooth transition of the threshold with DI and maintain physical reasonableness. Simultaneously, a reasonable constant threshold "barrier" is set at the high DI tail to limit the risk of missed detections caused by an infinitely increasing threshold. The final output is a piecewise threshold model: the threshold in the low-to-medium DI region is given by a quadratic function, and a validated upper limit of the constant is used after exceeding the boundary point, thus achieving robust coverage across the entire perturbation domain.

[0059] In terms of system integration, the real-time engine calculates the discrimination threshold (DI) online for each epoch and each satellite / frequency link, and reads the current threshold from the segmented model accordingly to make a threshold decision on the ΔGF. Once the threshold is exceeded, cycle slip marking and subsequent repair (such as ambiguity reset and phase segment splicing under the least squares or Kalman framework) are triggered; if the threshold is not exceeded, phase continuity is maintained and the filtering is carried out. Compared with threshold methods based on fixed thresholds or only on autoregression of historical residuals, this method explicitly introduces the "external disturbance intensity" through DI. Therefore, during the transition period of geomagnetic storms and increased flicker, the discrimination boundary can be adjusted in real time without long adaptive convergence, significantly reducing the alternating "tug-of-war" between instantaneous false alarms and missed detections. The algorithm only involves distribution statistics, percentile estimation, and low-order fitting, with low computational cost, making it suitable for real-time links embedded in receiver firmware or solution software.

[0060] From an industrial application perspective, the adaptive GF cycle slip threshold model can maintain the stability and consistency of cycle slip detection in scenarios such as surveying and mapping as-built layout, high-dynamic machinery control in mines and ports, millimeter-level deformation monitoring of bridges / dams, and power and communication time synchronization, even during periods of soaring Kp exponents and enhanced S4 scintillation. This reduces operational interruptions caused by abnormal triggers, shortens PPP convergence and reconvergence times, and suppresses coordinate drift caused by false detections. Its piecewise fitting structure is easily expandable with constellations and frequency bands, and also facilitates the formation of regionalized configurations of parameter re-estimation libraries based on regional ionospheric climate (high and low latitudes). Through this mechanism, terminals and platforms can achieve robustness to extreme space weather with lower maintenance costs, providing verifiable engineering assurance for the continuous operation of GNSS in critical industries.

[0061] like Figure 1 As shown, the adaptive GF cycle slip detection threshold model method based on DI provided in this embodiment of the invention includes the following steps:

[0062] Step S1: Obtain IGS tracking station data during the severe geomagnetic storm;

[0063] Step S2: Preprocess the data and calculate the inter-epoch difference GF value and DI value;

[0064] Step S3: Based on the distribution of the time-domain difference GF value and DI value, obtain the boundary point between the piecewise model constant threshold and the function threshold interval, and obtain the 99.7 percentile of the time-domain difference GF corresponding to DI in each step of the function threshold interval according to the 3σ criterion;

[0065] Step S4: Perform quadratic polynomial fitting on the percentiles and take a reasonable constant threshold to obtain the final piecewise threshold model;

[0066] Wherein, DI is the Doppler index, GF is the geometric cycle slip detection model, and IGS is the International GNSS Service Organization.

[0067] Step S1 provided in this embodiment of the invention, the process of obtaining IGS tracking station data during a severe geomagnetic storm, includes:

[0068] Raw observation data from IGS tracking stations with dual-frequency 30-second sampling rate throughout the day were selected during typical severe geomagnetic storms to support Doppler observations.

[0069] Step S2 in this embodiment of the invention, which preprocesses the data and calculates the inter-epoch difference GF value and DI value, includes:

[0070] The expression for dual-frequency GF combination (Geometry-Free) phase combination is:

[0071]

[0072] in, and They are frequencies and Phase measurement value in meters. and It is the phase wavelength at the corresponding frequency. and It is phase ambiguity; For ionospheric delay term, This is GF combined phase noise;

[0073] For consecutive epochs Perform the difference to obtain the time-domain difference. ,Right now The expression is:

[0074]

[0075] in, and These represent the GF phase combination values ​​for the next epoch and the previous epoch, respectively;

[0076] By calculating consecutive epochs The difference, without detecting cycle slips, eliminates the ambiguity term in the equation; therefore, the time-domain difference... The size, that is Dependence on time-domain differential residual ionospheric delay ,in It is a time-domain difference operator;

[0077] Under normal circumstances, The values ​​are very small, and the time intervals are short, such as 1 to 30 seconds. Therefore, when using GF combination to detect cycle slips, a tight empirical threshold of 0.05 m is usually set, which performs well during the ionospheric calm period. However, studies have shown that empirical thresholds cannot effectively handle ionospheric anomalies, leading to false cycle slip detections and generating a large number of unnecessary ambiguity initializations, thereby reducing positioning accuracy. Therefore, it is necessary to construct a GF cycle slip detection model that adaptively sets the threshold according to ionospheric changes.

[0078] The threshold model is fitted using the GNSS Doppler index (DI), which is defined as follows:

[0079]

[0080]

[0081]

[0082] in, Indicates frequency, These are Doppler measurements; This represents the rate of change of the geometric distance between the satellite and the receiver. At the speed of light, and These are the clock bias rates of the satellite and receiver, respectively. The ionospheric delay rate, The rate of change of tropospheric delay. For noise, This indicates the average calculation; the unit of DI is TECU / min, consistent with ROTI; in this study, a 5-minute sliding window was used to calculate DI, just like ROTI.

[0083] Step S3 of this embodiment of the invention, which uses the distribution of time-domain difference GF values ​​and DI values ​​to obtain the 99.7 percentile within each interval according to the 3σ criterion, includes:

[0084] Obtained from raw GNSS data of IGS stations selected during typical severe geomagnetic storms. Based on the distribution of DI, the boundary points between the constant threshold and the function threshold intervals should be reasonably set. For the function threshold model interval, a step size of 0.05 TECU / min should be used, and within each step size... The 99.7th percentile (3σ).

[0085] Step S4 provided in this embodiment of the invention involves performing a quadratic polynomial fitting on the percentiles to obtain the final threshold model;

[0086] For the function threshold interval, for DI within each step... A quadratic polynomial was fitted using the 99.7 percentile; for the constant threshold interval, a reasonable constant threshold was selected. This yielded the final piecewise threshold model.

[0087] like Figure 2 As shown, an embodiment of the present invention provides an adaptive GF cycle slip detection threshold model system based on DI, comprising:

[0088] The acquisition module is used to acquire data from IGS tracking stations during severe geomagnetic storms.

[0089] The preprocessing module is used to preprocess the data and calculate the interepoch difference GF value and DI value;

[0090] The distribution module is used to obtain the 99.7 percentile within each interval segment based on the 3σ criterion by using the distribution of time-domain difference GF and DI values.

[0091] The fitting module is used to perform quadratic polynomial fitting on percentiles to obtain the final threshold model.

[0092] This system uses a joint representation of "geometrically independent combination—epochal difference—Doppler index" as its core mechanism. The acquisition module extracts multi-constellation, multi-frequency carrier and Doppler data from the observation stream or RINEX files of the IGS tracking station during severe geomagnetic storms, and completes epoch synchronization and basic quality control. The preprocessing module constructs a geometry-free combination of carrier phases to eliminate geometric and clock bias terms, while highlighting ionospheric and ambiguity terms. Then, epochal difference is applied to the GF, suppressing integer ambiguity and revealing pulse-like characteristics of rapid ionospheric fluctuations and cycle slips. Simultaneously, a dimensionless Doppler Index (DI) is constructed using the carrier Doppler / phase rate as an instantaneous quantification of the "tracking loop disturbance intensity," forming a time series of <ΔGF, DI>, which provides a physical anchor for threshold adaptation.

[0093] The distribution module works by discretizing the DI axis into several step sizes, statistically analyzing the empirical distribution of the corresponding ΔGF amplitude within each step size, and extracting the 99.7 percentile using the 3σ criterion to obtain an empirical profile curve of "perturbation intensity → safe upper bound of ΔGF". Considering that DI has high sample density and significant nonlinearity in the low-to-medium perturbation region, while the sample size is sparse and the variance increases dramatically at the extreme perturbation tail, the system automatically determines the boundary between the function threshold region and the constant threshold region under data-driven conditions: the former retains the functional dependence on DI to delicately follow changes in ionospheric / scintillation intensity, while the latter uses a validated upper limit of the constant to suppress the risk of overfitting and missed detection caused by small samples with high DI. The fitting module performs quadratic polynomial fitting on the 99.7 percentile of each DI step size and applies continuous / monotonic constraints at the segment splicing points, outputting a smooth and stable piecewise threshold model T(DI).

[0094] During operation, the computer device's processor calculates DI and ΔGF online according to satellite-frequency-epoch, calls the coefficients and segmentation rules generated by the fitting module in memory, and provides T(DI) in real time. Cycle slip judgment and marking are performed using |ΔGF|>T(DI), triggering ambiguity reset or fragment splicing repair strategies within the solver. If no boundary is exceeded, phase continuity is maintained and participates in filtering. The above process mainly uses low-order statistics and polynomial calculations, resulting in low computational complexity. It is suitable for deployment in the firmware / software of information data processing terminals, and can also be executed on general-purpose CPUs / embedded SoCs by loading programs onto computer-readable storage media. Relying on the real-time quantification of external disturbances by DI and the adaptive scaling of T(DI), the system significantly reduces the trade-off between false alarms and missed detections during geomagnetic storm transitions, improving the availability and robustness of the PPP / RTK and deformation monitoring links.

[0095] Example 1:

[0096] This embodiment provides an adaptive GF cycle slip detection threshold model based on DI, the method including the following steps:

[0097] Step S1: Obtain IGS tracking station data during the severe geomagnetic storm;

[0098] Step S2: Preprocess the data and calculate the inter-epoch difference GF value and DI value;

[0099] Step S3: Based on the distribution of the time-domain difference GF value and DI value, obtain the boundary point between the piecewise model constant threshold and the function threshold interval, and obtain the 99.7 percentile of the time-domain difference GF corresponding to DI in each step of the function threshold interval according to the 3σ criterion;

[0100] Step S4: Perform quadratic polynomial fitting on the percentiles and take a reasonable constant threshold to obtain the final piecewise threshold model;

[0101] Wherein, DI is the Doppler index, GF is the geometric cycle slip detection model, and IGS is the International GNSS Service Organization.

[0102] In step S1, the process of acquiring IGS tracking station data during severe geomagnetic storms includes:

[0103] Raw observation data from 110 IGS tracking stations that support Doppler observations during the severe geomagnetic storm that occurred on May 11, 2024, were selected based on dual-frequency 30-second sampling rate data throughout the day.

[0104] An extreme geomagnetic storm event occurred on May 11, 2024. Dst, IMF-Bz, and Kp are three commonly used indicators reflecting ionospheric activity. Figure 3 This study presents the time series of the Dst, IMF-Bz, and Kp indices from May 10th to May 12th, 2024. The Dst index rose to 66nT around 17:00 on May 10th, 2024, and then rapidly declined around 18:00 on the same day. The IMF-Bz index also began to fluctuate significantly around this time. Furthermore, the Kp index began to rise after 13:00UT on May 10th, 2024, and rapidly increased to level 9+. Around 2:00 on May 11th, 2024, the Dst index reached its lowest point at -406nT. Therefore, this study uses this event as the basis for model construction and location validation.

[0105] On May 11, 2024, a total of 182 IGS stations supported dual-frequency Doppler observations. Among them, 110 stations (blue dots) were randomly selected for Geometry-Free Combined (GF) cycle slip threshold modeling, and the other 72 stations (red dots) used the established model to verify PPP performance on the same day.

[0106] In step S2, the process of preprocessing the data and calculating the inter-epoch difference GF value and DI value includes:

[0107] The expression for dual-frequency GF combination (Geometry-Free) phase combination is:

[0108]

[0109] in, and They are frequencies and Phase measurement value in meters. and It is the phase wavelength at the corresponding frequency. and It refers to phase ambiguity. For ionospheric delay term, This is the GF combined phase noise.

[0110] For consecutive epochs Perform the difference to obtain the time-domain difference. ,Right now The expression is:

[0111]

[0112] in, and These represent the GF phase combination values ​​for the subsequent epoch and the preceding epoch, respectively. This is achieved by calculating the phase combination values ​​between consecutive epochs. The difference, without detecting cycle slips, eliminates the ambiguity term in the equation. Therefore, time-domain difference... The size, that is Dependence on time-domain differential residual ionospheric delay ,in It is a time-domain difference operator.

[0113] Under normal circumstances, The values ​​are very small, and the time intervals are short, such as 1 to 30 seconds. Therefore, when using GF combination to detect cycle slips, a tight empirical threshold of 0.05 m is usually set, which performs well during periods of ionospheric calm. However, studies have shown that empirical thresholds cannot effectively handle ionospheric anomalies, leading to false cycle slip detections and generating a large number of unnecessary ambiguity initializations, thereby reducing positioning accuracy. Therefore, it is necessary to construct a GF cycle slip detection model that adaptively sets the threshold based on ionospheric changes.

[0114] The threshold model is fitted using the GNSS Doppler index (DI). The definition of DI is as follows:

[0115]

[0116]

[0117]

[0118] in, Indicates frequency, These are Doppler measurements. This represents the rate of change of the geometric distance between the satellite and the receiver. At the speed of light, and These represent the clock bias variation rates for the satellite and receiver, respectively. The ionospheric delay rate, The rate of change of tropospheric delay. For noise, This indicates the average calculation. The unit of DI is TECU / min, consistent with ROTI. In this study, a 5-minute sliding window was used to calculate DI, similar to ROTI.

[0119] In step S3, the process of obtaining the 99.7 percentile within each interval segment using the distribution of time-domain difference GF and DI values ​​according to the 3σ criterion includes:

[0120] Obtained from GNSS data of 110 IGS stations on May 11, 2024 The distribution plot of DI shows that when DI is between 0 and 1 TECU / min, The distribution is relatively concentrated. When 1 < DI ≤ 7 TECU / min, as DI increases, As the value increases, the trend stops when DI = 7 TECU / min. Therefore, a piecewise function is established with 1 and 7 as the dividing points. For DI in the interval of 1~7 TECU / min, a step size of 0.05 TECU / min is used, and values ​​within each step size are taken. The 99.7th percentile (3σ).

[0121] In step S4, a quadratic polynomial fitting is performed on the percentiles to obtain the final threshold model. The process includes:

[0122] For each step of DI within the range of 1 to 7 TECU / min A quadratic polynomial fit was performed using the 99.7 percentile, and constant thresholds were set for cases where DI was less than 1 TECU / min and greater than 7 TECU / min. The final result is as follows: Figure 4 The segmented threshold model shown:

[0123]

[0124] Example 2:

[0125] This embodiment provides corresponding experimental data for Embodiment 1, as follows:

[0126] Taking the PPP accuracy of 72 unmodeled IGS stations during the severe geomagnetic storm on May 11, 2024 as an example:

[0127] First, the dynamic PPP positioning error time series of three IGS stations located in different latitude zones (FRDN (45.93 °N, 66.66 °W), PRDS (50.87 °N, 114.29 °W), and HOFN (64.27 °N, 15.20 °W)) were analyzed. Second, the positioning accuracy of 72 stations was calculated to verify the model's performance. Two schemes were designed for comparison: Scheme 1 uses a traditional fixed GF threshold of 0.05 m, while Scheme 2 uses the adaptive GF threshold model proposed in this paper.

[0128] Figure 5 In the diagram, the blue line represents the result of Scheme 1, and the red line represents the result of Scheme 2. The RMS positioning error for both schemes is also marked in the upper right corner of the subplot. It is worth noting that Scheme 2 significantly improves positioning accuracy compared to Scheme 1. For these three stations at different latitudes, the positioning error time series in the three directions of Scheme 2 are smoother, with fewer outliers and spikes. This is because the adaptive GF threshold model of Scheme 2 filters out a large number of cycle slip misjudgments, thus avoiding a large amount of unnecessary ambiguity reinitialization and greatly improving positioning accuracy.

[0129] Based on the dynamic PPP results of 72 unmodeled stations using two schemes in the horizontal and vertical components as of May 11, 2024, it can be seen that globally, Scheme 2 significantly improves positioning accuracy in both the horizontal and vertical directions. Positioning errors are reduced at most stations, and the positioning accuracy at some stations improves from decimeter-level to centimeter-level.

[0130] The RMS values ​​of the dynamic PPP solution error in the horizontal and vertical components for all 72 non-modeling sites as of May 11, 2024 were calculated and summarized, as shown in Table 1. It is easy to see that the accuracy of Scheme 2 is significantly higher than that of Scheme 1. The model in this paper improves the accuracy of the horizontal and vertical components by 42.4% and 37.2% respectively compared with the traditional fixed threshold.

[0131] Table 1 shows the average root mean square values ​​of dynamic precise single-point positioning in the horizontal and vertical components using Scheme 1 and Scheme 2 as of May 11, 2024.

[0132]

[0133] Example 3:

[0134] This embodiment provides corresponding experimental data for Embodiment 1, as follows:

[0135] Taking the PPP accuracy of 215 IGS stations during the severe geomagnetic storm on October 11, 2024 as an example:

[0136] To verify the applicability of the model presented in this paper to generalized ionospheric activity, data from 215 IGS stations on October 11, 2024, were processed and analyzed. A significant geomagnetic storm also occurred on that day. Similarly, Figure 6 The PPP results for three IGS stations, FRDN (45.93 °N, 66.66 °W), PRDS (50.87 °N, 114.29 °W), and HOFN (64.27 °N, 15.20 °W), are shown. The convergence period of Scheme 1 fluctuates due to cycle slip detection errors, leading to frequent ambiguity resets. By applying the adaptive cycle slip detection threshold model in Scheme 2, greater tolerance is provided in cycle slip detection, largely avoiding false detections and maintaining high-precision positioning performance after convergence.

[0137] Based on the horizontal and vertical RMS values ​​of all experimental stations as of October 11, 2024, a similar conclusion can be drawn: Scheme 2 is superior to Scheme 1, with significant improvements in PPP accuracy for most stations in both the horizontal and vertical components. Table 2 presents the statistical averages of the horizontal and vertical components for Scheme 1 and Scheme 2 as of October 11, 2024. Again, it can be seen that Scheme 2 improves horizontal accuracy by 36.2% and vertical accuracy by 36.8% compared to Scheme 1.

[0138] Table 2 shows the average root mean square values ​​of dynamic precise single-point positioning in the horizontal and vertical components using Scheme 1 and Scheme 2 on October 11, 2024.

[0139]

[0140] It should be noted that embodiments of the present invention can be implemented in hardware, software, or a combination of both. The hardware portion can be implemented using dedicated logic; the software portion can be stored in memory and executed by a suitable instruction execution system, such as a microprocessor or dedicated-design hardware. Those skilled in the art will understand that the above-described devices and methods can be implemented using computer-executable instructions and / or included in processor control code, for example, such code provided on a carrier medium such as a disk, CD, or DVD-ROM, a programmable memory such as read-only memory (firmware), or a data carrier such as an optical or electronic signal carrier. The devices and modules of the present invention can be implemented by hardware circuitry such as very large-scale integrated circuits or gate arrays, semiconductors such as logic chips, transistors, or programmable hardware devices such as field-programmable gate arrays, programmable logic devices, etc., or by software executed by various types of processors, or by a combination of the above-described hardware circuitry and software, such as firmware.

[0141] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any modifications, equivalent substitutions, and improvements made by those skilled in the art within the scope of the technology disclosed in the present invention, and within the spirit and principles of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A method for an adaptive GF cycle slip detection threshold model based on DI, characterized in that, Includes the following steps: Step S1: Obtain IGS tracking station data during the severe geomagnetic storm; Step S2: Preprocess the data and calculate the inter-epoch difference GF value and DI value; Step S3: Using the distribution of time-domain difference GF and DI values, obtain the 99.7 percentile within each interval according to the 3σ criterion; Step S4: Perform quadratic polynomial fitting on the percentiles to obtain the final threshold model; Wherein, DI is the Doppler index, GF is the geometric cycle slip detection model, and IGS is the International GNSS Service Organization.

2. The adaptive GF cycle slip detection threshold model method based on DI as described in claim 1, characterized in that, Step S1, the process of acquiring IGS tracking station data during severe geomagnetic storms, includes: Raw observation data from IGS tracking stations with dual-frequency 30-second sampling rate throughout the day were selected during typical severe geomagnetic storms to support Doppler observations.

3. The adaptive GF cycle slip detection threshold model method based on DI as described in claim 1, characterized in that, Step S2, which involves preprocessing the data and calculating the inter-epoch difference GF and DI values, includes: The expression for dual-frequency GF combination (Geometry-Free) phase combination is: F GF =Φ1-Φ2=(λ1N1-λ2N2)-I GF +e Where Φ1 and Φ2 are the phase measurements in meters at frequencies f1 and f2, respectively, λ1 and λ2 are the phase wavelengths at the corresponding frequencies, and N1 and N2 are the phase ambiguities; ε represents the ionospheric delay term, and GF combined phase noise. For Φ of consecutive epochs GF Perform the difference to obtain the time-domain difference Φ GF That is, ΔΦ GF The expression is: in, and These represent the GF phase combination values ​​for the next epoch and the previous epoch, respectively; By calculating Φ between consecutive epochs GF The difference, without detecting cycle slips, eliminates the ambiguity term in the equation; therefore, the time-domain difference Φ GF The size, i.e., ΔΦ GF Depends on the time-domain differential residual ionospheric delay ΔI GF , where Δ is the time-domain difference operator; The threshold model is fitted using the GNSS Doppler index (DI). The definition of DI is as follows: Where i represents the frequency, and DI is the Doppler measurement value; The rate of change of the geometric distance between the satellite and the receiver is represented by c, where c is the speed of light. and These are the clock bias rates of the satellite and receiver, respectively. The ionospheric delay rate, ε is the rate of change of tropospheric delay. i For noise, <*> indicates average calculation; the unit of DI is TECU / min, consistent with ROTI; in this study, a 5-minute sliding window is used to calculate DI, just like ROTI.

4. The adaptive GF cycle slip detection threshold model method based on DI as described in claim 1, characterized in that, Step S3, which uses the distribution of time-domain difference GF and DI values ​​to obtain the 99.7 percentile within each interval according to the 3σ criterion, includes the following: ΔΦ was obtained from raw GNSS data of IGS stations selected during a typical severe geomagnetic storm. GF And DI, with a step size of 0.05TECU / min, take ΔΦ within each step size. GF The 99.7th percentile (3σ).

5. The adaptive GF cycle slip detection threshold model method based on DI as described in claim 1, characterized in that, In step S4, a quadratic polynomial fitting is performed on the percentiles to obtain the final threshold model. For ΔΦ in each step of DI GF The 99.7 percentile was fitted with a quadratic polynomial, according to ΔΦ. GF The distribution of DI is segmented and constant thresholds are set appropriately; finally, a segmented threshold model is obtained.

6. A DI-based adaptive GF cycle slip detection threshold model system implementing the DI-based adaptive GF cycle slip detection threshold model method as described in any one of claims 1-5, characterized in that, The DI-based adaptive GF cycle slip detection threshold model system includes: The acquisition module is used to acquire data from IGS tracking stations during severe geomagnetic storms. The preprocessing module is used to preprocess the data and calculate the interepoch difference GF value and DI value; The distribution module is used to obtain the 99.7 percentile within each interval segment based on the 3σ criterion by using the distribution of time-domain difference GF and DI values. The fitting module is used to perform quadratic polynomial fitting on percentiles to obtain the final threshold model.

7. A computer device, characterized in that, The computer device includes a memory and a processor. The memory stores a computer program, which, when executed by the processor, causes the processor to perform the steps of the DI-based adaptive GF cycle slip detection threshold model method as described in any one of claims 1-5.

8. A computer-readable storage medium storing a computer program, which, when executed by a processor, causes the processor to perform the steps of the DI-based adaptive GF cycle slip detection threshold model method as described in any one of claims 1-5.

9. An information data processing terminal, characterized in that, The information data processing terminal is used to implement the DI-based adaptive GF cycle slip detection threshold model system as described in claim 5.

Images