A high-precision ionospheric modeling method integrating multi-source ionospheric observation data

Abstract

The present invention belongs to the technical field of GNSS precision data processing, and specifically discloses a high-precision ionospheric modeling method that integrates multi-source ionospheric observation data, uses a corresponding method to extract multi-source measured ionospheric data at the ionospheric puncture point, then uses a corresponding ionospheric projection function to convert the ionospheric oblique delay into the ionospheric vertical delay, and uses a single-layer "normalization" method to "normalize" the multi-source ionospheric observation data to have the same detection range as the ground-based GNSS ionospheric data, and then uses the Helmert variance component estimation method to estimate the ionospheric modeling parameters of the multi-source ionospheric observation data after "normalization" and generate an ionospheric grid product. The present invention provides a method for integrating multi-system GNSS and multi-source ionospheric observation data to construct a high-precision ionospheric model, so as to improve the accuracy of the global ionospheric model in areas with sparse measurement stations and the positioning accuracy of GNSS single-frequency users.

Description

Technical Field

[0001] The present invention belongs to the technical field of GNSS precision data processing, and specifically relates to a high-precision ionospheric modeling method integrating multi-source ionospheric observation data. Background Art

[0002] At present, ground-based GNSS monitoring stations are unevenly distributed around the world, resulting in less ionospheric data in the ocean and areas with harsh natural conditions, which restricts the accuracy and reliability of the model in the corresponding areas. To this end, scholars use different methods to improve the accuracy and reliability of the ionospheric TEC model.

[0003] The first method is to use the empirical ionospheric model as the background field constraint, such as the International Reference Ionosphere (IRI) model, the final ionospheric grid product (Global IonosphereMap, GIM) provided by the International GNSS Service (IGS) the day before, the Klobuchar model and the NeQuick model. However, limited by the accuracy of the empirical model, the accuracy of the ionospheric TEC model constructed by this method is limited.

[0004] The second method is to construct an ionospheric TEC model by combining multi-source ionospheric data. Since low Earth Orbit (LEO) scientific satellites carry different types of ionospheric detection equipment, the ionospheric data they provide can fill the gap in ionospheric data in the ocean and areas with harsh natural conditions to a certain extent. Therefore, the ionospheric model can usually be constructed using ground-based GNSS ionospheric data, Jason ocean altimetry satellite ionospheric data, occultation ionospheric data, Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS) observation data, very long baseline interferometry data, and ship-borne GNSS and satellite-borne GNSS data, which improves the accuracy of the constructed ionospheric model in the ocean and areas with harsh natural conditions to a certain extent.

[0005] In summary, the problem with the prior art is that the existing research regards the TEC value of the undetected part of the ionosphere as a constant to achieve compensation. However, if the multi-source ionosphere data used is the LEO bottom ionosphere data, the difference between it and the ground-based GNSS ionosphere detection range is the ionosphere from the LEO satellite orbit altitude to the GNSS orbit altitude, that is, the plasma layer; if the multi-source ionosphere observation data used is the LEO top ionosphere data, the difference between it and the ground-based GNSS ionosphere detection range is the ionosphere from the bottom of the ionosphere to the LEO satellite orbit altitude. In the former case, since the plasma layer can extend to an altitude of 3 to 5 earth radii, its proportion at different times is different: it accounts for about 10% during the day and up to 60% at night; in the latter case, the undetected part includes the ionosphere D layer, E layer and F layer, and the proportion of the electron content of the undetected part is at least 40% regardless of daytime or nighttime. Therefore, the TEC value of the undetected part of the ionosphere cannot be ignored in terms of magnitude or variation.

[0006] Difficulty in solving the above technical problems: Due to the differences in the detection range of the ionosphere by different types of ionosphere observation data, simply fusing different types of ionosphere observation data during modeling will only result in limited improvement in the accuracy of the constructed ionosphere model. In addition, different types of ionosphere observation data are different, and only using the least squares method to estimate the modeling parameters will affect the accuracy of the parameters to be estimated. Summary of the invention

[0007] In view of the deficiencies in the prior art, the present invention provides a high-precision ionospheric modeling method that fuses multi-source ionospheric observation data to solve the problem that different types of ionospheric observation data have different detection ranges for the ionosphere, resulting in only simple fusion of different types of ionospheric observation data during modeling, which will result in limited improvement in the accuracy of the constructed ionospheric model. In addition, different types of ionospheric observation data are different, and only using the least squares method to estimate the modeling parameters will affect the accuracy of the parameters to be estimated.

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

[0009] A high-precision ionospheric modeling method integrating multi-source ionospheric observation data comprises the following steps:

[0010] In the first step, the ground-based GNSS observation data is preprocessed, and then the ionospheric slant delay in the line of sight is extracted;

[0011] The second step is to obtain the puncture point position and the corresponding zenith distance according to the satellite ephemeris and the ground-based GNSS station coordinates;

[0012] The third step is to obtain the puncture point position and zenith distance based on the multi-source observation data and extract the multi-source ionospheric slant delay; then, according to the type of ionospheric data, use the corresponding ionospheric projection function to convert it into ionospheric vertical delay;

[0013] The fourth step is to use the single-layer “normalization” method to “normalize” the multi-source ionospheric data so that it has the same detection range as the ground-based GNSS ionospheric data;

[0014] The fifth step is to construct the observation equation based on the converted ionospheric slant delay data of various types and the obtained puncture point position information, and then construct the normal equation based on the observation equation;

[0015] In the sixth step, the Helmert variance component estimation method is used to estimate the parameters to be estimated, and the differential code bias file and the global ionospheric grid file are generated.

[0016] Furthermore, the first step of data preprocessing includes: using the TurboEdit method to detect and repair the phase cycle slip; during processing, the satellite cutoff elevation angle is 10°.

[0017] Further, the method for obtaining the ionospheric slant delay observation value in the first step includes:

[0018] a) Combining the pseudorange and phase observations into geometric range combinations, we get the following formula:

[0019]

[0020] In the formula, and Denotes the geometry-free range combined observation values ​​of pseudorange and phase respectively; DPB r and DPB s Represent the differential phase deviation of the receiver and the satellite respectively; γ i =(f0 / f i ) 2 , where f0 = 10.23 MHz; and Respectively represent the pseudorange and phase observation values ​​at frequency i; represents the geometric distance from the receiver r to the satellite s; c represents the propagation speed of electromagnetic waves in a vacuum; dt r and dt s denote the clock difference between the receiver r and the satellite s respectively; Tro denotes the tropospheric delay on the observation path from the receiver r to the satellite s; Indicates the ionospheric delay on the observation path from receiver r to satellite s, unit: meter; B r,i and are the pseudorange hardware delays of the receiver and the satellite at frequency i; b r,i and Respectively represent the phase hardware delay of the receiver and the satellite at frequency i; λ i represents the wavelength of frequency i; represents the integer ambiguity of the phase observation at frequency i; represents the multipath effect of the jth satellite on frequency i; ε P and Φ denote the observation noise of pseudorange and phase observations respectively;

[0021] b) Since there is no cycle slip in an observation arc, the ambiguity values ​​of each epoch in the arc remain unchanged and DCB and DPB remain stable within a day, which can be estimated as constants; therefore, in a continuous arc, the difference between the phase geometry-free range observation value and the pseudorange geometry-free range observation value can be expressed as:

[0022]

[0023] Substituting it into the phase-free geometric distance combined observation value, we can get the following formula:

[0024]

[0025] Where T represents the time length of the continuous arc segment; <> represents the averaging operation; thus, is the ionospheric delay extracted using the phase smoothing pseudorange method, unit: m.

[0026] Furthermore, the sampling rate of the observation data selected in the second step is 30s; the height of the ionosphere single layer model used is 450km.

[0027] Furthermore, the multi-source ionospheric observation data of the third step includes: occultation ionospheric data, DORIS ionospheric data, top ionospheric data provided by LEO onboard GNSS receiver, Jason-2 VTEC and ground-based multi-system GNSS ionospheric data.

[0028] Further, the method for extracting multi-source ionospheric observation data in the third step includes:

[0029] a) Occultation ionospheric observation data: provided by the "ionPrf" product, which provides satellite upper and lower ionospheric observation data to obtain TEC values ​​for the entire ionospheric path:

[0030] VTEC=VTEC bellow +VTEC up

[0031] Where VTEC bellow Indicates the VTEC value below the COSMIC satellite; VTEC uprepresents the VTEC value above the COSMIC satellite, which can be obtained by extrapolation;

[0032] b) DORIS ionospheric observation data: The calculation steps include:

[0033] 1) Based on the position of the DORIS ground-based monitoring station and the LEO satellite position, the puncture point position and the zenith distance at the puncture point of the DORIS observation data are calculated, and the VTEC value at the puncture point is calculated using the ionospheric TEC model constructed by the ground-based GNSS;

[0034] 2) Calculate the mean deviation between the DORIS-containing STEC value and the STEC value converted by GIM in each continuous arc segment;

[0035] 3) Correct the obtained average deviation within the arc segment to the corresponding arc segment to obtain the corrected STEC value;

[0036] c) Jason2 VTEC is obtained based on the pseudo-range ionospheric delay correction of the Ku-band frequency of the satellite altimeter. The specific formula is as follows:

[0037]

[0038] Where, dR Ku is the Ku-band ionospheric distance correction number; f Ku is the Ku-band frequency in GHz. Due to the large noise of altimeter observation, a 15-25s window can usually be used in combination with the median filter method to extract ionospheric data. The accuracy of the smoothed ionospheric data is better than 1TECU.

[0039] d) The LEO top ionosphere observation data and ground-based GNSS ionosphere observation data can be extracted using the phase smoothed pseudorange method.

[0040] Further, the different types of ionospheric projection functions in the third step mainly include:

[0041] a) The LEO top ionosphere projection function can be expressed as follows:

[0042]

[0043] In the formula, R shell =R+h shell represents the distance between the upper boundary of the equivalent thick layer and the center of the earth, h shell is the height of the upper boundary of the thick layer; R orbit represents the orbital altitude of the LEO satellite, R = 6378.137 km, which is the average radius of the earth; z is the zenith angle between the onboard GNSS receiver and the GNSS satellite direction at the onboard GNSS receiver;

[0044] b) The ionospheric projection function used by other types of ionospheric observation data can be expressed as follows:

[0045]

[0046] In the formula, H opt =506.7; α=0.9782, which is the projection function adjustment coefficient; R=6378.137km, which is the average radius of the earth; H represents the height of the ionosphere single layer model, which is usually a fixed value between 350 and 450km; z is the zenith distance at the ground-based GNSS receiver.

[0047] Furthermore, the single-layer "normalization" method used in the fourth step is specifically as follows:

[0048] a) Use the NeQuick2 model to calculate the VTEC values ​​at the corresponding ionospheric puncture points in different detection ranges. The calculation formula is as follows:

[0049]

[0050] In the formula, H1 and H2 represent the detection ranges corresponding to different types of ionospheric data; h1 and h2 represent the detection ranges after “normalization”; represents the VTEC values ​​calculated at the puncture point for different types of ionospheric data using the NeQuick2 model; represents the “normalized” VTEC value calculated by the NeQuick2 model at the puncture point; ∫N e (s)ds represents the integral of the electron density on the corresponding observation path.

[0051] b) Calculate the scale factor between different detection ranges at time t, as shown in the following formula:

[0052]

[0053] c) The STEC extracted from different line of sight directions is converted into VTEC using the corresponding projection function. Assume that the observation data between the ground-based GNSS receiver and the GNSS satellite can be expressed as STEC GNSS The bottom ionosphere data provided by the Low Earth Orbit (LEO) satellite can be expressed as STEC LEO-BTM , the top ionosphere data provided by LEO satellites are expressed as STEC LEO-UP , which can be expressed as VTEC after projection transformation GNSS 、VTEC LEO-BTM and VTEC LEO-UP ; The specific conversion formula is as follows:

[0054]

[0055] Where mf1 and mf2 represent the ground-based ionosphere projection function and the satellite-based ionosphere projection function, respectively, corresponding to the MSLM projection function and the F&K projection function; it is worth noting that due to the LEO satellite orbit altitude h m Above F2, so the cutoff elevation angle of the observation data of the LEO satellite onboard GNSS receiver should be as large as possible to weaken the projection error;

[0056] d) According to the calculated scale factor, different types of ionospheric data are “normalized”; by “normalization”, the ionospheric data (VTEC) on the corresponding observation path can be normalized. Obs (t)) is converted into the ionospheric data VTEC on the path after “normalization” nomalized (t), the specific calculation formula is as follows:

[0057]

[0058] Further, the form of the method equation constructed in the fifth step is as follows:

[0059]

[0060] Wherein, the subscript GNSS represents the ground-based GNSS ionospheric data; the subscript Topside represents the top ionospheric data of the LEO onboard GNSS receiver; the subscript Altimeter represents the Jason-2 VTEC; the subscript Occultation represents the LEO occultation ionospheric data; the subscript DORIS represents the DORIS ionospheric data; N represents the normal equation matrix composed of different types of ionospheric data; B represents the construction matrix of different types of ionospheric data; and P represents the weights corresponding to different types of ionospheric data.

[0061] Furthermore, the specific steps of the sixth step of estimating the parameters to be estimated by the Helmert variance component estimation method are as follows:

[0062] a) Estimating the prior weights of different types of ionospheric data and determining the initial values ​​of the weights of different types of ionospheric data: P1, P2, ..., P n ;

[0063] b) Perform the first adjustment on all observations to obtain the different ionospheric data

[0064] c) Estimation of the first variance component is performed as follows:

[0065]

[0066] Obtaining a first estimate of the unit weight variance for different types of ionospheric data Then determine the weight according to the following formula:

[0067]

[0068] In the formula, c is A value in ;

[0069] The adjustment and variance component estimation are repeated until the unit weight variances of different ionospheric data are equal.

[0070] Another object of the present invention is to provide an information data processing terminal that applies the method of fusing multi-system GNSS and multi-source ionospheric observation data to construct a high-precision ionospheric model.

[0071] In summary, the advantages and positive effects of the present invention are as follows: the present invention provides a method for constructing a high-precision ionospheric model by integrating multi-system GNSS and multi-source ionospheric observation data, so as to improve the ionospheric modeling accuracy and the positioning accuracy of GNSS single-frequency users.

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

[0073] 1) Use the single-layer “normalization” method to “normalize” different types of ionospheric observation data so that they have the same detection range as the ground-based GNSS ionospheric data.

[0074] 2) The integration of multi-source ionospheric observation data can effectively improve the coverage of ionospheric observation data over the ocean.

[0075] 3) The Helmert variance component estimation method is used to estimate the parameters to be estimated, and the weights of different types of observation data are adaptively determined. BRIEF DESCRIPTION OF THE DRAWINGS

[0076] Figure 1 The present invention is a method flow chart of an embodiment of a high-precision ionospheric modeling method for fusing multi-source ionospheric observation data.

[0077] Figure 2 The present invention is a flowchart of an implementation of a high-precision ionospheric modeling method for fusing multi-source ionospheric observation data.

[0078] Figure 3 The present invention is a flowchart for implementing DORIS observation data correction provided by an embodiment of a high-precision ionospheric modeling method for fusing multi-source ionospheric observation data.

[0079] Figure 4 An embodiment of the present invention provides a schematic diagram of the detection range of different types of ionospheric data in accordance with a high-precision ionospheric modeling method integrating multi-source ionospheric observation data. DETAILED DESCRIPTION

[0080] 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.

[0081] In view of the problems existing in the prior art, the present invention provides a high-precision ionospheric modeling method integrating multi-source ionospheric observation data. The present invention is described in detail below in conjunction with the accompanying drawings.

[0082] The embodiment is basically as follows Figure 1 As shown, a high-precision ionospheric modeling method integrating multi-source ionospheric observation data, the following steps:

[0083] S101: Preprocess the ground-based GNSS observation data and then extract the ionospheric slant delay (Slant Total Electron Content, STEC) in the line of sight direction;

[0084] S102: Obtain the ionospheric pierce point (IPP) position and the corresponding zenith distance according to the satellite ephemeris and the ground-based GNSS station coordinates;

[0085] S103: Obtain the puncture point position and zenith distance based on multi-source observation data and extract multi-source STEC. Then, according to the type of ionospheric data, use the corresponding ionospheric projection function to convert it into ionospheric vertical delay (Vertical Total Electron Content, VTEC);

[0086] S104: Use the single-layer “normalization” method to “normalize” the multi-source ionospheric data so that it has the same detection range as the ground-based GNSS ionospheric data;

[0087] S105: constructing an observation equation based on various types of VTEC data obtained by conversion and combining the obtained puncture point position information, and constructing a normal equation based on the observation equation;

[0088] S106: Estimate the parameters to be estimated using the Helmert variance component estimation method, and generate a Differential Code Bias (DCB) file and a Global Ionosphere Map (GIM) file.

[0089] The technical solution of the present invention is further described below in conjunction with the accompanying drawings.

[0090] like Figure 2As shown, the high-precision ionospheric modeling method for fusing multi-source ionospheric observation data provided by the embodiment of the present invention includes the following steps:

[0091] Step 1: Download and organize different types of ionospheric observation data and preprocess them. The selected ground-based GNSS stations should be as evenly distributed as possible around the world and in large numbers to improve the accuracy of ionospheric modeling. Data preprocessing includes: using the TurboEdit method to detect and repair phase cycle slips; during processing, the satellite cutoff elevation angle is 10°.

[0092] Step 2: Use the phase smoothing pseudorange method to extract ionospheric observation data. The specific steps are as follows:

[0093] a) Combining the pseudorange and phase observations into geometric range combinations, we get the following formula:

[0094]

[0095] In the formula, and Denotes the geometry-free range combined observation values ​​of pseudorange and phase respectively; DPB r and DPB s Represent the differential phase deviation of the receiver and the satellite respectively; γ i =(f0 / f i ) 2 , where f0 = 10.23 MHz; and Respectively represent the pseudorange and phase observation values ​​at frequency i; represents the geometric distance from the receiver r to the satellite s; c represents the propagation speed of electromagnetic waves in a vacuum; dt r and dt s denote the clock difference between the receiver r and the satellite s respectively; Tro denotes the tropospheric delay on the observation path from the receiver r to the satellite s; Indicates the ionospheric delay on the observation path from receiver r to satellite s, unit: meter; B r,i and are the pseudorange hardware delays of the receiver and the satellite at frequency i; b r,i and Respectively represent the phase hardware delay of the receiver and the satellite at frequency i; λ i represents the wavelength of frequency i; represents the integer ambiguity of the phase observation at frequency i; represents the multipath effect of the jth satellite on frequency i; ε P and Φ denote the observation noise of pseudorange and phase observations respectively.

[0096] b) Since there is no cycle slip in an observation arc, the ambiguity values ​​of each epoch in the arc remain unchanged and DCB and DPB remain stable within a day, which can be estimated as constants. Therefore, in a continuous arc, the difference between the phase geometry-free range observation value and the pseudorange geometry-free range observation value can be expressed as:

[0097]

[0098] Substituting it into the phase-free geometric distance combined observation value, we can get the following formula:

[0099]

[0100] Where T represents the time length of the continuous arc segment; <> represents the averaging operation. Therefore, it is the ionospheric delay extracted using the phase smoothing pseudorange method, in meters.

[0101] Step 3: Calculate the IPP position and zenith distance information of the ground-based GNSS observation data. When calculating, the selected observation data sampling rate is 30s; the ionosphere single layer model height used is 450km.

[0102] Step 4: Calculate multi-source ionospheric observation data. The specific method mainly includes:

[0103] a) Occultation ionospheric observation data. Mainly provided by the "ionPrf" product, this data product provides satellite upper and lower ionospheric observation data to obtain the TEC value of the entire ionospheric path:

[0104] VTEC=VTEC bellow +VTEC up

[0105] Where VTEC bellow Indicates the VTEC value below the COSMIC satellite; VTEC up represents the VTEC value above the COSMIC satellite, which can be obtained by extrapolation.

[0106] b) DORIS ionospheric observation data. The calculation process is as follows: Figure 3 As shown, the specific steps mainly include:

[0107] 1) Based on the position of the DORIS ground-based monitoring station and the LEO satellite position, the puncture point position and the zenith distance at the puncture point of the DORIS observation data are calculated, and the VTEC value at the puncture point is calculated using the ionospheric TEC model constructed by ground-based GNSS.

[0108] 2) Calculate the mean deviation between the DORIS-derived STEC value and the STEC value converted by GIM in each continuous arc segment.

[0109] 3) Correct the obtained average deviation within the arc segment to the corresponding arc segment to obtain the corrected STEC value.

[0110] c) Usually, Jason2 VTEC can be obtained based on the pseudo-range ionospheric delay correction of the Ku-band frequency of the satellite altimeter. The specific formula is as follows:

[0111]

[0112] Where, dR Ku is the Ku-band ionospheric distance correction number; f Ku is the Ku-band frequency in GHz. Due to the large noise of altimeter observations, a 15-25s window can usually be used in combination with the median filter method to extract ionospheric data. The accuracy of the smoothed ionospheric data is better than 1TECU.

[0113] d) The LEO top ionosphere observation data and ground-based GNSS ionosphere observation data can be extracted using the phase smoothed pseudorange method.

[0114] Step 5: Use the corresponding ionospheric projection function to convert different types of STEC into VTEC. The specific ionospheric projection function mainly includes:

[0115] a) The LEO top ionosphere projection function can be expressed as follows:

[0116]

[0117] In the formula, R shell =R+h shell represents the distance between the upper boundary of the equivalent thick layer and the center of the earth, h shell is the height of the upper boundary of the thick layer; R orbit represents the orbital altitude of the LEO satellite, R = 6378.137 km, which is the average radius of the earth; z is the zenith angle between the onboard GNSS receiver and the GNSS satellite direction at the onboard GNSS receiver.

[0118] b) The ionospheric projection function used by other types of ionospheric observation data can be expressed as follows:

[0119]

[0120] In the formula, H opt =506.7; α=0.9782, which is the projection function adjustment coefficient; R=6378.137km, which is the average radius of the earth; H represents the height of the ionosphere single layer model, which is usually a fixed value between 350 and 450km; z is the zenith distance at the ground-based GNSS receiver.

[0121] Step 6: Figure 4 As shown, the specific process of using the single-layer "normalization" method is as follows:

[0122] a) Use the NeQuick2 model to calculate the VTEC values ​​at the corresponding ionospheric puncture points in different detection ranges. The calculation formula is as follows:

[0123]

[0124] In the formula, H1 and H2 represent the detection ranges corresponding to different types of ionospheric data; h1 and h2 represent the detection ranges after “normalization”; represents the VTEC values ​​calculated at the puncture point for different types of ionospheric data using the NeQuick2 model; represents the “normalized” VTEC value calculated by the NeQuick2 model at the puncture point; ∫N e (s)ds represents the integral of the electron density on the corresponding observation path.

[0125] b) Calculate the scale factor between different detection ranges at time t, as shown in the following formula:

[0126]

[0127] c) The STEC extracted from different line of sight directions is converted into VTEC using the corresponding projection function. Assume that the observation data between the ground-based GNSS receiver and the GNSS satellite can be expressed as STEC GNSS The bottom ionosphere data provided by the Low Earth Orbit (LEO) satellite can be expressed as STEC LEO-BTM , the top ionosphere data provided by LEO satellites are expressed as STEC LEO-UP , which can be expressed as VTEC after projection transformation GNSS 、VTEC LEO-BTM and VTEC LEO-UP . The specific conversion formula is as follows:

[0128]

[0129] Where mf1 and mf2 represent the ground-based ionosphere projection function and the satellite-based ionosphere projection function, respectively, corresponding to the MSLM projection function and the F&K projection function. It is worth noting that since the LEO satellite orbit altitude is h m Therefore, the cutoff elevation angle of the observation data of the LEO satellite onboard GNSS receiver should be as large as possible to weaken the projection error.

[0130] d) According to the calculated scale factor, different types of ionospheric data are “normalized”. Through “normalization”, the ionospheric data (VTEC) on the corresponding observation path can be normalized. Obs (t)) is converted into the ionospheric data VTEC on the path after “normalization” nomalized (t), the specific calculation formula is as follows:

[0131]

[0132] Step 7, the normal equation obtained when constructing the global ionosphere model is as follows:

[0133]

[0134] In the formula, the subscript GNSS represents the ground-based GNSS ionospheric data; the subscript Topside represents the top ionospheric data of the LEO onboard GNSS receiver; the subscript Altimeter represents Jason-2 VTEC; the subscript Occultation represents the LEO occultation ionospheric data; the subscript DORIS represents the DORIS ionospheric data; N represents the normal equation matrix composed of different types of ionospheric data; B represents the construction matrix of different types of ionospheric data; P represents the weights corresponding to different types of ionospheric data. Subsequently, the Helmert variance component estimation method is used to estimate the parameters to be estimated. The specific steps are as follows:

[0135] a) Estimating the prior weights of different types of ionospheric data and determining the initial values ​​of the weights of different types of ionospheric data: P1, P2, ..., P n ;

[0136] b) Perform the first adjustment on all observations to obtain the different ionospheric data

[0137] c) Estimation of the first variance component is performed as follows

[0138]

[0139] Obtaining a first estimate of the unit weight variance for different types of ionospheric data Then determine the weight according to the following formula:

[0140]

[0141] In the formula, c is A value in .

[0142] The adjustment and variance component estimation are repeated until the unit weight variances of different ionospheric data are equal.

[0143] 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.

[0144] 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 ionospheric modeling method integrating multi-source ionospheric observation data, characterized in that: The following steps are involved: In the first step, the ground-based GNSS observation data is preprocessed, and then the ionospheric slant delay in the line of sight is extracted; The method for obtaining the ionospheric slant delay observation value includes: a) Combining the pseudorange and phase observations into geometric range combinations, we get the following formula: In the formula, and Denotes the geometry-free range combined observation values ​​of pseudorange and phase respectively; DPB r and DPB s Represent the differential phase deviation of the receiver and the satellite respectively; γ i =(f0 / f i ) 2 , where f0 = 10.23 MHz; and Respectively represent the pseudorange and phase observation values ​​at frequency i; represents the geometric distance from the receiver r to the satellite s; c represents the propagation speed of electromagnetic waves in a vacuum; Indicates the ionospheric delay on the observation path from receiver r to satellite s, unit: meter; B r,i and are the pseudorange hardware delays of the receiver and the satellite at frequency i; b r,i and Respectively represent the phase hardware delay of the receiver and the satellite at frequency i; λ i represents the wavelength of frequency i; represents the integer ambiguity of the phase observation at frequency i; represents the multipath effect of the jth satellite on frequency i; ε P and Φ denote the observation noise of pseudorange and phase observations respectively; b) Since a cycle slip does not occur in an observation arc, the ambiguity values ​​of each epoch in the arc remain unchanged and DCB and DPB remain stable within a day, and are estimated as constants; therefore, in a continuous arc, the difference between the phase geometry-free range observation value and the pseudorange geometry-free range observation value is expressed as: Substituting it into the phase-free geometric distance combined observation value, we get the following formula: Where T represents the time length of the continuous arc segment; <> represents the averaging operation; thus, it is the ionospheric delay extracted using the phase smoothing pseudorange method, unit: m; The second step is to obtain the puncture point position and the corresponding zenith distance according to the satellite ephemeris and the ground-based GNSS station coordinates; The third step is to obtain the puncture point position and zenith distance based on the multi-source observation data and extract the multi-source ionospheric slant delay; then, according to the type of ionospheric data, use the corresponding ionospheric projection function to convert it into ionospheric vertical delay; The third step is to use different types of ionospheric projection functions including: a) The projection function of the LEO top ionosphere is expressed as follows: In the formula, R shell =R+h shell represents the distance between the upper boundary of the equivalent thick layer and the center of the earth, h shell is the height of the upper boundary of the thick layer; R orbit represents the orbital altitude of the LEO satellite, R = 6378.137 km, which is the average radius of the earth; z is the zenith angle between the onboard GNSS receiver and the GNSS satellite direction at the onboard GNSS receiver; b) The ionospheric projection function used by other types of ionospheric observation data is expressed as follows: In the formula, H opt =506.7; α=0.9782, which is the projection function adjustment coefficient; R=6378.137km, which is the average radius of the earth; H represents the height of the ionosphere single layer model, which is a fixed value between 350 and 450km; z is the zenith distance at the ground-based GNSS receiver; The fourth step is to use the single-layer normalization method to normalize the multi-source ionospheric data so that it has the same detection range as the ground-based GNSS ionospheric data; The single-layer normalization process is as follows: a) Use the NeQuick2 model to calculate the VTEC values ​​at the corresponding ionospheric puncture points in different detection ranges. The calculation formula is as follows: Where H1 and H2 represent the detection ranges corresponding to different types of ionospheric data; h1 and h2 represent the normalized detection ranges; represents the VTEC values ​​calculated at the puncture point for different types of ionospheric data using the NeQuick2 model; represents the normalized VTEC value calculated by the NeQuick2 model at the puncture point; ∫N e (s)ds represents the integral of the electron density on the corresponding observation path; b) Calculate the scale factor between different detection ranges at time t, as shown in the following formula: c) The STEC extracted from different sight lines is converted into VTEC using the corresponding projection function; the observation data between the ground-based GNSS receiver and the GNSS satellite is represented as STEC GNSS , the bottom ionosphere data provided by low-orbit satellites are expressed as STEC LEO-BTM , the top ionosphere data provided by LEO satellites are expressed as STEC LEO-UP , after projection transformation, they are expressed as VTEC GNSS 、VTEC LEO-BTM and VTEC LEO-UP ; The conversion formula is as follows: Where mf1 and mf2 represent the ground-based ionospheric projection function and the satellite-based ionospheric projection function, corresponding to the MSLM projection function and the F&K projection function, respectively; d) According to the calculated scale factor, different types of ionospheric data are normalized; the ionospheric data (VTEC) on the corresponding observation path are normalized. Obs (t)) is converted into the ionospheric data VTEC on the normalized path nomalized (t), the specific calculation formula is as follows: The fifth step is to construct the observation equation based on the converted ionospheric slant delay data of various types and the obtained puncture point position information, and then construct the normal equation based on the observation equation; The form of the law equation constructed in the fifth step is as follows: Wherein, the subscript GNSS represents the ground-based GNSS ionospheric data; the subscript Topside represents the top ionospheric data of the LEO onboard GNSS receiver; the subscript Altimeter represents Jason-2 VTEC; the subscript Occultation represents the LEO occultation ionospheric data; the subscript DORIS represents the DORIS ionospheric data; N represents the normal equation matrix composed of different types of ionospheric data; B represents the construction matrix of different types of ionospheric data; P represents the weights corresponding to different types of ionospheric data; Step 6: Use the Helmert variance component estimation method to estimate the parameters to be estimated, and generate differential code bias files and global ionospheric grid files; The specific steps of the sixth step Helmert variance component estimation method to estimate the parameters to be estimated are as follows: a) Estimating the prior weights of different types of ionospheric data and determining the initial values ​​of the weights of different types of ionospheric data: P1, P2, ..., P n ; b) Perform the first adjustment on all observations to obtain V for different ionospheric data i T P i V i ; c) Estimation of the first variance component is performed as follows: Obtaining a first estimate of the unit weight variance for different types of ionospheric data Then determine the weight according to the following formula: In the formula, c is A value in ; The adjustment and variance component estimation are repeated until the unit weight variances of different ionospheric data are equal.

2. The high-precision ionospheric modeling method for integrating multi-source ionospheric observation data according to claim 1, characterized in that: The first step of data preprocessing includes: using the TurboEdit method to detect and repair the phase cycle slip; during processing, the satellite cutoff elevation angle is 10°.

3. The high-precision ionospheric modeling method for integrating multi-source ionospheric observation data according to claim 1, characterized in that: The sampling rate of the observation data selected in the first step is 30s; the height of the ionosphere single layer model used is 450km.

4. The high-precision ionospheric modeling method for fusing multi-source ionospheric observation data according to claim 1, characterized in that: The multi-source ionospheric observation data in the third step include: occultation ionospheric data, DORIS ionospheric data, top ionospheric data provided by LEO onboard GNSS receiver, Jason-2 VTEC and ground-based multi-system GNSS ionospheric data.

5. The high-precision ionospheric modeling method for fusing multi-source ionospheric observation data according to claim 1, characterized in that: The method for extracting multi-source ionospheric observation data in the third step comprises: a) Occultation ionospheric observation data: provided by the "ionPrf" product, which provides satellite upper and lower ionospheric observation data to obtain TEC values ​​for the entire ionospheric path: VTEC=VTEC bellow +VTEC up Where VTEC bellow Indicates the VTEC value below the COSMIC satellite; VTEC up represents the VTEC value above the COSMIC satellite, obtained by extrapolation; b) DORIS ionospheric observation data: The calculation steps include: 1) Based on the position of the DORIS ground-based monitoring station and the LEO satellite position, the puncture point position and the zenith distance at the puncture point of the DORIS observation data are calculated, and the VTEC value at the puncture point is calculated using the ionospheric TEC model constructed by the ground-based GNSS; 2) Calculate the mean deviation between the DORIS-containing STEC value and the STEC value converted by GIM in each continuous arc segment; 3) Correct the obtained average deviation within the arc segment to the corresponding arc segment, that is, obtain the corrected STEC value; c) Jason2 VTEC is obtained based on the pseudo-range ionospheric delay correction of the Ku-band frequency of the satellite altimeter. The specific formula is as follows: Where, dR Ku is the Ku-band ionospheric distance correction; f Ku is the Ku-band frequency in GHz; d) LEO top ionosphere observation data and ground-based GNSS ionosphere observation data are extracted using the phase smoothed pseudorange method.

Images