Global ionosphere TEC inversion method based on Galileo HAS service

Abstract

The invention relates to a global ionosphere TEC inversion method based on a Galileo HAS service, and the method comprises the steps: obtaining correction information through the HAS service, carrying out the preprocessing of observation data of a reference station according to the correction information, and obtaining a preprocessed carrier phase observation value; calculating carrier phase ambiguity and a purified carrier phase observation value according to the preprocessed carrier phase observation value; calculating the total electron content STEC of an inclined path according to the carrier phase ambiguity and the purified carrier phase observation value; constructing a double-layer ionosphere network model according to the geographic coordinates of the ionosphere puncture point and STEC; and carrying out inversion on the vertical total electron content of the ionized layer by using Kalman filtering to obtain the vertical total electron content of the ionized layer and carrying out compression. The method breaks through the network dependence of traditional foundation enhancement, and can still provide continuous and reliable ionized layer delay correction under the scenes of ocean, desert and the like.

Description

Technical Field

[0001] The present invention relates to the field of satellite navigation enhancement technology, and in particular to a global ionospheric TEC inversion method based on the Galileo HAS service. Background Art

[0002] Ionospheric delay is one of the main sources of error in high-precision Global Navigation Satellite System (GNSS) applications. Its temporal and spatial variations make real-time, accurate correction a technical challenge. Traditional ionospheric total electron content (TEC) inversion methods rely heavily on precise ephemeris data from Ground-Based Augmentation Systems (GBAS) or post-process internet access. However, existing technologies struggle to meet the real-time, global, and robust requirements in areas without network coverage, such as the polar regions and open oceans, or in emergency scenarios such as disaster-induced communication outages.

[0003] In recent years, the launch of the Galileo system's High Accuracy Service (HAS) has provided new ideas for global real-time ionospheric monitoring. HAS directly broadcasts the precise orbit and clock corrections of GPS / Galileo satellites through the E6-B signal, which theoretically can get rid of the dependence on ground communication networks. Current research focuses on the performance verification of HAS in the field of precise positioning, while existing ionospheric modeling methods face two major technical bottlenecks: First, the error of regional ionospheric models increases significantly at low latitudes or during magnetic storms, and the TEC prediction deviation in the equatorial region can reach 20TECU (Total Electron Content Unit, 1TECU=10 16 electrons / m 2 ) and above; secondly, while the Global Ionospheric Grid has wide coverage, it suffers from insufficient spatiotemporal resolution and poor timeliness. Although international GNSS services provide high-frequency TEC data in real-time, they still rely on internet transmission and cannot serve scenarios with limited communication.

[0004] The core technical architecture of the HAS service relies on the Galileo system's high-stability space segment and high-precision ground segment. The space segment utilizes Medium Earth Orbit (MEO) satellites capable of broadcasting in the E6 frequency band, creating a global enhanced signal broadcast platform. The ground segment comprises a closed-loop system consisting of multi-frequency, multi-mode monitoring stations, high-precision atomic clocks, and a data processing center. This system utilizes multi-frequency observation data combined with a Kalman filter algorithm to generate high-precision corrections, such as precise satellite orbits and clock errors, in real time. These corrections are broadcast with low latency via the E6-B channel, supported by a data assurance mechanism to ensure reliable transmission.

[0005] Galileo HAS, with its high-precision orbit and clock correction capabilities, has broad application prospects in many fields such as drone navigation, surveying and mapping, and disaster monitoring. Current ionospheric delay correction methods mainly rely on ground-based augmentation networks or post-precision ephemeris, which makes it difficult to cover areas without networks such as oceans and deserts. In addition, the global ionospheric model has insufficient temporal and spatial resolution and cannot adapt to the dynamic changes of the ionosphere in real time. Especially under complex disturbances such as magnetic storms or equatorial anomalies, the error increases significantly, resulting in a decrease in the stability of high-precision positioning. Although the existing Galileo HAS service provides satellite-based precision corrections, it lacks the ability to adapt to ionospheric disturbances in real time, which limits its application in dynamic scenarios such as drone navigation and autonomous driving. Summary of the Invention

[0006] The purpose of this invention is to provide a global ionospheric TEC inversion method based on the Galileo HAS service, breaking through the network dependence of traditional ground-based augmentation, and still providing continuous and reliable ionospheric delay correction in scenarios such as oceans and deserts, providing technical support for high-precision applications such as drone navigation and autonomous driving.

[0007] To achieve the above object, the present invention provides the following solutions:

[0008] A global ionospheric TEC inversion method based on the Galileo HAS service, including:

[0009] Acquire correction information through the HAS service, preprocess the observation data of the reference station according to the correction information, and obtain preprocessed carrier phase observation values;

[0010] Calculating carrier phase ambiguity and cleaned carrier phase observation values ​​based on the preprocessed carrier phase observation values;

[0011] Calculating the slant-path ionospheric total electron content (STEC) based on the carrier phase ambiguity and the cleaned carrier phase observation value;

[0012] A double-layer ionospheric network model is constructed based on the geographical coordinates of the ionospheric puncture point and the ionospheric total electron content (STEC) of the slant path.

[0013] According to the double-layer ionosphere network model, the vertical total electron content of the ionosphere is inverted using Kalman filtering to obtain the vertical total electron content of the ionosphere and compress it.

[0014] Optionally, preprocessing the observation data of the reference station according to the correction information includes:

[0015] Calculating the orbit correction and clock correction at the current moment based on the correction information;

[0016] Calculating the real-time precise clock error according to the clock error correction number;

[0017] The observation data of the reference station is preprocessed according to the orbit correction number and the real-time precise clock error.

[0018] Optionally, preprocessing the observation data of the reference station according to the orbit correction number and the real-time precise clock error includes:

[0019]

[0020] Among them, Φ' is the carrier phase observation value after preprocessing, Φ is the carrier phase observation value, f is the carrier frequency, is the real-time precise clock error after HAS correction, C is the speed of light, Δr ECEF is the orbit correction number under ECEF, and u is the unit vector from the satellite to the receiver.

[0021] Optionally, computing carrier phase ambiguities and cleaned carrier phase observations includes:

[0022] Calculating double-difference observations and ionospheric free combination observations based on the preprocessed carrier phase observations;

[0023] Using the double-difference observations and the ionospheric free combination observations, fixed wide-lane ambiguities and floating-point ambiguities are obtained respectively;

[0024] According to the fixed wide lane ambiguity and the floating point ambiguity, the geometric free ambiguity is calculated and resolved to obtain the carrier phase ambiguity and the purified carrier phase observation value.

[0025] Optionally, computing double-difference observations includes:

[0026]

[0027] in, is the double-difference observation, is the wide-lane ambiguity, Φ'1 and Φ'2 represent the carrier phase observations of the two frequency bands that have been corrected by the orbit and clock corrections of the HAS; j and k represent the inter-station and inter-satellite differences, respectively; A and B represent two different receivers or satellites, and f1 and f2 are the carrier frequencies of the two different frequency bands of the Galileo satellite navigation system.

[0028] Optionally, building a two-layer ionospheric network model includes:

[0029] Processing the slant path ionospheric total electron content (STEC), calculating a regional STEC spatial gradient, stratifying the ionosphere according to the regional STEC spatial gradient, and dividing the network resolution;

[0030] Based on the geographical coordinates of the ionospheric puncture point and STEC, a linear relationship between STEC and the vertical total electron content of the grid point corresponding to the ionospheric puncture point is established;

[0031] According to the linear relationship, the vertical total electron content is filled into the corresponding grid points as the initial value to construct the double-layer ionosphere network model.

[0032] Optionally, calculation of regional STEC spatial gradients includes:

[0033]

[0034] in, is the regional STEC spatial gradient, To correct the magnetic latitude, λ LT For local time.

[0035] Optionally, the inversion of the ionospheric vertical total electron content based on Kalman filtering includes:

[0036] S1. Constructing an initialization state vector according to VTEC in the double-layer ionospheric network model;

[0037] S2, perform time update and predict the prior state vector and prior covariance matrix at the current moment;

[0038] S3, based on the current moment’s prior state vector and prior covariance matrix, combined with the STEC observation value, constructs the observation equation and calculates the Kalman gain;

[0039] S4. Update the state estimate and the posterior covariance matrix according to the Kalman gain, output the inverted real-time VTEC value, add one at each moment, and return to S2.

[0040] Optionally, before obtaining the ionospheric vertical total electron content for compression, the method includes: extracting the VTEC error standard deviation of each grid point from the covariance matrix of the Kalman filter, and performing error evaluation on the inversion result.

[0041] Optionally, obtain the ionospheric vertical total electron content for compression including:

[0042] L total =L bottom +L top +L GIVEI ≤26·B

[0043] Among them, L total is the total message length, L bottom is the length of the compressed underlying grid data, L top is the top grid data length, L GIVEI is the length of GIVEI data, and B is the number of bytes that each page of telegram can accommodate.

[0044] The beneficial effects of the present invention are: (1) Satellite-based broadcasting and global coverage: No ground base stations or the Internet are required, and global seamless coverage is achieved by relying on the Galileo MEO satellite constellation, solving the coverage blind spot problem of traditional RTK / SBAS in polar regions, oceans and other regions.

[0045] (2) High Precision and Real-Time: High-precision inversion of ionospheric TEC is achieved by utilizing HAS centimeter-level orbits, sub-nanosecond clock errors, and real-time Kalman filtering. A dynamic update mechanism ensures that the inversion results can respond to instantaneous changes in the ionosphere in real time, meeting the real-time correction requirements in dynamic scenarios.

[0046] (3) Dynamic construction and optimization of the two-layer ionospheric grid model: Based on the real-time STEC spatial gradient and MODIP-LT coordinate system, the ionospheric layer height and grid resolution are dynamically adjusted. By dynamically adapting the layer height, the vertical heights of the bottom and top layers are adjusted in real time to solve the problem of ionospheric structure misalignment at low latitudes or during magnetic storms. By setting the grid resolution differently, the grid resolution is dynamically refined to improve the ability to capture changes in ionospheric details. Finally, the dynamic puncture point correction is performed to ensure that it matches the actual electron density distribution and reduce mapping errors.

[0047] (4) Anti-interference and robustness: Through real-time integrity monitoring (GIVEI) and dynamic error assessment mechanisms, we can effectively cope with complex ionospheric disturbances and ensure the stability and reliability of the inversion results. Dynamic adjustment of error weights in sparse areas and extreme space weather conditions further improves the adaptability of the model. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0049] Figure 1 Schematic diagram of a flow chart of a global ionospheric TEC inversion method based on the Galileo HAS service according to an embodiment of the present invention;

[0050] Figure 2 A flowchart for constructing a double-layer ionospheric grid model according to an embodiment of the present invention;

[0051] Figure 3 This is a Kalman filter VTEC inversion flow chart of an embodiment of the present invention. DETAILED DESCRIPTION

[0052] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0053] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0054] A global ionospheric TEC inversion method based on the Galileo HAS service, including:

[0055] Correction information is obtained through the HAS service, and the observation data of the reference station is preprocessed according to the correction information to obtain the preprocessed carrier phase observation value;

[0056] Calculating carrier phase ambiguity and cleaned carrier phase observation values ​​based on the preprocessed carrier phase observation values;

[0057] The slant path ionospheric total electron content (STEC) is calculated based on the carrier phase ambiguity and the cleaned carrier phase observations.

[0058] A double-layer ionospheric network model is constructed based on the geographical coordinates of the ionospheric puncture point and the ionospheric total electron content (STEC) of the slant path.

[0059] According to the double-layer ionosphere network model, the vertical total electron content of the ionosphere is inverted using Kalman filtering to obtain the vertical total electron content of the ionosphere and compress it.

[0060] Furthermore, preprocessing of the observation data of the reference station according to the correction information includes:

[0061] Calculate the orbit correction and clock correction at the current moment based on the correction information;

[0062] Calculate the real-time precise clock error based on the clock correction number;

[0063] The observation data of the reference station are preprocessed according to the orbit correction number and real-time precise clock error.

[0064] Furthermore, the observation data of the reference station is preprocessed according to the orbit correction number and the real-time precise clock error, including:

[0065]

[0066] Among them, Φ' is the carrier phase observation value after preprocessing, Φ is the carrier phase observation value, f is the carrier frequency, is the real-time precise clock error after HAS correction, C is the speed of light, Δr ECEF is the orbit correction number under ECEF, and u is the unit vector from the satellite to the receiver.

[0067] Furthermore, the calculation of the carrier phase ambiguity and the cleaned carrier phase observation value includes:

[0068] Calculate double difference observations and ionospheric free combination observations based on preprocessed carrier phase observations;

[0069] Double-difference observations and ionospheric free combination observations are used to obtain fixed wide-lane ambiguities and floating-point ambiguities respectively.

[0070] According to the fixed wide-lane ambiguity and floating-point ambiguity, the geometric free ambiguity is calculated and resolved to obtain the carrier phase ambiguity and the purified carrier phase observation value.

[0071] Furthermore, calculating the double-difference observations includes:

[0072]

[0073] in, is the double-difference observation, is the wide-lane ambiguity, Φ'1 and Φ'2 represent the carrier phase observations of the two frequency bands that have been corrected by the orbit and clock corrections of the HAS; j and k represent the inter-station and inter-satellite differences, respectively; A and B represent two different receivers or satellites, and f1 and f2 are the carrier frequencies of the two different frequency bands of the Galileo satellite navigation system.

[0074] Furthermore, constructing a double-layer ionospheric network model includes:

[0075] The total electron content (STEC) of the ionospheric slant path is processed to calculate the regional STEC spatial gradient. The ionosphere is layered according to the regional STEC spatial gradient and the network resolution is divided.

[0076] Based on the geographical coordinates of the ionospheric puncture point and STEC, a linear relationship between STEC and the vertical total electron content of the grid point corresponding to the ionospheric puncture point is established;

[0077] According to the linear relationship, the vertical total electron content is filled into the corresponding grid points as the initial value to construct a double-layer ionospheric network model.

[0078] Furthermore, the calculation of regional STEC spatial gradients includes:

[0079]

[0080] in, is the regional STEC spatial gradient, To correct the magnetic latitude, λ LT For local time.

[0081] Furthermore, the inversion of the vertical total electron content of the ionosphere based on the Kalman filter includes:

[0082] S1. Construct an initialization state vector based on VTEC in the double-layer ionospheric network model;

[0083] S2, perform time update and predict the prior state vector and prior covariance matrix at the current moment;

[0084] S3, based on the current moment’s prior state vector and prior covariance matrix, combined with the STEC observation value, constructs the observation equation and calculates the Kalman gain;

[0085] S4. Update the state estimate and the posterior covariance matrix according to the Kalman gain, output the inverted real-time VTEC value, add one at each moment, and return to S2.

[0086] Furthermore, before obtaining the ionospheric vertical total electron content for compression, the method includes extracting the VTEC error standard deviation of each grid point from the covariance matrix of the Kalman filter and performing error evaluation on the inversion result.

[0087] Furthermore, obtaining the vertical total electron content of the ionosphere for compression includes:

[0088] L total =L bottom +L top +L GIVEI ≤26·B

[0089] Among them, Ltotal is the total message length, L bottom is the length of the compressed underlying grid data, L top is the top grid data length, L GIVEI The length of the GIVEI data.

[0090] The present embodiment will be further described below with reference to the accompanying drawings:

[0091] like Figure 1 A global ionospheric TEC inversion method based on the Galileo HAS service of this embodiment is shown, and the method includes the following steps:

[0092] Step 1: Real-time acquisition and preprocessing of HAS correction information;

[0093] Step 2: Carrier phase ambiguity layer resolution and data purification;

[0094] Step 3: Calculate the ionospheric IPP and STEC in the MODIP-LT coordinate system;

[0095] Step 4: Construction and initialization of the double-layer ionosphere grid model;

[0096] Step 5: Real-time inversion of VTEC based on Kalman filtering;

[0097] Step 6: Error evaluation of inversion results and generation of GIVEI;

[0098] Step 7: VTEC grid data differential compression encoding;

[0099] Step 8: User-side VTEC interpolation.

[0100] Specific steps: 1. Use the precise orbit and clock correction information provided by HAS to pre-process the reference station observation data to eliminate the influence of satellite orbit error and clock error on the observation value, and ensure the benchmark accuracy of subsequent ambiguity resolution and TEC inversion. First, the receiver receives the HAS correction information broadcast by the Galileo satellite via the E6-B band, which contains the orbit correction (Δ r Radial, Δ a Tangential, Δ c Normal) and corresponding velocity correction Information such as the clock correction polynomial coefficients (a0, a1, a2) and the reference time (t0).

[0101] (1) In order to convert the real-time orbit correction information of HAS into a precise orbit suitable for TEC inversion, the following steps are performed:

[0102] Since the orbit correction number is based on the satellite body fixed system (SBF), and the reference station observation data is collected in the Earth-Centered Earth-Fixed (ECEF) coordinate system, first use the formula Calculate the orbit correction number in the satellite-fixed system at the current time t, where is the orbit correction vector at reference time t0, is the corresponding velocity correction vector. Then through the rotation matrix:

[0103]

[0104] where r broadcast and The position and velocity vector of the satellite in the ECEF coordinate system calculated by the broadcast ephemeris are The orbit correction number in the SBF coordinate system is converted to ECEF. The real-time orbit correction number after the coordinate system conversion is superimposed on the satellite orbit calculated by the broadcast ephemeris to obtain the real-time precise orbit r of the satellite at time t. precise =r broadcast +Δr ECEF .

[0105] (2) Calculate the real-time clock correction value at the current time t from the clock correction information: For clock correction Follow the steps below to get the real-time precise clock error: In the formula and They are respectively the satellite clock error calculated by broadcast ephemeris and the real-time precise clock error after HAS correction, C is the speed of light, δC s The relativistic effect can be calculated based on the position and speed of the satellite. The calculation method is:

[0106] Since the orbit error calculated by traditional broadcast ephemeris is about decimeter level and the clock error can reach meter level, while after HAS correction, they can reach centimeter level and sub-nanosecond level respectively, the real-time satellite orbit error Δr obtained by HAS data preprocessing is ECEF and satellite clock error Correct the observations to eliminate satellite errors.

[0107] For the carrier phase observation value Φ, the preprocessed carrier phase observation value Φ' can be expressed as:

[0108]

[0109] Where f is the carrier frequency, u is the unit vector from the satellite to the receiver, and its direction points to the satellite.

[0110] Step 2: Because traditional PPP-RTK or RTK ambiguity resolution also relies on a ground reference station network or post-precision ephemeris, it cannot be implemented in areas without network coverage, such as oceans and deserts. Moreover, a single ambiguity resolution method is prone to failure when the ionosphere is active (such as in low latitudes or during magnetic storms). Therefore, this embodiment adopts a layered and progressive ambiguity resolution strategy to gradually improve accuracy. The specific steps are as follows:

[0111] (1) Using the long wavelength characteristic of the wide-lane (HMW) combined observation value, the initial ambiguity is quickly fixed, providing a basis for subsequent accurate solution. Based on the preprocessed carrier phase observation value Φ', the wide-lane combined observation value L is calculated using dual frequency. WL for:

[0112]

[0113] Among them, f1 and f2 are the carrier frequencies of two different frequency bands of the Galileo satellite navigation system; Φ'1 and Φ'2 are the preprocessed carrier phase observation values ​​of the receiver f1 and f2 frequency bands respectively.

[0114] The double difference technique is used to eliminate the common error and the wide lane ambiguity N is reduced within a certain period of time. WL Fixed to an integer, for the wide lane combination observation value L WL Inter-station and inter-satellite differences are performed to significantly improve the success rate and accuracy of wide-lane ambiguity fixation. Since residual orbit errors (especially the orbit error of the broadcast ephemeris) still affect ambiguity fixation, the precise orbit correction number Δr pre-processed by the HAS correction information is directly used when constructing the double-difference observation equation. ECEF and clock correction number Replacing the broadcast ephemeris parameters makes the double difference residual smaller, improves the success rate and accuracy of wide lane ambiguity fixation, and thus more accurately reflects information such as ionospheric delay. Expressed as:

[0115]

[0116] By constantly adjusting The integer property is satisfied, and the wide-lane ambiguity is finally fixed. Where Φ'1 and Φ'2 represent the carrier phase observations of the two frequency bands that have been corrected by the orbit and clock corrections of the HAS; j and k represent the inter-station and inter-satellite differences, respectively; and A and B represent two different receivers or satellites.

[0117] (2) After completing the wide-lane ambiguity fixation, in order to eliminate the influence of the ionospheric first-order delay on the ambiguity resolution, it is necessary to construct an observation value combination without ionospheric error. The integer characteristics of are fixed, which can constrain L IF First, calculate the ionospheric free combination observation value L IF , for dual-frequency signals, its expression is: Where Φ'1 and Φ'2 represent the carrier phase observations of the two frequency bands that have been corrected by the orbit and clock corrections of the HAS; r and s represent the receiver and satellite identifiers, respectively, which can be represented by different receivers A and B or satellites j and k.

[0118] Although the wavelength of this combination is short, it completely eliminates the first-order ionospheric delay, laying the foundation for high-precision floating-point ambiguity resolution. rec , use the PPP-RTK model to solve the floating point ambiguity N IF , estimate the parameter dt by the least squares method rec and The observation equation is: where ρ = || r precise -r rec ||.

[0119] ρ r,s : The geometric distance from the receiver to the satellite;

[0120] Satellite precise clock corrections provided by HAS;

[0121] λ IF : wavelength of the ionospheric-free combination;

[0122] dt rec : receiver clock error;

[0123] ε: observation noise

[0124] The floating point blur to resolve.

[0125] This step is provided by HAS correction The broadcast clock error is replaced by the satellite clock error, significantly reducing the satellite clock error residual in the equation. This makes the receiver clock error estimate closer to the true value, thereby improving the accuracy of floating-point ambiguity resolution.

[0126] (3) In order to eliminate the influence of frequency difference, N IF and a fixed N WL , calculate the unambiguous geometrically free ambiguity N GF , whose expression is:

[0127]

[0128] After the ambiguity resolution is completed, the fixed carrier phase ambiguities N1 and N2 of the frequency bands f1 and f2 and the purified carrier phase observation values ​​Φ1” and Φ2” are output for use in calculating the slant path total electron content (STEC) in step 4.

[0129]

[0130]

[0131] Where, is the carrier wavelength.

[0132] Step 3: Calculate the geographical coordinates of the ionospheric piercing point (IPP) (λ IPP , ), for a satellite and reference station observation pair, the IPP The calculation formula is:

[0133]

[0134] Among them, R E is the radius of the Earth, θ is the satellite zenith distance, and h is the ionospheric grid height (in step 4, the bottom layer h1 or the top layer h2).

[0135] IPP Lambda IPP It can be calculated based on information such as satellite azimuth and reference station longitude.

[0136] Since the MODIP-LT (Modified Dip Latitude-Local Time) coordinate system more accurately characterizes the magnetic field and solar radiation dependence of the ionospheric electron density by correcting the magnetic dip latitude and local time conversion, it is significantly better than the geographic coordinate system in the equator and polar regions, and can significantly reduce modeling errors. After converting the IPP geographic coordinates to MODIP-LT coordinates, the STEC is mapped to the corresponding grid points. Perform the following conversion

[0137]

[0138] Where Δλ LT When making local amendments; For the geomagnetic latitude of IPP, the International Geomagnetic Reference Field (IGRF) model needs to be called t is the date;

[0139] Based on the fixed single-frequency carrier phase ambiguities N1 and N2 in step 2 and the cleaned carrier phase observations Φ1″ and Φ2″, the slant path ionospheric total electron content (STEC) is calculated using the following formula:

[0140] Step 4: Construct a two-layer ionospheric grid model, dynamically adjusting the layer height based on the real-time STEC spatial gradient and the MODIP-LT coordinate system. This avoids the misalignment between the IPP and the actual electron density concentration area caused by fixed altitude, which can introduce mapping errors. By breaking through the limitations of fixed altitude stratification, the two-layer grid height is dynamically adapted to the ionospheric activity level (STEC gradient) and spatial distribution (MODIP partitioning), resolving the misalignment between the puncture point and the actual ionospheric structure in low-latitude or magnetic storm scenarios, and improving the model's coverage of global ionospheric differences.

[0141] like Figure 2 As shown in the figure, the STEC observations of different satellites at the same time are first processed in the MODIP-LT coordinate system to extract the regional STEC spatial gradient in real time. To correct the magnetic latitude, λ LT When it is local time, it is obtained by converting the coordinates of the ionospheric piercing point in step 3.

[0142] Layered highly dynamic adaptation mechanism:

[0143] Real-time STEC spatial gradients The coordinated judgment with the Modified Magnetic Latitude (MODIP) realizes the intelligent adjustment of the ionospheric vertical layer height:

[0144] (1) Dynamic formula adjustment of active areas:

[0145] when When the ionosphere is determined to be active (such as the equatorial anomaly belt, magnetic storm events), the layer height is dynamically calculated according to the partition characteristics:

[0146] Low magnetic latitude area (MODIP<30°): bottom layer height Top floor height Capture dramatic changes in electron density by extending vertical coverage;

[0147] In the mid-magnetic latitudes (30°≤MODIP≤60°): the default altitude will be temporarily adjusted by ±20km to accommodate sudden disturbances.

[0148] High magnetic latitudes (MODIP>60°): Response to ionospheric compression caused by geomagnetic activity.

[0149] Transition zone resolution priority optimization:

[0150] when When the disturbance is moderate, the default layer height is maintained (low latitude h1 = 260 km / h2 = 1700 km, mid-latitude h1 = 270 km / h2 = 1600 km, high latitude h1 = 280 km / h2 = 1500 km), and only the horizontal grid resolution is encrypted to improve the ability to capture details.

[0151] Calm zone fixed parameters:

[0152] when When using the preset fixed layer height, it reflects the stable distribution characteristics of electron density (such as the quiet ionosphere at mid-latitude night).

[0153] (2) Rules for setting grid resolution differentiation:

[0154] Based on MODIP partition and Dual gradient criteria, implementing three-level resolution dynamic control:

[0155] Base resolution:

[0156] The default grid is 0.5°×0.5°, which is suitable for mid-magnetic latitudes and The peaceful scene.

[0157] Gradient-triggered encryption:

[0158] Low magnetic latitude area: default resolution 0.3°×0.3°, when The density is increased to 0.2°×0.2°. The transition is encrypted to 0.25°×0.25°;

[0159] Medium magnetic latitude: Default 0.5°×0.5°, only When encrypted to 0.3°×0.3°;

[0160] High magnetic latitude area: default 0.4°×0.4°, When the gradient is medium, it is encrypted to 0.25°×0.25°, and when the gradient is medium, it is 0.3°×0.3°.

[0161] Special event adaptation:

[0162] During magnetic storms, a 0.3°×0.3° encrypted grid is forcibly used in mid-latitude areas, and during geomagnetic substorms, the high-latitude areas switch to a 0.25°×0.25° mode to ensure modeling accuracy under extreme events.

[0163] Dynamic puncture point (IPP) coordinate correction:

[0164] Recalculate the geographical coordinates of the ionospheric puncture point based on the adjusted layer heights h1' and h'2: where h' iThe dynamically adjusted bottom or top altitude ensures that the IPP coordinates match the actual ionospheric vertical structure.

[0165] The relationship between STEC and the vertical total electron content (VTEC) of the grid point where the corresponding IPP is located is established, and VTEC is filled into the corresponding grid point as the initial value. Considering the projection function mf(Z), we have: STEC = VTEC init mf(θ);

[0166] The calculation formula of mf(θ) is: θ is the zenith distance at IPP; e is the altitude angle of the satellite, h2 and h1 are the heights of the upper and lower layers of the double-layer ionosphere model, respectively.

[0167] Step 5. After completing the construction of the initial two-layer ionosphere network model, the VTEC of each grid point is dynamically updated based on the Kalman filter. Through the "prediction-correction" iterative mechanism, the static error of the spatial interpolation bias of the initial grid model of the ionosphere model and the dynamic error of the rapid change of electron density caused by magnetic storms, equatorial anomalies, etc. are eliminated. Based on the new STEC observation values, the VTEC grid data is continuously optimized to ensure that the inversion results keep up with the instantaneous state of the ionosphere. Figure 3 The specific steps are as follows:

[0168] According to the VTEC in the double-layer ionospheric network model, the state equation is constructed with the initial VTEC value VTEC of each grid point in the double-layer grid model. init As the core, let’s initialize the state vector:

[0169]

[0170] Among them, VTEC i,k is the VTEC of the i-th grid point at time k init Value, δ r,0 It is the initial parameters of other states such as satellite / receiver phase deviation.

[0171] At the same time, the covariance matrix P0 is initialized as a diagonal matrix, the diagonal elements correspond to the initial variances of VTEC and each deviation parameter, and the variances of the remaining parameters are set according to the noise characteristics.

[0172] The state equation can be expressed as:

[0173] Among them, Φ k,k-1 is the state transfer matrix (describing the change of VTEC from time k-1 to time k), is the process noise (obeying Gaussian distribution N(0,Q k-1 ), reflecting the error not captured by the model).

[0174] The static initial VTEC of the double-layer grid model is converted into a mathematical model that can describe the dynamic changes over time, providing a basic framework for subsequent "prediction-correction". At the same time, interference parameters such as phase deviation are incorporated to reduce their impact on VTEC inversion.

[0175] (2) Predict the current state vector based on the state equation and propagate the covariance matrix:

[0176] Based on the state vector of the previous moment k-1 and the covariance matrix P k-1 , predict the prior state vector of the current moment k and the prior covariance matrix

[0177]

[0178] Among them, Q k-1 is the process noise covariance matrix. In the absence of new observation data, the current VTEC and uncertainty (covariance) are estimated based on historical states and time-varying patterns, providing a “priori reference” for subsequent corrections using observations, ensuring dynamic and real-time performance.

[0179] (3) According to the observation equation, calculate the Kalman gain observation update and Kalman gain calculation:

[0180] With the prior state of (2) and the prior covariance Based on the observation equation, the STEC observations were combined to construct Calculate the Kalman gain K k :

[0181]

[0182] in, is the observation vector, which contains the STEC observation value at the corresponding time calculated in step 3; it is mapped to the VTEC grid through the projection function mf(Z), H k is the observation matrix, which connects the state vector with the observation value; is the observation noise, which obeys the Gaussian distribution N(0,R k );R k is the observation noise covariance matrix.

[0183] The Kalman gain is a weighting factor that balances prediction error (prior covariance) with observation error (observation noise). A larger value indicates a stronger influence of the observation on the correction. This step links the static grid model's predictions with the dynamic STEC observations, providing a quantitative weight for subsequent corrections.

[0184] (4) Update the state estimate and the posterior covariance matrix according to the Kalman gain, output the inverted real-time VTEC value, and add one at each time to return to (2):

[0185] Use the Kalman gain K k Correct the prior state to obtain the posterior state (optimal VTEC)

[0186] Update the posterior covariance matrix (reflecting the uncertainty of the corrected VTEC) P k :

[0187] Output the current real-time VTEC value, add 1 to the time k, and return to step S2 to enter the next iteration.

[0188] The residuals between the observed and predicted values The prior VTEC is corrected to eliminate prediction bias and achieve more accurate real-time VTEC. The posterior covariance matrix quantifies the uncertainty after correction, providing a basis for subsequent GIVEI generation. This iterative process ensures that VTEC is dynamically updated over time, meeting "real-time" requirements.

[0189] Step 6: To ensure the reliability of the VTEC inversion results, the covariance matrix P of the Kalman filter in step 5 is used. k Extract the VTEC error standard deviation σ of each grid point VTEC , generates the Grid Ionospheric Vertical Error Indicator (GIVEI), which provides the user end with a basis for integrity monitoring. The calculation formula is:

[0190] GIVEI=3·σ VTEC ;

[0191] In areas with sparse observation data such as the ocean and polar regions, the inversion error is relatively large due to the small number of observations. The GIVEI value can be automatically increased by increasing the weight coefficient β of the error assessment:

[0192] GIVEI sparse =β·3·σ VTEC ;

[0193] During special space weather events such as magnetic storms, the abnormal term ΔVTEC, which is the change in ionospheric electron density during magnetic storms, is considered. storm , recalculate the error standard deviation Then we get the GIVEI value applicable to magnetic storms:

[0194]

[0195]

[0196] By dynamically adjusting the error weights in sparse areas and during magnetic storms, the problem of underestimation of errors in traditional models under extreme space weather conditions is solved.

[0197] Step 7: To adapt to the message length limit of Galileo HAS protocol, block differential compression coding is implemented on the VTEC grid and GIVEI data generated in step 6, giving priority to ensuring the accuracy of the underlying grid. i and VTEC i+1 , transmit the difference ΔVTEC i =VTEC i+1 -VTEC i In this way, a higher compression ratio is achieved to meet the real-time broadcast requirements. After adding the top-level grid VTEC data and GIVEI data, ensure that the total message length L total To meet the 26-page message limit stipulated by the HAS protocol, let the number of bytes that each page of message can accommodate be B, then:

[0198] L total =L bottom +L top +L GIVEI ≤26·B;

[0199] Among them, L bottom is the length of the compressed underlying grid data, L top is the top grid data length, L GIVEI The length of the GIVEI data.

[0200] Step 8. Based on the compressed VTEC grid data broadcasted in step 7, the user receiver can use bilinear interpolation to calculate the VTEC value corresponding to its own position according to the VTEC values ​​of the neighboring grid points, thereby achieving high-precision ionospheric delay correction. The coordinates of the user position in the grid are (x, y), and the coordinates and VTEC values ​​of the four neighboring grid points around it are (x1, y1, VTEC1), (x1, y2, VTEC2), (x2, y1, VTEC3), and (x2, y2, VTEC4). The VTEC value corresponding to the user position is user The calculation formula is:

[0201]

[0202] This embodiment utilizes precise orbit, clock correction and other parameters broadcast by HAS via the E6-B signal, combined with multi-frequency observation data and the Kalman filter algorithm, to achieve real-time ionospheric vertical TEC inversion in areas without global network coverage. Through dual-layer ionospheric grid modeling, ionospheric puncture point mapping and error assessment technology, high-precision VTEC products are generated and broadcast globally. The user end can obtain the ionospheric delay correction value for a specific location through interpolation calculation based on its own location, thereby serving high-precision applications. This embodiment breaks through the traditional ground-based augmentation network dependence and can still provide continuous and reliable ionospheric delay correction in scenarios such as oceans and deserts, providing technical support for high-precision applications such as drone navigation and autonomous driving.

[0203] The embodiments described above are merely descriptions of preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Without departing from the spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by persons skilled in the art should fall within the scope of protection defined by the claims of the present invention.

Claims

1. A global ionospheric TEC inversion method based on the Galileo HAS service, characterized in that: include: Acquire correction information through the HAS service, preprocess the observation data of the reference station according to the correction information, and obtain preprocessed carrier phase observation values; Calculating carrier phase ambiguity and cleaned carrier phase observation values ​​based on the preprocessed carrier phase observation values; Calculating the slant-path ionospheric total electron content (STEC) based on the carrier phase ambiguity and the cleaned carrier phase observation value; A double-layer ionospheric network model is constructed based on the geographical coordinates of the ionospheric puncture point and the ionospheric total electron content (STEC) of the slant path. According to the double-layer ionosphere network model, the vertical total electron content of the ionosphere is inverted using Kalman filtering to obtain the vertical total electron content of the ionosphere and compress it.

2. The global ionospheric TEC inversion method based on the Galileo HAS service according to claim 1, characterized in that: Preprocessing the observation data of the reference station according to the correction information includes: Calculating the orbit correction and clock correction at the current moment based on the correction information; Calculating the real-time precise clock error according to the clock error correction number; The observation data of the reference station is preprocessed according to the orbit correction number and the real-time precise clock error.

3. The global ionospheric TEC inversion method based on the Galileo HAS service according to claim 2, characterized in that: Preprocessing the observation data of the reference station according to the orbit correction number and the real-time precise clock error includes: Among them, Φ' is the carrier phase observation value after preprocessing, Φ is the carrier phase observation value, f is the carrier frequency, is the real-time precise clock error after HAS correction, C is the speed of light, Δr ECEF is the orbit correction number under ECEF, and u is the unit vector from the satellite to the receiver.

4. The global ionospheric TEC inversion method based on the Galileo HAS service according to claim 3, characterized in that: Calculating the carrier phase ambiguity and the cleaned carrier phase observations involves: Calculating double-difference observations and ionospheric free combination observations based on the preprocessed carrier phase observations; Using the double-difference observations and the ionospheric free combination observations, fixed wide-lane ambiguities and floating-point ambiguities are obtained respectively; According to the fixed wide lane ambiguity and the floating point ambiguity, the geometric free ambiguity is calculated and resolved to obtain the carrier phase ambiguity and the purified carrier phase observation value.

5. The global ionospheric TEC inversion method based on the Galileo HAS service according to claim 4, characterized in that: Calculating double-difference observations involves: in, is the double-difference observation, is the wide-lane ambiguity, Φ'1 and Φ'2 represent the carrier phase observations of the two frequency bands that have been corrected by the orbit and clock corrections of the HAS; j and k represent the inter-station and inter-satellite differences, respectively; A and B represent two different receivers or satellites, and f1 and f2 are the carrier frequencies of the two different frequency bands of the Galileo satellite navigation system.

6. The global ionospheric TEC inversion method based on the Galileo HAS service according to claim 1, characterized in that: Building a two-layer ionospheric network model includes: Processing the slant path ionospheric total electron content (STEC), calculating a regional STEC spatial gradient, stratifying the ionosphere according to the regional STEC spatial gradient, and dividing the network resolution; Based on the geographical coordinates of the ionospheric puncture point and STEC, a linear relationship between STEC and the vertical total electron content of the grid point corresponding to the ionospheric puncture point is established; According to the linear relationship, the vertical total electron content is filled into the corresponding grid points as the initial value to construct the double-layer ionosphere network model.

7. The global ionospheric TEC inversion method based on the Galileo HAS service according to claim 6, characterized in that: Calculation of regional STEC spatial gradients includes: Where ▽STEC is the regional STEC spatial gradient, To correct the magnetic latitude, λ LT For local time.

8. The global ionospheric TEC inversion method based on the Galileo HAS service according to claim 1, characterized in that: The inversion of the vertical total electron content of the ionosphere based on the Kalman filter includes: S1. Constructing an initialization state vector according to VTEC in the double-layer ionospheric network model; S2, perform time update and predict the prior state vector and prior covariance matrix at the current moment; S3, based on the current moment’s prior state vector and prior covariance matrix, combined with the STEC observation value, constructs the observation equation and calculates the Kalman gain; S4. Update the state estimate and the posterior covariance matrix according to the Kalman gain, output the inverted real-time VTEC value, add one at each moment, and return to S2.

9. The global ionospheric TEC inversion method based on the Galileo HAS service according to claim 1, characterized in that: Before obtaining the ionospheric vertical total electron content for compression, the following steps are involved: extracting the VTEC error standard deviation of each grid point from the covariance matrix of the Kalman filter and performing error evaluation on the inversion results.

10. The global ionospheric TEC inversion method based on the Galileo HAS service according to claim 1, characterized in that: Obtaining the vertical total electron content of the ionosphere for compression includes: L total =L bottom +L top +L GIVEI ≤26 B Among them, L total is the total message length, L bottom is the length of the compressed underlying grid data, L top is the top grid data length, L GIVEI is the length of GIVEI data, and B is the number of bytes that each page of telegram can accommodate.

Images