Global broadcast ionospheric time delay correction method and device based on occultation observation
By constructing a global broadcast ionospheric model using an ionospheric delay correction method based on occultation observations, the problem of uneven coverage and insufficient robustness caused by reliance on ground-based GNSS observation networks in existing technologies is solved, and efficient and autonomous ionospheric delay correction is achieved.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- INNOVATION ACAD FOR PRECISION MEASUREMENT SCI & TECH CAS
- Filing Date
- 2025-12-22
- Publication Date
- 2026-05-12
AI Technical Summary
Existing broadcast ionospheric models rely on ground-based GNSS observation networks, resulting in uneven global coverage, low correction accuracy in sparse data areas, and insufficient robustness during ionospheric disturbances, making it difficult to achieve efficient, globally covered ionospheric time delay correction.
A global broadcast ionospheric delay correction method is constructed based on occultation observation data. By acquiring electron density profile data obtained from occultation inversion and performing quality control, a single-layer global ionospheric spherical harmonic function model is constructed, and the zero-order spherical harmonic coefficient is used to correct the ionospheric delay, thus eliminating the dependence on ground-based GNSS observation networks.
It achieves efficient, simple, and autonomous ionospheric delay correction globally, enhancing the autonomy and robustness of broadcast ionospheric services, and is suitable for data-sparse regions and ionospheric disturbance environments.
Smart Images
Figure CN122017882A_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of ionospheric time delay correction technology, and particularly relates to a global broadcast ionospheric time delay correction method and device based on occultation observation. Background Technology
[0002] The ionosphere is one of the most challenging sources of error affecting the performance of global navigation satellite systems (GNSS). Broadcasting global ionospheric correction parameters in satellite navigation ephemeris is the primary technique used by various GNSS systems to achieve ionospheric delay correction. To meet the real-time correction needs of global single-frequency users, major GNSS systems such as the US GPS, the EU Galileo, and China's BeiDou all broadcast ionospheric model parameters through navigation ephemeris. GPS uses the 8-parameter Klobuchar model, Galileo uses the 3-parameter NeQuick-G model, and the BeiDou-3 system uses the 9-parameter BeiDou Global Broadcast Ionospheric Model (BDGIM). These broadcast ionospheric models have played a crucial role in improving single-frequency positioning accuracy. However, the parameter estimation of existing broadcast ionospheric models is highly dependent on ground-based GNSS observation networks. This dependence faces numerous challenges in practical applications, limiting the model's performance consistency and dynamic adaptability globally. First, the spatial distribution of global ground-based GNSS observation stations is severely uneven. Existing monitoring stations are mostly concentrated on land and in economically developed regions, while observation stations are scarce or even absent in vast oceans, polar regions, and many developing countries. This uneven spatial coverage directly affects the global correction accuracy of the model. Secondly, in actual operation, due to limitations in data sharing mechanisms, operational management policies, or communication conditions, the collection of global observation data may be delayed or incomplete. For example, the NeQuick-G and NTCM-G models require daily assimilation of global GNSS observation data to estimate ionization level parameters. If data from key areas cannot be accessed in a timely manner, it will directly affect the accuracy and timeliness of the model parameters. To address the problem of insufficient data in poorly observed areas, some systems have introduced climatological background models based on long-term statistics as a supplementary means. For example, while BDGIM primarily uses ground-based GNSS data from China and surrounding areas for parameter estimation, it uses background field models to provide TEC estimates in other regions. However, these models are built based on historical averages and lack the ability to respond to short-term changes in the ionosphere. During periods of high solar activity or under disturbances such as geomagnetic storms, the ionospheric state changes drastically, and climatological models struggle to capture its dynamic characteristics, leading to a significant decrease in correction effectiveness and insufficient model robustness. Therefore, the core problem facing current broadcast ionospheric models is: how to construct a new type of modified model with simple parameters, high computational efficiency, global coverage, and good dynamic adaptability without relying on global or regional ground-based GNSS observation networks, especially to maintain stable and reliable performance in data-sparse regions and during ionospheric disturbances.
[0003] In recent years, with the rapid development of low Earth orbit satellite constellations such as FY-3, Tianmu, Yunyao, and COSMIC-2, ionospheric observations based on GNSS radio occultation technology have demonstrated unique advantages. Occultation observations are characterized by near-global coverage, lack of geographical limitations, and high vertical resolution, effectively filling the observation gaps of ground-based GNSS in regions such as the ocean and polar regions, and providing a new data foundation for constructing global ionospheric models independent of ground networks.
[0004] Current research has attempted to refine global ionospheric spherical harmonic function (SHHF) models using occultation data, but significant shortcomings remain. Firstly, existing studies do not adequately consider the differences in spatial altitude coverage between occultation-derived TECs and ground-based GNSS TECs, making it difficult to directly serve single-frequency GNSS users. Secondly, SHHF-based models typically require a large number of parameters (e.g., at least 16 parameters for orders ≥3), far exceeding the parameter capacity of current broadcast models, which is detrimental to real-time broadcasting. Therefore, effectively utilizing multi-source occultation observation data, achieving high-quality data preprocessing, driving lightweight ionospheric models, and constructing a global broadcast correction system independent of ground-based GNSS data remain key technical challenges for improving the autonomy, continuity, and global consistency of GNSS services. Summary of the Invention
[0005] The purpose of this invention is to overcome the above-mentioned defects and problems in the prior art and provide a method and device for global broadcast ionospheric delay correction based on occultation observation. The method is free from dependence on ground-based GNSS observation networks and has the advantages of simple model structure, high computational efficiency, few broadcast parameters, and strong global coverage, which significantly enhances the autonomy and robustness of broadcast ionospheric services.
[0006] To achieve the above objectives, the technical solution of the present invention is:
[0007] In a first aspect, the present invention provides a global broadcast ionospheric delay correction method based on occultation observations, comprising:
[0008] Obtain electron density profile data retrieved from occultation observations and perform quality control to obtain qualified electron density profile data;
[0009] Ionospheric TEC data is obtained by integrating the electron density in qualified electron density profile data and then converting the ionospheric TEC data to the height range of ground-based GNSS observations to obtain normalized ionospheric TEC data.
[0010] Based on normalized ionospheric TEC data, a single-layer global ionospheric harmonic function model is constructed, and the zero-order spherical harmonic coefficient is obtained by inversion estimation using this model;
[0011] A global broadcast ionospheric function model was constructed, and the global broadcast ionospheric function model was solved based on the zero-order spherical harmonic coefficients to obtain the fixed parameters of the global broadcast ionospheric function model;
[0012] The zero-order spherical harmonic coefficients are broadcast as ionospheric delay correction parameters, and the ionospheric delay correction is calculated using the broadcast ionospheric delay correction parameters and the fixed parameters of the global broadcast ionospheric function model.
[0013] Preferably, the spherical harmonic function model of the single-layer global ionosphere is:
[0014] ;
[0015] In the formula, For vertical ionospheric delay; and These are the latitude and longitude of the puncture point in the Sun-fixed geomagnetic coordinate system; for Spend The regularized Legendre function of order 1; and These are the coefficients of the global ionospheric harmonic function model to be estimated.
[0016] Preferably, the step of obtaining qualified electron density profile data through quality control includes:
[0017] Calculate the relative rate of change of electron density between adjacent height points for each effective electron density profile;
[0018] The height axis is divided into several intervals. For each effective electron density profile, a robust percentile index of the relative electron density change rate is statistically analyzed in each height interval to generate noise upper limit threshold, sawtooth amplitude threshold and extreme mutation threshold.
[0019] Based on the noise upper limit threshold, extreme mutation threshold, and sawtooth amplitude threshold, a three-level progressive judgment process is established, including high noise background criteria, extreme single-point jump criteria, and valid sawtooth structure criteria. If any criterion is triggered, the electron density profile is determined to be invalid and discarded.
[0020] Preferably, the method for calculating the relative electron density change rate is as follows:
[0021] For each valid electron density profile, calculate the absolute value of electron density between adjacent height points:
[0022] ;
[0023] In the formula, This represents the absolute difference in electron density. For height point electron density; For height point electron density;
[0024] Normalizing based on the smaller electron density value between two adjacent points, we obtain the relative rate of change of electron density:
[0025] ;
[0026] In the formula, This represents the rate of change in relative electron density.
[0027] Preferably, the step of generating a robust percentile index for the rate of change of relative electron density within each height interval, and generating a noise upper limit threshold, a sawtooth amplitude threshold, and an extreme abrupt change threshold, includes:
[0028] Within each height interval, the 75th percentile, 95th percentile, and maximum value of the relative electron density change rate are statistically analyzed to obtain the first, second, and third statistics, respectively. The first, second, and third statistics of all electron density profiles within the same height interval are summarized. Then, the 75th percentile of the first statistic is taken to obtain the sawtooth amplitude threshold, the 75th percentile of the second statistic is taken to obtain the noise upper limit threshold, and the 75th percentile of the third statistic is taken and multiplied by 1.2 to obtain the extreme mutation threshold.
[0029] Preferably, the high-noise background criterion is:
[0030] The proportion of the relative electron density change rate exceeding the upper limit threshold of the noise level. :
[0031] ;
[0032] In the formula, It represents the rate of change of relative electron density; For height point The upper limit threshold for noise; This represents the total number of data points in the electron density profile data.
[0033] If the proportion exceeds the preset threshold, it indicates that the electron density profile is contaminated by high-frequency noise, and the electron density profile is deemed invalid.
[0034] Preferably, the extreme single-point jump criterion is:
[0035] If there are two or more height points within the valid height range that satisfy:
[0036] ;
[0037] In the formula, It represents the rate of change of relative electron density; For height point The extreme mutation threshold;
[0038] This indicates that the electron density profile contains a non-physical abrupt change, and the electron density profile is deemed invalid.
[0039] Preferably, the criterion for the effective serrated structure is:
[0040] Calculate the sign change of the electron density difference sequence to locate potential sawtooth vertices;
[0041] For each sign flip point, if at least one of the two adjacent relative electron density change rates is not lower than the sawtooth amplitude threshold, then it is a valid sawtooth vertex.
[0042] The maximum continuous length of the effective sawtooth vertices is counted. If the maximum continuous length is greater than or equal to the set threshold, it indicates that there is a systematic, non-physical periodic oscillation in the electron density profile, and the electron density profile is judged to have sawtooth anomalies.
[0043] Preferably, the normalized ionospheric TEC data for:
[0044] ;
[0045] In the formula, For ionospheric TEC data; TEC data for the altitude range from the ground to GNSS satellite orbit; This is TEC data for the altitude range from ground level to low Earth orbit satellite orbit.
[0046] Secondly, the present invention provides a global broadcast ionospheric delay correction device based on occultation observation, the device being used to implement the aforementioned global broadcast ionospheric delay correction method based on occultation observation, the device comprising:
[0047] The electron density profile data acquisition module is used to acquire electron density profile data retrieved from occultation observations and to perform quality control to obtain qualified electron density profile data.
[0048] The ionospheric TEC data acquisition module is used to integrate the electron density in qualified electron density profile data to obtain ionospheric TEC data, and convert the ionospheric TEC data to the height range of ground-based GNSS observations to obtain normalized ionospheric TEC data.
[0049] The zero-order spherical harmonic coefficient acquisition module is used to construct a single-layer global ionospheric harmonic function model based on normalized ionospheric TEC data, and to obtain the zero-order spherical harmonic coefficients by inversion estimation through this model.
[0050] The solidified parameter acquisition module is used to construct a global broadcast ionospheric function model and solve the global broadcast ionospheric function model based on the zero-order spherical harmonic coefficients to obtain the solidified parameters of the global broadcast ionospheric function model;
[0051] The ionospheric delay correction acquisition module is used to broadcast the zero-order spherical harmonic coefficient as the ionospheric delay correction parameter, and to calculate the ionospheric delay correction number using the broadcast ionospheric delay correction parameter and the fixed parameters of the global broadcast ionospheric function model.
[0052] Compared with the prior art, the beneficial effects of the present invention are as follows:
[0053] This invention discloses a method and apparatus for global broadcast ionospheric delay correction based on occultation observations. The method eliminates reliance on ground-based GNSS observation networks, utilizing globally covered GNSS radio occultation observation data to construct an ionospheric correction model. Simultaneously, by introducing a rigorous data quality control mechanism, particularly a novel method for identifying and eliminating sawtooth electron density profiles, it effectively solves the problem of traditional quality control methods failing to detect non-physical oscillation anomalies. Furthermore, the use of TEC normalization processing technology addresses the technical challenge of inconsistent spatial altitude coverage between occultation observation TEC and ground-based GNSS observation TEC. The model constructed by this invention possesses advantages such as simple structure, high computational efficiency, fewer broadcast parameters, and strong global coverage, significantly improving the autonomy, continuity, and robustness of broadcast ionospheric correction services. It is particularly suitable for data-sparse regions and complex space environments such as ionospheric disturbances. Attached Figure Description
[0054] Figure 1 This is a flowchart of a global broadcast ionospheric delay correction method based on occultation observations proposed in an embodiment of the present invention.
[0055] Figure 2 This is a flowchart illustrating the quality control of electron density profile data as proposed in an embodiment of the present invention.
[0056] Figure 3 This is a schematic diagram illustrating the analysis of Case 1, which shows the qualified electron density profile proposed in the embodiments of the present invention.
[0057] Figure 4 This is a schematic diagram illustrating the analysis of Case 2, which shows the qualified electron density profile proposed in the embodiments of the present invention.
[0058] Figure 5 This is a schematic diagram illustrating the analysis of Case 1, a non-compliant electron density profile proposed in an embodiment of the present invention.
[0059] Figure 6 This is a schematic diagram illustrating the analysis of Case 2, which shows an electron density profile of substandard quality, as presented in an embodiment of the present invention.
[0060] Figure 7 This is a schematic diagram illustrating the analysis of Case 3, a non-compliant electron density profile proposed in an embodiment of the present invention.
[0061] Figure 8 This is a structural block diagram of a global broadcast ionospheric delay correction device based on occultation observation proposed in an embodiment of the present invention.
[0062] Figure 9 This is a structural block diagram of a global broadcast ionospheric delay correction device based on occultation observation proposed in an embodiment of the present invention. Detailed Implementation
[0063] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.
[0064] See Figure 1 This invention provides a global broadcast ionospheric delay correction method based on occultation observations, comprising:
[0065] S1. Obtain electron density profile data retrieved from occultation observations and perform quality control to obtain qualified electron density profile data; specifically, obtain electron density profile data retrieved from radio occultation (RO) observations by low Earth orbit (LEO) satellites within a preset time window before the target time. The electron density profile data includes the time of occurrence of the occultation event, the peak height of the maximum electron density, the maximum electron density, and electron density values at different altitudes.
[0066] S2. Based on the three-dimensional electron density empirical model Nequick, the electron density in the qualified electron density profile data is integrated to obtain ionospheric TEC data, and the ionospheric TEC data is converted to the height range of ground-based GNSS observation to obtain normalized ionospheric TEC data.
[0067] S3. Based on normalized ionospheric TEC data, a single-layer global ionospheric harmonic function model is constructed, and the zero-order spherical harmonic coefficients are obtained by inversion estimation through this model.
[0068] S4. Construct a global broadcast ionospheric function model and solve the global broadcast ionospheric function model based on the zero-order spherical harmonic coefficients to obtain the fixed parameters of the global broadcast ionospheric function model.
[0069] S5. The zero-order spherical harmonic coefficients are broadcast as ionospheric delay correction parameters. The broadcast ionospheric delay correction parameters and the fixed parameters of the global broadcast ionospheric function model are used to calculate the ionospheric delay correction number.
[0070] Currently, the broadcast ionospheric models of major satellite navigation systems all rely on ground-based GNSS observation data for parameter fitting. Therefore, under special conditions where ground-based observation data is missing or unavailable (such as communication interruptions, station failures, or areas with no coverage), the models struggle to maintain continuous and reliable correction capabilities. This invention addresses the problem of the heavy dependence of existing broadcast ionospheric models on global or regional ground-based GNSS observation networks during parameter estimation by proposing a global broadcast ionospheric delay correction method based on occultation observations. This invention only requires broadcasting a single broadcast coefficient, significantly reducing the communication load and parameter update frequency of the navigation satellite system. Simultaneously, it relies entirely on space-based radio occultation observation data, without depending on ground monitoring networks, effectively solving the performance degradation problem of traditional ground-dependent models under conditions of regional data shortage. This provides a highly reliable and robust ionospheric delay correction solution for future global navigation systems operating in wide-area conditions without ground coverage.
[0071] Furthermore, under certain conditions, the vertical electron density profile obtained from occultation inversion may exhibit non-physical sawtooth fluctuations. These anomalies often stem from numerical instability, noise amplification, or multipath effects in the inversion algorithm, and do not represent the true ionospheric structure. Failure to identify and remove these anomalies will directly impact the accuracy and stability of subsequent modeling. Current occultation data processing workflows generally lack specific checks for such sawtooth profiles, leading to the misuse of anomalous data and weakening the reliability of the model.
[0072] See Figure 2 Quality control is performed to obtain qualified electron density profile data, including:
[0073] S101. Calculate the relative rate of change of electron density between adjacent height points of each effective electron density profile.
[0074] S102. Divide the height axis into several intervals. For each effective electron density profile, calculate the robust percentile index of the relative electron density change rate in each height interval to generate the noise upper limit threshold, sawtooth amplitude threshold and extreme mutation threshold.
[0075] S103. Based on the noise upper limit threshold, extreme mutation threshold, and sawtooth amplitude threshold, a three-level progressive judgment process is formulated, including high noise background criteria, extreme single-point jump criteria, and valid sawtooth structure criteria. If any criterion is triggered, the electron density profile is determined to be invalid and discarded.
[0076] To overcome the problem that a globally fixed threshold cannot adapt to the non-uniformity of ionospheric height, this invention constructs a height-stratified dynamic threshold system based on historical high-quality electron profile data (such as long-term observations from COSMIC, Tianmu-1, Fengyun-3, the Yunyao constellation, and ionospheric altimeters). This invention is applicable to the automatic quality screening of multi-source observation data such as radio occultation (RO), ionospheric altimeter, and GNSS inversion, and is particularly useful for identifying non-physical structural anomalies such as high-frequency sawtooth oscillations and abrupt changes.
[0077] The core idea of this invention is that local changes in the real ionospheric profile exhibit high dependence and statistical regularity, while non-physical anomalies manifest as random or periodic, drastic fluctuations unrelated to altitude. Based on this, this invention proposes a detection strategy combining segmented dynamic thresholds and multiple criteria fusion. Specifically, it includes: constructing a highly adaptive dynamic threshold system: utilizing a large amount of historical "normal" profile data, calculating robust percentile indices of relative electron density change rates for each altitude range to generate three types of thresholds—noise upper limit threshold, sawtooth amplitude threshold, and extreme abrupt change threshold; designing a three-level progressive discrimination process: sequentially screening for high-noise background, extreme single-point jumps, and valid sawtooth structures; triggering any condition results in an invalid profile. Compared to traditional methods, this invention has the advantages of high adaptability, low false alarm rate, and high detection rate, making it particularly suitable for quality control of multi-source ionospheric data in complex space environments.
[0078] Furthermore, the original electron density profile data is first processed using conventional quality control procedures to remove obviously invalid data. Then, through manual visual screening, abnormal profiles with non-physical oscillation characteristics such as "sawtooth" are further removed, ultimately forming a pure "normal profile" training set, i.e., an effective electron density profile dataset.
[0079] Furthermore, the method for calculating the relative rate of change of electron density is as follows:
[0080] For each valid electron density profile, calculate the absolute value of electron density between adjacent height points:
[0081] ;
[0082] In the formula, This represents the absolute difference in electron density. For height point electron density; For height point electron density;
[0083] Normalizing based on the smaller electron density value between two adjacent points, we obtain the relative rate of change of electron density:
[0084] ;
[0085] In the formula, It represents the rate of change of relative electron density; This indicates the smaller electron density value between two adjacent points.
[0086] Furthermore, robust percentile indices of relative electron density change rate are statistically analyzed within each height interval to generate noise upper limit threshold, sawtooth amplitude threshold, and extreme abrupt change threshold, including:
[0087] The rate of change of relative electron density was statistically analyzed within each altitude interval (e.g., altitude intervals were divided into intervals of 50 km along the altitude axis). The 75th percentile (P75), 95th percentile (P95), and maximum value (Max) are used to obtain the first, second, and third statistics, respectively. These statistics are then combined across all electron density profiles within the same height range. The 75th percentile of the first statistic is used to obtain the sawtooth amplitude threshold (P75_of_P75), the 75th percentile of the second statistic is used to obtain the noise upper limit threshold (P75_of_P95), and the 75th percentile of the third statistic is multiplied by 1.2 to obtain the extreme mutation threshold (1.2×P75_of_Max). This results in a dynamic threshold table that varies with height, used for anomaly detection in subsequent profiles.
[0088] This invention, through hierarchical statistics and a robust percentile strategy, accurately identifies non-physical jagged anomalies while preserving normal physical fluctuations, significantly improving the adaptability and reliability of quality control.
[0089] Based on the aforementioned dynamic thresholds, this invention designs a three-level progressive criterion, whereby the profile is deemed invalid if any criterion is triggered. First, the upper limit threshold for noise is used to assess whether the profile is in an overall high-noise state. Then, the extreme mutation threshold is used to check whether there are drastic single-point jumps in the profile. Finally, the sawtooth amplitude threshold is used to identify more deceptive "regular sawtooth" anomalies (manifested as alternating rising / falling oscillations).
[0090] Furthermore, the high-noise background criterion (overall fluctuation intensity screening) is as follows:
[0091] The proportion of the relative electron density change rate exceeding the upper limit threshold of the noise level. :
[0092] ;
[0093] In the formula, It represents the rate of change of relative electron density; For height point The upper limit threshold for noise; This represents the total number of data points in the electron density profile data.
[0094] If the proportion exceeds a preset threshold (e.g., 20%), it indicates that the electron density profile is contaminated by high-frequency noise and lacks physical reliability, thus the electron density profile is deemed invalid.
[0095] Furthermore, the extreme single-point transition criterion (non-physical transition screening) is as follows:
[0096] Extreme mutation thresholds are used to identify points with abnormally high relative electron density change rates. If two or more height points exist within the effective height range that satisfy the following conditions:
[0097] ;
[0098] In the formula, It represents the rate of change of relative electron density; For height point The extreme mutation threshold;
[0099] This indicates that the electron density profile contains a non-physical abrupt change, and the electron density profile is deemed invalid.
[0100] Furthermore, the effective sawtooth structure criterion (periodic oscillation identification) adopts a joint discrimination strategy of structure + amplitude, specifically as follows:
[0101] (1) Detect sign reversal. Calculate the sign change of the electron density difference sequence to locate potential sawtooth vertices (i.e., locate potential local maxima / minimum points);
[0102] (2) Verify that the amplitude meets the standard. For each sign flip point, if at least one of the two adjacent relative electron density change rates is not lower than the sawtooth amplitude threshold zigzagThrVec, then it is a valid sawtooth vertex;
[0103] (3) Continuity assessment. The maximum continuous length of the effective sawtooth vertices is statistically analyzed. If the maximum continuous length is greater than or equal to the set threshold (the set threshold is 5), it indicates that there is a systematic, non-physical periodic oscillation in the electron density profile, and the electron density profile is judged to have sawtooth anomalies.
[0104] Furthermore, the essence of sawtooth anomalies is that electron density alternately rises and falls with height, forming a continuous sequence of local extrema. This invention accurately locates potential sawtooth vertices through the following steps; the method for determining the potential sawtooth vertices is as follows:
[0105] (1) Calculate the first-order difference (reflecting the local trend of change):
[0106] Given electron density profile and corresponding height The first-order difference (i.e., change) between adjacent height points is calculated as follows:
[0107] ;
[0108] In the formula, This represents the change in electron density between two adjacent altitude points; For height point electron density; For height point electron density;
[0109] (2) Extract the difference symbol sequence :
[0110] ;
[0111] (3) Check if the symbol has changed:
[0112] Calculate the difference of symbol sequences :
[0113] ;
[0114] like This indicates that at the altitude point At this point, the trend of electron density change reverses:
[0115] If the difference sign sequence changes from +1 to -1, it indicates the height point. This is a local maximum (peak);
[0116] If the difference sign sequence changes from -1 to +1, it indicates the height point. This is a local minimum point (valley);
[0117] The locations of local maxima and minima are potential sawtooth vertices. By combining amplitude and continuity constraints, the true ionospheric structure and pseudo-oscillations can be effectively distinguished.
[0118] This invention, through the fusion of segmented dynamic thresholds and multi-level criteria, can efficiently and accurately identify non-physical sawtooth oscillations and extreme jumps, significantly improving the quality control capability of data before ionospheric modeling; it has the following advantages: (1) Strong adaptability: the threshold is dynamically adjusted with height, taking into account both the stability of the E layer and the high gradient characteristics of the F layer; (2) Low false alarm rate: through structural features (sign flip) and amplitude dual verification, it avoids misjudging the real gradient as an anomaly; (3) High detection rate: it covers a variety of anomaly types such as high noise, extreme jumps, and periodic oscillations; (4) Strong engineering applicability: the algorithm is simple and easy to integrate into the existing ionospheric data processing flow.
[0119] The conventional electron density profile quality control method uses the following four quality control indicators to screen the raw electron density profile data.
[0120] (1) Average relative deviation:
[0121] ;
[0122] In the formula, The average relative deviation; This represents the total number of electron density samples in a single profile; For the first Electron density inversion value at each sampling point; This is the background value obtained by filtering with a 7-point moving average.
[0123] If a certain cross-section If the value is greater than 0.25, it is judged as an abnormal profile and is removed.
[0124] (2) Noise factor:
[0125] ;
[0126] In the formula, Noise factor; This refers to the number of sampling points within an altitude range above 300km.
[0127] like If the profile is deemed to have excessive fluctuations in the height region, it is marked as a problem profile.
[0128] (3) Top electron density gradient:
[0129] ;
[0130] In the formula, The top electron density gradient; and The values are electron density at altitudes of 490 km and 420 km, respectively, with an altitude interval of 70 km.
[0131] like If the density at the top decreases too slowly, it may reflect an inversion bias, and the profile is marked as a problem profile.
[0132] (4) and Reasonable range:
[0133] Based on existing research experience It is usually located between 180 and 450 km, while The reasonable range is to If any parameter exceeds the above range, the entire profile will be considered abnormal and discarded.
[0134] Based on the above four criteria, any electron density profile that meets any of the rejection conditions is marked as a problem profile.
[0135] Figure 3 and Figure 4 The application effect of the abnormal profile quality discrimination method proposed in this invention is demonstrated by taking the profiles ionPrf_TM01.2024.215.10.48.E24_0001.0001.nc (quality qualified electron density profile case 1) and ionPrf_TM01.2024.215.00.01.G03_0001.0001.nc (quality qualified electron density profile case 2) as examples. Figure 3 (a) and Figure 4 (a) is a schematic diagram of electron density distribution with height, from Figure 3 (a) and Figure 4 (a) It can be seen that both electron density profiles exhibit smooth and continuous variation characteristics, without obvious sawtooth oscillations; Figure 3 (b) and Figure 4 (b) is a schematic diagram showing the distribution of electron density change and dynamic threshold with height. Figure 3 (b) and Figure 4 (b) It can be seen that the corresponding relative electron density change rate does not exceed the three dynamic thresholds (sawtooth amplitude threshold, noise upper limit threshold, and extreme mutation threshold) set by the present invention throughout the entire height range. The results show that the present invention can accurately identify and retain qualified profiles that conform to physical properties.
[0136] Figure 5 Taking the typical sawtooth-shaped abnormal electron density profile ionPrf_TM01.2024.215.00.57.R24_0001.0001.nc (Case 1 of substandard electron density profile) as an example, Figure 6 Taking the typical sawtooth-shaped abnormal electron density profile ionPrf_TM02.2024.215.00.57.R20_0001.0001.nc (Case 2 of substandard electron density profile) as an example, Figure 7 Taking the typical sawtooth-shaped abnormal electron density profile ionPrf_TM01.2024.215.00.35.G05_0001.0001.nc (Case 3 of unqualified electron density profile) as an example, the application effect of the abnormal profile quality discrimination method proposed in this invention is demonstrated. Figure 5 (a) Figure 6 (a) and Figure 7 (a) is a schematic diagram of electron density distribution with height. Figure 5 (b) Figure 6 (b) and Figure 7 (b) is a schematic diagram showing the distribution of electron density change and dynamic threshold with height. Figure 5 (a) Figure 6 (a) and Figure 7(a) It can be seen that all three cross sections exhibit obvious non-physical sawtooth-like undulation characteristics, among which Figure 5 (a) and Figure 6 The oscillations in the cross-section shown in (a) are particularly violent; by Figure 5 (b) Figure 6 (b) and Figure 7 (b) It is evident that the corresponding relative electron density change rate significantly exceeds the three dynamic thresholds set by this invention in multiple height ranges—including the sawtooth amplitude threshold (P75_of_P75), the noise upper limit threshold (P75_of_P95), and the extreme mutation threshold (1.2×P75_of_Max). The results show that this invention can effectively identify such high-frequency oscillation anomalies and mark them as unqualified profiles.
[0137] Table 1 lists the traditional quality control indicators corresponding to the three typical sawtooth-shaped anomalous electron density profiles mentioned above. Among them, the R24 profile, with a maximum height of only 489 km, does not fully cover the 420-490 km interval, therefore the top gradient cannot be calculated; the remaining indicators are as follows: the mean relative deviation (MRD) are 0.07, 0.08, and 0.03, respectively, all far below the rejection threshold of 0.25; the upper-level noise factor Δ is respectively... , and Both are less than 0.01; the top gradients of R20 and G05 are respectively and Both meet the "rapid descent" requirement; hmF2 (294.95-350.00 km) and NmF2 ( All three sections are within the physically reasonable range. Therefore, under existing industry standards, these three sections, which clearly exhibit jagged anomalies, did not trigger any traditional quality control criteria and were incorrectly judged as "qualified" sections.
[0138] Table 1. Statistics of quality discrimination indicators for traditional methods
[0139]
[0140] This comparison clearly demonstrates that traditional methods are ineffective in identifying jagged structural anomalies, while the present invention significantly improves the detection capability of such hidden anomalies by fusing segmented dynamic thresholds with structural features.
[0141] Furthermore, normalized ionospheric TEC data for:
[0142] ;
[0143] In the formula, TEC data of the ionosphere obtained by surface integration of electron density profile from occultation; TEC data for the ground-to-GNSS satellite orbital altitude range calculated based on the NeQuick background model; This is TEC data for the orbital altitude range from ground to low Earth orbit (LEO) satellites, calculated based on the NeQuick background model.
[0144] By employing TEC normalization processing technology, the technical challenge of inconsistent spatial altitude coverage between occultation observation TEC and ground-based GNSS observation TEC was resolved.
[0145] Furthermore, based on the ionospheric delay observations extracted from occultation data, a single-layer global ionospheric spherical harmonic function model with a 2-hour resolution is constructed daily. The mathematical expression of this model is:
[0146] ;
[0147] In the formula, For vertical ionospheric delay; and These are the latitude and longitude of the puncture point in the Sun-fixed geomagnetic coordinate system; for Spend The regularized Legendre function of order 1; and These are the coefficients of the global ionospheric harmonic function model to be estimated.
[0148] The model is solved using the least squares algorithm to obtain the spherical harmonic coefficients of each order. The first coefficient is the zeroth-order spherical harmonic coefficient. This invention performs high-order modeling (5*5) on occultation data to solve for the coefficients of the spherical harmonic function model, and uses the zeroth-order spherical harmonic coefficient (SHCO) as the broadcast coefficients of the MNTCM broadcast model, thus achieving a balance between limiting the broadcast coefficients and ensuring model accuracy.
[0149] Furthermore, the mathematical expression for the global broadcast ionospheric function model is:
[0150] ;
[0151] ;
[0152] ;
[0153] ;
[0154] In the formula, and These represent the variations of the ionosphere with local time and season, respectively. and These are, respectively, the variation of the ionosphere with geomagnetic latitude and the low-latitude anomalies of the ionosphere; The impact of solar activity on the ionosphere; The latitude of the puncture point in the solar-fixed geomagnetic coordinate system; Local time; It is accumulated over a year; These are the broadcast parameters for the broadcast ionospheric model; , , , , , , , , , , , These are the 12 solidification parameters for the broadcast ionosphere model; The angle of the sun's deflection; This refers to geomagnetic latitude. and This refers to the latitude and precision of the geomagnetic north pole.
[0155] Furthermore, the zero-order spherical harmonic coefficients obtained from the occultation data of the most recent day are used as broadcast parameters for the broadcast ionospheric delay correction method of the global satellite navigation system for message encoding and transmission to each GNSS satellite, with an update cycle of 1 day.
[0156] Furthermore, the ionospheric delay correction is calculated using the broadcast ionospheric delay correction parameters and the fixed parameters of the global broadcast ionospheric function model, including:
[0157] The S501 GNSS receiver receives broadcast ephemeris data from various satellites, including satellite orbits, clock biases, and ionospheric delay correction parameters sent to the satellites.
[0158] S502. Calculate the positions of GNSS satellites and receivers based on GNSS satellite orbit and clock bias information.
[0159] S503. Calculate the vertical total electron content of the ionosphere at the puncture point using the ionospheric delay correction parameters and the solidification parameters of the global broadcast ionospheric function model.
[0160] S504. Convert the vertical total electron content to the tilted total electron content to determine the ionospheric delay correction number on the signal propagation path.
[0161] Specifically, the vertical total electron content is converted to the tilted total electron content according to the following formula:
[0162] ;
[0163] In the formula, This is the projection function of the ionosphere; The radius of the Earth; The height of the ionosphere thin layer; This is the satellite's elevation angle.
[0164] Referring to Table 2, Klobuchar represents the broadcast ionospheric delay correction technique used by GPS; RO-MNTCM-SHC0 represents the global broadcast ionospheric model established in this invention; BDGIM represents the ionospheric delay correction technique used by the BeiDou-3 system; and NTCM-G represents the ionospheric delay correction technique used by Galileo. Among these, RO-MNTCM-SHC0 (this invention), Klobuchar, BDGIM, and NTCM-G broadcast 1, 8, 9, and 3 broadcast coefficients, respectively. As can be seen from Table 2, compared with other ionospheric delay correction techniques, this invention significantly reduces the number of broadcast parameters. Furthermore, the ionospheric delay accuracy of this invention is roughly equivalent to the NTCM-G model during years of high solar activity and is significantly superior to the Klobuchar and BDGIM models. The reason for the higher accuracy of this invention compared to other ionospheric delay correction techniques is that it adopts zero-order spherical harmonic function model coefficients calculated based on occultation data as input parameters, which more accurately reflects the overall trend of ionospheric variation. Furthermore, this invention, based on historical average data of global ionospheric total electron content covering one or more solar activity cycles, embeds parameters characterizing the temporal and spatial variations of global ionospheric total electron content into the receiver. This makes it possible to directly calculate the ionospheric delay at any time and location globally by inputting only a single broadcast parameter. This invention fundamentally eliminates the dependence of existing broadcast models on ground-based GNSS data and can significantly reduce the number of broadcast parameters, thereby alleviating the communication burden on satellites.
[0165] Table 2. Comparison of Ionospheric Delay Correction Accuracy of Different Broadcast Ionospheric Delay Correction Techniques
[0166]
[0167] See Figure 8 The present invention also provides a global broadcast ionospheric delay correction device based on occultation observations. The device is used to implement the aforementioned global broadcast ionospheric delay correction method based on occultation observations. The device includes:
[0168] The electron density profile data acquisition module is used to acquire electron density profile data retrieved from occultation observations and to perform quality control to obtain qualified electron density profile data.
[0169] The ionospheric TEC data acquisition module is used to integrate the electron density in qualified electron density profile data to obtain ionospheric TEC data, and convert the ionospheric TEC data to the height range of ground-based GNSS observations to obtain normalized ionospheric TEC data.
[0170] The zero-order spherical harmonic coefficient acquisition module is used to construct a single-layer global ionospheric harmonic function model based on normalized ionospheric TEC data, and to obtain the zero-order spherical harmonic coefficients by inversion estimation through this model.
[0171] The solidified parameter acquisition module is used to construct a global broadcast ionospheric function model and solve the global broadcast ionospheric function model based on the zero-order spherical harmonic coefficients to obtain the solidified parameters of the global broadcast ionospheric function model;
[0172] The ionospheric delay correction acquisition module is used to broadcast the zero-order spherical harmonic coefficient as the ionospheric delay correction parameter, and to calculate the ionospheric delay correction number using the broadcast ionospheric delay correction parameter and the fixed parameters of the global broadcast ionospheric function model.
[0173] See Figure 9 The present invention also provides a global broadcast ionospheric delay correction device based on occultation observation, including a memory and a processor;
[0174] The memory is used to store computer program code and transmit the computer program code to the processor;
[0175] The processor is configured to execute, according to instructions in the computer program code, a global broadcast ionospheric delay correction method based on occultation observations as described above.
[0176] The present invention also provides a computer-readable storage medium storing a computer program, which, when executed by a processor, implements the global broadcast ionospheric delay correction method based on occultation observation as described above.
[0177] Generally, the computer instructions for implementing the method of the present invention can be carried on any combination of one or more computer-readable storage media. Non-transitory computer-readable storage media can include any computer-readable medium except for the signal itself, which is temporarily propagating.
[0178] Computer-readable storage media can be, for example, but not limited to, electrical, magnetic, optical, electromagnetic, infrared, or semiconductor systems, apparatuses, or any combination thereof. More specific examples (a non-exhaustive list) of computer-readable storage media include: electrical connections having one or more wires, portable computer disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EKROM or flash memory), optical fibers, portable compact disk read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination thereof. In this invention, a computer-readable storage medium can be any tangible medium containing or storing a program that can be used by or in conjunction with an instruction execution system, apparatus, or device.
[0179] Computer program code for performing the operations of this invention can be written in one or more programming languages or a combination thereof. These programming languages include object-oriented programming languages—such as Java, Smalltalk, and C++—as well as conventional procedural programming languages—such as the "C" language or similar programming languages. In particular, Python, suitable for neural network computation, and platform frameworks such as TensorFlow and PyTorch can be used. The program code can be executed entirely on the user's computer, partially on the user's computer, as a standalone software package, partially on the user's computer and partially on a remote computer, or entirely on a remote computer or server. In cases involving remote computers, the remote computer can be connected to the user's computer or to an external computer (e.g., via the Internet using an Internet service provider) through any type of network, including a local area network (LAN) or a wide area network (WAN).
[0180] For details regarding the aforementioned equipment and non-transitory computer-readable storage media, please refer to the specific description of a global broadcast ionospheric delay correction method based on occultation observations and its beneficial effects, which will not be repeated here.
[0181] Although embodiments of the present invention have been shown and described above, it should be understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those skilled in the art can make changes, modifications, substitutions and variations to the above embodiments within the scope of the present invention.
Claims
1. A global broadcast ionospheric delay correction method based on occultation observations, characterized in that, include: Obtain electron density profile data retrieved from occultation observations and perform quality control to obtain qualified electron density profile data; Ionospheric TEC data is obtained by integrating the electron density in qualified electron density profile data and then converting the ionospheric TEC data to the height range of ground-based GNSS observations to obtain normalized ionospheric TEC data. Based on normalized ionospheric TEC data, a single-layer global ionospheric harmonic function model is constructed, and the zero-order spherical harmonic coefficient is obtained by inversion estimation using this model; A global broadcast ionospheric function model was constructed, and the global broadcast ionospheric function model was solved based on the zero-order spherical harmonic coefficients to obtain the fixed parameters of the global broadcast ionospheric function model; The zero-order spherical harmonic coefficients are broadcast as ionospheric delay correction parameters, and the ionospheric delay correction is calculated using the broadcast ionospheric delay correction parameters and the fixed parameters of the global broadcast ionospheric function model.
2. The global broadcast ionospheric delay correction method based on occultation observations according to claim 1, characterized in that, The spherical harmonic function model of the single-layer global ionosphere is as follows: ； In the formula, For vertical ionospheric delay; and These are the latitude and longitude of the puncture point in the Sun-fixed geomagnetic coordinate system; for Spend The regularized Legendre function of order 1; and These are the coefficients of the global ionospheric harmonic function model to be estimated.
3. The global broadcast ionospheric delay correction method based on occultation observations according to claim 1, characterized in that, The process of obtaining qualified electron density profile data through quality control includes: Calculate the relative rate of change of electron density between adjacent height points for each effective electron density profile; The height axis is divided into several intervals. For each effective electron density profile, a robust percentile index of the relative electron density change rate is statistically analyzed in each height interval to generate noise upper limit threshold, sawtooth amplitude threshold and extreme mutation threshold. Based on the noise upper limit threshold, extreme mutation threshold, and sawtooth amplitude threshold, a three-level progressive judgment process is established, including high noise background criteria, extreme single-point jump criteria, and valid sawtooth structure criteria. If any criterion is triggered, the electron density profile is determined to be invalid and discarded.
4. The global broadcast ionospheric delay correction method based on occultation observations according to claim 3, characterized in that, The method for calculating the relative electron density change rate is as follows: For each valid electron density profile, calculate the absolute value of electron density between adjacent height points: ； In the formula, This represents the absolute difference in electron density. For height point electron density; For height point electron density; Normalizing based on the smaller electron density value between two adjacent points, we obtain the relative rate of change of electron density: ； In the formula, This represents the rate of change in relative electron density.
5. The global broadcast ionospheric delay correction method based on occultation observations according to claim 3, characterized in that, The robust percentile index of the relative electron density change rate within each height interval is used to generate noise upper limit threshold, sawtooth amplitude threshold, and extreme abrupt change threshold, including: Within each height interval, the 75th percentile, 95th percentile, and maximum value of the relative electron density change rate are statistically analyzed to obtain the first, second, and third statistics, respectively. The first, second, and third statistics of all electron density profiles within the same height interval are summarized. Then, the 75th percentile of the first statistic is taken to obtain the sawtooth amplitude threshold, the 75th percentile of the second statistic is taken to obtain the noise upper limit threshold, and the 75th percentile of the third statistic is taken and multiplied by 1.2 to obtain the extreme mutation threshold.
6. The global broadcast ionospheric delay correction method based on occultation observations according to claim 3, characterized in that, The high-noise background criterion is: The proportion of the relative electron density change rate exceeding the upper limit threshold of the noise level. : ； In the formula, It represents the rate of change of relative electron density; For height point The upper limit threshold for noise; This represents the total number of data points in the electron density profile data. If the proportion exceeds the preset threshold, it indicates that the electron density profile is contaminated by high-frequency noise, and the electron density profile is deemed invalid.
7. The global broadcast ionospheric delay correction method based on occultation observations according to claim 3, characterized in that, The extreme single-point jump criterion is: If there are two or more height points within the valid height range that satisfy: ； In the formula, It represents the rate of change of relative electron density; For height point The extreme mutation threshold; This indicates that the electron density profile contains a non-physical abrupt change, and the electron density profile is deemed invalid.
8. The global broadcast ionospheric delay correction method based on occultation observations according to claim 3, characterized in that, The criterion for an effective sawtooth structure is: Calculate the sign change of the electron density difference sequence to locate potential sawtooth vertices; For each sign flip point, if at least one of the two adjacent relative electron density change rates is not lower than the sawtooth amplitude threshold, then it is a valid sawtooth vertex. The maximum continuous length of the effective sawtooth vertices is counted. If the maximum continuous length is greater than or equal to the set threshold, it indicates that there is a systematic, non-physical periodic oscillation in the electron density profile, and the electron density profile is judged to have sawtooth anomalies.
9. The global broadcast ionospheric delay correction method based on occultation observations according to claim 1, characterized in that, The normalized ionospheric TEC data for: ； In the formula, For ionospheric TEC data; TEC data for the altitude range from the ground to GNSS satellite orbit; This is TEC data for the altitude range from ground level to low Earth orbit satellite orbit.
10. A global broadcast ionospheric delay correction device based on occultation observations, characterized in that, The apparatus is used to implement the method according to any one of claims 1 to 9, the apparatus comprising: The electron density profile data acquisition module is used to acquire electron density profile data retrieved from occultation observations and to perform quality control to obtain qualified electron density profile data. The ionospheric TEC data acquisition module is used to integrate the electron density in qualified electron density profile data to obtain ionospheric TEC data, and convert the ionospheric TEC data to the height range of ground-based GNSS observations to obtain normalized ionospheric TEC data. The zero-order spherical harmonic coefficient acquisition module is used to construct a single-layer global ionospheric harmonic function model based on normalized ionospheric TEC data, and to obtain the zero-order spherical harmonic coefficients by inversion estimation through this model. The solidified parameter acquisition module is used to construct a global broadcast ionospheric function model and solve the global broadcast ionospheric function model based on the zero-order spherical harmonic coefficients to obtain the solidified parameters of the global broadcast ionospheric function model; The ionospheric delay correction acquisition module is used to broadcast the zero-order spherical harmonic coefficient as the ionospheric delay correction parameter, and to calculate the ionospheric delay correction number using the broadcast ionospheric delay correction parameter and the fixed parameters of the global broadcast ionospheric function model.