A global ionospheric TEC inversion method based on Galileo HAS service

Abstract

The application relates to a global ionospheric TEC inversion method based on a Galileo HAS service, which comprises the following steps: obtaining correction information through the HAS service, pre-processing observation data of a reference station according to the correction information, and obtaining carrier phase observation values after preprocessing; calculating carrier phase ambiguity and purified carrier phase observation values according to the carrier phase observation values after preprocessing; calculating slant path total electron content (STEC) according to the carrier phase ambiguity and the purified carrier phase observation values; constructing a double-layer ionospheric network model according to ionospheric piercing point geographical coordinates and the STEC; and inversing ionospheric vertical total electron content by using Kalman filtering, obtaining the ionospheric vertical total electron content and compressing the ionospheric vertical total electron content. The application breaks through the traditional network dependence of ground-based enhancement and can still provide continuous and reliable ionospheric delay correction in marine, desert and other scenes.

Description

Technical Field

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

[0002] Ionospheric delay is one of the main error sources in high-precision applications of Global Navigation Satellite System (GNSS), and its spatiotemporal variation characteristics make real-time accurate correction a technical challenge. Traditional methods for retrieving total ionospheric electron content (TEC) heavily rely on ground-based augmentation systems (GBAS) or precise ephemeris data acquired post-hoc via the internet. However, in areas without network coverage, such as polar regions and open oceans, or in emergency scenarios such as communication disruptions caused by disasters, existing technologies struggle to meet the requirements for real-time performance, global coverage, and robustness.

[0003] In recent years, the introduction of the Galileo High Accuracy Service (HAS) has provided a new approach for global real-time ionospheric monitoring. HAS directly broadcasts precise orbits and clock corrections for GPS / Galileo satellites via the E6-B signal, theoretically eliminating reliance on terrestrial communication networks. Current research largely focuses on verifying the performance of HAS in precise positioning, while existing ionospheric modeling methods face two major technical bottlenecks: First, regional ionospheric models exhibit significantly increased errors at low latitudes or during geomagnetic storms, with TEC prediction biases in the equatorial region reaching up to 20 TECUs (Total Electron Content Units, 1 TECU = 10^6 TECUs). 16 electrons / m 2 Secondly, while the global ionospheric grid has a wide coverage, it suffers from insufficient spatiotemporal resolution and poor timeliness. Although international GNSS services provide high-frequency TEC data through real-time services, they still rely on internet transmission and cannot serve communication-constrained scenarios.

[0004] The core technical architecture of HAS services relies on the highly stable space segment and high-precision ground segment of the Galileo system. The space segment uses Medium Earth Orbit (MEO) satellites with E6 band broadcasting capabilities to build a global enhanced signal broadcasting platform. The ground segment consists of a closed-loop system composed of multi-frequency, multi-mode monitoring stations, high-precision atomic clocks, and a data processing center. It uses multi-frequency observation data combined with Kalman filtering algorithms to generate high-precision corrections such as satellite precise orbits and clock bias corrections in real time, and broadcasts them with low latency through the E6-B channel, supplemented by a data protection mechanism to ensure reliable transmission.

[0005] Galileo HAS, with its high-precision orbit and clock error correction capabilities, has broad application prospects in various fields such as UAV navigation, surveying and mapping, and disaster monitoring. Current ionospheric delay correction methods mainly rely on ground-based augmentation networks or post-hoc precise ephemeris, which are insufficient to cover network-free areas such as oceans and deserts. Furthermore, the spatiotemporal resolution of global ionospheric models is insufficient to adapt to dynamic changes in the ionosphere in real time, especially under complex disturbances such as geomagnetic storms or equatorial anomalies, where errors increase significantly, leading to a decrease in the stability of high-precision positioning. While existing Galileo HAS services provide satellite-based precise correction, they lack real-time adaptive modeling capabilities for ionospheric disturbances, limiting their application in dynamic scenarios such as UAV 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 Galileo HAS services, which breaks through the network dependence of traditional ground-based augmentation 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 UAV navigation and autonomous driving.

[0007] To achieve the above objectives, the present invention provides the following solution:

[0008] A global ionospheric TEC inversion method based on Galileo HAS services includes:

[0009] Correction information is obtained through HAS service, and the observation data of the reference station is preprocessed based on the correction information to obtain the preprocessed carrier phase observation value.

[0010] Based on the preprocessed carrier phase observations, calculate the carrier phase ambiguity and the cleaned carrier phase observations;

[0011] The total electron content (STEC) of the tilted path ionosphere is calculated based on the carrier phase ambiguity and the purified carrier phase observation.

[0012] A two-layer ionospheric network model was constructed based on the geographical coordinates of the ionospheric puncture point and the total electron content (STEC) of the ionospheric ionosphere along the inclined path.

[0013] Based on the aforementioned two-layer ionospheric network model, the total vertical electron content of the ionosphere is inverted using Kalman filtering to obtain and compress the total vertical electron content of the ionosphere.

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

[0015] Calculate the orbital correction and clock error correction for the current moment based on the correction information;

[0016] Calculate the real-time precision clock error based on the clock error correction value;

[0017] The observation data of the reference station are preprocessed based on the orbital correction and the real-time precision clock difference.

[0018] Optionally, preprocessing the observation data of the reference station based on the orbital correction and the real-time precision clock difference includes:

[0019]

[0020] Where Φ' is the preprocessed carrier phase observation, Φ is the carrier phase observation, and f is the carrier frequency. The real-time precision clock error is HAS corrected, where C is the speed of light, and Δr is the speed of light. ECEF Here, represents the orbital correction under ECEF, and u is the unit vector from the satellite to the receiver.

[0021] Optionally, the calculation of carrier phase ambiguity and cleaned carrier phase observations includes:

[0022] Based on the preprocessed carrier phase observations, calculate the double-difference observations and the ionospheric free combination observations;

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

[0024] Based on the fixed wide-lane ambiguity and the floating-point ambiguity, the geometric free ambiguity is calculated, and the carrier phase ambiguity and the purified carrier phase observation value are obtained by solving the calculation.

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

[0026]

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

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

[0029] The total electron content (STEC) of the ionosphere along the inclined path is processed, the regional STEC spatial gradient is calculated, and the ionosphere is layered and the network resolution is divided according to the regional STEC spatial gradient.

[0030] Based on the geographic coordinates of the ionospheric puncture point and STEC, a linear relationship is established between STEC and the vertical total electron content of the corresponding grid point where the ionospheric puncture point is located.

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

[0032] Optionally, calculating the regional STEC spatial gradient includes:

[0033]

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

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

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

[0037] S2. Perform time updates and predict the prior state vector and prior covariance matrix at the current time.

[0038] S3. Based on the prior state vector and prior covariance matrix at the current moment, construct the observation equation and calculate the Kalman gain using the STEC observations.

[0039] S4. Update the state estimate and posterior covariance matrix based on the Kalman gain, output the inverted real-time VTEC value, increment by one at time step S2.

[0040] Optionally, before compressing the total vertical electron content of the ionosphere, the following steps are taken: extracting the standard deviation of the VTEC error for each grid point from the covariance matrix of the Kalman filter, and evaluating the error of the inversion results.

[0041] Optionally, compressing the vertical total electron content of the ionosphere includes:

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

[0043] Among them, L total L represents the total message length. bottom L represents the length of the compressed underlying grid data. top L represents the length of the top-level grid data. GIVEI B is the length of the GIVEI data, and B is the number of bytes that can be held per page of the message.

[0044] The beneficial effects of the present invention are: (1) Satellite-based broadcasting, global coverage: without the need for ground base stations or the Internet, 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 and ocean regions.

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

[0046] (3) Dynamic construction and optimization of the double-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 height of the bottom and top layers is adjusted in real time to solve the problem of ionospheric structural misalignment during low latitudes or geomagnetic storms; by setting differentiated grid resolution, the grid resolution is dynamically refined to improve the ability to capture detailed changes in the ionosphere; finally, dynamic puncture point correction ensures matching with the actual electron density distribution and reduces mapping errors.

[0047] (4) Anti-interference and robustness: Through real-time integrity monitoring (GIVEI) and dynamic error assessment mechanism, the model effectively copes with complex ionospheric disturbances, ensuring the stability and reliability of the inversion results. Dynamic adjustment of error weights under sparse regions and extreme space weather further enhances the model's adaptability. Attached Figure Description

[0048] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

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

[0050] Figure 2 This is a flowchart illustrating the construction process of a double-layer ionosphere mesh model according to an embodiment of the present invention.

[0051] Figure 3 This is a flowchart of the Kalman filter VTEC inversion process according to an embodiment of the present invention. Detailed Implementation

[0052] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0053] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0054] A global ionospheric TEC inversion method based on Galileo HAS services includes:

[0055] Correction information is obtained through HAS service, and the observation data of the reference station is preprocessed based on the correction information to obtain the preprocessed carrier phase observation values.

[0056] Based on the preprocessed carrier phase observations, calculate the carrier phase ambiguity and the cleaned carrier phase observations;

[0057] The total electron content (STEC) of the tilted path ionosphere is calculated based on the carrier phase ambiguity and the purified carrier phase observation.

[0058] A two-layer ionospheric network model was constructed based on the geographical coordinates of the ionospheric puncture point and the total electron content (STEC) of the ionospheric ionosphere along the inclined path.

[0059] Based on the two-layer ionospheric network model, the total vertical electron content of the ionosphere is inverted using Kalman filtering to obtain and compress the total vertical electron content of the ionosphere.

[0060] Furthermore, the preprocessing of the reference station's observation data based on the correction information includes:

[0061] Calculate the orbital correction and clock error correction for the current moment based on the correction information;

[0062] Calculate the real-time precision clock error based on the clock error correction.

[0063] The observation data from the reference station are preprocessed based on orbital corrections and real-time precision clock errors.

[0064] Furthermore, the preprocessing of the reference station's observation data based on orbital corrections and real-time precision clock errors includes:

[0065]

[0066] Where Φ' is the preprocessed carrier phase observation, Φ is the carrier phase observation, and f is the carrier frequency. The real-time precision clock error is HAS corrected, where C is the speed of light, and Δr is the speed of light. ECEF Here, represents the orbital correction under ECEF, and u is the unit vector from the satellite to the receiver.

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

[0068] Based on the preprocessed carrier phase observations, calculate the double-difference observations and the ionospheric free combination observations;

[0069] By using double-difference observations and ionospheric free combination observations, fixed wide-lane ambiguity and floating-point ambiguity are obtained respectively;

[0070] Based on the fixed wide-lane ambiguity and floating-point ambiguity, the geometric free ambiguity is calculated, and the carrier phase ambiguity and the purified carrier phase observation value are obtained by solving the problem.

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

[0072]

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

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

[0075] The total electron content (STEC) of the ionosphere along the inclined path is processed, the regional STEC spatial gradient is calculated, and the ionosphere is layered according to the regional STEC spatial gradient and the network resolution is divided.

[0076] Based on the geographic coordinates of the ionospheric puncture point and STEC, a linear relationship is established between STEC and the vertical total electron content of the corresponding grid point where the ionospheric puncture point is located.

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

[0078] Furthermore, the calculation of the STEC spatial gradient in the computational region includes:

[0079]

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

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

[0082] S1. Construct an initial state vector based on VTEC in the two-layer ionosphere network model;

[0083] S2. Perform time updates and predict the prior state vector and prior covariance matrix at the current time.

[0084] S3. Based on the prior state vector and prior covariance matrix at the current moment, construct the observation equation and calculate the Kalman gain using the STEC observations.

[0085] S4. Update the state estimate and posterior covariance matrix based on the Kalman gain, output the inverted real-time VTEC value, increment by one at time step S2.

[0086] Furthermore, before compressing the vertical total electron content of the ionosphere, the following steps are taken: extracting the standard deviation of the VTEC error for each grid point from the covariance matrix of the Kalman filter, and evaluating the error of the inversion results.

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

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

[0089] Among them, Ltotal L represents the total message length. bottom L represents the length of the compressed underlying mesh data. top L represents the length of the top-level grid data. GIVEI This is the length of the GIVEI data.

[0090] The following description, in conjunction with the accompanying drawings, further illustrates this embodiment:

[0091] like Figure 1 The figure shown is a global ionospheric TEC inversion method based on Galileo HAS service in this embodiment. The method includes the following steps:

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

[0093] Step 2: Layered resolution of carrier phase ambiguity and data cleansing;

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

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

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

[0097] Step Six: Error Assessment of Inversion Results and Generation of GIVEI;

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

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

[0100] The specific steps are as follows: First, preprocess the reference station observation data using the precise orbit and clock correction information provided by HAS to eliminate the influence of satellite orbit errors and clock errors on the observations, ensuring the accuracy of the reference for subsequent ambiguity resolution and TEC inversion. The receiver first receives the HAS correction information broadcast by the Galileo satellite via the E6-B band, which includes orbit correction (Δ... r Radial, Δ a Tangential, Δ c (Normal) and corresponding velocity correction Information such as clock error correction polynomial coefficients (a0, a1, a2) and 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 orbital corrections are expressed in a Spacecraft Body Fixed System (SBF) coordinate system, while the reference station observation data are collected in an Earth-Centered Earth-Fixed (ECEF) coordinate system, we first apply the formula... Calculate the orbital correction in the star-fixed system at time t, where The orbital correction vector at reference time t0. Correct the corresponding velocity vector. Then apply the rotation matrix:

[0103]

[0104] Where r broadcast and The position and velocity vectors of the satellite in the ECEF coordinate system calculated for broadcast ephemeris are obtained through the formula... This process transforms orbital corrections from the SBF coordinate system to the ECEF coordinate system. The transformed real-time orbital corrections are then superimposed onto the satellite orbit calculated from the broadcast ephemeris to obtain the satellite's real-time precise orbit r at time t. precise =r broadcast +Δr ECEF .

[0105] (2) Calculate the real-time clock error correction for the current time t from the clock error correction information: For clock error correction The real-time precision clock error is obtained by processing it using the following steps: In the formula and These are the satellite clock bias calculated from the broadcast ephemeris and the real-time precise clock bias after HAS correction, respectively, where C is the speed of light and δC is the speed of light. s This is a correction for relativistic effects. Relativistic effects can be calculated based on the satellite's position and velocity, using the following method:

[0106] Because traditional broadcast ephemeris calculations result in orbital errors at the decimeter level and clock errors at the meter level, while HAS corrections achieve centimeter-level and sub-nanosecond-level errors respectively, the real-time satellite orbital error Δr obtained from HAS data preprocessing is further utilized. ECEF and satellite clock bias The observations are corrected to eliminate satellite-side errors.

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

[0108]

[0109] In the formula, f is the carrier frequency, and u is the unit vector from the satellite to the receiver, pointing towards the satellite.

[0110] Step 2: Since traditional PPP-RTK or RTK ambiguity resolution also relies on ground-based reference station networks or post-hoc precise ephemeris analysis, it cannot be implemented in areas without network coverage, such as oceans and deserts. Furthermore, single ambiguity resolution methods are prone to failure when the ionosphere is active (e.g., in low-latitude regions or during geomagnetic storms). Therefore, this embodiment adopts a layered, progressive ambiguity resolution strategy to gradually improve accuracy. The specific steps are as follows:

[0111] (1) Utilizing the longer wavelength of wide-lane (HMW) combined observations, the initial ambiguity is quickly fixed, providing a foundation for subsequent accurate resolution. Based on the preprocessed carrier phase observation Φ', the wide-lane combined observation L is calculated using dual-frequency calculation. WL for:

[0112]

[0113] Where 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's f1 and f2 frequency bands, respectively.

[0114] The double-difference technique is used to eliminate common errors and reduce the wide-lane ambiguity N within a certain time. WL Fixed as an integer, for the wide lane combination observation value L WL Inter-station and inter-satellite differencing is performed, significantly improving the success rate and accuracy of wide-lane ambiguity fixation. Since residual orbital errors (especially those from broadcast ephemeris) still affect ambiguity fixation, the precise orbital correction Δr preprocessed from HAS correction information is directly used when constructing the double-difference observation equations. ECEF Clock error correction Replacing broadcast ephemeris parameters results in smaller double-difference residuals, improving the success rate and accuracy of wide-lane ambiguity fixation, and thus more accurately reflecting information such as ionospheric delay. At this point, the double-difference observations... Represented as:

[0115]

[0116] Through continuous adjustments Satisfying its integer characteristics, the wide-lane ambiguity is finally fixed. In the formula, Φ'1 and Φ'2 represent the two frequency band carrier phase observations that have been corrected by the orbit and clock bias corrections of HAS; j and k represent inter-station and inter-satellite differentials, respectively; A and B represent two different receivers or satellites.

[0117] (2) After fixing the wide-lane ambiguity, in order to eliminate the influence of the first-order ionospheric delay on the ambiguity resolution, it is necessary to construct an observation set without ionospheric error, and due to the wide-lane ambiguity... The integer properties of L are fixed and can be constrained. IF Ambiguity resolution in the ionosphere. First, calculate the ionospheric free combination observation value L. IF For a dual-frequency signal, its expression is: In the formula, Φ'1 and Φ'2 represent the two frequency band carrier phase observations that have been corrected by the orbit and clock bias corrections of HAS; r and s represent the receiver and satellite identifiers, respectively, which can be represented as different receivers A and B or satellites j and k.

[0118] Although the wavelength of this combination is relatively short, it completely eliminates first-order ionospheric delay, laying the foundation for high-precision floating-point ambiguity resolution. Combined with the precisely known coordinates r of the reference station... rec The floating-point ambiguity N is solved using the PPP-RTK model. IF The parameter dt is estimated using the least squares method. rec and The observation equation is: Where ρ=||r precise -r rec ||.

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

[0120] HAS provides precise satellite clock correction;

[0121] λ IF : Wavelength without ionospheric assemblies;

[0122] dt rec Receiver clock bias;

[0123] ε: Observation noise

[0124] The floating-point ambiguity to be solved.

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

[0126] (3) To eliminate the influence of frequency differences, through N IF and the fixed N WL Calculate the unambiguous geometric free fuzziness N GF Its expression is:

[0127]

[0128] After completing the ambiguity resolution, the fixed carrier phase ambiguities N1 and N2 of frequency bands f1 and f2 and the purified carrier phase observations Φ1” and Φ2” are output for use in step four to calculate the total electron content (STEC) along the slant path.

[0129]

[0130]

[0131] In the formula, The carrier wavelength.

[0132] Step 3: Calculate the geographic coordinates (λ) of the ionospheric puncture point (IPP). IPP , For a given satellite and reference station observation pair, IPP's... The calculation formula is:

[0133]

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

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

[0136] Because the MODIP-LT (Modified Dip Latitude-Local Time) coordinate system, through correction of magnetic tilt latitude and local time transformation, more accurately characterizes the magnetic field and solar radiation dependence of ionospheric electron density, especially significantly superior to geographic coordinate systems in the equator and polar regions, it can significantly reduce modeling errors. After converting the IPP geographic coordinates to MODIP-LT coordinates, the STEC is then mapped to the corresponding grid points. The following transformation is performed.

[0137]

[0138] In the formula, Δλ LT When making local corrections; To determine the geomagnetic latitude of the IPP, the International Geomagnetic Reference Field (IGRF) model needs to be called. t represents the date;

[0139] Based on the fixed single-frequency carrier phase ambiguities N1 and N2 from step two and the purified carrier phase observations Φ1” and Φ2”, the slant total electron content (STEC) of the ionosphere is calculated using the following formula:

[0140] Step 4: Construct a two-layer ionospheric grid model, dynamically adjusting the layer height based on real-time STEC spatial gradients and the MODIP-LT coordinate system. This avoids the misalignment between the IPP and the actual electron density concentration areas caused by a fixed height, which introduces mapping errors. It overcomes the limitations of fixed-height layering, allowing the two-layer grid height to dynamically adapt to the ionospheric activity level (STEC gradient) and spatial distribution (MODIP partitions), solving the problem of misalignment between the puncture point and the actual ionospheric structure in low-latitude or geomagnetic storm scenarios, and improving the model's ability to cover global ionospheric differences.

[0141] like Figure 2 As shown, STEC observations from 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 for magnetic latitude, λ LT When the location is determined, the coordinates of the ionospheric puncture point are obtained through step three.

[0142] Layered height dynamic adaptation mechanism:

[0143] Real-time STEC spatial gradient In conjunction with the Modified Magnetochore Indices (MODIP) criterion, intelligent adjustment of the vertical segmentation height of the ionosphere is achieved:

[0144] (1) Dynamic formulaic adjustment of the active area:

[0145] when When an active ionospheric state is identified (such as the equatorial anomaly zone or geomagnetic storm events), the stratification height is dynamically calculated based on the zonal characteristics.

[0146] Low magnetic latitude zone (MODIP < 30°): Bottom layer height Top floor height Capture dramatic changes in electron density by expanding the vertical coverage;

[0147] Mid-magnetic latitude region (30°≤MODIP≤60°): Temporarily fine-tuned by ±20km based on the default altitude to adapt to sudden disturbances;

[0148] High magnetic latitude region (MODIP > 60°): In response to ionospheric compression caused by geomagnetic activity.

[0149] Prioritize resolution optimization in the transition region:

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

[0151] Calm zone fixed parameters:

[0152] when At that time, a preset fixed stratification height is used to reflect the stable distribution characteristics of electron density (such as the calm ionosphere at night in mid-latitudes).

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

[0154] Based on MODIP partitioning and The gradient uses a dual criterion to implement dynamic adjustment at three levels of resolution:

[0155] Base resolution:

[0156] The default grid size is 0.5° × 0.5°, suitable for the mid-magnetic latitude region. A peaceful scene.

[0157] Gradient-triggered encryption:

[0158] Low magnetic latitude region: Default resolution 0.3° × 0.3°, when Encryption is applied to 0.2°×0.2°. The transition encryption is reduced to 0.25°×0.25°.

[0159] Mid-magnetic latitude zone: Default 0.5° × 0.5°, only applicable to... Encryption is applied up to 0.3° × 0.3°.

[0160] High magnetic latitude region: default 0.4° × 0.4°. For medium gradients, the encryption is set to 0.25°×0.25°. For medium gradients, the encryption is set to 0.3°×0.3°.

[0161] Special event adaptation:

[0162] During geomagnetic storms, a 0.3°×0.3° dense grid is forcibly activated in the mid-latitude region, while during geomagnetic substorms, the high-latitude region switches to a 0.25°×0.25° mode to ensure modeling accuracy under extreme events.

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

[0164] Based on the adjusted stratification heights h1' and h'2, the geographic coordinates of the ionospheric puncture point were recalculated: Where h' iThe height of the bottom or top layer is dynamically adjusted to ensure that the IPP coordinates match the actual vertical structure of the ionosphere.

[0165] Establish the relationship between STEC and the vertical total electron content (VTEC) of the corresponding IPP grid point, and fill the corresponding grid point with VTEC as the initial value. Considering the projection function mf(Z), we have: STEC = VTEC init ·mf(θ);

[0166] The formula for calculating mf(θ) is: θ is the zenith distance at IPP; e represents the satellite's elevation angle; h2 and h1 represent the heights of the upper and lower layers of the double ionospheric model, respectively.

[0167] Step 5: After constructing the initial two-layer ionospheric network model, the VTEC of each grid point is dynamically updated based on Kalman filtering. Through a "prediction-correction" iterative mechanism, static errors from spatial interpolation biases in the initial ionospheric grid model and dynamic errors from rapid electron density changes caused by geomagnetic storms, equatorial anomalies, etc., are eliminated. Based on new STEC observations, the VTEC grid data is continuously optimized to ensure the inversion results closely reflect the instantaneous state of the ionosphere. For example... Figure 3 The specific steps are as follows:

[0168] Based on the VTEC in the two-layer ionospheric network model, a state equation is constructed using the initial VTEC value of each grid point in the two-layer grid model. init As the core, let the initial state vector be:

[0169]

[0170] Among them, VTEC i,k For the VTEC of the i-th grid point at time k init Value, δ r,0 These are initial parameters for satellite / receiver phase deviation and other states.

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

[0172] The state equation can be expressed as:

[0173] Where, Φ k,k-1 This is the state transition matrix (describing the change of VTEC from time k-1 to time k). The process noise (following a Gaussian distribution N(0,Q)) k-1 (This reflects errors that the model did not capture).

[0174] The static initial VTEC of the two-layer grid model is transformed 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 time k-1 The covariance matrix P k-1 Predict the prior state vector at time k. and prior covariance matrix

[0177]

[0178] Among them, Q k-1 This is the process noise covariance matrix. In the absence of new observation data, based on historical states and time-varying patterns, the VTEC and uncertainty (covariance) at the current moment are estimated, providing a "prior reference" for subsequent correction using observed values, ensuring dynamism and real-time performance.

[0179] (3) Calculate the Kalman gain observation update and Kalman gain calculation based on the observation equation:

[0180] With the prior state of (2) and prior covariance Based on this, observation equations are constructed using STEC observations. Calculate the Kalman gain K k :

[0181]

[0182] in, The observation vector contains the STEC observation values ​​for the corresponding time calculated in step three; it is mapped to the VTEC grid through the projection function mf(Z), H k For the observation matrix, establish the relationship between the state vector and the observed values; To observe the noise, we assume it follows a Gaussian distribution N(0,R). k ); R k To observe the noise covariance matrix.

[0183] Kalman gain is a weighting coefficient that balances prediction error (prior covariance) and observation error (observation noise). The larger the value, the stronger the influence of the observations on the correction. This step correlates the predicted values ​​of the static grid model with the dynamic STEC observations, providing quantified weights for subsequent corrections.

[0184] (4) Update the state estimate and posterior covariance matrix based on the Kalman gain, output the inverted real-time VTEC value, and increment it by one at time step (2):

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

[0186] Updated post-hoc covariance matrix (reflecting the uncertainty of the corrected VTEC) P k :

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

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

[0189] Step Six: To ensure the reliability of the VTEC inversion results, the covariance matrix P based on the Kalman filter from Step Five is used... k Extract the standard deviation σ of VTEC error for each grid point VTEC This generates a gridded ionospheric vertical error indicator (GIVEI), providing a basis for integrity monitoring at the user end. The calculation formula is as follows:

[0190] GIVEI = 3·σ VTEC ;

[0191] In sparsely observed regions of the ocean and polar regions, the inversion error is relatively large due to the limited number of observations. The GIVEI value can be automatically increased by increasing the weighting coefficient β for error assessment.

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

[0193] During special space weather events such as geomagnetic storms, the anomaly ΔVTEC of ionospheric electron density changes during geomagnetic storms is considered. storm Recalculate the standard deviation of the error This leads to the GIVEI value applicable during geomagnetic storms:

[0194]

[0195]

[0196] By dynamically adjusting the error weights during sparse regions and geomagnetic storms, the problem of underestimating errors in traditional models under extreme space weather conditions is addressed.

[0197] Step 7: To adapt to the message length limitations of the Galileo HAS protocol, block differential compression encoding is applied to the VTEC grid and GIVEI data generated in Step 6, prioritizing the accuracy of the underlying grid. For the VTEC values ​​of adjacent grid points within a certain partition of the underlying grid... i and VTEC i+1 Transmit its difference ΔVTEC i =VTEC i+1 -VTEC i This method achieves a high compression ratio, meeting the requirements of real-time broadcasting. After appending the top-level mesh VTEC data and GIVEI data, the total message length L is ensured. total To satisfy the 26-page message limit stipulated by the HAS protocol, let B be the number of bytes that each page of message can hold, then:

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

[0199] Among them, L bottom L represents the length of the compressed underlying grid data. top L represents the length of the top-level grid data. GIVEI This is the length of the GIVEI data.

[0200] Step 8: Based on the compressed VTEC grid data broadcast 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 neighboring grid points, achieving high-precision ionospheric delay correction. The user's position coordinates in the grid are (x, y), and the coordinates and VTEC values ​​of its four neighboring grid points are (x1, y1, VTEC1), (x1, y2, VTEC2), (x2, y1, VTEC3), and (x2, y2, VTEC4), respectively. Therefore, the VTEC value corresponding to the user's position is... user The calculation formula is:

[0201]

[0202] This embodiment utilizes precise orbit and clock correction parameters broadcast by HAS via the E6-B signal, combined with multi-frequency observation data and a 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 techniques, a high-precision VTEC product is generated and broadcast globally. Users can obtain ionospheric delay correction values ​​for specific locations through interpolation calculations, thus serving high-precision applications. This embodiment overcomes the network dependence of traditional ground-based augmentation systems, providing continuous and reliable ionospheric delay correction even in marine and desert environments, providing technical support for high-precision applications such as UAV navigation and autonomous driving.

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

Claims

1. A global ionospheric TEC inversion method based on Galileo HAS services, characterized in that, include: Correction information is obtained through the HAS service. Based on this correction information, the observation data from the reference station is preprocessed to obtain preprocessed carrier phase observations, including: Calculate the orbital correction and clock error correction for the current moment based on the correction information; Calculate the real-time precision clock error based on the clock error correction value; The observation data of the reference station are preprocessed based on the orbital correction and the real-time precision clock error: ； in, These are the preprocessed carrier phase observations. For carrier phase observations, For carrier frequency, This is the real-time precision clock bias after HAS correction. At the speed of light, For orbital corrections under ECEF, This is the unit vector from the satellite to the receiver; Based on the preprocessed carrier phase observations, calculate the carrier phase ambiguity and the cleaned carrier phase observations, including: Based on the preprocessed carrier phase observations, calculate the double-difference observations and the ionospheric free combination observations; Using the double-difference observations and the ionospheric free combination observations, fixed wide-lane ambiguity and floating-point ambiguity are obtained respectively; Based on the fixed wide-lane ambiguity and the floating-point ambiguity, calculate the geometric free ambiguity, and solve it to obtain the carrier phase ambiguity and the cleaned carrier phase observation value; Calculating double-difference observations includes: ； in, These are double-difference observations. For the ambiguity of the wide alley, , This represents the carrier phase observations of the two frequency bands that have been corrected for orbit and clock bias using HAS. , These represent inter-station and inter-satellite differential communication, respectively; A and B represent two different receivers or satellites. , These are the carrier frequencies of two different frequency bands of the Galileo satellite navigation system; The total electron content (STEC) of the tilted path ionosphere is calculated based on the carrier phase ambiguity and the purified carrier phase observation. A two-layer ionospheric network model was constructed based on the geographical coordinates of the ionospheric puncture point and the total electron content (STEC) of the ionospheric ionosphere along the inclined path. Based on the aforementioned two-layer ionospheric network model, the total vertical electron content of the ionosphere is inverted using Kalman filtering to obtain and compress the total vertical electron content of the ionosphere.

2. The global ionospheric TEC inversion method based on Galileo HAS service according to claim 1, characterized in that, Constructing a two-layer ionospheric network model includes: The total electron content (STEC) of the ionosphere along the inclined path is processed, the regional STEC spatial gradient is calculated, and the ionosphere is layered and the network resolution is divided according to the regional STEC spatial gradient. Based on the geographic coordinates of the ionospheric puncture point and STEC, a linear relationship is established between STEC and the vertical total electron content of the corresponding grid point where the ionospheric puncture point is located. Based on the linear relationship, the vertical total electron content is filled to the corresponding grid points as the initial value to construct the double-layer ionosphere network model.

3. The global ionospheric TEC inversion method based on Galileo HAS service according to claim 2, characterized in that, The calculation of the STEC spatial gradient in the region includes: ； in, For the regional STEC spatial gradient, To correct the magnetic latitude, Local time.

4. The global ionospheric TEC inversion method based on Galileo HAS service according to claim 1, characterized in that, The inversion of the vertical total electron content of the ionosphere based on Kalman filtering includes: S1. Construct an initialization state vector based on VTEC in the double-layer ionosphere network model; S2. Perform time updates and predict the prior state vector and prior covariance matrix at the current time. S3. Based on the prior state vector and prior covariance matrix at the current moment, construct the observation equation and calculate the Kalman gain using the STEC observations. S4. Update the state estimate and posterior covariance matrix based on the Kalman gain, output the inverted real-time VTEC value, increment by one at time step S2.

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

6. The global ionospheric TEC inversion method based on Galileo HAS service according to claim 1, characterized in that, Compressing the vertical total electron content of the ionosphere includes: in, Total message length The length of the compressed underlying grid data. The length of the top-level grid data. B is the length of the GIVEI data, and B is the number of bytes that can be held per page of the message.

Images