A global ionosphere modeling method based on adaptive stochastic model

Abstract

The present invention relates to a global ionosphere modeling method based on an adaptive random model. First, the phase-smoothed pseudorange observation value after removing the cycle slip is used, and then the random model is adaptively adjusted using the multi-system GNSS observation value, and the observation equation is obtained by combining the spherical harmonic function model. Secondly, the time correlation constraint is added to the estimated parameter, the vertical ionospheric total electron content (VTEC) estimation is added with a non-negative constraint, and the satellite DCB is added with a zero mean constraint, and the algorithm equation is obtained and solved, and finally the global ionosphere map and the code pseudorange deviation (DCB) estimation are obtained. The present invention is designed for the observation data of the four major global navigation satellite systems (GNSS); it is suitable for the post-processing and solving of the ionosphere model of the global GNSS station; the random model determined based on the observation data has an adaptive characteristic; based on the additional constraint of the time correlation between the model parameters, the model has a higher degree of fitting; based on the additional non-negative constraint of the physical characteristics of the ionosphere itself, the model can better reflect the real ionosphere change law.

Description

Technical Field

[0001] The present invention relates to the research field of satellite positioning and navigation GNSS ionosphere model, and in particular to an ionosphere modeling method for global ionosphere monitoring and research. Background Art

[0002] With the development and popularization of electronic communication technology, people's daily life has become inseparable from electronic devices. However, due to their own characteristics, electronic communication signals will inevitably be interfered by the ionosphere during the transmission between the satellite and the ground. This has a great impact on people's production and life. Therefore, studying the impact of the ionosphere on the propagation of electronic signals has become an important research field in the study of the ionosphere and even the entire earth's space environment. Due to the dispersion characteristics of the ionosphere itself, electronic signals passing through the ionosphere will produce delays related to their own frequencies. Based on this feature, multi-frequency GNSS observation data can be used to extract ionospheric delays. With the continuous improvement of the global navigation satellite system, ionospheric monitoring technology based on GNSS observation data has become a new means. Compared with other traditional observation technologies, the GNSS method has quickly become a research hotspot at home and abroad due to its advantages such as high inversion accuracy, wide coverage, and all-weather continuous observation.

[0003] Although the research on ionospheric modeling at home and abroad is relatively mature, there are still some shortcomings. For example, in terms of random model determination, how to determine the accuracy of observations between different systems and give a reasonable and reliable prior weight ratio is still a problem worth studying. At present, various analysis centers at home and abroad generally use GPS+GLONASS dual system observation data to generate GIM. When determining the weights of observations, the differences in the accuracy of observations from different systems are not considered, resulting in poor fitting of GIM. Summary of the invention

[0004] The present invention proposes: using the single difference method between GNSS stations to eliminate satellite-related errors, using the receiver DCB as a constant, and using the peak-to-peak value (difference between the maximum and minimum values) of single difference observations of different systems as an indicator to measure the accuracy of observations, thereby determining the weighted ratio relationship between systems. The present invention uses observation data from the four major global navigation satellite systems (GPS+GLONASS+Galileo+BeiDou) and adopts the single difference method between stations to estimate the accuracy of observations of different systems, thereby obtaining a reasonable and reliable weighted ratio relationship between systems.

[0005] In summary, the present invention proposes a method for constructing a global ionospheric model based on an adaptive random model, which analyzes the accuracy of observation values ​​based on multi-system observation data and determines the weight relationship between systems, thereby establishing an accurate and reasonable random model with better adaptive characteristics.

[0006] To achieve the above object, the present invention provides a global ionosphere modeling method based on an adaptive stochastic model, comprising:

[0007] Step L1: Acquire GNSS observation data, and acquire carrier phase smoothed pseudorange observation values ​​according to the GNSS observation data;

[0008] Step L2: According to the observation values ​​described in step L1, obtain the weighted ratio relationship between the observation quantities of different systems;

[0009] Step L3: Obtain the precise satellite ephemeris and the precise position of the ground station, and obtain the GNSS observation equation based on the relationship between the observed value and the parameter to be estimated;

[0010] Step L4: According to the observation equation described in step L3, obtain the normal equation after adding constraints and solve it;

[0011] Step L5: According to the solution of the normal equation described in L4, obtain the estimated parameters and the VTEC estimation at the global ionospheric grid point;

[0012] Step L6: Based on the VTEC estimation described in L5, obtain the parameters and observation value estimates solved after adding the inequality constraints;

[0013] Step L7: According to the observed value estimation described in L6, obtain the observed value after residual test, and repeat steps L2 to L6;

[0014] Step L8: Obtain the GIM file and DCB product based on the parameter estimation described in step L7.

[0015] Optionally, in the above-mentioned global ionospheric TEC model construction method, the step L1 comprises:

[0016] Obtain GNSS observation data of the reference station of the object for the whole day and two hours before and after;

[0017] Performing cycle slip detection and elimination processing on the GNSS observation data;

[0018] Carrier phase smoothing pseudorange processing is performed on the data to obtain piecewise smoothed GNSS observation values.

[0019] Optionally, in the above-mentioned global ionospheric TEC model construction method, the step L2 includes:

[0020] Obtain the station pairs with a distance of less than 200m and containing four-system observation values ​​among the global IGS stations;

[0021] Performing antenna phase center correction on the observation data of the station pair to obtain a corrected carrier phase smoothed pseudorange observation value;

[0022] The observation values ​​are subjected to single-difference processing between measuring stations, the peak-to-peak value of the single-difference observation values ​​is calculated, and the weight ratio relationship between systems is determined.

[0023] Optionally, in the above-mentioned global ionospheric TEC model construction method, the step L3 includes:

[0024] Calculate the coordinates of the puncture point in the solar geomagnetic reference system based on the data;

[0025] Calculating the observation altitude angle according to the satellite and receiver position coordinates, and determining the weight of each observation using an altitude angle weighted formula;

[0026] According to the spherical harmonic function model and the piecewise linear model of the ionosphere, the coefficients of the parameters to be estimated at the puncture point position of each epoch are calculated;

[0027] Construct GNSS observation equations.

[0028] Optionally, in the above-mentioned global ionospheric TEC model construction method, the step L4 includes:

[0029] According to the observation equation, constructing a normal equation according to the least squares criterion;

[0030] Add zero mean constraints to satellite DCB parameters;

[0031] Add random walk constraints to the spherical harmonic coefficients.

[0032] Optionally, in the above-mentioned global ionospheric TEC model construction method, step L6 includes:

[0033] According to the VTEC estimation, obtaining the longitude and latitude of the grid point where the VTEC estimation is negative;

[0034] Based on the grid points, construct the inequality constraint equation

[0035] Iteratively solve the inequality constraint equations to obtain parameter estimates.

[0036] Optionally, in the above-mentioned global ionospheric TEC model construction method, step L7 includes:

[0037] Calculating an observation residual based on the observation estimate;

[0038] According to the Laida criterion, the residuals of the observed values ​​are tested successively;

[0039] If the test passes, the observation is retained;

[0040] If the test fails, the observation is eliminated.

[0041] In summary, the present invention adopts a spherical harmonic function model with multiple constraints, adds multiple constraints to parameters according to the global variation law of TEC, and improves the fitting degree and accuracy of the global ionosphere model.

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

[0043] The existing global ionosphere model is generally built using GPS+GLONASS dual-system observations, and does not fully utilize multi-system observation data. The global ionosphere modeling method based on GNSS four-system observations fully utilizes the observation data of each system, with a higher degree of model fitting and more realistic and reliable results.

[0044] In the existing global ionosphere modeling methods, the weight ratio relationship between systems is generally given by empirical methods, which cannot maintain long-term stability in the long-term ionosphere model construction, and has a certain impact on the reliability of long-term ionosphere monitoring. The present invention uses the observation data of the day of modeling to determine the accuracy of the observation values ​​of each system, which has strong real-time and applicability. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 It is a schematic flow chart of a global ionospheric TEC modeling method in a preferred embodiment of the present invention;

[0046] Figure 2 for Figure 1 Specific flow diagram of step S2;

[0047] Figure 3 for Figure 1 Specific flow diagram of step S3;

[0048] Figure 4 for Figure 1 Specific flow diagram of step S6; DETAILED DESCRIPTION

[0049] The following figures will be combined with schematic diagrams to describe the specific embodiments of the present invention in more detail. The advantages and features of the present invention will become more apparent from the following description and claims. It should be noted that the drawings are in very simplified form and are not in exact proportion, and are only used for convenience and clarity to assist in explaining the purpose of the embodiments of the present invention.

[0050] Related terms:

[0051] GNSS is the Global Navigation Satellite System

[0052] VTEC is the vertical ionospheric total electron content

[0053] IGS is the International GNSS Service

[0054] GIM is the Global Ionospheric Map

[0055] DCB is the differential code deviation

[0056] refer to Figure 1 In a preferred embodiment of the present invention, a global ionospheric TEC modeling method includes:

[0057] Step S1: Acquire GNSS observation data, and acquire carrier phase smoothed pseudorange observation values ​​according to the GNSS observation data;

[0058] Step S11: Obtain dual-frequency GNSS pseudorange phase observation data of globally distributed IGS reference stations.

[0059] Step S12: using the GF-MW method to detect and eliminate cycle slips from the GNSS observation data.

[0060] The GF-MW test method is as follows:

[0061] The GF combination of the carrier phase observation value obtained by subtracting the GNSS dual-frequency phase observation value is as follows:

[0062]

[0063] in, represents the GF test quantity, the subscript j represents the jth frequency, and the corresponding wavelength is λ j , frequency is f j The carrier phase observation value is a j is the frequency f j where I is the ionospheric delay on the first frequency and ε is the random error.

[0064] The GF combination is subtracted between epochs and the ionospheric delay is ignored, and the GF verification quantity is obtained as follows:

[0065] Δa GF =λ1Δa1-λ2Δa2 (2)

[0066] According to the integer characteristics of cycle slips, the minimum value of the ionospheric residual test quantity |Δa GF |≈0.2833. Therefore, comparing the absolute value of the GF test quantity with the threshold value of 0.2833, if |Δa GF |>0.2833, a cycle slip is considered to have occurred at that epoch, otherwise it is considered that no cycle slip has occurred at that epoch.

[0067] The MW combined observation values ​​are obtained by subtracting the GNSS wide lane combined observation values ​​from the pseudorange combined observation values ​​as follows:

[0068]

[0069] in, c is the speed of light, a MW =a1-a2 is the wide lane ambiguity, and the cycle slip test quantity can be expressed as:

[0070]

[0071] The mean and variance of the wide lane ambiguity are calculated epoch by epoch according to the following recursive formula:

[0072]

[0073]

[0074] Determine whether the following two widelane ambiguity-related conditions are met at the same time:

[0075]

[0076]

[0077] If yes, it is considered that no cycle slip occurred in the kth epoch;

[0078] If not, it is considered that a cycle slip occurs in the kth epoch.

[0079] Step S13: performing carrier phase smoothing pseudorange processing on the data after cycle slip removal.

[0080] Specifically, the carrier phase smoothing pseudorange formula is:

[0081]

[0082] Where i, j represent the epoch number, P4(i) = P1(i) - P2(i), L4(i) = L1(i) - L2(i), P1(i), P2(i) represent the dual-frequency pseudorange observations of the i-th epoch, L1(i), L2(i) represent the dual-frequency carrier phase observations of the i-th epoch, is the carrier phase smoothed pseudorange of the ith epoch, VTEC is the ionospheric delay on the vertical path, f1, f2 are signal frequencies, MF(z) is the projection function, DCB s , DCB r Represent the DCB parameters of the satellite and receiver respectively.

[0083] Step S2: According to the observation values ​​in step S1, the weighted ratio relationship between the observation quantities of different systems is obtained;

[0084] Specifically, refer to Figure 2 , step S2 comprises:

[0085] Step S21: Obtain the global IGS station position coordinates, and select base station pairs within 200m of each other and whose observation files contain four-system observation values ​​according to the position coordinates.

[0086] Step S22: performing antenna phase center correction on the observation data according to the base station pair in S21 to obtain a corrected carrier phase smoothed pseudorange observation value;

[0087] Step S23: Perform inter-station single difference on the carrier phase smoothed pseudorange observation values ​​respectively, and obtain the difference between the maximum and minimum single difference observation values ​​of each system as the observation value accuracy of the system.

[0088] Step S24: Obtain the weighted relationship of each system observation value according to the multi-system observation value accuracy.

[0089] Step S3: Obtain the precise satellite ephemeris and the precise position of the ground station, and obtain the GNSS observation equation according to the file.

[0090] Specifically, refer to Figure 3 , step S3 comprises:

[0091] Step S31: Establish a spherical harmonic function model to fit the expression of the total electron content (VTEC) on the vertical path:

[0092]

[0093] in, is the unnormalized associated Legendre function, are the spherical harmonic coefficients, N is the model order, λ, are the longitude and latitude of the puncture point in the solar geomagnetic coordinate system.

[0094] Step S32: In order to accurately describe the global diurnal variation characteristics of TEC, a day is divided into multiple periods, and only the spherical harmonic function parameters at the beginning and end of each period are estimated. The ionospheric delay parameters at other times in the period are represented by the linear combination of the spherical harmonic function coefficients at these two specific times. The piecewise linear representation method is as follows:

[0095]

[0096] Among them, t represents any time in the time period, T t ,T t `represent the start and end time of the period respectively, and X(*) represents the spherical harmonic function coefficient at time *.

[0097] Step S33: According to the precise satellite ephemeris and the precise coordinates of the ground station, the satellite altitude angle of the observation epoch is obtained and the altitude angle is weighted.

[0098] Preferably, the altitude angle weighted formula used is:

[0099]

[0100] Wherein, the subscript i represents the i-th satellite, θ represents its altitude angle, σ0 is the a priori unit weighted mean error, and σ is the a priori standard deviation of the current observation.

[0101] Step S34: Simultaneously establish the GNSS observation equations using equations (9), (10), and (11).

[0102] Step S4: According to the observation equation described in step S3, obtain the normal equation after adding constraints and solve it.

[0103] Specifically, step S4 includes:

[0104] Step S41: adding a random walk constraint to the spherical harmonic coefficient parameters in the observation equation, that is, assuming that the variation of the short-term ionospheric delay coefficient to be estimated over time is a random walk process, and its expression is as follows:

[0105]

[0106] Δσ is the accuracy of the spherical harmonic coefficient increment to be estimated, is a constant variance rate (TECU 2 / h), Δt is the time interval between two sets of parameters to be estimated.

[0107] Step S42: Add a zero mean constraint to the satellite DCB parameters in the observation equation, that is, consider that the sum of the satellite DCBs of each system is zero.

[0108] Preferably, the present invention uses the least square method to solve, and the present invention does not impose any limitation on this, and its normal equation is:

[0109]

[0110] Among them, N RW , N DCB are the normal matrices corresponding to the random walk and DCB zero mean constraint, R and S are the DCB coefficient matrices corresponding to the satellite and receiver, P is the observation weight matrix, and l is the observation matrix.

[0111] Step S5: According to the solution of the normal equation described in S4, the estimated parameters to be solved and the VTEC estimation at the global ionospheric grid point are obtained.

[0112] Step S6: According to the VTEC estimation described in S5, the parameters and observation value estimates solved after adding the inequality constraints are obtained.

[0113] Specifically, refer to Figure 4 , step S6 comprises:

[0114] Step S61: According to the VTEC estimation of each grid point, the longitude and latitude grid coordinates of the area where the VTEC estimation is negative are obtained.

[0115] Step S62: Add the following inequality constraint to the VTEC estimate in the observation equation:

[0116] VTEC=C coef X coef ≥0 (14)

[0117] Among them, X coef is the spherical harmonic coefficient parameter, C coef is the corresponding coefficient matrix. Since only equations can be used as constraints in the least squares estimation, a non-negative vector h≥0 is introduced to transform equation (14) into the following equation:

[0118] C coef X coef -h=0 (15)

[0119] According to the least squares criterion and combined with the Lagrange multiplier method, the following equation is obtained:

[0120]

[0121] in, C=(C coef ,C DCB ) T , C DCB They represent the DCB parameters and their corresponding coefficient matrices, P is the observation value weight matrix, q is the Lagrange coefficient, and q≥0,q T h=0. Solving equation (16), we get the expression of parameter solution:

[0122]

[0123] Where N = C T PC,

[0124] Step S63: Substitute equation (17) into equation (15), thereby transforming the least squares estimation with additional inequality constraints into a linear compensation problem:

[0125] h=Mq+w (18) where h and q are the parameters to be determined.

[0126] Step S64: traverse w and determine whether it is always positive;

[0127] If so, the solution to the equation is h = w, q = 0;

[0128] If not, find the first negative element w in w i , by linear transformation, exchanging q i and h i , obtain a new equation and repeat step S64.

[0129] Step S65: Solve the spherical harmonic function parameter estimation:

[0130] X=N -1 U+N -1 G T q (19)

[0131] Step S7: According to the observed value estimation described in S6, obtain the observed value after residual test, and repeat steps S2 to S6.

[0132] Specifically, the steps of residual test are as follows:

[0133] Step S71: Calculate the observation value residual according to the observation value estimation;

[0134] Step S72: Test the residuals of the observations in turn according to the Laida criterion. If the test passes, the observation is retained; if the test fails, the observation is discarded.

[0135] Step S8: Obtain the GIM file and the DCB file according to the parameter estimation described in step S7.

[0136] Specifically, step S8 includes:

[0137] Step S81: Divide the global area into 5°×2.5° grids, calculate the VTEC estimates and their corresponding accuracy at each grid point, and write them into the GIM file;

[0138] Step S82: extract the DCB estimates of the satellites and receivers in each system from the parameter estimates in turn, calculate the accuracy and write it into the DCB file.

[0139] In summary, the present invention adopts a spherical harmonic function model with multiple constraints, adds multiple constraints to parameters according to the global variation law of TEC, and improves the fitting degree and accuracy of the global ionosphere model.

[0140] Specifically, compared with the prior art, the present invention has the following advantages:

[0141] The existing global ionosphere model is generally built using GPS+GLONASS dual-system observations, and does not fully utilize multi-system observation data. The global ionosphere modeling method based on GNSS four-system observations fully utilizes the observation data of each system, with a higher degree of model fitting and more realistic and reliable results.

[0142] In the existing global ionosphere modeling methods, the weight ratio relationship between systems is generally given by empirical methods, which cannot maintain long-term stability in the long-term ionosphere model construction, and has a certain impact on the reliability of long-term ionosphere monitoring. The present invention uses the observation data of the day of modeling to determine the accuracy of the observation values ​​of each system, which has strong real-time and applicability.

[0143] The above is only a preferred embodiment of the present invention and does not limit the present invention in any way. Any technician in the relevant technical field, without departing from the scope of the technical solution of the present invention, makes any form of equivalent replacement or modification to the technical solution and technical content disclosed in the present invention, which does not depart from the content of the technical solution of the present invention and still falls within the protection scope of the present invention.

Claims

1. A global ionosphere modeling method based on an adaptive stochastic model, characterized in that: include: Step L1: Acquire GNSS observation data, and acquire carrier phase smoothed pseudorange observation values ​​according to the GNSS observation data; Step L2: According to the observation values ​​described in step L1, obtain the weighted ratio relationship between the observation quantities of different systems; Step L3: Obtain the precise satellite ephemeris and the precise position of the ground station, and obtain the GNSS observation equation based on the relationship between the observed value and the parameter to be estimated; Step L4: According to the observation equation described in step L3, obtain the normal equation after adding constraints and solve it; Step L5: According to the solution of the normal equation described in L4, obtain the estimated parameters and the VTEC estimation at the global ionospheric grid point; Step L6: Based on the VTEC estimation described in L5, obtain the parameters and observation value estimates solved after adding the inequality constraints; Step L7: According to the observed value estimation described in L6, obtain the observed value after residual test, and repeat steps L2 to L6; Step L8: Obtain the GIM file and DCB product according to the parameter estimation described in step L7; The step L1 comprises: Obtain GNSS observation data of the base station for the whole day and two hours before and after; Performing cycle slip detection and elimination processing on the GNSS observation data; Performing carrier phase smoothing pseudorange processing on the data to obtain piecewise smoothed GNSS observation values; The step L2 comprises: Obtain the station pairs with a distance of less than 200m and containing four-system observation values ​​among the global IGS stations; Performing antenna phase center correction on the observation data of the station pair to obtain a corrected carrier phase smoothed pseudorange observation value; The observation values ​​are subjected to single-difference processing between measuring stations, the peak-to-peak value of the single-difference observation values ​​is calculated, and the weight ratio relationship between systems is determined.

2. The global ionosphere modeling method according to claim 1, characterized in that: The step L3 comprises: Calculate the coordinates of the puncture point in the solar geomagnetic reference system based on the data; Calculating the observation altitude angle according to the satellite and receiver position coordinates, and determining the weight of each observation using an altitude angle weighted formula; According to the spherical harmonic function model and the piecewise linear model of the ionosphere, the coefficients of the parameters to be estimated at the puncture point position of each epoch are calculated; Construct GNSS observation equations.

3. The global ionosphere modeling method according to claim 1, characterized in that: The step L4 comprises: According to the observation equation, constructing a normal equation according to the least squares criterion; Add zero mean constraints to satellite DCB parameters; Add random walk constraints to the spherical harmonic coefficients.

4. The global ionosphere modeling method according to claim 1, characterized in that: The step L6 comprises: According to the VTEC estimation, obtaining the longitude and latitude of the grid point where the VTEC estimation is negative; constructing inequality constraint equations based on the grid points; Iteratively solve the inequality constraint equations to obtain parameter estimates.

5. The global ionosphere modeling method according to claim 1, characterized in that: The step L7 comprises: Calculating an observation residual based on the observation estimate; According to the Laida criterion, the residuals of the observed values ​​are tested successively; If the test passes, the observation is retained; If the test fails, the observation is eliminated.

6. The global ionosphere modeling method according to claim 2, characterized in that: The step L3 comprises: Step L31: Establish a spherical harmonic function model to fit the expression of the total electron content VTEC on the vertical path: in, is the unnormalized associated Legendre function, are the spherical harmonic coefficients, N is the model order, λ, are the longitude and latitude of the puncture point in the solar geomagnetic coordinate system; Step L32: In order to accurately describe the global diurnal variation characteristics of TEC, a day is divided into multiple periods, and only the spherical harmonic function parameters at the beginning and end of each period are estimated. The ionospheric delay parameters at other times in the period are represented by the linear combination of the spherical harmonic function coefficients at these two specific times; the piecewise linear representation method is as follows: Among them, t represents any time in the time period, T t ,T t ` respectively represent the start and end time of the period, and X(*) represents the spherical harmonic function coefficient at time *; Step L33: According to the precise satellite ephemeris and the precise coordinates of the ground station, the satellite altitude angle of the observation epoch is obtained and the altitude angle is weighted; The altitude angle weighting formula used is: Where, the subscript i represents the i-th satellite, θ represents its altitude angle, σ0 is the prior unit weighted mean error, and σ is the prior standard deviation of the current observation; Step L34: Simultaneously establish GNSS observation equations using equations (10) and (11).

7. The global ionosphere modeling method according to claim 4, characterized in that: The step L6 comprises: Step L61: According to the VTEC estimation of each grid point, obtaining the longitude and latitude grid coordinates of the area where the VTEC estimation is negative; Step L62: Add the following inequality constraint to the VTEC estimate in the observation equation: VTEC=C coef X coef ≥0 (14) Among them, X coef is the spherical harmonic coefficient parameter, C coef is the corresponding coefficient matrix; since only equations can be used as constraints in the least squares estimation, a non-negative vector h≥0 is introduced to transform equation (14) into the following equation: C coef X coef -h=0 (15) According to the least squares criterion and combined with the Lagrange multiplier method, the following equation is obtained: in, C=(C coef ,C DCB ) T , C DCB They represent the DCB parameters and their corresponding coefficient matrices, P is the observation value weight matrix, q is the Lagrange coefficient, and q≥0,q T h=0; Solve equation (16) to obtain the expression of parameter solution: Where N = C T PC, Step L63: Substitute equation (17) into equation (15), thereby transforming the least squares estimation with additional inequality constraints into a linear compensation problem: h=Mq+w (18) Among them, h and q are the parameters to be determined. Step L64: traverse w and determine whether it is always positive; If so, the solution to the equation is h = w, q = 0; If not, find the first negative element w in w i , by linear transformation, exchanging q i and h i , get the new equation, and repeat step L64; Step L65: Solve for the spherical harmonics parameter estimates: X=N -1 U+N -1 G T q (19)。

Images