A high-precision real-time ionosphere modeling method

Abstract

The present invention belongs to the technical field of GNSS precision data processing, and specifically discloses a method for constructing a high-precision real-time ionospheric model by combining Kalman filtering; firstly, a high-precision real-time ionospheric slant delay (Slant Electron Content, STEC) value at the ionospheric pierce point (IPP) is extracted by using the undifferenced ambiguity integer solution of the Global Navigation Satellite System (GNSS); an ionospheric prediction model is established by using a semi-parametric model to construct an ex post ionospheric model based on multi-source ionospheric observation data; hardware delays at the satellite end and the receiver end are extracted from the ex post ionospheric model constructed based on the undifferenced ambiguity integer solution; and a real-time ionospheric model is constructed by using Kalman filtering and spherical harmonic functions. The present invention provides a method for constructing a real-time ionospheric model by combining the undifferenced ambiguity integer solution and the Kalman filter, so as to improve the modeling accuracy of the real-time ionospheric model and the positioning accuracy of the GNSS single-frequency user, and provide a basis for the monitoring and early warning of space weather.

Description

Technical Field

[0001] The invention belongs to the technical field of GNSS precision data processing, and specifically relates to a high-precision ionosphere real-time modeling method. Background Art

[0002] The Global Navigation Satellite System (GNSS) has been widely used in the construction and monitoring of ionospheric models due to its advantages of high accuracy, wide coverage, and all-weather continuous monitoring. At the same time, with the rapid development of global navigation satellite systems such as GPS, GLONASS, BDS, and Galileo and the increasing number of ground-based GNSS tracking stations, the accuracy and reliability of the global ionospheric model have been continuously improved. At present, the two-dimensional ionospheric grid model constructed based on the thin ionosphere assumption has become an important data source for current GNSS navigation and positioning users and ionospheric application research. Although the post-event global ionospheric grid products (Global Ionosphere Map, GIM) provided by various ionospheric analysis centers have a high accuracy of about 2 to 8 TECU, there is a delay of at least 1 to 2 days, which not only cannot achieve real-time monitoring of the ionosphere, but also cannot provide real-time ionospheric delay correction for navigation and positioning users.

[0003] To this end, navigation and positioning users often use the following two methods to correct the real-time ionospheric delay problem. One is to use the broadcast ionospheric model, such as Klobuchar and its refined model, NeQuick model, BDS ionospheric correction model, NTCM and its refined model. However, these methods have limited accuracy in ionospheric delay correction, and can neither meet the real-time positioning needs of high-precision GNSS users nor reflect the real-time change characteristics of the ionosphere. Another method is to use the real-time ionospheric TEC model constructed by ground-based GNSS real-time observation data to correct the ionospheric delay, such as: the satellite-based augmentation system has constructed a global real-time ionospheric model, the German Aerospace Center has constructed a real-time ionospheric model for the European region, the National Oceanic and Atmospheric Administration of the United States has constructed a real-time ionospheric model for the United States, the Royal Observatory of Belgium has constructed a real-time ionospheric model for the European region, the French National Center for Space Studies has constructed a global real-time ionospheric model and broadcast it to global users, the National University of La Plata has constructed a real-time ionospheric model for South America, and the University of Catalonia and the Chinese Academy of Sciences have constructed a global ionospheric model.

[0004] In summary, the problems existing in the prior art are:

[0005] (1) Since the number of IGS monitoring stations that can provide real-time data streams is relatively small (about 180), the amount of available ionospheric data at each moment or in a short period of time is small and unevenly distributed, resulting in limited accuracy of the estimated parameters.

[0006] (2) When constructing a real-time ionospheric model in the future by combining the real-time data provided by massive navigation-augmented LEO satellites and the global regional CORS network, this method will reuse the data within the window to construct the normal equation, resulting in data redundancy and affecting the efficiency of parameter estimation.

[0007] Difficulty in solving the above technical problems: At present, the ionospheric observation data extracted using phase-smoothed pseudorange has a maximum deviation of about ±8TECU, which is poor in accuracy and affects the improvement of the accuracy of the constructed real-time ionospheric model. In addition, the number of monitoring stations that can provide real-time data is small and their distribution is limited, which limits the accuracy of the real-time ionospheric model in areas with sparse monitoring stations. Summary of the invention

[0008] In view of the shortcomings in the prior art, the present invention provides a high-precision real-time ionosphere modeling method to solve the problem that the ionosphere observation data extracted using phase smoothed pseudorange at this stage has a maximum deviation of about ±8TECU, which is poor in accuracy and affects the improvement of the accuracy of the constructed real-time ionosphere model. In addition, the number of monitoring stations that can provide real-time data is small and their distribution is limited, which limits the accuracy of the real-time ionosphere model in areas with sparse monitoring stations.

[0009] To achieve the above object, the present invention adopts the following technical solutions:

[0010] A high-precision ionosphere real-time modeling method comprises the following steps:

[0011] Step 1, calculating the ionospheric penetration point position and the zenith distance information of the penetration point according to the received real-time precise ephemeris and real-time station coordinates;

[0012] Step 2, extracting the ionospheric slant delay in real time according to the received real-time precise ephemeris, real-time observation data and real-time correction data;

[0013] Step 3: Use the undifferenced ambiguity integer solution to construct the post-event ionospheric TEC model to obtain the hardware delay parameters, and correct the hardware delay parameters in the real-time ionospheric slant delay; at the same time, use the improved ionospheric projection function to convert the ionospheric slant delay into the ionospheric vertical delay;

[0014] Step 4, using the post-event ionospheric total electron content model constructed using multi-source ionospheric data as historical data, and using a semi-parametric model to construct an ionospheric TEC prediction model;

[0015] Step 5: Use spherical harmonics to fit the ionospheric vertical delay at each puncture point and the ionospheric vertical delay provided by the ionospheric TEC prediction model; at the same time, use Kalman filtering to estimate the model parameters to be estimated every 5 minutes, and generate a real-time global ionospheric grid product.

[0016] Furthermore, the sampling rate of the precise ephemeris selected in step 1 is set to 30s; and the satellite cutoff elevation angle is set to 10°.

[0017] Further, the method for obtaining the high-precision ionospheric slant delay observation value in step 2 includes:

[0018] S1: Perform precise single-point positioning calculations of standard ionosphere-free combinations for each GNSS monitoring station around the world that can provide real-time observation data, and obtain the ionosphere-free combination ambiguity of each station and each satellite;

[0019] S2: Based on the original observation data of the GNSS monitoring station that provides real-time observation data, the MW combined observation value is obtained, and the wide lane integer ambiguity is obtained and fixed by directly taking the average and then rounding;

[0020] S3: The narrow lane real ambiguity is calculated using the ionospheric-free combined ambiguity obtained in S1 and the wide lane integer ambiguity obtained in S2, and the narrow lane fractional deviation estimation equation is established to solve the narrow lane fractional deviation value, and the narrow lane integer ambiguity is fixed using the LAMBDA method;

[0021] S4: Substitute the fixed wide lane integer ambiguity and the fixed narrow lane integer ambiguity into the geometry-free range-phase combination observation value to extract high-precision ionospheric TEC information.

[0022] Furthermore, the multi-source ionospheric observation data in step 3 includes: occultation ionospheric data, DORIS ionospheric data, top ionospheric data provided by LEO onboard GNSS receivers, Jason-2 VTEC and ground-based multi-system GNSS ionospheric data.

[0023] Further, the formula for determining the ionospheric vertical delay value at the puncture point in step 3 is:

[0024]

[0025] In the formula, H opt =506.7; α=0.9782, is the projection function adjustment coefficient; R=6378.137km, is the average radius of the earth; z is the zenith distance at the ground-based GNSS receiver; VTEC is the ionospheric vertical delay; STEC is the ionospheric slant delay.

[0026] Further, the formula of the post-event ionospheric model constructed based on multi-source data in step 3 is:

[0027]

[0028] In the formula, n and m represent the order and degree of spherical harmonics respectively, n max represents the maximum expansion order of spherical harmonics; represents the regularized Legendre function; Indicates the geomagnetic / geographic latitude at IPP; s = λ-λ 0 represents the day-solid longitude at the IPP; λ and λ 0 They represent the accuracy and solar longitude at the puncture point respectively; and represents the spherical harmonic coefficients.

[0029] Furthermore, the method for constructing the ionospheric TEC prediction model in step 4 is: using the ex post ionospheric model constructed by multi-source ionospheric observation data as historical observation data, and expressing it in the time domain in the following form:

[0030] L(t)=ψ(t)+S(t)+Δ

[0031] Where L(t) represents the ionospheric TEC value; S(t) represents the non-parametric part of the ionospheric TEC value; Δ represents the observation noise; ψ(t) is the parameter part. The above formula is expressed in the following form using Fourier series:

[0032]

[0033] In the formula, w k =2π / T k , T k is the period of the ionospheric TEC value, i.e., 1 / k (k = 1 / 1, 2, 3, 4, 5, 6) days; f represents the number of periods of the ionospheric TEC value; A k and B k They represent the parameters to be estimated respectively;

[0034] The parameters to be estimated are estimated by using the semi-parametric kernel estimation method to solve the equation. The specific steps are as follows:

[0035] First, assume that at time t p is relative to the reference time t 0 The prediction time; therefore, the kernel weight function W i (t p ) can be defined as follows:

[0036]

[0037] Wherein, i, j = 1, 2, ..., q; p = 1, 2, ..., l; q is the number of data of known ionospheric TEC values; l is the number of data of ionospheric TEC values ​​to be predicted; K(g) represents the selected kernel function;

[0038] Therefore, the residual of the ionospheric TEC value can be expressed as:

[0039]

[0040] At this time, assuming that X is known, let C = [c 1 , c 2 , …, c q ] T ,based on Then define the non-parametric part S(t p ) has the following kernel estimate:

[0041]

[0042] Let W = (W i (t p )) q×q , Y=LA 0 , then the forecast residual can be expressed as:

[0043] V=(EW)(CX-Y)

[0044] Where E represents the unit matrix;

[0045] Therefore, according to the least squares criterion, the estimate of X is It can be expressed as:

[0046]

[0047] So, we get the non-parametric part The valuation can be expressed as follows:

[0048]

[0049] So, according to The estimated value of the parameter part can be obtained; then, the estimated value of the parameter part and the estimated value of the non-parameter part are substituted back and extrapolated in the time domain to realize the prediction of the ionospheric TEC value.

[0050] Further, the parameter estimation in step 5 specifically includes the following steps:

[0051] a) One-step state prediction:

[0052]

[0053] b) Covariance one-step forecast:

[0054]

[0055] c) Calculation of filter gain:

[0056]

[0057] d) State estimation:

[0058]

[0059] e) Covariance Variance Estimation:

[0060]

[0061] In the formula, x k is the n-dimensional state vector of the ionospheric observation at time k; H k represents the n-dimensional observation matrix; v k represents the n-dimensional observation noise sequence; R k represents the covariance matrix of the observation noise; the superscript “^” represents the estimated value of the state parameter; the superscript “-” represents the prior value used in the kth iteration; P represents the symmetric positive definite variance covariance matrix, which is used to evaluate the state vector x k The estimation accuracy of ; K represents the Kalman gain matrix.

[0062] Another object of the present invention is to provide an information data processing terminal that applies the modeling method of constructing a real-time global ionosphere model based on GNSS real-time ionosphere observations in combination with Kalman filtering.

[0063] In summary, the advantages and positive effects of the present invention are as follows: the present invention provides a modeling method for constructing a real-time global ionospheric model based on GNSS real-time ionospheric observations and combined with Kalman filtering to improve the modeling accuracy of the real-time ionospheric model and the positioning accuracy of GNSS single-frequency users, and provide a basis for monitoring and early warning of space weather.

[0064] Compared with the prior art, the present invention has the following advantages:

[0065] 1) Using undifferenced ambiguity integer solutions to extract ionospheric observations in real time can improve the accuracy of the extracted ionospheric observations.

[0066] 2) The ionospheric TEC forecast model is established using a semi-parametric model for the ex post ionospheric TEC model constructed based on the undifferenced ambiguity integer solution, which can achieve constraints on data-sparse areas.

[0067] 3) Building a real-time ionospheric model based on Kalman filtering can improve the efficiency of parameter estimation and reduce data redundancy. BRIEF DESCRIPTION OF THE DRAWINGS

[0068] Figure 1 The present invention is a flow chart of a modeling method of a high-precision ionosphere real-time modeling method embodiment.

[0069] Figure 2 It is a flow chart for implementing a modeling method of an embodiment of a high-precision ionosphere real-time modeling method of the present invention. DETAILED DESCRIPTION

[0070] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0071] In view of the problems existing in the prior art, the present invention provides a modeling method for constructing a real-time global ionospheric model based on GNSS real-time ionospheric observations combined with Kalman filtering. The present invention is described in detail below in conjunction with the accompanying drawings.

[0072] The embodiment is basically as follows Figure 1 As shown, a high-precision real-time ionosphere modeling method is provided, and the modeling method for constructing a real-time global ionosphere model based on GNSS real-time ionosphere observations combined with Kalman filtering includes the following steps:

[0073] S101: Calculate the ionospheric pierce point (IPP) position and IPP zenith distance information based on the received real-time precise ephemeris and real-time station coordinates.

[0074] S102: Extracting ionospheric slant delay (Slant Electron Content, STEC) in real time according to the received real-time precise ephemeris, real-time observation data and real-time correction data.

[0075] S103: Use the undifferenced ambiguity integer solution to construct the post-event ionospheric TEC model to obtain the hardware delay parameters, and correct the hardware delay parameters in the real-time STEC. At the same time, use the modified ionospheric projection function (Modified SingleLayer Model, MSLM) to convert STEC into ionospheric vertical delay (Vertical Electron Content, VTEC).

[0076] S104: Using the post-event ionospheric total electron content (TEC) model constructed using multi-source ionospheric data as historical data, and using a semi-parametric model to construct an ionospheric TEC prediction model.

[0077] S105: Use spherical harmonics to fit the VTEC at each IPP position and the VTEC provided by the ionospheric TEC forecast model. At the same time, use Kalman filtering to estimate the model parameters every 5 minutes and generate a real-time global ionospheric grid product (GIM).

[0078] The following is combined with Figure 2 The technical solution of the present invention is further described.

[0079] like Figure 2 As shown, the modeling method for constructing a real-time global ionosphere model based on GNSS real-time ionosphere observations and combined with Kalman filtering provided in an embodiment of the present invention includes the following steps:

[0080] Step 1: Calculate the IPP position and IPP zenith distance information based on the received real-time precise ephemeris and real-time station coordinates. The sampling rate of the selected precise ephemeris is set to 30s; the satellite cutoff elevation angle is set to 10°.

[0081] Step 2: Extract STEC in real time based on the received real-time precise ephemeris, real-time observation data and real-time correction data. This is mainly done by forming wide lanes and narrow lanes in the ambiguity domain, using the overall solution method to estimate the phase fractional cycle deviation, and substituting the generated fractional cycle deviation product back into the overall solution equation to fix the integer ambiguity. The ionospheric delay observations extracted based on the GNSS non-differenced ambiguity integer solution contain non-differenced phase deviations. Based on the ionospheric modeling method, spherical harmonic functions are used to describe VTEC, separate the phase deviation part in the ionospheric observations, and obtain absolute high-precision STEC observations. The specific steps are as follows:

[0082] a) Perform PPP positioning solution of standard ionosphere-free combination for each GNSS monitoring station around the world that can provide real-time observation data, and obtain the ionosphere-free combination ambiguity of each station and each satellite;

[0083] b) obtaining MW combined observation values ​​based on the original observation data of the GNSS monitoring station providing the real-time observation data, and obtaining and fixing the wide lane integer ambiguity by directly averaging and then rounding;

[0084] c) calculating the narrow lane real ambiguity using the ionospheric-free combined ambiguity obtained in step a) and the wide lane integer ambiguity obtained in step b), establishing a narrow lane fractional deviation estimation equation, solving the narrow lane fractional deviation value, and fixing the narrow lane integer ambiguity using the LAMBDA method;

[0085] d) Substituting the fixed wide lane integer ambiguity and the fixed narrow lane integer ambiguity into the geometry-free range-phase combination observation value to extract high-precision ionospheric TEC information.

[0086] Step 3: Extract ionospheric data using occultation ionospheric data, DORIS ionospheric data, top ionospheric data provided by LEO spaceborne GNSS receivers, Jason-2 VTEC, and ground-based multi-system GNSS data.

[0087] Step 4: At the IPP, convert STEC to VTEC. The specific conversion method is as follows:

[0088]

[0089] In the formula, H opt =506.7; α=0.9782, is the projection function adjustment coefficient; R=6378.137km, is the average radius of the earth; z is the zenith distance at the ground-based GNSS receiver.

[0090] Step 5: The method for constructing the post-event ionospheric model based on multi-source data is as follows:

[0091]

[0092] In the formula, n and m represent the order and degree of spherical harmonics respectively, n max represents the maximum expansion order of spherical harmonics; represents the regularized Legendre function; Indicates the geomagnetic / geographic latitude at IPP; s = λ-λ 0 represents the day-solid longitude at the IPP; λ and λ 0 They represent the accuracy and solar longitude at IPP respectively; and represents the spherical harmonic coefficients.

[0093] Step 6: Use the post-event ionospheric model constructed from multi-source ionospheric observation data as historical observation data and express it in the time domain as follows:

[0094] L(t)=ψ(t)+S(t)+Δ

[0095] Where L(t) represents the ionospheric TEC value; S(t) represents the non-parametric part of the ionospheric TEC value; Δ represents the observation noise; ψ(t) is the parameter part, which can usually be expressed in the following form using Fourier series:

[0096]

[0097] In the formula, w k =2π / T k , T kis the period of the ionospheric TEC value, that is, 1 / k (k = 1 / 1, 2, 3, 4, 5, 6) days; f represents the number of periods of the ionospheric TEC value; A k and B k They represent the parameters to be estimated.

[0098] Then, the parameters to be estimated are estimated by solving the equation using the semi-parametric kernel estimation method. The specific steps are as follows:

[0099] First, assume that at time t p is relative to the reference time t 0 The prediction time. Therefore, the kernel weight function W i (t p ) can be defined as follows:

[0100]

[0101] Wherein, i, j = 1, 2, ..., q; p = 1, 2, ..., l; q is the number of data of known ionospheric TEC values; l is the number of data of ionospheric TEC values ​​to be predicted; K(g) represents the selected kernel function.

[0102] Therefore, the residual of the ionospheric TEC value can be expressed as:

[0103]

[0104] At this time, assuming that X is known, let C = [c 1 ,c 2 ,…,c q ] T ,based on Then define the non-parametric part S(t p ) has the following kernel estimate:

[0105]

[0106] Let W = (W i (t p )) q×q , Y=LA 0 , then the forecast residual can be expressed as:

[0107] V=(EW)(CX-Y)

[0108] Where E represents the unit matrix.

[0109] Therefore, according to the least squares criterion, the estimate of X is It can be expressed as:

[0110]

[0111] So, we get the non-parametric part The valuation can be expressed as follows:

[0112]

[0113] So, according to Then, the estimated value of the parameter part can be obtained. After that, the estimated value of the parameter part and the estimated value of the non-parameter part are substituted back and extrapolated in the time domain to realize the prediction of the ionospheric TEC value.

[0114] Step 7: Combine Kalman filtering and spherical harmonics to estimate the configuration parameters of real-time modeling. The specific estimation method is as follows:

[0115] a) One-step state prediction:

[0116]

[0117] b) Covariance one-step forecast:

[0118]

[0119] c) Calculation of filter gain:

[0120]

[0121] d) State estimation:

[0122]

[0123] e) Covariance Variance Estimation:

[0124]

[0125] In the formula, x k is the n-dimensional state vector of the ionospheric observation at time k; H k represents the n-dimensional observation matrix; v k represents the n-dimensional observation noise sequence; R k represents the covariance matrix of the observation noise; the superscript “^” represents the estimated value of the state parameter; the superscript “-” represents the prior value used in the kth iteration; P represents the symmetric positive definite variance covariance matrix, which is used to evaluate the state vector x k The estimation accuracy of ; K represents the Kalman gain matrix.

[0126] It should be noted that, in this article, the terms "comprises", "includes" or any other variations thereof are intended to cover non-exclusive inclusion, so that a process, method, article or apparatus that includes a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, article or apparatus.

[0127] Finally, it should be noted that the above embodiments are only used to illustrate the technical solution of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solution of the present invention can be modified or replaced by equivalents without departing from the purpose and scope of the technical solution of the present invention, which should be included in the scope of the claims of the present invention.

Claims

1. A high-precision real-time ionosphere modeling method, Features: The following steps are involved: Step 1, calculating the ionospheric penetration point position and the zenith distance information of the penetration point according to the received real-time precise ephemeris and real-time station coordinates; Step 2, extracting the ionospheric slant delay in real time according to the received real-time precise ephemeris, real-time observation data and real-time correction data; Step 3: Use the undifferenced ambiguity integer solution to construct the post-event ionospheric TEC model to obtain the hardware delay parameters, and correct the hardware delay parameters in the real-time ionospheric slant delay; at the same time, use the improved ionospheric projection function to convert the ionospheric slant delay into the ionospheric vertical delay; Step 4, using the post-event ionospheric total electron content model constructed using multi-source ionospheric data as historical data, and using a semi-parametric model to construct an ionospheric TEC prediction model; Step 5: Use spherical harmonics to fit the ionospheric vertical delay at each puncture point and the ionospheric vertical delay provided by the ionospheric TEC prediction model; at the same time, use Kalman filtering to estimate the model parameters to be estimated every 5 minutes, and generate a real-time global ionospheric grid product.

2. A high-precision ionosphere real-time modeling method according to claim 1, Features: The sampling rate of the precise ephemeris selected in step 1 is set to 30s; the satellite cutoff elevation angle is set to 10°.

3. A high-precision ionosphere real-time modeling method according to claim 1, Features: The method for obtaining the high-precision ionospheric slant delay observation value in step 2 includes: S1: Perform precise single-point positioning calculations of standard ionosphere-free combinations for each GNSS monitoring station around the world that can provide real-time observation data, and obtain the ionosphere-free combination ambiguity of each station and each satellite; S2: Based on the original observation data of the GNSS monitoring station that provides real-time observation data, the MW combined observation value is obtained, and the wide lane integer ambiguity is obtained and fixed by directly taking the average and then rounding; S3: The narrow lane real ambiguity is calculated using the ionospheric-free combined ambiguity obtained in S1 and the wide lane integer ambiguity obtained in S2, and the narrow lane fractional deviation estimation equation is established to solve the narrow lane fractional deviation value, and the narrow lane integer ambiguity is fixed using the LAMBDA method; S4: Substitute the fixed wide lane integer ambiguity and the fixed narrow lane integer ambiguity into the geometry-free range-phase combination observation value to extract high-precision ionospheric TEC information.

4. An information data processing terminal that applies the high-precision real-time ionosphere modeling method described in any one of claims 1 to 3 and combines it with a Kalman filter to construct a real-time global ionosphere model.

Images