A Spatiotemporal Kriging-Based Ionospheric Grid and Threat Modeling Method

Abstract

This invention discloses an ionospheric grid and threat modeling method based on spatiotemporal Kriging, belonging to the field of satellite navigation ionospheric modeling technology. The method includes the following steps: S1, calibrating the observation station location; S2, calibrating instrument bias; S3, constructing an IPP sample dataset; S4, establishing an ionospheric grid model; S5, optimizing relevant model parameters; S6, estimating ionospheric delay and error; S7, extending the ring-shaped data stripping scheme; and S8, establishing an ionospheric threat model. This invention employs the aforementioned spatiotemporal Kriging-based ionospheric grid and threat modeling method to establish a high-precision real-time ionospheric grid model based on spatiotemporal Kriging. Considering the spatial undersampling threat in four-dimensional spatiotemporal coordinates, it establishes a highly intact near-real-time ionospheric threat model based on spatiotemporal Kriging. The proposed data stripping scheme, combined with a sample grouping strategy, provides a solution for real-time ionospheric threat model establishment, contributing to the advancement of ionospheric research, application, and development.

Description

Technical Field

[0001] This invention relates to the field of satellite navigation ionospheric modeling technology, and in particular to an ionospheric grid and threat modeling method based on spatiotemporal Kriging. Background Technology

[0002] Global Navigation Satellite Systems (GNSS) primarily consist of four core satellite navigation constellations, regional satellite navigation systems, and satellite navigation augmentation systems. Ionospheric delay, one of the most significant sources of error in GNSS, remains a critical and long-standing issue requiring close attention. How to establish an effective ionospheric model based on theoretical research and observational data to correct ionospheric delay along the propagation path of satellite navigation signals has always been a pressing problem in the field of satellite navigation. The ionosphere, located approximately 60–1000 km above the Earth's surface, is an important component of the Sun-Earth space environment. It possesses a complex multi-layered physical structure and spatiotemporal variation patterns, significantly impacting modern radio engineering systems and human space activities. The propagation path of GNSS signals in the ionosphere is curved, and the propagation speed changes accordingly, resulting in path delay, also known as ionospheric delay or ionospheric refraction error. This ultimately leads to deviations in the satellite-user distance measurements and positioning results obtained by the user receiver, reducing the accuracy and reliability of satellite navigation.

[0003] Commonly used single-frequency ionospheric correction methods include ionospheric broadcast models and grid models. Broadcast models transmit ionospheric empirical model parameters to single-frequency users via GNSS broadcast ephemeris data. Examples include GPS broadcast Klobuchar model parameters, Galileo broadcast NeQuick model parameters, and the BeiDou-3 Global Navigation Satellite System broadcast BDGIM model parameters. Users can calculate the ionospheric delay at any time and location based on the model parameters and the corresponding empirical model, but the accuracy of ionospheric error correction is limited. Grid models, based on the thin ionospheric layer assumption, divide the global or regional ionospheric layer into numerous grids.

[0004] Based on integrity requirements, GIVE is the upper confidence (99.9% confidence) limit for the estimated delay residuals of the ionospheric grid model. It is calculated using the nominal error standard deviation and the spatial undersampling threat. To mitigate integrity risks caused by ionospheric irregularities not observed by ground stations, especially in the relatively more ionospheric Asia-Pacific region, a safety margin needs to be added to the overall error standard deviation in addition to the IGP nominal error standard deviation (or model uncertainty). This safety margin is typically described using the spatial undersampling threat. However, current spatial undersampling threat models developed by scholars both domestically and internationally are often fixed constants based on long-term observation data from multiple stations, representing the maximum undersampling risk at any time in the entire service area. These models lack location- and time-varying characteristics and can only satisfy the needs of relatively calm high-latitude regions. They cannot provide a stable and reliable ionospheric safety margin in the mid-to-low latitude Asia-Pacific region, such as China, where the ionosphere is more active. Therefore, in addition to establishing a global or regional ionospheric grid model, it is also necessary to establish an ionospheric threat model that can generate spatial undersampled threats, while simultaneously meeting the accuracy and integrity requirements of global or regional GIM products, or the accuracy and integrity requirements of SBAS grid augmentation information and the positioning accuracy, integrity, availability, and continuity requirements of the SBAS system.

[0005] Currently, scholars and institutions both domestically and internationally widely employ spatial Kriging methods for real-time ionospheric modeling. Spatial Kriging, a least mean square estimation method, provides a smooth description of the spatially distributed IPP delay, which can match the stochastic structure of the ionosphere. However, spatial Kriging only utilizes spatial information under the thin ionospheric assumption, without simultaneously considering the temporal trend and spatial correlation of IPP data, thus ignoring the important temporal information contained in the spatiotemporal dataset. The less-discussed spatiotemporal Kriging method, on the other hand, uses both the sampling time and spatial location of the data as independent variables, transforming the spatial interpolation problem into a high-dimensional spatiotemporal interpolation problem, and has the potential to solve complex spatiotemporal anisotropic interpolation modeling problems such as the ionosphere. Although the ionospheric broadcast model has better real-time performance and higher correction accuracy than the ionospheric grid model, and neither model improves navigation and positioning as much as the dual-frequency ionospheric-free combination method, given the booming development of GNSS applications and the fact that single-frequency users remain the primary service targets, and considering the advantages and disadvantages of the ionospheric broadcast model (lacking integrity information) versus the ionospheric grid model (facilitating the generation of integrity parameters), continuing to develop high-precision real-time ionospheric grid modeling methods and high-integrity near-real-time ionospheric threat modeling methods is essential and has significant application value for both global / regional GIM products and regional SBAS grid augmentation information. Summary of the Invention

[0006] The purpose of this invention is to provide an ionospheric grid and threat modeling method based on spatiotemporal Kriging, and to establish a high-precision real-time ionospheric grid model based on spatiotemporal Kriging. Considering the spatial undersampling threat in four-dimensional spatiotemporal coordinates, a highly intact near-real-time ionospheric threat model based on spatiotemporal Kriging is established. The proposed data stripping scheme combined with the sample grouping strategy can provide a solution for the real-time establishment of ionospheric threat models, and help promote the research, application and development of the ionosphere.

[0007] To achieve the above objectives, this invention provides a spatiotemporal Kriging-based ionospheric grid and threat modeling method, comprising the following steps:

[0008] S1, Location of the calibration observation station

[0009] Using data from multiple ground observation stations collected by the master control station, and employing real-time dynamic technology or precise single-point positioning technology based on carrier phase observations, the precise three-dimensional coordinates of the ground observation stations are calculated through post-processing, and the positions of the observation stations are calibrated daily or monthly.

[0010] S2, Calibration Instrument Deviation

[0011] Using data from multiple ground observation stations collected by the master control station and the calibrated station location coordinates, the differential code deviation of GNSS satellites and the differential code deviation of ground observation station receivers that meet constant constraints are calibrated daily using ionospheric mathematical modeling methods or precise single-point positioning technology.

[0012] S3. Construct the IPP sample dataset

[0013] Using ground-based observation stations distributed globally or regionally, calculate the three-dimensional coordinates, observation elevation angle, and vertical ionospheric delay of the IPP between the observation station and the GNSS core constellation and / or regional system satellites, and construct an IPP sample dataset for ionospheric modeling together with the observation time.

[0014] S4. Establish an ionospheric grid model

[0015] Based on elevation angle, observation time, spatial distance, and number of IPPs, the IPP sample dataset was filtered in sequence. Using the linearly normalized IPP dataset, an ionospheric grid model based on spatiotemporal Kriging was established.

[0016] S5. Optimize relevant model parameters

[0017] We select the correlation function corresponding to the normalized IPP observation time and three-dimensional location, design a four-dimensional spatiotemporal correlation model of the stochastic process in the spatiotemporal Kriging model, construct the objective function using the correlation function and process variance, iteratively update the correlation model parameters to minimize the objective function, and then obtain the optimal correlation model parameters.

[0018] S6. Estimating Ionospheric Delay and Error

[0019] Using the ionospheric grid model and related model parameters, the ionospheric delay and nominal error variance of the IGP or IPP to be estimated are estimated, and then the estimation results are inversely normalized.

[0020] S7, Extended Circular Data Preemption Scheme

[0021] The planar two-dimensional annular data expropriation scheme is improved into a spatial three-dimensional annular data expropriation scheme, which uses the spatial three-dimensional coordinates of IPP to group the available IPP sample points;

[0022] S8. Establish an ionospheric threat model

[0023] Based on the fitting radius and inner loop radius, the IPP sample points are grouped in a loop, the sample IPPs are spatiotemporally Kriging modeled, and the undersampled IPPs are threat detected. After the loop ends, all threat samples are counted, and the maximum value of the sample is taken as the undersampled threat in the ionospheric space.

[0024] Preferably, in step S3, the method for calculating the three-dimensional coordinates of IPP is as follows:

[0025] IPP Latitude The calculation formula is:

[0026]

[0027] in, It's a user dimension. It is the geocentric angle between the user and the IPP sub-satellite point. and These are the azimuth and elevation angles of the GNSS satellites observed by the user. It is the approximate radius of the Earth's ellipse. It is the height of the ionosphere thin layer;

[0028] The formula for calculating IPP longitude is:

[0029]

[0030] The user's longitude is known to be If the user dimension meets And the user's latitude, geocentric angle, and azimuth satisfy... Or the user dimension meets And the user's latitude, geocentric angle, and azimuth satisfy... The formula for calculating IPP longitude is:

[0031] ;

[0032] Vertical ionospheric delay calculations include:

[0033] After cycle slip detection and repair, the length of the continuous arc segment also needs to be checked. If the length of the continuous arc segment is... Less than the predetermined minimum value If so, the TEC for that epoch will not be solved; if Then, using dual-frequency carrier phase observations without geometric combination to smooth dual-frequency code pseudorange observations without geometric combination, the code smoothing formula is:

[0034]

[0035] in, Dual-frequency code pseudorange has no geometric combination of observations. and These are satellites in and The pseudorange observations on the code; These are dual-frequency carrier phase observations without geometric combination. and These are satellites in and Carrier phase observations on;

[0036] By combining the differential code bias of the observation station receiver and the differential code bias of the GNSS satellite in the calibration dataset, as well as the code smoothing results, the vertical total electron content (VTEC) of the ionosphere is obtained:

[0037]

[0038] in, and These are the frequencies of the GNSS dual-frequency signal. It is the projection function at the ionospheric intersection. These are dual-frequency code pseudorange observations without geometric combination after code smoothing. It's the speed of light. and These are the station receiver differential code bias and the GNSS satellite differential code bias in the calibration dataset, respectively.

[0039] Finally, the vertical delay in the IPP sample dataset That is, the ionospheric VTEC at frequency The ionospheric delay caused by the above is calculated using the following formula:

[0040]

[0041] Among them, frequency It depends on the GNSS navigation signal selected by the user for positioning.

[0042] Preferably, in step S4, establishing the ionospheric grid model includes the following steps:

[0043] S41. Select the IGP to be estimated and calculate its estimated time. and three-dimensional coordinates ;

[0044] S42. Select a single IPP from the IPP sample dataset and extract the observation time of that single IPP. 3D coordinates Elevation angle observed by ground observation station and vertical ionospheric delay ;

[0045] S43. If the elevation angle observed by the ground observation station is higher than the lowest elevation angle ,Right now If yes, proceed to the next step; otherwise, select the next IPP and re-execute step S42.

[0046] S44, According to the sampling interval and maximum number of time samples The observation time of a single IPP is checked, and if it meets the requirements... If yes, proceed to the next step; otherwise, select the next IPP and re-execute step S42.

[0047] S45. Calculate the spatial distance between IGP and IPP based on their three-dimensional coordinates. , and the set fitting radius Compare, if satisfied If yes, proceed to the next step; otherwise, select the next IPP and re-execute step S42.

[0048] S46. Select samples from IPP points that meet the time and distance requirements. If the maximum fitted radius... Number of IPPs within Exceeding the predetermined upper limit ,Right now Then gradually decrease the fitting radius. To the minimum fitting radius ,until ;

[0049] If the number of IPPs within the maximum fitting radius is lower than the predetermined upper limit but not lower than the predetermined lower limit ,Right now If the number of IPPs within the maximum fitting radius is lower than the predetermined lower limit, then proceed to the next step; If so, the IGP to be estimated will not be modeled;

[0050] S47. Normalize the IPP that satisfies steps S43, S44, S45, and S46, i.e., use... The average of individual IPPs and standard deviation Construct a normalized IPP dataset (7);

[0051] S48. Based on the spatiotemporal Kriging method, a grid model of the ionosphere of the IGP to be estimated is established using the normalized IPP dataset.

[0052] Preferably, in step S47, the normalization method is as follows:

[0053]

[0054] in, It is the mean of the time, location, and delay of all IPP samples. It is the variance of the time, location, and delay of all sample IPPs.

[0055] Preferably, in step S48, establishing the ionospheric grid model of the IGP to be estimated includes the following steps:

[0056] Define the normalized vertical ionospheric delay at the i-th IGP or IPP. for:

[0057]

[0058] in, It is the normalization moment. It is the normalized three-dimensional position of IGP or IPP. , , , and It is the regression coefficient. Represents a random process. Indicates measurement noise;

[0059] Assume a steady-state stochastic process The mean of is zero, and the non-zero correlation covariance satisfies:

[0060]

[0061]

[0062] in, It is process variance. It has parameters Models related to stochastic processes;

[0063] Assume the mean and non-zero correlation covariance of the measurement noise satisfy the following:

[0064]

[0065]

[0066] in, The tilt factor, derived from the vertical ionospheric delay, is the tilt factor for obtaining the tilted ionospheric delay and is defined as the satellite elevation angle. The function; The variance of the measurement error is obtained by considering the dual-frequency code carrier smoothing CCL of the receiver DCB and the satellite DCB:

[0067]

[0068] in, It is a smooth length or arc segment count. It refers to the accuracy of pseudorange measurement. It refers to the accuracy of carrier phase measurement. It is the standard deviation of DCB at the ground station. It is the standard deviation of satellite DCB;

[0069] A linear combination of the vertical delays at the IPP of the samples within the fitting radius is performed to obtain the result at the normalized time. and normalized position Spatiotemporal Kriging estimates of IGP or IPP at the location:

[0070]

[0071] in, It is the number of IPP samples. It is a weight coefficient vector. It is the ionospheric observation vector; the unbiased constraint for the spatiotemporal Kriging estimate is:

[0072]

[0073]

[0074]

[0075] The weighting coefficients satisfying the unbiased constraints are obtained using the Lagrange multiplier method with equality constraints, and the result is as follows:

[0076]

[0077] in, It is the weighted matrix of IPP for all samples, derived from the stochastic process. covariance matrix and measuring noise covariance matrix constitute:

[0078]

[0079]

[0080]

[0081]

[0082]

[0083]

[0084] in, It is the covariance vector between all sample IPPs and the estimated IGP or IPP; assuming that the stochastic processes and measurement noise in the ionospheric model are uncorrelated, the nominal error variance of the spatiotemporal Kriging estimate is:

[0085]

[0086] By treating time as another independent variable besides IPP location, a spatiotemporal Kriging correlation model is constructed.

[0087] Preferably, in step S6, the final delay estimate at the IGP to be estimated is calculated using inverse normalization. and final error standard deviation This includes the following steps:

[0088] To ensure consistency of the vertical ionospheric delay dimensions between the sample IPP and the IGP to be estimated, it is also necessary to... and Perform the inverse normalization operation to obtain the final estimate. and standard deviation of error ,Right now:

[0089]

[0090] .

[0091] Preferably, in step S7, the planar two-dimensional annular data expropriation scheme is transformed into a spatial three-dimensional annular data expropriation scheme, and the available IPP sample points are grouped using the spatial three-dimensional coordinates of the IPP, including the following steps:

[0092] S701, Select the estimation time for the IGP to be estimated. Calculate its three-dimensional coordinates ;

[0093] S702. Select a single IPP from the available IPP sample points and extract its observation time. 3D coordinates Elevation angle observed by ground observation station ;

[0094] S703, If the elevation angle observed by the ground observation station is higher than the minimum elevation angle ,Right now If yes, proceed to the next step; otherwise, select the next IPP and re-execute step S702.

[0095] S704, According to the sampling interval and maximum number of time samples The observation time of a single IPP is checked, and if it meets the requirements... If yes, proceed to the next step; otherwise, select the next IPP and re-execute step S702.

[0096] S705, Set the fitting radius Minimum Fitting Radius ,Right now ;

[0097] S706. Calculate the spatial distance between IGP and IPP based on their three-dimensional coordinates. , and the set fitting radius Compare, if satisfied If the IPP is included in the threat IPP sample dataset, it is included in the dataset; otherwise, it is discarded.

[0098] S707. Select the next IPP and continue with steps S702-S706 until all IPPs have been filtered.

[0099] S708, Set Inner Ring Radius It is 0, that is ;

[0100] S709. In the threat IPP sample dataset, select the fitting radius. For IPPs located inside the ring but outside the ring, construct modeled IPP samples; select IPPs inside the ring to construct undersampled IPP samples.

[0101] S710, inner ring radius increased ,Right now Repeat step S709 until the outer ring radius equals the fitted radius, i.e. ;

[0102] S711, Fitting radius increased ,Right now Repeat steps S706~S710 until the fitted radius equals the maximum fitted radius, i.e. ;

[0103] S712. Select the next IGP to be estimated and repeat steps S701 to S711 until all IGPs to be estimated have been sampled.

[0104] Preferably, in step S8, for a single IGP to be estimated, in a single loop with a fixed fitting radius and inner loop radius, when the estimation error of the k-th undersampled IGP is greater than... times At that time, there is a threat, namely:

[0105]

[0106] The above formula is equivalent to:

[0107]

[0108] in, It is the confidence quantile corresponding to the ionospheric integrity risk; the ionospheric undersampling threat sample of the kth undersampled IPP is defined as:

[0109]

[0110] Based on the data expropriation scheme, undersampled threats corresponding to different fitting radii are identified by iteratively building a spatiotemporal Kriging model and searching for the maximum undersampled threat sample.

[0111]

[0112] Finally, the overall error variance, including the inflated nominal error variance and the threat of spatial undersampling, is:

[0113]

[0114] in, It is the nominal error inflation factor, calculated using the chi-square distribution of the number of IPPs minus the number of regression coefficients to be estimated. Calculations show that It is the probability of a false alarm. This is the probability of missed detection; correspondingly, the error envelope of the IGP to be estimated, i.e., the grid ionospheric vertical error GIVE, is:

[0115]

[0116] in, It is the confidence quantile corresponding to the 99.9% error envelope probability.

[0117] Therefore, the present invention employs the aforementioned spatiotemporal Kriging-based ionospheric grid and threat modeling method, the technical effects of which are as follows:

[0118] (1) This invention is the first to introduce the spatiotemporal Kriging method into ionospheric modeling;

[0119] (2) This invention designs a correlation model of the stochastic process and a covariance matrix of measurement noise in the spatiotemporal Kriging model. The weighted matrix generated by the two can improve the accuracy of ionospheric modeling.

[0120] (3) This invention provides a more accurate real-time ionospheric grid modeling method for global / regional GIM products and SBAS grid augmentation information;

[0121] (4) This invention extends the ring-shaped data stripping scheme for modeling undersampling threats in ionospheric space and provides an improvement direction for planar two-dimensional data stripping schemes such as the three-quadrant data stripping scheme;

[0122] (5) This invention provides a near-real-time ionospheric threat modeling method with higher integrity for global / regional GIM products and SBAS grid augmentation information;

[0123] (6) The data stripping scheme proposed in this invention, combined with the sample grouping strategy, can provide a solution for establishing an ionospheric threat model in real time, which will help promote the research, application and development of the ionosphere.

[0124] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0125] Figure 1 This is a schematic diagram of the overall implementation scheme of the ionospheric grid and threat modeling method based on spatiotemporal Kriging of the present invention;

[0126] Figure 2 A schematic diagram illustrating ionospheric grid and threat modeling;

[0127] Figure 3 This is a schematic diagram of GNSS observation data preprocessing;

[0128] Figure 4 Schematic diagram of ionospheric grid modeling;

[0129] Figure 5 A schematic diagram illustrating the optimization of relevant model parameters;

[0130] Figure 6 Schematic diagram for modeling ionospheric threats;

[0131] Figure 7 This is a schematic diagram of a three-dimensional annular data stripping scheme. Detailed Implementation

[0132] The technical solution of the present invention will be further described below with reference to the accompanying drawings and embodiments.

[0133] Unless otherwise defined, the technical or scientific terms used in this invention shall have the ordinary meaning as understood by one of ordinary skill in the art to which this invention pertains.

[0134] Example 1

[0135] (I) Overall Implementation Plan for Ionospheric Modeling

[0136] A comprehensive implementation plan for ionospheric grid and threat modeling based on spatiotemporal Kriging is as follows: Figure 1 As shown, the system mainly comprises three modules: GNSS observation data preprocessing, ionospheric grid modeling, and ionospheric threat modeling. The GNSS observation data preprocessing module provides necessary data support for ionospheric modeling (including the ionospheric grid modeling module and the ionospheric threat modeling module), including pre-calibrated observation station locations, satellite differential code bias (DCB), observation station differential code bias (DCB), and the location and total electron content (TEC) of sample IPPs. The ionospheric grid modeling module establishes a spatiotemporal Kriging model of the ionosphere, and after optimization of relevant model parameters, estimates the ionospheric delay of the IGP and its nominal error variance. The ionospheric threat modeling module uses a data stripping scheme to group the sample IPP data, then performs spatiotemporal Kriging modeling on the modeled IPPs, uses undersampled IPPs for threat detection, and finally establishes a spatial undersampled threat model.

[0137] By utilizing the calibrated observation station locations, receiver differential code offsets, and satellite differential code offsets, combined with the navigation ephemeris and dual-frequency observation data of the GNSS core constellation, the three-dimensional coordinates and vertical total electron content of the IPP point can be calculated. This constructs an IPP sample dataset containing the observation time, 3D coordinates, and vertical ionospheric delay for all IPPs. The specific method is described in the GNSS observation data preprocessing module. Using the IPP sample dataset, individual IGPs can be modeled one by one until the ionospheric grid and threat modeling for all IGPs is completed. The modeling scheme is as follows: Figure 2 As shown, the specific modeling steps are as follows:

[0138] (1) Select the IGP to be estimated and calculate its estimated time. and three-dimensional coordinates ;

[0139] (2) Select a single IPP in the IPP sample dataset and extract its observation time. 3D coordinates Elevation angle observed by ground observation station and vertical ionospheric delay ;

[0140] (3) If the elevation angle observed by the ground observation station is higher than the lowest elevation angle ,Right now If yes, proceed to the next step; otherwise, select the next IPP and re-execute step (2).

[0141] (4) Based on the sampling interval and maximum number of time samples The observation time of a single IPP is checked, and if it meets the requirements... If yes, proceed to the next step; otherwise, select the next IPP and re-execute step (2).

[0142] (5) Calculate the spatial distance between IGP and IPP based on their three-dimensional coordinates. , and the set fitting radius Compare, if satisfied If yes, proceed to the next step; otherwise, select the next IPP and re-execute step (2).

[0143] (6) Screen the IPP points that meet the time and distance requirements. If the maximum fitting radius is... Number of IPPs within Exceeding the predetermined upper limit ,Right now Then gradually decrease the fitting radius. To the minimum fitting radius ,until If the number of IPPs within the maximum fitting radius is lower than the predetermined upper limit but not lower than the predetermined lower limit. ,Right now If the number of IPPs within the maximum fitting radius is lower than the predetermined lower limit, then proceed to the next step; If so, the IGP to be estimated will not be modeled;

[0144] (7) Normalize the IPP that satisfies steps (3), (4), (5) and (6), that is, use The average of individual IPPs and standard deviation Construct a normalized IPP dataset. The specific normalization method is shown in formula (7).

[0145] (8) Based on the spatiotemporal Kriging method, a grid model of the ionosphere of the IGP to be estimated is established using the normalized IPP dataset, the relevant model parameters are optimized, and the vertical ionospheric delay at the IGP to be estimated is estimated. and nominal standard deviation of error ;

[0146] (9) Calculate the final delay estimate at the IGP to be estimated by inverse normalization. and final error standard deviation ;

[0147] (10) Based on the designed ring data stripping scheme, the IPP sample dataset is grouped and a spatiotemporal Kriging model is iteratively established to build an ionospheric threat model and obtain the spatial undersampling threat at the IGP to be estimated. ;

[0148] (11) Based on the ionospheric integrity risk and error envelope probability requirements, combined with the final error standard deviation of step (8) and the undersampling threat in step (9) Calculate the overall error standard deviation of the IGP to be estimated. And error envelope (grid ionospheric vertical error, GIVE);

[0149] (12) Select the next IGP to be estimated and repeat steps (1) to (10) until the overall error standard deviation of all IGPs to be estimated is calculated. And error envelope (grid ionospheric vertical error, GIVE).

[0150] (II) Preprocessing of GNSS observation data

[0151] The GNSS observation data preprocessing module takes navigation ephemeris and dual-frequency observations from GPS, BDS, Galileo, and GLONASS core constellation satellites as input, and outputs calibration datasets, IPP positions, and TECs, such as... Figure 3 As shown. In order to establish a high-precision ionospheric grid and threat model, this invention first needs to continuously calibrate the position of the ground observation station, the receiver differential code deviation, and the satellite differential code deviation. Utilizing the short-term stability and invariance of this calibration dataset, the calculation results of the previous day or month are regarded as true values ​​and substituted into the elevation angle calculation formula and the code smoothing equation, thereby calculating the position coordinates of the IPP and the total vertical electron content of the ionosphere in real time.

[0152] 2.1 Data Calibration

[0153] GNSS ground stations mainly consist of observation stations, master control stations, and injection stations. Numerous widely distributed observation stations at known locations continuously track and monitor downlink space signals from constellations such as GPS, BDS, Galileo, GLONASS, QZSS, and IRNSS, obtaining raw observation data (ephemeris, pseudorange, carrier phase, etc.) and transmitting it to the master control station. The master control station collects GNSS observation data from the observation stations, corrects satellite orbits and clock errors, and generates ephemeris information for each satellite, sending it to the injection stations. The injection stations are responsible for uploading the ephemeris information to the GNSS satellites, completing data transmission between the ground and satellites. The ionospheric grid and threat modeling method proposed in this invention belongs to the master control station algorithm. Based on GNSS observation data transmitted to the master control station from multiple ground stations (ideally more than 20 evenly distributed), it processes the data in real time to generate global or regional ionospheric grid models and ionospheric threat models. The applicability of the models, the location and number of IGPs depend on the distribution of the selected ground observation stations.

[0154] Because the coverage or applicable area of ​​ionospheric grid models and threat models is global or local, it cannot be achieved through a single observation station or a few observation stations. Moreover, the specific selection of ground observation stations and the applicable scope of the model need to be designed based on the actual station construction situation. This invention requires multiple ground observation stations within the service area to achieve its purpose.

[0155] The master control station receives raw observation data transmitted from the observation station in real time, but it does not need to calibrate the station's position, receiver differential code deviation, and GNSS satellite differential code deviation in real time. This is because the position has long-term stability and invariance, and the differential code deviation also has short-term stability and invariance. Therefore, it is only necessary to update the calibration dataset according to a fixed time frequency and using data within a fixed measurement segment. This invention suggests using the raw data from the previous day to calibrate the station position, receiver differential code deviation, and satellite differential code deviation once a day. The calibration dataset is used to separate instrument deviations and calculate ionospheric delay in real time. For the station position, mature real-time dynamic (RTK) technology based on carrier phase observations or precise point positioning (PPP) technology can be used to calculate the three-dimensional coordinates of the ground observation station through post-processing. For the receiver differential code deviation and GNSS satellite differential code deviation, it can be assumed that they are constant within a day, and the model parameters and instrument deviation can be estimated simultaneously using an ionospheric mathematical model. Commonly used ionospheric mathematical models include spherical harmonic functions, polynomial functions, generalized trigonometric series functions, and spherical cap harmonic functions. In addition to the step-by-step update method for the calibration dataset described above, the observation station location, the differential code offset of the observation station receiver, and the differential code offset of the GNSS satellite can also be used simultaneously as parameters to be estimated, and the calibration dataset can be directly estimated and updated in one step using PPP technology. This invention focuses on extracting ionospheric information and establishing an ionospheric grid and threat model. The calibration dataset only provides the data foundation for this invention and is not its focus; therefore, the already mature RTK, PPP, and ionospheric mathematical modeling techniques will not be described in detail.

[0156] 2.2 IPP Location Calculation

[0157] Based on the thin ionospheric layer hypothesis, the free electrons in the ionosphere along the GNSS signal propagation path are concentrated on an infinitely thin sphere at a specific height. The choice of the thin layer height varies depending on the industry context and application purpose, but is usually selected between the peak ionospheric density, i.e., 350km to 450km. For example, the thin layer height chosen for IGS ionospheric products is 450km, while the thin layer height chosen for SBAS enhanced ionospheric products is 350km. The ionospheric puncture point (IPP) is defined as the intersection of the line connecting the satellite to the ground receiver and the thin ionospheric layer. It is assumed that all the total electron content (STEC) of the ionosphere along the line of sight (tilted) is compressed onto the IPP.

[0158] IPP Latitude The calculation formula is:

[0159] (1)

[0160] in, It is the latitude of the user (ground observation station receiver). It is the geocentric angle between the user and the IPP sub-satellite point. and These are the azimuth and elevation angles of the GNSS satellites observed by the user. It is the approximate radius of the Earth's ellipse. It is the height of the ionosphere.

[0161] The formula for calculating IPP longitude is:

[0162] (2)

[0163] Specifically, if the user dimension meets the requirements And the user's latitude, geocentric angle, and azimuth satisfy... Or the user dimension meets And the user's latitude, geocentric angle, and azimuth satisfy... The formula for calculating IPP longitude is:

[0164] (3)

[0165] This invention uses the three-dimensional coordinates of the IPP in the geocentric geofixed (ECEF) coordinate system to index the IPP. Therefore, after calculating the longitude and latitude of the IPP, it is also necessary to combine the selected thin layer height to perform coordinate transformation on the IPP location. The coordinate-transformed IPP location is used to construct the IPP sample dataset.

[0166] 2.3 TEC Solution

[0167] To accurately extract ionospheric TEC, in addition to pre-calibrating the differential code bias of the ground station receiver and the satellite, it is also necessary to detect and correct any cycle slips that may exist in the carrier phase observations. Obstacles blocking the satellite signal, external interference, or adverse receiver conditions can all cause a brief loss of signal lock, resulting in interruptions in integer counting and cycle slips. Methods for cycle slip detection and correction are well-established. For example, cycle slip detection can be performed by subtracting epochs of ionospheric residual combination observations or by subtracting epochs of integer ambiguity from wide-lane observations; cycle slip correction can be performed using double-difference or triangular observations. However, cycle slip detection and correction only ensure that there are no systemic biases caused by cycle slips in the carrier phase observations between consecutive epochs. To improve the accuracy of determining ionospheric TEC using dual-frequency observations, it is also necessary to determine the continuity of epochs without cycle slips (including those after correction), and to use dual-frequency observations within continuous arcs to determine the magnitude of the combination of integer ambiguity and instrument bias, thereby obtaining the absolute original ionospheric TEC observation information. Therefore, after cycle slip detection and repair, it is also necessary to check the length of the continuous arc segment. If the length of the continuous arc segment is... Less than the predetermined minimum value ,Right now If so, the TEC for that epoch will not be solved; if That is, continuous If no cycle slip occurs in the dual-frequency carrier phase observations within a given epoch, then the dual-frequency code pseudorange observations without geometric combinations are smoothed using the dual-frequency carrier phase observations without geometric combinations. The code smoothing formula is as follows:

[0168] (4)

[0169] in, Dual-frequency code pseudorange has no geometric combination of observations. and These are satellites in and The pseudorange observations on the code; These are dual-frequency carrier phase observations without geometric combination. and These are satellites in and The carrier phase observation values.

[0170] By combining the differential code bias of the observation station receiver and the differential code bias of the GNSS satellite in the calibration dataset, as well as the code smoothing results, the vertical total electron content (VTEC) of the ionosphere can be obtained:

[0171] (5)

[0172] in, and These are the frequencies of the GNSS dual-frequency signal. It is the projection function at the ionospheric intersection. These are dual-frequency code pseudorange observations without geometric combination after code smoothing. It's the speed of light. and These are the station receiver differential code bias and the GNSS satellite differential code bias in the calibration dataset, respectively.

[0173] Finally, the vertical delay in the IPP sample dataset That is, the ionospheric VTEC at frequency The ionospheric delay caused by the above is calculated using the following formula:

[0174] (6)

[0175] Among them, frequency Depending on the GNSS navigation signal selected for user positioning, for example, SBAS users use GPS L1C / A code signals for positioning, and their frequency... .

[0176] (III) Ionospheric grid modeling

[0177] To simplify calculations, this invention transforms dimensional time, location, and delay into dimensionless parameters. Specifically, within each ionospheric information update interval, after constructing the IPP sample dataset using the GNSS observation data preprocessing module, it is necessary to use the sample mean and sample standard deviation to calculate the observation time of all IPPs. (in seconds within a week), 3D coordinates (Geocentric Earth-Fixed Coordinate System) and Vertical Delayed Observations of the Ionosphere (After removing DCB) Perform linear normalization:

[0178] (7)

[0179] in, It is the mean of the time, location, and delay of all IPP samples. It is the variance of the time, location, and delay of all sample IPPs.

[0180] The ionospheric grid modeling module takes the normalized IPP dataset processed as described above as input and outputs the estimated delay and nominal error at the IGP. This invention uses spatiotemporal Kriging as its core. First, a local ionospheric grid model is established with the spatiotemporal center of the sample IPP sample as the origin. Then, the correlation model of the stochastic process in the model and the covariance matrix of the measurement noise are designed. Finally, the optimized correlation model parameters are used to estimate the IGP information.

[0181] 3.1 Spatiotemporal Kriging Modeling

[0182] This invention, based on the traditional spatial Kriging method that uses only three-dimensional spatial coordinates as independent variables, introduces a fourth independent variable—observation time—to jointly construct an ionospheric model based on four-dimensional spatiotemporal Kriging. The modeling process is as follows: Figure 4 As shown.

[0183] Define the normalized vertical ionospheric delay at the i-th IGP or IPP. for:

[0184] (8)

[0185] in, It is the normalization moment. It is the normalized three-dimensional position of IGP or IPP. , , , and It is the regression coefficient. Represents a random process. This represents measurement noise. The normalization operation is equivalent to establishing the origin of the local ionospheric model's coordinates at the center of the sample IPP, rather than at the IGP.

[0186] Assume a steady-state stochastic process The mean of is zero, and the non-zero correlation covariance satisfies:

[0187] (9)

[0188] (10)

[0189] in, It is process variance. It has parameters The stochastic process correlation model. Similarly, assume that the mean and non-zero correlation covariance of the measurement noise satisfy:

[0190] (11)

[0191] (12)

[0192] in, The tilt factor, derived from the vertical ionospheric delay, is the tilt factor for obtaining the tilted ionospheric delay and is defined as the satellite elevation angle. The function. It is the variance of the measurement error, which is determined by considering the dual frequencies of the receiver DCB and the satellite DCB. and Code loading smoothing (CCL) yields:

[0193] (13)

[0194] in, It is a smooth length or arc segment count. It refers to the accuracy of pseudorange measurement. It refers to the accuracy of carrier phase measurement. It is the standard deviation of DCB at the ground station. It is the standard deviation of satellite DCB.

[0195] By linearly combining the vertical delays at the IPP of the samples within the fitting radius, we can obtain the result at the normalized time. and normalized position Spatiotemporal Kriging estimates of IGP or IPP at the location:

[0196] (14)

[0197] in, It is the number of IPP samples. It is a weight coefficient vector. This is the ionospheric observation vector. The unbiased constraint for this spatiotemporal Kriging estimate is:

[0198] (15)

[0199] (16)

[0200] (17)

[0201] The weighting coefficients satisfying the unbiased constraints can be obtained using the Lagrange multiplier method with equality constraints, and the result is:

[0202] (18)

[0203] in, It is the weighted matrix of IPP for all samples, derived from the stochastic process. covariance matrix and measuring noise covariance matrix constitute:

[0204] (19)

[0205] (20)

[0206] (twenty one)

[0207] (twenty two)

[0208] (twenty three)

[0209] (twenty four)

[0210] in, This is the covariance vector between all sample IPPs and the IGP (or IPP) to be estimated. Assuming that stochastic processes and measurement noise in the ionospheric model are uncorrelated, the nominal error variance of the spatiotemporal Kriging estimate is:

[0211] (25)

[0212] By treating time as another independent variable besides IPP location, a spatiotemporal Kriging correlation model can be constructed, and the parameters of this model can be obtained through iterative optimization.

[0213] 3.2 Optimization of Relevant Model Parameters

[0214] This invention first designs correlation functions for the normalized time and three-dimensional position, then multiplies the four correlation functions to obtain a multiplicative overall correlation model. Common correlation functions include Gaussian functions, exponential functions, linear functions, etc. For example, four Gaussian functions can be used to establish a correlation model for a stochastic process:

[0215] (26)

[0216] in, , , and They are time The parameters to be determined for the relevant function, and their positions Parameters and positions of related functions Parameters and positions of related functions The parameters of the correlation function are derived to obtain the optimal model parameters. , , and That is, the maximum likelihood estimate, using the correlation function. and process variance Construct the objective function:

[0217] (27)

[0218] in, It is a related function The determinant, It is the maximum likelihood estimate of the process variance. The regression parameters are obtained based on the general least squares estimation method. The model parameters are iteratively updated using optimization methods such as steepest descent, conjugate gradient, and Newton's method until the objective function reaches its minimum value. The model parameters at this minimum value are defined as the optimal model parameters. The optimization process is as follows: Figure 5 As shown.

[0219] 3.3 Delay and Variance Estimation

[0220] according to Figure 4 The ionospheric grid modeling process shown utilizes a stochastic process-related model with four obtained Gaussian function parameters. Calculate the correlation function between sample IPPs Covariance Matrix And the covariance vector between all sample IPPs and the estimated IGP (or IPP). Combined with the measurement noise covariance matrix Construct the final weighted matrix. Finally, the weight coefficient vector is calculated. Estimate of the IGP (or IPP) to be estimated and nominal error variance .

[0221] To ensure consistency of the vertical ionospheric delay dimensions between the sample IPP and the IGP to be estimated, it is also necessary to... and Performing the inverse normalization operation will yield the final estimate. and standard deviation of error ,Right now:

[0222] (28)

[0223] (29)

[0224] Dimensions restored and This data can be provided to GNSS users as ionospheric grid data. Users can then use interpolation algorithms to calculate the ionospheric delay and error standard deviation at their IGP (Integrated Point of Presence), thereby meeting their positioning and integrity requirements. However, to mitigate the risk of spatial undersampling caused by ionospheric irregularities, an ionospheric threat model needs to be established to detect potential ionospheric spatial undersampling threats and to expand the standard deviation of the IGP ionospheric error estimated by spatiotemporal Kriging, thereby improving integrity.

[0225] (iv) Ionospheric threat modeling

[0226] Traditional methods for modeling ionospheric threats using undersampled data often require observational data from at least one solar cycle (11 years), resulting in an ionospheric threat model with a fitted radius. and relative center of mass A two-dimensional function, when using a fixed... and When used as a modeling constraint, the ionospheric threat model represents the maximum threat value over the entire observation period and area, and does not change with different observation times and locations. For SBAS, a larger ionospheric threat increases the overall error standard deviation (or uncertainty) of IGP, thereby increasing the overall uncertainty of the GNSS navigation signal. This results in a higher user protection level, increasing the probability of false alarms and reducing system integrity and availability. Therefore, traditional methods establish overly conservative ionospheric threat models, which cannot meet the integrity and availability requirements of real-time navigation augmentation systems like SBAS.

[0227] This invention, based on a three-dimensional annular data expulsion scheme, groups the sample IPPs into two groups: one group is used for spatiotemporal Kriging modeling, and the other group is used to detect spatial undersampling threats. This is repeated at different fitting radii. and inner ring radius The following grid modeling and threat detection operations are performed to construct threat samples and establish an undersampled threat model for IGP. The specific process is as follows: Figure 6 As shown. When there are too many IGPs, establishing an ionospheric threat model for all IGPs requires a large computational cost, which increases the computational load and complexity of the master control station and may not meet the real-time requirements. Therefore, under the premise that the dynamic changes of ionospheric characteristics are small in a short period of time (such as one hour), this invention proposes a strategy of establishing an ionospheric threat model with observation data at fixed time intervals (such as one hour), and then recalculating the overall error standard deviation of the ionosphere at the same time interval to generate a near real-time ionospheric product, that is, the overall error standard deviation at all IGPs.

[0228] 4.1 Grouping of Sample Data

[0229] To simulate potential spatial undersampling environments, all available IPP sample points within the fitting radius need to be grouped. One group is used for ionospheric grid modeling, and the other group is used to find undersampling threats (i.e., ionospheric threat models). Traditional annular data expropriation schemes are built in an east-north two-dimensional plane, using the two-dimensional coordinates of the IPP in the east and north directions under the northeast-north sky coordinate system to represent the IPP's position, and using these coordinates to calculate the distance between the IPP and IGP. This allows for grouping of available IPP sample points and building an ionospheric threat model based on spatial Kriging. However, traditional planar two-dimensional annular data expropriation schemes neglect the spatial variation characteristics of the ionospheric thin layer in the sky direction (altitude direction), and traditional spatial Kriging modeling methods also ignore the temporal variation characteristics of the ionosphere. Therefore, this invention extends and designs a spatial three-dimensional annular data expropriation scheme, such as... Figure 7 As shown, by combining the spatiotemporal Kriging method that takes into account both the temporal and spatial characteristics of the ionosphere, a near real-time ionospheric threat model is generated, providing dynamic undersampled threats for real-time navigation augmentation systems such as SBAS.

[0230] This invention directly uses the three-dimensional coordinates of the IPP in a geocentric-fixed coordinate system (such as WGS84) without performing coordinate transformation, to calculate the spatial distance between all available IPPs and the IGP to be estimated within a certain time interval. ,according to The IPP sample points are grouped as follows:

[0231] (1) Select the estimation time of the IGP to be estimated Calculate its three-dimensional coordinates ;

[0232] (2) Select a single IPP from the available IPP sample points and extract its observation time. 3D coordinates Elevation angle observed by ground observation station ;

[0233] (3) If the elevation angle observed by the ground observation station is higher than the lowest elevation angle ,Right now If yes, proceed to the next step; otherwise, select the next IPP and re-execute step (2).

[0234] (4) Based on the sampling interval and maximum number of time samples The observation time of a single IPP is checked, and if it meets the requirements... If yes, proceed to the next step; otherwise, select the next IPP and re-execute step (2).

[0235] (5) Set the fitting radius Minimum Fitting Radius ,Right now ;

[0236] (6) Calculate the spatial distance between IGP and IPP based on their three-dimensional coordinates. , and the set fitting radius Compare, if satisfied If the IPP is included in the threat IPP sample dataset, it is included in the dataset; otherwise, it is discarded.

[0237] (7) Select the next IPP and continue with steps (2) to (6) until all IPPs have been filtered;

[0238] (8) Set the inner ring radius It is 0, that is ;

[0239] (9) In the threat IPP sample dataset, select the fitting radius Inside, but located in a ring (inner ring radius is...) The outer ring radius is Width is ) outside (i.e., satisfy) or Construct modeling IPP samples from the IPPs within the annulus; select samples within the annulus that satisfy the following conditions. Construct undersampled IPP samples from the IPP of ) ;

[0240] (10) Increasing the inner ring radius ,Right now Repeat step (9) until the outer ring radius equals the fitted radius, i.e. ;

[0241] (11) The fitting radius increases ,Right now Repeat steps (6) to (10) until the fitted radius equals the maximum fitted radius, i.e. ;

[0242] (12) Select the next IGP to be estimated and repeat steps (1) to (11) until all IGPs to be estimated have been sampled.

[0243] 4.2 Undersampling Threat Modeling

[0244] When the number of modeling IPP samples in grouping step (9) of section 4.1 meets the requirements of spatiotemporal Kriging modeling, i.e., meets the requirements of ionospheric modeling step (6) in section (I), and there are undersampled IPP samples, a spatiotemporal Kriging model is established using the modeling IPP samples, and the vertical delay and error standard deviation of all undersampled IPPs are estimated, thereby detecting and estimating potential ionospheric threats. Fitting radius Need to start from the minimum fitting radius Gradually increase to the maximum fitting radius The inner ring radius needs to be gradually increased from 0 to the fitted radius. minus The undersampled threat modeling of the IGP to be estimated can only be completed after the above two rounds of iteration. When there is more than one IGP to be estimated, each IGP needs to go through two rounds of iteration: fitting radius and inner loop radius. The modeling and threat detection are carried out cyclically. In each single loop, i.e., for each fixed pair of fitting radius and inner loop radius, the spatiotemporal Kriging model established by the sample IPP is used to detect the threats that may exist in all undersampled IPPs. The detected threats are used as ionospheric threat samples. After all loops are completed, all threat samples are counted, and the maximum value among them is taken as the undersampled threat of the IGP to be estimated.

[0245] Ionospheric threat is typically defined as the standard deviation of a Gaussian distribution, which utilizes the nominal error standard deviation. and estimation error Establishment, of which and These are the model-estimated delay and the actual delay, respectively. The actual delay is the ionospheric delay caused by VTEC, as output by the GNSS observation data preprocessing module. For a single IGP to be estimated, in a single loop with a fixed fitting radius and inner loop radius, when the estimation error of the k-th undersampled IGP is greater than... times At that time, there is a threat, namely:

[0246] (30)

[0247] The above formula is equivalent to:

[0248] (31)

[0249] in, This is the confidence quantile corresponding to the ionospheric integrity risk. Therefore, the ionospheric undersampling threat sample for the kth undersampled IPP is defined as:

[0250] (32)

[0251] Based on the data expropriation scheme, undersampled threats corresponding to different fitting radii are identified by iteratively building a spatiotemporal Kriging model and searching for the maximum undersampled threat sample.

[0252] (33)

[0253] Finally, the overall error variance, including the inflated nominal error variance and the threat of spatial undersampling, is:

[0254] (34)

[0255] in, It is the nominal error inflation factor, calculated using the chi-square distribution of the number of IPPs minus the number of regression coefficients to be estimated. Calculations show that It is the probability of a false alarm. This is the probability of missed detection. Correspondingly, the error envelope of the IGP to be estimated, i.e., the grid ionospheric vertical error (GIVE), is...

[0256] (35)

[0257] in, It is the confidence quantile corresponding to the 99.9% error envelope probability.

[0258] Therefore, this invention employs the aforementioned spatiotemporal Kriging-based ionospheric grid and threat modeling method to establish a high-precision real-time ionospheric grid model based on spatiotemporal Kriging. Considering the spatial undersampling threat in four-dimensional spatiotemporal coordinates, a highly intact near-real-time ionospheric threat model based on spatiotemporal Kriging is established. The proposed data stripping scheme, combined with a sample grouping strategy, can provide a solution for establishing a real-time ionospheric threat model, and contribute to promoting the research, application, and development of the ionosphere.

[0259] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A method for ionospheric grid and threat modeling based on spatiotemporal Kriging, characterized in that, Includes the following steps: S1, Location of the calibration observation station Using data from multiple ground observation stations collected by the master control station, and employing real-time dynamic technology or precise single-point positioning technology based on carrier phase observations, the precise three-dimensional coordinates of the ground observation stations are calculated through post-processing, and the positions of the observation stations are calibrated daily or monthly. S2, Calibration Instrument Deviation Using data from multiple ground observation stations collected by the master control station and the calibrated station location coordinates, the differential code deviation of GNSS satellites and the differential code deviation of ground observation station receivers that meet constant constraints are calibrated daily using ionospheric mathematical modeling methods or precise single-point positioning technology. S3. Construct the IPP sample dataset Using ground-based observation stations distributed globally or regionally, calculate the three-dimensional coordinates, observation elevation angle, and vertical ionospheric delay of the IPP between the observation station and the GNSS core constellation and / or regional system satellites, and construct an IPP sample dataset for ionospheric modeling together with the observation time. S4. Establish an ionospheric grid model Based on elevation angle, observation time, spatial distance, and number of IPPs, the IPP sample dataset was filtered in sequence. Using the linearly normalized IPP dataset, an ionospheric grid model based on spatiotemporal Kriging was established. S5. Optimize relevant model parameters We select the correlation function corresponding to the normalized IPP observation time and three-dimensional location, design a four-dimensional spatiotemporal correlation model of the stochastic process in the spatiotemporal Kriging model, construct the objective function using the correlation function and process variance, iteratively update the correlation model parameters to minimize the objective function, and then obtain the optimal correlation model parameters. S6. Estimating Ionospheric Delay and Error Using the ionospheric grid model and related model parameters, the ionospheric delay and nominal error variance of the IGP or IPP to be estimated are estimated, and then the estimation results are inversely normalized. S7, Extended Circular Data Preemption Scheme The planar two-dimensional annular data expropriation scheme is improved into a spatial three-dimensional annular data expropriation scheme, which uses the spatial three-dimensional coordinates of IPP to group the available IPP sample points; S8. Establish an ionospheric threat model Based on the fitting radius and inner loop radius, the IPP sample points are grouped in a loop, the sample IPPs are spatiotemporally Kriging modeled, and the undersampled IPPs are threat detected. After the loop ends, all threat samples are counted, and the maximum value of the sample is taken as the undersampled threat in the ionospheric space.

2. The ionospheric grid and threat modeling method based on spatiotemporal Kriging according to claim 1, characterized in that, In step S3, the method for calculating the three-dimensional coordinates of IPP is as follows: IPP Latitude The calculation formula is: in, It's a user dimension. It is the geocentric angle between the user and the IPP sub-satellite point. and These are the azimuth and elevation angles of the GNSS satellites observed by the user. It is the approximate radius of the Earth's ellipse. It is the height of the ionosphere thin layer; The formula for calculating IPP longitude is: The user's longitude is known to be If the user dimension meets And the user's latitude, geocentric angle, and azimuth satisfy... Or the user dimension meets And the user's latitude, geocentric angle, and azimuth satisfy... The formula for calculating IPP longitude is: Vertical ionospheric delay calculations include: After cycle slip detection and repair, the length of the continuous arc segment also needs to be checked. If the length of the continuous arc segment is... Less than the predetermined minimum value If so, the TEC for that epoch will not be solved; if Then, using dual-frequency carrier phase observations without geometric combination to smooth dual-frequency code pseudorange observations without geometric combination, the code smoothing formula is: in, Dual-frequency code pseudorange has no geometric combination of observations. and These are satellites in and The pseudorange observations on the code; These are dual-frequency carrier phase observations without geometric combination. and These are satellites in and Carrier phase observations on; By combining the differential code bias of the observation station receiver and the differential code bias of the GNSS satellite in the calibration dataset, as well as the code smoothing results, the vertical total electron content (VTEC) of the ionosphere is obtained: in, and These are the frequencies of the GNSS dual-frequency signal. It is the projection function at the ionospheric intersection. These are dual-frequency code pseudorange observations without geometric combination after code smoothing. It's the speed of light. and These are the station receiver differential code bias and the GNSS satellite differential code bias in the calibration dataset, respectively. Finally, the vertical delay in the IPP sample dataset That is, the ionospheric VTEC at frequency The ionospheric delay caused by the above is calculated using the following formula: Among them, frequency It depends on the GNSS navigation signal selected by the user for positioning.

3. The ionospheric grid and threat modeling method based on spatiotemporal Kriging according to claim 1, characterized in that, In step S4, the ionospheric grid model is established, including the following steps: S41. Select the IGP to be estimated and calculate its estimated time. and three-dimensional coordinates ; S42. Select a single IPP from the IPP sample dataset and extract the observation time of that single IPP. 3D coordinates Elevation angle observed by ground observation station and vertical ionospheric delay ; S43. If the elevation angle observed by the ground observation station is higher than the lowest elevation angle ,Right now If yes, proceed to the next step; otherwise, select the next IPP and re-execute step S42. S44, According to the sampling interval and maximum number of time samples The observation time of a single IPP is checked, and if it meets the requirements... If yes, proceed to the next step; otherwise, select the next IPP and re-execute step S42. S45. Calculate the spatial distance between IGP and IPP based on their three-dimensional coordinates. , and the set fitting radius Compare, if satisfied If yes, proceed to the next step; otherwise, select the next IPP and re-execute step S42. S46. Select samples from IPP points that meet the time and distance requirements. If the maximum fitted radius... Number of IPPs within Exceeding the predetermined upper limit ,Right now Then gradually decrease the fitting radius. To the minimum fitting radius ,until ; If the number of IPPs within the maximum fitting radius is lower than the predetermined upper limit but not lower than the predetermined lower limit ,Right now If the number of IPPs within the maximum fitting radius is lower than the predetermined lower limit, then proceed to the next step; If so, the IGP to be estimated will not be modeled; S47. Normalize the IPP that satisfies steps S43, S44, S45, and S46, i.e., use... The average of individual IPPs and standard deviation Construct a normalized IPP dataset; S48. Based on the spatiotemporal Kriging method, a grid model of the ionosphere of the IGP to be estimated is established using the normalized IPP dataset.

4. The ionospheric grid and threat modeling method based on spatiotemporal Kriging according to claim 3, characterized in that, In step S47, the normalization method is as follows: in, It is the mean of the time, location, and delay of all IPP samples. It is the variance of the time, location, and delay of all sample IPPs.

5. The ionospheric grid and threat modeling method based on spatiotemporal Kriging according to claim 4, characterized in that, In step S48, a grid model of the ionosphere to be estimated (IGP) is established, including the following steps: Define the normalized vertical ionospheric delay at the i-th IGP or IPP. for: in, It is the normalization moment. It is the normalized three-dimensional position of IGP or IPP. , , , and It is the regression coefficient. Represents a random process. Indicates measurement noise; Assume a steady-state stochastic process The mean of is zero, and the non-zero correlation covariance satisfies: in, It is process variance. It has parameters Models related to stochastic processes; Assume the mean and non-zero correlation covariance of the measurement noise satisfy the following: in, The tilt factor, derived from the vertical ionospheric delay, is the tilt factor for obtaining the tilted ionospheric delay and is defined as the satellite elevation angle. The function; The variance of the measurement error is obtained by considering the dual-frequency code carrier smoothing CCL of the receiver DCB and the satellite DCB: in, It is a smooth length or arc segment count. It refers to the accuracy of pseudorange measurement. It refers to the accuracy of carrier phase measurement. It is the standard deviation of DCB at the ground station. It is the standard deviation of satellite DCB; A linear combination of the vertical delays at the IPP of the samples within the fitting radius is performed to obtain the result at the normalized time. and normalized position Spatiotemporal Kriging estimates of IGP or IPP at the location: in, It is the number of IPP samples. It is a weight coefficient vector. It is the ionospheric observation vector; the unbiased constraint for the spatiotemporal Kriging estimate is: The weighting coefficients satisfying the unbiased constraints are obtained using the Lagrange multiplier method with equality constraints, and the result is as follows: in, It is the weighted matrix of IPP for all samples, derived from the stochastic process. covariance matrix and measuring noise covariance matrix constitute: in, It is the covariance vector between all sample IPPs and the estimated IGP or IPP; assuming that the stochastic processes and measurement noise in the ionospheric model are uncorrelated, the nominal error variance of the spatiotemporal Kriging estimate is: By treating time as another independent variable besides IPP location, a spatiotemporal Kriging correlation model is constructed.

6. The ionospheric grid and threat modeling method based on spatiotemporal Kriging according to claim 5, characterized in that, In step S6, the final delay estimate at the IGP to be estimated is calculated using inverse normalization. and final error standard deviation This includes the following steps: To ensure consistency of the vertical ionospheric delay dimensions between the sample IPP and the IGP to be estimated, it is also necessary to... and Perform the inverse normalization operation to obtain the final estimate. and standard deviation of error ,Right now: 。 7. The ionospheric grid and threat modeling method based on spatiotemporal Kriging according to claim 6, characterized in that, In step S7, the planar two-dimensional annular data expropriation scheme is transformed into a spatial three-dimensional annular data expropriation scheme. The available IPP sample points are grouped using the spatial three-dimensional coordinates of the IPP, including the following steps: S701, Select the estimation time for the IGP to be estimated. Calculate its three-dimensional coordinates ; S702. Select a single IPP from the available IPP sample points and extract its observation time. 3D coordinates Elevation angle observed by ground observation station ; S703, If the elevation angle observed by the ground observation station is higher than the minimum elevation angle ,Right now If yes, proceed to the next step; otherwise, select the next IPP and re-execute step S702. S704, According to the sampling interval and maximum number of time samples The observation time of a single IPP is checked, and if it meets the requirements... If yes, proceed to the next step; otherwise, select the next IPP and re-execute step S702. S705, Set the fitting radius Minimum Fitting Radius ,Right now ; S706. Calculate the spatial distance between IGP and IPP based on their three-dimensional coordinates. , and the set fitting radius Compare, if satisfied If the IPP is included in the threat IPP sample dataset, it is included in the dataset; otherwise, it is discarded. S707. Select the next IPP and continue with steps S702-S706 until all IPPs have been filtered. S708, Set Inner Ring Radius It is 0, that is ; S709. In the threat IPP sample dataset, select the fitting radius. For IPPs located inside the ring but outside the ring, construct modeled IPP samples; select IPPs inside the ring to construct undersampled IPP samples. S710, inner ring radius increased ,Right now Repeat step S709 until the outer ring radius equals the fitted radius, i.e. ; S711, Fitting radius increased ,Right now Repeat steps S706~S710 until the fitted radius equals the maximum fitted radius, i.e. ; S712. Select the next IGP to be estimated and repeat steps S701 to S711 until all IGPs to be estimated have been sampled.

8. The ionospheric grid and threat modeling method based on spatiotemporal Kriging according to claim 7, characterized in that, In step S8, for a single IGP to be estimated, in a single loop with a fixed fitting radius and inner loop radius, when the estimation error of the k-th undersampled IGP is greater than... times At that time, there is a threat, namely: in, and These are the model-estimated latency and the actual latency, respectively. The above formula is equivalent to: in, It is the confidence quantile corresponding to the ionospheric integrity risk; the ionospheric undersampling threat sample of the kth undersampled IPP is defined as: Based on the data expropriation scheme, undersampled threats corresponding to different fitting radii are identified by iteratively building a spatiotemporal Kriging model and searching for the maximum undersampled threat sample. in, The fitting radius is... The inner ring radius is... Finally, the overall error variance, including the inflated nominal error variance and the threat of spatial undersampling, is: in, It is the nominal error inflation factor, calculated using the chi-square distribution of the number of IPPs minus the number of regression coefficients to be estimated. Calculations show that It is the probability of a false alarm. This is the probability of missed detection; correspondingly, the error envelope of the IGP to be estimated, i.e., the grid ionospheric vertical error GIVE, is: in, It is the confidence quantile corresponding to the 99.9% error envelope probability.

Images