A dual-frequency TEC inversion optimization method based on real-time height data

Abstract

This invention belongs to the field of satellite navigation and space environment monitoring technology, and relates to a dual-frequency TEC inversion optimization method based on real-time altimetry data. Based on GNSS ground-based dual-frequency observations, it integrates ionospheric altimetry data and absolute signal deviation products to achieve high-precision conversion of STEC (Slant Total Electron Content) to VTEC (Vertical Total Electron Content). The method includes data preprocessing, TEC inversion, cycle slip detection, deviation correction, and mapping correction, and generates ROTI (Rate of TEC Index, i.e., ionospheric scintillation characterization parameter) products. By introducing a single-layer model for correcting the effective ionospheric height, the ability to characterize ionospheric disturbances is improved.

Description

Technical Field

[0001] This invention belongs to the field of satellite navigation and space environment monitoring technology, specifically relating to a method for ionospheric TEC inversion and correction based on GNSS data. Background Technology

[0002] TEC (Special Emission Terrain) is a crucial physical parameter describing the state and spatiotemporal variations of the ionosphere, directly impacting the delay effect of satellite navigation signals during propagation. With the widespread application of GNSS in high-precision positioning, navigation and timing, and space weather monitoring, accurately acquiring ionospheric TEC information is of significant practical importance for improving navigation and positioning accuracy, analyzing ionospheric disturbance characteristics, and ensuring stable system operation. Based on the frequency-dependent characteristics of ionospheric delay in dual-frequency GNSS observation data, the STEC along the satellite signal propagation path can be inverted without relying on external empirical models. This method offers advantages such as continuous observation, wide coverage, and low cost, and has become one of the important technical means for current ionospheric monitoring.

[0003] In existing TEC inversion studies, methods and processing procedures based on the GPS system are relatively mature, while studies based on BDS (BeiDou Navigation Satellite System) are relatively few. The BeiDou system, with its multi-constellation configuration and multi-frequency observation signals, offers a richer variety of observation types and frequency combinations, demonstrating significant potential for TEC inversion. However, this also places higher demands on data processing methods and model adaptability. Therefore, developing a stable and applicable TEC inversion procedure for the BeiDou system remains of considerable research and application value.

[0004] In existing dual-frequency TEC retrieval methods, a single-layer ionospheric model is typically used, assuming that the ionospheric electron density is concentrated on a thin shell at a fixed height. A geometric mapping function is then used to convert STEC to VTEC. However, the actual ionospheric electron density distribution exhibits significant spatiotemporal variations, and its equivalent height changes with latitude, local time, and spatial environmental conditions. While the fixed-height assumption simplifies the calculation process to some extent, it can easily introduce systematic errors when the ionospheric structure is complex or highly perturbed, thus limiting further improvements in retrieval accuracy.

[0005] With the development of ionospheric detection technology, observation methods such as ionospheric altimeters can provide data reflecting the characteristics of ionospheric height changes, providing conditions for characterizing the dynamic changes of the equivalent ionospheric height. Introducing ionospheric altimeter data into the VTEC inversion process to correct the mapping relationship in the traditional single-layer model is theoretically feasible, helping to reduce errors caused by the fixed height assumption and improve the reliability of VTEC inversion results.

[0006] It should be noted that even after considering variations in ionospheric height and improving the mapping model, TEC inversion is still inevitably affected by various factors. These include the spatiotemporal resolution limitations of ionospheric altimetry data, consistency issues between different observation data, and noise, biases, and anomalous observations inherent in the GNSS observation data itself. All of these factors can impact the final results. Therefore, how to further optimize the TEC inversion process while ensuring model stability and engineering feasibility remains a key issue in current research. Summary of the Invention

[0007] The purpose of this invention is to provide a dual-frequency TEC inversion optimization method based on real-time altimetry data, in order to solve the problems of traditional dual-frequency GNSS-based TEC inversion process, which does not consider ionospheric height changes, has limited mapping model accuracy, and lacks stability of inversion results.

[0008] The technical solution of the present invention is as follows:

[0009] A dual-frequency TEC inversion optimization method based on real-time altimetry data includes the following steps:

[0010] S1: Set the daily data missing rate requirement, download the preliminary data, read and preprocess the multi-source data;

[0011] S2: To address the potential mismatch between GNSS observation data and GNSS navigation data, this invention further designs and introduces a satellite data validity screening mechanism. It performs initial inversion calculations of ionospheric TEC based on dual-frequency GNSS observation data, and calculates the satellite's spatial position and clock bias information at a given observation epoch using conventional broadcast ephemeris calculation methods, based on the read navigation message data and observation time information.

[0012] S3: After obtaining the satellite spatial position and receiver spatial coordinates corresponding to the observation epoch, to further describe the spatial geometric relationship between the satellite and the receiver, it is necessary to calculate the elevation and azimuth angles of available satellites relative to the receiver at each observation epoch. To obtain higher quality observation data, a satellite elevation angle cutoff threshold is set to remove low-quality observation data. After completing the calculation of the satellite elevation and azimuth angles relative to the receiver, to reduce the impact of multipath effects introduced by low-elevation satellites and tropospheric and ionospheric mapping errors on the ionospheric total electron content inversion results, this invention performs a screening process on the observation data based on the elevation angle threshold.

[0013] S4: Cycle slip error detection and correction in carrier observation, and on this basis, smoothing and calibration processing of pseudorange and carrier observation STEC;

[0014] S5: To further improve the accuracy of ionospheric TEC inversion results, this invention uses OSB data to correct satellite-side errors in the smoothed results;

[0015] S6: At this stage, the present invention introduces an automatic identification and elimination mechanism for STEC abnormal segments, and further implements integrity constraint checks on this basis;

[0016] S7: Introduce observation data from the digital altimeter of the ionosphere and use the key ionospheric height parameter hmF2 to adaptively correct the mapping function from STEC to VTEC;

[0017] S8: The receiver DCB estimation method based on minimizing VTEC dispersion is used to invert the differential code deviation at the receiver end;

[0018] S9: The final TEC result undergoes zero-baseline adjustment, which means that the minimum TEC value throughout the entire time period is used as a reference, and the TEC data is shifted as a whole. The shift amount is the absolute value of the current satellite's minimum value, thereby avoiding negative values ​​or inconsistencies in the baseline caused by residual system bias. Attached Figure Description

[0019] The dates of the following figures are all July 7, 2024. To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the figures used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the figures described below are only some embodiments of the present invention. For those skilled in the art, other figures can be obtained based on these figures without creative effort.

[0020] Figure 1 This is a flowchart illustrating the optimization method of the traditional dual-frequency TEC inversion based on altimeter data according to the present invention.

[0021] Figure 2 This is a schematic diagram of the initial inversion results of dual-frequency observation TEC in a specific embodiment;

[0022] Figure 3 This is a schematic diagram of the TEC smoothing calibration results in a specific embodiment;

[0023] Figure 4 This is a schematic diagram of TEC after removing satellite system errors in a specific embodiment;

[0024] Figure 5 This is a schematic diagram of satellite system error statistics in a specific embodiment;

[0025] Figure 6 This is a schematic diagram of receiver system error results in a specific embodiment;

[0026] Figure 7 This is a schematic diagram of the TEC zero baseline adjustment result in a specific embodiment;

[0027] Figure 8 This is a schematic diagram of the ROTI results for a 5-minute time window in a specific embodiment. Detailed Implementation

[0028] The present invention will be further described below with reference to the accompanying drawings and embodiments. It should be noted that, unless otherwise specified, the embodiments and technical features described in this application can be combined with each other. It should also be pointed out that, unless otherwise indicated, all technical and scientific terms used in this application have the same meaning as commonly understood by one of ordinary skill in the art to which this application pertains. The terms "comprising" or "including" and similar words used in this invention refer to elements or objects preceding the word that encompass the elements or objects listed following the word and their equivalents, without excluding other elements or objects.

[0029] like Figure 1 As shown, this invention provides a dual-frequency TEC inversion optimization method based on real-time altimetry data, comprising the following steps:

[0030] S1: Obtain RINEX format observation data from the ionospheric scintillation monitors of domestic observation stations via the Meridian Project Phase II Data Center, and simultaneously acquire the corresponding RINEX format navigation message data; concurrently, download ionospheric parameters from digital altimeters within or adjacent areas of the observation region, used to extract and construct improved model data for the key ionospheric height parameter hmF2. Obtain OSB data within the corresponding observation time range via the IGS Ionospheric Analysis Center of Wuhan University. The functions for monitoring various missing observation data and supplementing missing data are explained in detail below:

[0031] The system checks for instances where a 5-minute data segment is missing. If such an instance occurs, all missing data is assigned the value NaN, and the missing timestamps and other data are continuously added to ensure the normal operation of subsequent calculations. It also checks for instances where all or some satellite data is lost within a very short period of time. If such an instance occurs, the data is automatically assigned the value NaN.

[0032] In the RINEX 3.03 observation data version, erroneous satellite data will be set to 0. Therefore, it is necessary to check for similar erroneous data. If it occurs, the corresponding satellite data at that epoch needs to be assigned the value NaN.

[0033] In traditional models for TEC inversion based on dual-frequency GNSS data, we need to obtain the receiver's geocentric and geofixed coordinates from the observation data and convert them into geographic coordinates including longitude, latitude, and altitude. However, in high-resolution observation files, 288 raw coordinates are recorded in a single day. To conform to the model's calculation method, this invention takes the average of several raw coordinates read and uses them as the receiver coordinates in the model.

[0034] The navigation data also follows the RINEX 3.03 data format. The invention uses a self-written function to implement the initial reading function, and nests it with the observation data reading function to achieve automatic reading.

[0035] S2: In a specific example, the matching of observation and navigation data for different stations or at different times varies. In this example, satellite 14 corresponding to the ZAZP station has navigation data but no corresponding observation data, resulting in the ability to determine only the satellite-to-ground geometric parameters but not to perform preliminary TEC inversion. Similarly, satellite 62 only has observation data but no corresponding navigation data, meaning only preliminary TEC inversion can be performed, but further parameter correction is not possible. Therefore, in the data preprocessing stage, this invention matches the satellite numbers in the observation and navigation data, retaining only satellites with both valid observation and corresponding navigation data within the same observation epoch, and automatically discarding satellite data that does not meet the above conditions, thereby avoiding invalid satellites from participating in ionospheric parameter inversion calculations. Considering the signal system characteristics of the BeiDou satellite system and the completeness and stability of existing observation data, channels 2 and 6 are selected from multiple available observation channels as a dual-frequency combination, using C2I and C6I pseudorange observations and L2I and L6I carrier observations for calculation. The reason for choosing the above channel combination is that the signal frequency intervals corresponding to channels 2 and 6 are relatively large, which can effectively enhance the sensitivity of the dual-frequency ionospheric delay differential to TEC; at the same time, the availability and temporal continuity of the observation data corresponding to these two types of channels are better, which is conducive to the unified correction of subsequent systematic errors. This step includes the following formula:

[0036]

[0037] Where STECp and STECl are the inversion results of pseudorange observation and carrier observation, respectively; k is the conversion coefficient from dual-frequency difference to total electron content of the oblique ionosphere; C6I, C2I, L2I, and L6I are the observed values ​​at each frequency point; f2 is the frequency of the second channel of the Beidou satellite, which is taken as 1561.098 MHz here; f6 is the frequency of the sixth channel of the Beidou satellite, which is taken as 1268.52 MHz here; and 40.3 is the ionospheric refractive constant.

[0038] By performing dual-frequency combined calculations using pseudorange and carrier observations respectively, two types of initial STEC results with different error characteristics can be obtained, such as... Figure 2 As shown in the figure, STECp is unaffected by integer ambiguity but has a relatively high noise level; STECl has high measurement accuracy and time continuity but contains unknown integer ambiguity terms. These two types of STEC results provide fundamental data support for subsequent carrier smoothing pseudorange results and bias correction.

[0039] This method, based on the broadcast ephemeris parameters corresponding to the ZAZP station, uses the Kepler orbital model to calculate the satellite's mean perihelion angle and obtains the off-perihelion angle by iteratively solving the Kepler equation. It further combines orbital perturbation correction parameters to calculate the satellite's position in the orbital plane. A rotation matrix is ​​constructed using the right ascension of the ascending node, orbital inclination, and argument of perihelion to transform the satellite coordinates from the orbital coordinate system to the geocentric-fixed coordinate system. Simultaneously, a satellite clock bias polynomial model and relativistic effect corrections are introduced to obtain the high-precision satellite position and satellite clock bias corresponding to the station. The formula for calculating the relativistic clock bias correction term is as follows:

[0040]

[0041] Where t r This is the relativistic clock error correction value, F is the relativistic constant, e is the satellite orbital eccentricity, a is the orbital semi-major axis, and E is the anomalous angle.

[0042] S3: In a specific embodiment, to obtain the observation geometric relationship parameters between the satellite and the receiver, the elevation and azimuth angles of the satellite relative to the receiver are calculated based on the position coordinates of the satellite and the receiver in the ECEF (Earth-Centered, Earth-Fixed coordinate system). First, the geocentric coordinates of the receiver and the satellite are converted into geodetic coordinates to obtain the longitude and latitude information of the receiver. Then, a local station-centered coordinate system is constructed with the receiver's location as the origin, and a rotational transformation relationship from the geocentric coordinate system to the station-centered coordinate system is established based on the receiver's longitude and latitude. By performing coordinate transformation on the relative position vector between the satellite and the receiver, the east, north, and zenith components of the satellite in the station-centered coordinate system are obtained. Then, the satellite elevation angle is calculated based on the geometric relationship between the zenith component and the horizontal component, and the satellite azimuth angle is calculated based on the projection relationship of the horizontal component. Through the above steps, the elevation and azimuth information of the satellite can be accurately obtained at each observation epoch, providing the necessary geometric parameter support for subsequent elevation threshold screening, ionospheric path modeling, and STEC to VTEC mapping calculation.

[0043] Furthermore, to mitigate the adverse effects of ionospheric horizontal gradient enhancement, tropospheric delay amplification, and multipath effects on low-elevation satellite signals during propagation, this invention sets a satellite elevation angle screening condition for satellite observation data participating in the ionospheric total electron content inversion, setting the satellite elevation angle cutoff threshold to 15°. Specifically, after completing the satellite elevation angle calculation, the elevation angle information of each observation epoch and the corresponding satellite is judged. When the satellite elevation angle is less than 15°, it is considered that the slant path of the observation path through the ionosphere is too long, the observation error is significantly amplified, and the reliability of the corresponding pseudorange and carrier observation values ​​is low. Therefore, this type of low-elevation satellite data is removed from the subsequent TEC inversion processing. When the satellite elevation angle is greater than or equal to 15°, its observation data is retained for subsequent calculations. When the satellite elevation angle is lower than the above-mentioned extreme elevation angle threshold of 15°, the reliability of the corresponding observation data is considered low; for satellite observation data whose elevation angle meets the threshold condition, it is retained for subsequent processing. By constructing an elevation angle mask matrix, observation epochs that do not meet the elevation angle conditions and the corresponding satellite observation data are marked as invalid values, and this is simultaneously applied to the STECp obtained from pseudorange observations, the STECl obtained from carrier observations, and the corresponding elevation angle, azimuth angle, and time information.

[0044] By using the unified screening mechanism based on elevation angle, low-quality observation data can be effectively eliminated without changing the data time structure, ensuring the consistency of different physical quantities in time and satellite dimensions, thereby improving the stability and reliability of subsequent ionospheric parameter inversion and mapping results.

[0045] S4: In a specific embodiment, the corresponding STECl time series is first extracted for each valid satellite by its number, and invalid observation epochs with empty time labels are removed. For each satellite, the difference in STECl between adjacent epochs is calculated, and the standard deviation of this difference sequence is used as an adaptive threshold reference. The cycle slip determination threshold is set within the range of 1-2 TECu (Total Electron Content unit). When the absolute value of the STECl difference between adjacent epochs exceeds the above threshold, a carrier cycle slip is determined to have occurred at that location.

[0046] Meanwhile, to distinguish between the discontinuity caused by normal continuous observations and long-term observation interruptions, this invention further introduces a time-time-based criterion. When the time interval between adjacent observation epochs exceeds a preset threshold (set to 3 hours in this embodiment), the STECl sequence is considered to have a time break, and the entire data segment is divided into several time-continuous subsequences, which are then subjected to cycle slip detection and correction processing respectively.

[0047] During cycle slip correction, for each detected cycle slip position, the STECl subsequences before and after the cycle slip are treated as two consecutive segments. By calculating the offset of the STECl values ​​at the boundary between the two segments, the preceding STECl segment is shifted as a whole to ensure continuity with the initial value of the following STECl segment at the cycle slip position. This eliminates the step error caused by abrupt changes in carrier integer ambiguity to the STECl sequence. For isolated missing data points that occur during observation, STECl values ​​from adjacent epochs are used for interpolation to maintain sequence continuity.

[0048] After completing cycle slip detection and correction for carrier STECl, to further eliminate residual integer ambiguity constants in carrier observations, the observation epoch with the largest elevation angle is selected as the reference point in the continuous observation sequence of each satellite. The cycle slip-corrected carrier STECl sequence is then shifted as a whole to return to the original observation level while maintaining temporal continuity.

[0049] The smoothing process is as follows: First, using a single satellite as the processing object, the entire-day STECl time series is filtered for valid observation epochs, and the data is judged to have long-term interruptions based on the time interval between adjacent valid epochs. When the time interval between two adjacent epochs exceeds a preset threshold (e.g., 3 hours), the satellite observation sequence is considered to be discontinuous in the time dimension and is divided into multiple independent time segments; if there is no time interval exceeding the threshold, the entire-day observation is treated as a continuous time window for processing.

[0050] Subsequently, within each continuous time window, a difference sequence between STECp and STECl is constructed, reflecting the systematic shift of STECp relative to STECl. To address the potential outlier problem in the difference sequence, this embodiment employs a robust statistical bias estimation strategy: first, the median of the difference sequence is calculated as a central reference value, and then the median absolute deviation is calculated to describe the discrete characteristics of the data; based on this, an outlier discrimination threshold is constructed to remove outlier difference samples that significantly deviate from the central value, retaining only reliable observation data for bias estimation. For cases with excessive outliers or degraded difference sequences, an adaptive degradation strategy is employed to ensure the stability of the bias estimation.

[0051] After obtaining the robust bias estimate within the current time window, this bias is applied as a translation amount to the STECl sequence within the corresponding time period, achieving overall translational alignment of STECl to STECp. The results are as follows: Figure 3 As shown.

[0052] This processing method effectively eliminates the absolute deviation caused by the initial ambiguity while maintaining the high precision characteristics of the carrier wave, so that the smoothed STEC has both good absoluteness and time continuity.

[0053] By integrating carrier cycle slip detection, correction, and pseudorange-carrier STEC smoothing calibration processes, the stability and reliability of STEC time series can be significantly improved, providing high-quality input data for subsequent VTEC mapping and ionospheric parameter inversion.

[0054] S5: In a specific embodiment, firstly, based on OSB data, the corresponding frequency point deviations are combined to calculate the required relative deviation values ​​between the C6I and C2I frequency points. Then, the obtained frequency point deviations are uniformly converted from nanosecond-level to time delay quantities. Furthermore, considering the physical relationship that ionospheric delay is inversely proportional to the square of the electromagnetic wave frequency, the speed of light constant and ionospheric scaling factor are introduced to perform a scaling transformation on the time deviation, thereby converting the differential code deviation at the satellite end from the time domain to the TEC unit domain. The results are as follows: Figure 5 As shown, and the formula is as follows:

[0055]

[0056] Where Bias_Time is the satellite-end bias time, Bias_TEC is the bias TEC caused by the satellite-end bias time, OSB_C6I and OSB_C2I are the absolute signal bias values ​​at the corresponding frequency points, in nanoseconds (ns), so they need to be converted to seconds (s) later; c is the speed of light; constant A is a proportionality coefficient of 40.3 derived from the refractive properties of ionospheric plasma; f1 and f2 are the frequencies corresponding to the frequency points, respectively, taken as 1561.098 MHz and 1268.52 MHz.

[0057] After completing the unit and scale conversions described above, the STEC bias term for each satellite is constructed, and a corresponding satellite-end differential code bias sequence is generated for each valid satellite. Finally, this satellite differential code bias is subtracted epoch-by-epoch from the STEC observations that already include satellite bias, effectively correcting for systematic errors at the satellite end and obtaining STEC data free from the influence of satellite bias. Specific results are as follows: Figure 4 As shown.

[0058] S6: In a specific embodiment, firstly, statistical analysis is performed on the STEC time series of each satellite, and the steady-state background level of STEC is constructed by calculating the epoch median reference sequence of the full satellite observation data.

[0059] Based on this, for the STEC sequences of a single satellite, the data is segmented according to temporal continuity: when the time interval between adjacent valid observation epochs exceeds a preset threshold (e.g., 1 hour), it is considered as a different continuous observation segment and judged separately. For continuous observation sequences without obvious time breaks, the overall deviation of the satellite's STEC sequence from the median reference is compared. If its average deviation exceeds a preset threshold, it is determined that there is a systematic anomaly in the satellite's observations for that day, and the corresponding STEC sequence is removed entirely.

[0060] For STEC sequences with time breaks, anomaly detection is performed on each consecutive observation segment. Specifically, the deviation between the STEC and the median reference within each consecutive observation segment is counted. When the number of epochs with deviations exceeding a set threshold accounts for a certain proportion of the valid observation epochs in that segment, the consecutive observation segment is determined to be an anomaly segment, and the entire data in that segment is removed.

[0061] Finally, the number of valid observation epochs for each satellite at the end of each day is counted. If the number of valid observations is lower than a preset threshold, it is considered that the satellite has too few observation data in that period and the data reliability is low. In this case, the corresponding STEC data will be removed as a whole.

[0062] S7: In a specific embodiment, ionospheric digital altimeter observation data is introduced, and the mapping function from STEC to VTEC is adaptively corrected using the key ionospheric height parameter hmF2:

[0063]

[0064] Where S(E) is the mapping factor obliquely to VTEC, R is the average radius of the Earth, h is the peak height of the F2 layer, and E is the satellite elevation angle.

[0065] S8: In a specific embodiment, firstly, the receiver DCB search range is set to -30 to 30 nanoseconds, with a step size of 0.1 nanoseconds, and a systematic scan of the original STEC is performed. For each candidate DCB value, it is converted into a deviation in TEC units, and this deviation is subtracted from the TEC data after removing the satellite DCBs to obtain the corrected STEC. Subsequently, the STEC is converted to VTEC using a mapping factor, and the cumulative standard deviation of VTEC for all satellites within the same time period is calculated. The DCB corresponding to the minimum standard deviation is used as the optimal estimate of the receiver deviation, and the search step size is gradually reduced through multiple iterations to achieve high-precision inversion of the receiver DCB. The results are as follows. Figure 6 As shown.

[0066] After the receiver DCB estimation is completed, it is applied to the absolute TEC calculation. Specifically, the estimated receiver DCB value is subtracted from the smoothed calibrated TEC containing receiver bias to obtain the corrected STEC; the corrected VTEC is then calculated by combining the mapping factor.

[0067] S9: In a specific embodiment, firstly, using the minimum STEC value after correction for the satellite-to-ground system error as a reference, its absolute value is used as an offset and added to both STEC and VTEC, thereby obtaining non-negative absolute STEC and absolute VTEC. The result is as follows. Figure 7 As shown.

[0068] Based on this, the ROTI (Representation of Ionospheric Scintillation Indicator) was calculated. In the ROTI calculation process, firstly, short-time data gaps in the STEC time series were linearly interpolated to ensure the continuity of the time series; then, STEC was downsampled according to a preset time window, and the rate of change of TEC between adjacent epochs was calculated; finally, the standard deviation of the rate of change series within the sliding time window was calculated to obtain the ROTI value for the corresponding epoch. The results are as follows: Figure 8 As shown.

[0069] In summary, this invention addresses the needs of high-resolution GNSS ionospheric monitoring by constructing a method for ionospheric TEC inversion and scintillation characterization that integrates multi-source observation data, based on the traditional dual-frequency TEC inversion method. RINEX format observation data and corresponding navigation messages from domestic GNSS observation stations are acquired through the Meridian Project Phase II Data Center. Simultaneously, observation parameters from ionospheric digital altimeters are incorporated to construct an improved mapping model including the peak height hmF2 of the F2 layer. Combined with the OSB product provided by the IGS Ionospheric Analysis Center of Wuhan University, accurate correction of satellite-side system errors is achieved. This invention also addresses the common problems of missing data, erroneous values, and temporal discontinuities in high temporal resolution observation data by designing a multi-level data integrity detection and repair mechanism. While maintaining temporal structure consistency, it effectively eliminates low-quality observations, providing a stable and reliable data foundation for subsequent inversion.

[0070] In the TEC inversion process, this invention combines the signal characteristics of the BeiDou satellite system, selecting C2I / C6I pseudorange observations and L2I / L6I carrier observations to construct a dual-frequency combination, obtaining initial STEC values ​​with different error characteristics respectively. Satellite position and clock errors are then accurately calculated using broadcast ephemeris and orbital dynamics models, establishing a high-precision satellite-to-ground geometric relationship. Based on this, the system error introduced by low-elevation angle observations is effectively reduced through elevation and azimuth angle calculations and elevation angle masking constraints. Addressing the cycle slip and integer ambiguity issues in carrier observations, this invention proposes a cycle slip detection and correction strategy that integrates time continuity criteria and robust statistical methods. Through pseudorange-carrier STEC smoothing calibration, a STEC time series with both high precision and absolute consistency is achieved.

[0071] Furthermore, this invention introduces an adaptive oblique-to-VTEC mapping model based on the hmF2 parameter, improving the physical rationality of VTEC projection under the traditional single-layer ionosphere assumption. Regarding receiver-side system error processing, a receiver DCB calculation method based on minimizing VTEC dispersion is adopted. Multiple iterative searches are performed within the range of -30 to 30 nanoseconds to achieve high-precision estimation of the receiver DCB, which is then uniformly applied to the correction of STEC and VTEC. Subsequently, a non-negative absolute TEC sequence is constructed through a minimum zero-point shifting strategy, ensuring the continuity and comparability of the TEC time series.

[0072] Building upon this foundation, the present invention further calculates the ionospheric scintillation characterization index ROTI. By interpolating, downsampling, and performing sliding window statistics on the corrected STEC time series, the time-varying characteristics of small-scale irregularities in the ionosphere are successfully characterized. Experimental results show that the present invention outperforms existing technologies in terms of TEC inversion accuracy, temporal continuity, systematic error suppression, and ionospheric scintillation characterization capability, demonstrating significant technological advancement and engineering application value.

[0073] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some modifications or alterations to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent changes, and alterations made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the scope of the present invention.

Claims

1. A dual-frequency TEC inversion optimization method based on real-time altimetry data, comprising the following steps: S1: Set the daily data missing rate requirement, download the preliminary data, read and preprocess the multi-source data; S2: To address the potential mismatch between GNSS observation data and GNSS navigation data, this invention further designs and introduces a satellite data validity screening mechanism. It performs initial inversion calculations of ionospheric TEC based on dual-frequency GNSS observation data, and calculates the satellite's spatial position and clock bias information at a given observation epoch using conventional broadcast ephemeris calculation methods, based on the read navigation message data and observation time information. S3: After obtaining the satellite spatial position and receiver spatial coordinates corresponding to the observation epoch, to further describe the spatial geometric relationship between the satellite and the receiver, it is necessary to calculate the elevation and azimuth angles of available satellites relative to the receiver at each observation epoch. To obtain higher quality observation data, a satellite elevation angle cutoff threshold is set to remove low-quality observation data. After completing the calculation of the satellite elevation and azimuth angles relative to the receiver, to reduce the impact of multipath effects introduced by low-elevation satellites and tropospheric and ionospheric mapping errors on the ionospheric total electron content inversion results, this invention performs a screening process on the observation data based on the elevation angle threshold. S4: Cycle slip error detection and correction in carrier observation, and on this basis, smoothing and calibration processing of pseudorange STEC and carrier STEC; S5: To further improve the accuracy of ionospheric TEC inversion results, this invention uses OSB (Orbital Signal Bias) data to correct satellite-end errors in the smoothed results; S6: At this stage, the present invention introduces an automatic identification and elimination mechanism for STEC abnormal segments, and further implements integrity constraint checks on this basis; S7: Introduce observation data from the digital altimeter of the ionosphere, and use the key height parameter hmF2 (Height of maximum electron density in the F2 layer, i.e., peak height of the F2 layer) to adaptively correct the mapping function S(E) from STEC to VTEC. S8: The differential code bias at the receiver end is calculated using the DCB (Differential Code Bias) estimation method based on minimizing VTEC dispersion. S9: The final TEC result undergoes zero-baseline adjustment, which means that the minimum TEC value throughout the entire time period is used as a reference, and the TEC data is shifted as a whole. The shift amount is the absolute value of the current satellite's minimum value, thereby avoiding negative values ​​or inconsistencies in the baseline caused by residual system bias.

2. The dual-frequency TEC inversion optimization method based on real-time altimetry data according to claim 1, characterized in that, Step S1: In the embodiments of the present invention, the GNSS observation data and navigation data used are from Dongfang Station, the time resolution of the observation data is 1 second, and the time resolution of the navigation data is 20 minutes; the ionospheric digital altimeter data used are from Danzhou Station, and the time resolution is 15 minutes. The downloaded satellite code offset data includes relevant data from all satellites corresponding to the observed data. The final resolution is referenced to the highest resolution in the data, meaning the final TEC product has a time resolution of 1 second. The observed data follows the RINEX 3.03 (Receiver Independent Exchange Format (Version 3.03) data format. Matlab 2021 does not have a latest built-in function for reading observed data; therefore, this invention uses a custom-written function. Since the observed data is recorded in 5-minute single-file formats, a batch data reading method is required. Furthermore, this invention includes functions for monitoring various missing observation data and supplementing missing data. This invention is not limited to the aforementioned sites; data from other sites with similar or identical observation capabilities can also be applied to the method described herein.

3. The dual-frequency TEC inversion optimization method based on real-time altimetry data according to claim 1, characterized in that, In step S2, regarding the satellite data validity screening mechanism, specifically, during actual observation, some satellites may have complete or partial observation data in the observation files, but lack navigation message data within the corresponding time range. Since TEC's inversion calculation relies on precise satellite orbit and clock bias information, satellites with only observation data but lacking navigation data cannot participate in subsequent calculations and are considered invalid satellite data. Therefore, invalid satellite data needs to be removed at this stage. Furthermore, considering the signal system characteristics of the BeiDou satellite system and the completeness and stability of existing observation data, this invention selects channels 2 and 6 from multiple available observation channels as a dual-frequency combination, i.e., using C2I and C6I pseudorange observations and L2I and L6I carrier observations, replacing the traditional C / A code and P code combination method in GPS (Global Positioning System) for ionospheric parameter inversion calculation. The reason for choosing the above channel combination is that the signal frequency intervals corresponding to the second and sixth channels are relatively large, which can effectively enhance the sensitivity of the dual-frequency ionospheric delay differential to TEC; at the same time, the availability and temporal continuity of the observation data corresponding to these two types of channels are better, which is conducive to the unified correction of subsequent system errors. Where STECp and STECl are the inversion results of pseudorange observation and carrier observation, respectively; k is the conversion coefficient from dual-frequency difference to total electron content of oblique ionosphere; C6I, C2I, L2I, and L6I are the observed values ​​at each frequency point; f2 is the frequency of the second channel of Beidou satellite, which is taken as 1561.098 MHz here; f6 is the frequency of the sixth channel of Beidou satellite, which is taken as 1268.52 MHz here. For the RINEX 3.03 hybrid navigation message format used, the navigation data reading and parsing process was adapted to ensure that ephemeris parameters from different satellite systems can be used in calculations under a unified time reference. This guarantees the consistency and availability of the multi-system satellite position and clock bias calculation results. The satellite position and clock bias information obtained through this processing serves as the basic input data for subsequent ionospheric parameter calculations and spatial mapping analysis based on satellite-receiver geometry.

4. The dual-frequency TEC inversion optimization method based on real-time altimetry data according to claim 1, characterized in that, In step S3, the results provide necessary geometric constraints for subsequent ionospheric parameter mapping based on satellite-receiver geometry, low-quality observation removal, and spatial distribution analysis. In this invention, the satellite elevation angle cutoff threshold is set to 15°. By setting a predetermined satellite elevation angle cutoff threshold, when the satellite elevation angle is below this threshold, the reliability of the corresponding observation data is considered low; satellite observation data whose elevation angle meets the threshold condition are retained for subsequent processing.

5. The dual-frequency TEC inversion optimization method based on real-time altimetry data according to claim 1, characterized in that, In step S4, although carrier observation theoretically has high measurement accuracy, it is susceptible to cycle slips due to factors such as signal interruption, obstruction, and receiver state changes, resulting in discontinuities in the phase observation time series. To address this issue, this invention constructs a time continuity criterion to identify anomalous changes in carrier observation and corrects or removes identified cycle slip epochs, thereby ensuring the continuity of the STEC1 sequence in the time dimension. After completing carrier cycle slip processing, utilizing the consistency between pseudorange observation and carrier observation in ionospheric delay sensitivity, STECp is used as an absolute reference to smoothly calibrate STEC1, obtaining the STEC to be processed in the next step.

6. The dual-frequency TEC inversion optimization method based on real-time altimetry data according to claim 1, characterized in that, In step S5, specifically: the satellite absolute signal deviation mainly originates from the hardware delay differences of signals at different frequencies in the satellite transmission link. This deviation is relatively stable in time, but varies between different satellites and different frequency channels. If not corrected, it will directly affect the absolute accuracy of the dual-frequency combined ionospheric delay. In this step, based on the selected dual-frequency observation channel combination, the satellite absolute signal deviation product at the corresponding frequency is read and introduced into the STEC calculation results to systematically correct the deviation of STEC. Through this processing, the inter-frequency deviation introduced by the satellite-end hardware characteristics can be effectively eliminated, ensuring that the STEC results maintain a consistent reference benchmark across different satellites.

7. The dual-frequency TEC inversion optimization method based on real-time altimetry data according to claim 1, characterized in that, In step S6, after completing the absolute signal bias correction for the satellites, considering the influence of factors such as observation noise, signal obstruction, multipath effects, and incomplete local correction errors, some satellite STEC time series may still contain anomalous data segments that significantly deviate from the actual characteristics of ionospheric changes. To address this issue, it is necessary to automatically identify and process anomalous segments in the STEC data that has undergone systematic bias correction. This process constructs a statistical reference quantity based on multi-satellite observation results and analyzes the STEC time series of individual satellites segment by segment. By comparing the degree of deviation between STEC and the reference statistic within a continuous time period, it is determined whether there are significant anomalies in that time period. When anomalous samples occupy a major proportion within a certain continuous time period and the deviation exceeds a preset threshold, the entire time period is determined to be an anomalous segment, and the corresponding STEC data is removed, thereby eliminating the abrupt changes or drift effects caused by non-ionospheric physical changes. After removing the aforementioned outlier segments, this invention further performs an integrity check on the number of valid observations of each satellite's STEC data within the critical time window. Specifically, the number of valid STEC observation points for each satellite within a specified time range is counted. When the number of valid observations is lower than a preset threshold, the data continuity and reliability of that satellite within that time period are considered insufficient, making it unsuitable for subsequent ionospheric parameter inversion calculations. The STEC data of the corresponding satellite within that time range are then completely removed.

8. The dual-frequency TEC inversion optimization method based on real-time altimetry data according to claim 1, characterized in that, In step S7, after completing the quality control and outlier data removal for STEC, it is necessary to convert STEC to VTEC to characterize the true electron density distribution characteristics of the ionosphere above the observation area. Traditional STEC-to-VTEC mapping methods typically employ a single-layer thin-shell model at a fixed height and assume a constant ionospheric height. This assumption is difficult to satisfy in actual ionospheric environments, especially when ionospheric activity is enhanced or regional differences are significant, easily introducing systematic mapping errors. To address the aforementioned issues, this invention introduces observation data from a digital altimeter for the ionosphere and utilizes the key ionospheric height parameter hmF2 to adaptively correct the mapping function from STEC to VTEC. By introducing a time-varying effective ionospheric height parameter, the mapping process can reflect the true spatial distribution characteristics of the ionospheric main electron density layer, thereby improving the physical rationality and accuracy of the VTEC inversion results. In the specific implementation process, based on the Earth's average radius and satellite elevation information, and combined with the effective ionospheric height hmF2 provided by ionospheric altimetry data, an improved oblique factor mapping function is constructed. This mapping function describes the geometric relationship between the oblique propagation path and the vertical direction, and its expression is as follows: Where S(E) is the mapping factor to VTEC obliquely; R is the Earth's average radius; h is the peak height of the F2 layer; and E is the satellite elevation angle. VTEC is then obtained using the following formula: By incorporating time-varying ionospheric height information, this mapping method can effectively reduce the errors caused by the fixed height assumption, improve the stability of VTEC mapping under low elevation angle conditions, and enhance the responsiveness of the mapping results to the actual changes in the ionosphere. This improved mapping method based on ionospheric altimetry data enables the organic integration of the geometric model and the actual ionospheric structure during the conversion from STEC to VTEC, providing a more reliable technical means for the generation of high-precision ionospheric TEC products and their spatiotemporal evolution analysis.

9. The dual-frequency TEC inversion optimization method based on real-time altimetry data according to claim 1, characterized in that, In step S8, after completing satellite-side systematic error correction, STEC anomaly removal, and mapping of STEC to VTEC based on ionospheric altimetry data, the resulting TEC data still contains systematic biases introduced by receiver hardware delays. This bias manifests as a fixed delay difference between different frequency channels, commonly referred to as receiver differential code bias. If not estimated and corrected, it will directly affect the absolute accuracy of the TEC and the comparability between different observation periods. In this step, a receiver DCB estimation method based on minimizing VTEC dispersion is used to calculate the differential code deviation at the receiver end. The basic idea of ​​this method is: within a reasonable receiver deviation search range, different candidate deviation values ​​are introduced to correct the STEC data, and the corrected STEC is mapped to VTEC; when the receiver deviation value is reasonable, the VTEC mapped from STEC at the same time should have good consistency, and its overall dispersion should be minimized. In the specific implementation process, the initial search interval and step size for receiver bias are first set, and the candidate bias values ​​are iterated one by one. For each candidate bias value, it is subtracted from the STEC data after removing satellite bias, and the corrected STEC is converted into VTEC using the aforementioned improved oblique mapping factor. Subsequently, the statistical dispersion index is calculated for the VTEC results corresponding to multiple satellites within the same epoch, and accumulated in the time dimension to characterize the overall consistency level of VTEC under the candidate bias value. By comparing the VTEC dispersion corresponding to different receiver biases, the bias value that minimizes the dispersion is selected as the receiver differential code bias estimation result for that observation station. To improve estimation accuracy, this invention further employs a stepwise narrowing of the search interval based on the initial search, performing multiple rounds of refined estimation of the receiver bias to obtain a stable and reliable receiver bias value. After completing the receiver differential code bias estimation, the estimated receiver bias is uniformly incorporated into the STEC data for correction, resulting in the absolute STEC after removing systematic biases between the satellite end and the receiver end. Simultaneously, based on the improved mapping function, the corrected STEC is converted into VTEC, obtaining the final high-precision VTEC mapping result.

10. The dual-frequency TEC inversion optimization method based on real-time altimetry data according to claim 1, characterized in that, In step S9, after completing the receiver differential code bias estimation and correction, the STEC and its corresponding VTEC inversion results are obtained after removing the systematic biases between the satellite end and the receiver end. To unify the reference benchmark of the TEC products and ensure their physical rationality, this invention performs zero-reference adjustment processing on the final TEC results. That is, using the minimum value of the TEC over the entire time period as a reference, the TEC data is shifted as a whole to make its minimum value zero, thereby avoiding negative values ​​or benchmark inconsistencies caused by residual systematic biases. After the above zero-baseline adjustment, the final STEC and VTEC products are obtained respectively. The TEC products remain consistent with the original high-resolution observation data in the time dimension and can be used to finely characterize the rapid changes in the total electron content of the ionosphere. Building upon this, the present invention further calculates the ionospheric scintillation characterization index ROTI based on STEC data obtained through carrier inversion. ROTI reflects the fluctuation intensity of TEC over short timescales and is an important parameter for characterizing ionospheric irregularity activity and space weather disturbances. In the ROTI calculation process, firstly, linear interpolation is performed on the short-time data discontinuities in the STEC time series to ensure the continuity of the time series; then, STEC is downsampled according to a preset time window, and the rate of change of TEC between adjacent epochs is calculated; finally, the standard deviation of the rate of change series is calculated within the sliding time window to obtain the ROTI value of the corresponding epoch.

Images