Global ionosphere inversion method based on BP (Back Propagation) neural network fused with multi-source data

Abstract

The invention relates to the technical field of satellite positioning navigation enhancement, in particular to a global ionosphere inversion method based on BP neural network fusion multi-source data, and the method comprises the steps: inverting STEC through ground-based GNSS observation data, then employing a projection function to project the STEC to a vertical direction of an ionosphere puncture point, and obtaining VTEC; a BP neural network and a Kriging interpolation method are utilized to fuse foundation observation data into an IRI-2016 ionosphere model to solve the problem that GNSS tracking stations are distributed unevenly in the world, so that the inversion precision of the ionosphere is improved, a spherical harmonic function is adopted to express VTEC space distribution, and then the ocean altimetry satellite observation data is utilized to achieve the purpose of solving the problem that the GNSS tracking stations are distributed unevenly in the world. According to the method, the problem of deficiency of ground-based GNSS observation data in an ocean area is solved by using ground-based GNSS observation data, DORIS observation data and COSMIC occultation data, a global ionosphere observation equation is established, a normal equation is constructed, and finally global ionosphere grid data VTEC is generated through Helmert variance estimation precise weighting calculation, so that the global ionosphere observation precision is improved.

Description

Technical Field

[0001] The invention relates to the technical field of satellite positioning and navigation enhancement, and in particular to a global ionosphere inversion method based on BP neural network fusion of multi-source data. Background Art

[0002] The ionosphere is an important part of near-Earth space. Research on the ionosphere can promote people's understanding of the Sun-Earth space system and thus better serve human space activities. Research on the physical characteristics of the ionosphere is of great significance. The requirements for improving the accuracy of ionospheric inversion are getting higher and higher. It is also particularly important to obtain high-precision global ionospheric VTEC. In practical applications, dual-frequency receiver users can use dual-frequency corrections to eliminate the main ionospheric errors, but for single-frequency receiver users, high-precision ionospheric models are extremely important for ionospheric error correction.

[0003] At present, the global ionospheric TEC products of the ionospheric analysis centers of various international GNSS service organizations (International GNSS Service, IGS) still only use ground-based GNSS observation data. Since the spherical harmonic function used to establish the global ionosphere model mathematically requires that the globally distributed GNSS observation stations must be uniform in order to obtain the ideal effect, but the GNSS tracking stations are not evenly distributed around the world, especially in the vast ocean areas, there is a phenomenon of blank tracking stations, resulting in low model accuracy in the ocean area. Therefore, researchers use ionospheric empirical models (IRI model, NeQuick model, Klobuchar model, etc.) as virtual observation values ​​to fill the blank observation areas, and use the VTEC values ​​calculated by the model as virtual observation values ​​to constrain blank areas such as the equator or the southern hemisphere, which can improve the problem of uneven distribution of GNSS tracking stations. However, due to the accuracy of the empirical model, the accuracy of the fused model is limited. To address this problem, researchers have also used a series of algorithms for smoothing and constraints, such as Kalman filtering, sequential least squares, and piecewise linearization, but the global VTEC accuracy still needs to be improved. On the other hand, scholars have integrated multi-source data with ground-based GNSS observation data (ocean altimetry satellites, DORIS, COSMIC, etc.) to estimate the global ionospheric VTEC, thereby providing the possibility of improving the accuracy of the global ionospheric model in the ocean area to a certain extent. However, due to the different design principles, satellite altitudes and tracking methods of various satellite systems, the accuracy of ionospheric data of various systems is different, and there are systematic deviations between systems. It is necessary to accurately weight various types of data. If the observation data of each system are weighted according to the classical adjustment method, it may cause problems such as inaccurate weighting, which will affect its accuracy and make it difficult to meet the accuracy requirements. Summary of the invention

[0004] The purpose of the present invention is to provide a global ionospheric inversion method based on BP neural network fusion of multi-source data, aiming to solve the problem of limited accuracy of the existing global ionospheric TEC fusion based on the ionospheric empirical model as the background model combined with the ground-based GNSS data inversion, accurately weight the data of different satellite systems, and improve the accuracy of the global ionospheric grid data.

[0005] To achieve the above object, the present invention provides a global ionosphere inversion method based on BP neural network fusion multi-source data, comprising the following steps:

[0006] Step 1: Preprocess the multi-source and multi-system GNSS observation data, invert the oblique ionospheric delay of each system through the observation data, calculate the coordinates of the ionospheric piercing point of each system through the ephemeris data, and use the projection function to project to the vertical direction of the ionospheric piercing point to obtain the ionospheric VTEC;

[0007] Step 2: Using the international reference ionosphere IRI-2016 model as the background model, the BP neural network is used to fuse the IRI-2016 model and the ionospheric VTEC obtained by ground-based inversion to obtain ΔVTEC;

[0008] Step 3: Obtain ΔVTEC at the location to be determined by Kriging interpolation method, and use the IRI-2016 model at the location to be determined as the background model to obtain VTEC at the location to be determined. Combined with the VTEC obtained by ground foundation inversion, high-precision VTEC can be obtained.

[0009] Step 4: Invert the VTEC data of ocean satellite track points through the Ku- and C-band dual-frequency data of ocean altimetry satellites, and invert the absolute VTEC data through the DORIS dual-frequency observation data and the VTEC data obtained by the COSMIC occultation system;

[0010] Step 5: Use spherical harmonics combined with other obtained VTEC data to establish the global ionosphere estimation observation equation and construct the normal equation in combination with the observation equation;

[0011] Step 6: Estimate the weights of observation data from different sources and systems through the Helmert variance component estimation method, and solve the parameters to be estimated in the model to obtain the global ionospheric VTEC grid data.

[0012] Optionally, the execution process of step 1 includes the following steps:

[0013] Step 1.1: Preprocessing of ground-based multi-system GNSS observation data: cycle slip detection and repair of observation data by combining MW combination method and ionospheric residual method; setting sampling rate to 30s and cutoff elevation angle to 20°;

[0014] Step 1.2: Preprocessing of Jason-2 / 3 data: Obtain ionospheric correction data of Jason-2 / 3 satellite altimeter through GDR product files, obtain TEC, remove data on ocean and ice surface, and then perform a 25s window smoothing.

[0015] Step 1.3: The preprocessing method of DORIS system observation data is the same as that of ground-based multi-system GNSS observation data;

[0016] Step 1.4: Obtain ground-based dual-frequency pseudorange and carrier phase observation data, without considering the influence of multipath effect and observation noise, subtract the pseudorange observation values ​​and carrier phase observation values ​​on the two frequencies to form an observation value without geometric distance combination, and use the carrier phase smoothed pseudorange algorithm to introduce the smoothing amount to solve the high-precision STEC;

[0017] Step 1.5: Select the ionospheric single layer model height as 450km, select the commonly used SLM projection function as the projection function, and project the ionospheric VTEC.

[0018] Optionally, the execution process of step 2 is specifically that the BP neural network generates a predicted ΔVTEC by using the VTEC value provided by the IRI model and the VTEC value obtained by the foundation inversion measurement as a forward transmission signal, and then calculates the prediction error, and then adjusts the network weights by back propagation of the error, so that the prediction result is closer to the expected output ΔVTEC.

[0019] Optionally, in step 3, the weights are determined by first performing optimal and unbiased estimation on the values ​​of the variables using the Kriging interpolation method, and then solving for the electron content difference at the location to be determined.

[0020] Optionally, the data used by the ocean altimetry satellite in step 4 are Jason-2 and Jason-3 ocean altimetry satellite data, and the data is obtained from the GDR product of CNSE.

[0021] Optionally, the DORIS data acquisition method in step 4 includes the following steps:

[0022] Step 4.1: Invert the GIM VTEC at the DORIS ionospheric penetration point using ground-based GNSS observation data;

[0023] Step 4.2: Project the GIM VTEC onto the signal propagation path using the SLM projection function to obtain the GIM STEC, and then calculate the average value of the difference between the GIM STEC and the relative DORIS STEC in each DORIS observation arc;

[0024] Step 4.3: Add this average value to the DORIS STEC value to get the absolute DORIS STEC, and finally project the DORIS STEC using the SLM projection function to get the DORIS VTEC;

[0025] Step 4.4: Use COSMIC occultation data to obtain VTEC, where the VTEC data source is extracted from the ionPrf product provided by UCAR, the Universities Consortium for Atmospheric Research in the United States.

[0026] Optionally, the observation equation establishment process in step 5 is specifically to form a normal equation for observation data from different sources by superposition of normal equations, and the variance of each type of observation data is calculated based on the Helmert variance component estimation.

[0027] Optionally, the process of determining the weights of observation data from different sources and systems using the Helmert variance component estimation method in step 6 includes the following steps:

[0028] Step 6.1: Determine the prior variance of GNSS data, DORIS data, Jason-2 / 3 data, and COSMIC data as

[0029] Step 6.2: Perform the first adjustment to obtain the corresponding P i V i ;

[0030] Step 6.3: Use the simplified formula for variance component estimation to calculate the posterior unit weight variance of each type of observation:

[0031]

[0032] Step 6.4: Determine the post-test variance of observation data from different sources based on the post-test unit weight variance

[0033] Step 6.5: Repeat the adjustment, variance component estimation, and posterior variance calculation until the posterior unit weight variances of all types of observations are equal, that is,

[0034] The invention provides a global ionosphere inversion method based on BP neural network fusion of multi-source data. The oblique ionosphere TEC is inverted through ground-based GNSS observation data, and then a projection function is used to project the oblique ionosphere TEC to the vertical direction of the ionosphere puncture point to obtain the vertical ionosphere TEC. The ground-based observation data is fused into the IRI-2016 ionosphere model by using the BP neural network and the Kriging interpolation method to solve the problem of uneven global distribution of GNSS tracking stations, thereby improving the inversion accuracy of the ionosphere. Spherical harmonic functions are used to describe the VTEC spatial distribution, and then the lack of ground-based GNSS observation data in ocean areas is improved by using ocean altimetry satellite observation data, DORIS observation data and COSMIC occultation data. A global ionosphere observation equation is established, a normal equation is constructed, and finally the global ionosphere grid data VTEC is generated by accurately weighting calculation through Helmert variance estimation, so as to improve the global ionosphere observation accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0036] Figure 1 It is a schematic flow chart of the steps of the global ionosphere inversion method based on BP neural network fusion of multi-source data of the present invention.

[0037] Figure 2 It is a specific solution flow chart of the global ionosphere inversion method based on BP neural network fusion of multi-source data of the present invention.

[0038] Figure 3 It is an algorithm flow chart of the BP neural network of the global ionosphere inversion method based on the BP neural network fusion of multi-source data of the present invention.

[0039] Figure 4 The present invention is a flow chart of VTEC inversion by the DORIS system based on the global ionosphere inversion method of the BP neural network fusion multi-source data. DETAILED DESCRIPTION

[0040] Embodiments of the present invention are described in detail below, examples of which are shown in the accompanying drawings, wherein the same or similar reference numerals throughout represent the same or similar elements or elements having the same or similar functions. The embodiments described below with reference to the accompanying drawings are exemplary and are intended to be used to explain the present invention, and should not be construed as limiting the present invention.

[0041] The Chinese and English meanings of the English abbreviations appearing in this application are as follows:

[0042] BP: backpropagation, back propagation algorithm;

[0043] VTEC: Vertical Total Electron Content, ionosphere vertical total electron content;

[0044] STEC: Slant Total Electron Content, the total electron content in the ionosphere line of sight;

[0045] ΔVTEC, difference in ionospheric electron content;

[0046] DORIS: Doppler Orbitography and Radio-positioning Integrated by Satellite, abbreviated as Doris system or DORIS system;

[0047] COSMIC: the Constellation Observing System for Meteorology, Ionosphere and Climate

[0048] MW: Melbourne-Wubbena, Melbourne-Wubbena combination;

[0049] CNSE: Centre National d'Etudes Spatiales, French National Space Research Center;

[0050] See also Figure 1 The present invention provides a global ionosphere inversion method based on BP neural network fusion multi-source data, comprising the following steps:

[0051] S1: Preprocess the multi-source and multi-system GNSS observation data, invert the oblique ionospheric delay of each system through the observation data, calculate the coordinates of the ionospheric piercing point of each system through the ephemeris data, and use the projection function to project to the vertical direction of the ionospheric piercing point to obtain the ionospheric VTEC;

[0052] S2: The international reference ionosphere IRI-2016 model is used as the background model, and the ionospheric VTEC obtained by ground-based inversion is fused with the IRI-2016 model using a BP neural network to obtain ΔVTEC;

[0053] S3: The ΔVTEC of the location to be determined is obtained by the Kriging interpolation method, and the VTEC of the location to be determined is obtained by using the IRI-2016 model of the location to be determined as the background model. The high-precision VTEC can be obtained by combining the VTEC obtained by the foundation inversion.

[0054] S4: VTEC data of ocean satellite track points are inverted from Ku- and C-band dual-frequency data of ocean altimetry satellites, and absolute VTEC data are inverted from DORIS dual-frequency observation data and VTEC data are obtained from COSMIC occultation system;

[0055] S5: Use spherical harmonics combined with other obtained VTEC data to establish the global ionosphere estimation observation equation and construct the normal equation by combining the observation equation;

[0056] S6: The weights of observation data from different sources and systems are estimated through the Helmert variance component estimation method, and the parameters to be estimated in the model are solved to obtain the global ionospheric VTEC grid data.

[0057] The specific ionospheric inversion method process is as follows: Figure 2 As shown, the following is further described in conjunction with specific embodiments and execution steps:

[0058] The execution process of step S1 can be divided into the following steps:

[0059] (1.1) Preprocessing of ground-based multi-system GNSS observation data: The MW combination method and the ionospheric residual method are used to detect and repair cycle slips in the observation data. The sampling rate is set to 30s and the cutoff elevation angle is set to 20°.

[0060] (1.2) VTEC value processing of Jason-2 / 3 ocean altimeter satellite data: The ionospheric correction data of the Jason3 satellite altimeter is obtained through the GDR product file. After VTEC acquisition, the data on the ocean and ice surface are removed, and then a 25 s window smoothing is performed.

[0061] (1.3) DORIS system observation data preprocessing: the same as the ground-based multi-system GNSS observation data preprocessing method.

[0062] (1.4) Using pseudorange and carrier phase observation data, the total electron content STEC on the propagation path is obtained by using the carrier phase smoothed pseudorange method. The specific steps are as follows:

[0063] Obtain ground-based dual-frequency pseudorange and carrier phase observation data, and its observation equation is:

[0064]

[0065] In the formula, f represents the signal frequency of each system; r represents the station receiver identifier; s represents the satellite identifier; represents the pseudorange observation value at frequency f; represents the carrier phase observation value at frequency f; represents the geometric distance between the receiver and the satellite; c represents the speed of light; dt r Indicates the receiver clock error; dt s represents the satellite clock error; represents the ionospheric delay on the propagation path from receiver r to satellite s; represents the tropospheric delay on the propagation path from receiver r to satellite s; b r,f and represents the ranging code hardware delay between satellite s and receiver r on frequency f; B r,f and represents the hardware delay in the carrier phase of the satellite s and the receiver r at frequency f; Δρ represents the error caused by the phase center error of the receiver and satellite antenna, the relativistic effect, the rotation of the earth, the solid tide and the ocean load tide; λ f represents the carrier phase wavelength at frequency f; represents the carrier phase integer ambiguity at frequency f; represents the pseudorange observation noise of satellite s and receiver r at frequency f; represents the carrier phase observation noise of satellite s and receiver r at frequency f.

[0066] Without considering the influence of multipath effect and observation noise, the pseudo-range observation value and carrier phase observation value on the two frequencies are subtracted to form the observation value without geometric distance combination. The observation equation is as follows:

[0067]

[0068] In the formula, and They represent the pseudorange and phase geometry-free range combined observation values ​​respectively;

[0069] The carrier phase smoothing pseudorange algorithm is used to introduce the smoothing value D rs Solve the high-precision STEC, the expression is:

[0070]

[0071] Among them, N TECp is the STEC calculated from the dual-frequency pseudorange observations, N TECl is the STEC calculated from the dual-frequency carrier phase observation value, and N is the effective sampling number. Then the high-precision STEC can be calculated as: STEC = N TECl +D rs

[0072] It should be noted that the data preprocessing can select pseudorange and carrier phase observation data obtained by different satellite systems, and use the wide-lane MW combination method and the ionospheric residual method to detect and repair the observation data. The sampling rate is set to 30s and the cut-off elevation angle is 20°. The preprocessing of Jason-2 / 3 ocean altimetry satellite data is to remove the data on the ocean and ice surface from the VTEC value calculated by the ocean altimetry data, and then perform a 25s window smoothing. The preprocessing method for the dual-frequency observation data of the DORIS system is the same as the preprocessing method for the GNSS ground-based observation data. Since the N TECp The value is the absolute total electron content, but the pseudorange has low accuracy and large noise, while the N obtained by carrier phase observation data is TECl The value is relative to the total electron content. Although the carrier phase has high accuracy, it can only reflect the relative change of the total electron content due to the influence of the integer ambiguity. To solve the above problem, the present invention adopts the carrier phase smoothing pseudorange method to solve STEC, which combines the advantages of pseudorange and carrier phase and obtains the absolute change of the total electron content with high accuracy.

[0073] (1.5) The specific method for calculating the vertical ionospheric delay is as follows: select the ionospheric single layer model height as 450km, and select the commonly used SLM projection function as the projection function. The SLM projection function expression is as follows:

[0074]

[0075] Among them, z is the zenith distance at IPP, E is the satellite altitude angle, R is the radius of the earth, and H is the height of the single-layer model 450km;

[0076] The execution process of step S2 includes the following steps:

[0077] (2.1) The IRI model is an empirical model obtained based on a large amount of ionospheric detection data. The IRI-2016 model data is downloaded from the official website of the International Reference Ionosphere. It can calculate the electron density in the range of 50 to 2000 km above sea level. The electron density in the range of 2000 km to the satellite orbit altitude can be obtained by extrapolation. The equation for calculating the electron density below 2000 km by the IRI model is:

[0078]

[0079] Where f(h) is the electron density corresponding to height h, y i is the height h i The corresponding electron density, k is the coefficient to be determined. After solving the coefficient k by least squares, the electron density above 2000km altitude can be extrapolated, and finally the electron density within the range of 50km to the satellite orbit altitude can be accumulated to obtain the total electron content.

[0080] (2.2) BP neural network is selected to fuse the VTEC obtained by IRI model and ground-based inversion to obtain ΔVTEC. The core part of BP neural network is the process of continuously correcting the weights of hidden layer and output layer during continuous training. The training process of neural network is as follows: given the VTEC value generated by IRI model as input vector and the VTEC value obtained by ground-based inversion as expected output, the error between expected output and actual output is calculated, and then the error gradient is calculated. The weights are adjusted to continuously train and correct the weights of hidden layer and output layer, and the error ΔVTEC of the area to be calculated is output. The error ΔVTEC is fused with IRI background mode to obtain the VTEC value of the selected ionospheric inversion area, thereby solving the problem of low VTEC inversion accuracy in some areas due to blank tracking stations.

[0081] The execution process of step S3 includes the following steps:

[0082] (3.1) The method for obtaining the ΔVTEC at the desired location by Kriging interpolation is as follows:

[0083]

[0084] Where, ΔVTEC IPP,i is the vertical ionospheric electron content difference at the puncture point, W i is the inverse of the variation function of the electron content data.

[0085] The execution process of step S4 includes the following steps:

[0086] (4.1) The data used by the ocean altimetry satellite is the Jason-2 / 3 ocean altimetry satellite data, which is obtained by extracting the GDR product of CNSE (Centre National d'Etudes Spatiales). The VTEC below the satellite altitude is inverted using the Jason-2 / 3 altimetry data. Jas2 / 3 , and its calculation formula is:

[0087]

[0088] Among them, f Ku =13.6GHz, f C =5.3GHz, h Ku 、h C are the distances from the satellite to the ocean surface measured by the Ku and C bands, respectively. Ku , D C are the sea state deviation corrections for Ku and C bands, B Ku , B C These are the instrument bias corrections for the Ku and C bands, respectively.

[0089] It should be noted that the orbital altitude of Jason-2 / 3 series satellites is 1336km, and the obtained VTEC does not include the contribution of the ionosphere above the satellite orbit; this product provides the distance from the satellite to the ocean surface measured in different bands and the corresponding profile bias correction value; among them, the distance from the satellite to the ocean surface has included the instrument error correction including DCB;

[0090] (4.2) DORIS data acquisition is to obtain the ranging code and phase observation values ​​in RINEX format through the DORIS receiver, and then use the DORIS satellite dual-frequency data to invert STEC, and its expression is:

[0091]

[0092] Where STEC is relative STEC, f1=2036.25MHz, f2=401.25MHz, λ1 and λ2 are the corresponding wavelengths, φ1 and φ2 are phase observations, and ΔD=D1-D2 includes the integer ambiguity in the two-frequency phase observations and the dual-frequency phase observation deviation.

[0093] In order to obtain absolute TEC, it is necessary to use the GIM VTEC of the ground-based inversion to realize the conversion from relative TEC to absolute TEC. The calculation flow chart is as follows: Figure 4 As shown, the specific steps are: invert the GIM VTEC where the DORIS ionospheric piercing point is located through ground-based GNSS observation data; project the GIM VTEC to the signal propagation path through the SLM projection function to obtain GIM STEC, and then calculate the average value of the difference between GIM STEC and the relative DORIS STEC in each DORIS observation arc; then add this average value to the DORIS STEC value to obtain the absolute DORIS STEC, and finally project the DORIS STEC through the SLM projection function to obtain the DORIS VTEC. It should be noted that since the average value already includes the TEC part above the DORIS satellite orbit altitude, the final DORIS VTEC covers the entire ionospheric altitude. Specifically, the DORIS system VTEC inversion flow chart is shown below Figure 4 shown.

[0094] (4.3) VTEC is obtained using COSMIC occultation data. The VTEC data is extracted from the ionPrf product provided by the University Corporation for Atmospheric Research (UCAR).

[0095] The execution process of step S5 includes the following steps:

[0096] (5.1) The global ionosphere model adopts a spherical harmonic model, where the order of the spherical harmonic model is 15×15, the time resolution is 2 hours, the interval in the latitude direction is 2.5°, and the interval in the longitude direction is 5°. The expression of the spherical harmonic function is:

[0097]

[0098] Where β is the geomagnetic latitude at the IPP, n and m are the orders of the spherical harmonics, and n max To expand the maximum order, are the model parameters to be determined, is the normalized Langrangian function of order n and m, MC(n,m) is the regularization function, δ om represents the Kronecker function. s is the IPP longitude in the solar-solid system, and the calculation formula is: s = λ-λ0≈UT+λ-π, where UT is the universal time, λ is the IPP geographical longitude, and λ0 represents the solar longitude.

[0099] (5.2) The specific observation equation is as follows:

[0100]

[0101] Where VTEC(β,s) is expressed by spherical harmonics; VTEC BP-IRI ,VTEC Jas2 / 3 ,VTEC DORIS ,VTEC COSMIC Respectively, they represent BP-IRI2016 VTEC, Jason-2 / 3 altimeter VTEC, DORIS system VTEC, and COSMIC occultation system VTEC; DCB r ,DCB s They represent the hardware delay of the receiver and the hardware delay of the satellite respectively; bias represents the systematic deviation of various types of observation data.

[0102] (5.3) For observation data from different sources, the normal equation is formed by superposition of normal equations, which can be expressed as:

[0103]

[0104] Where m represents the m-type observation value, including GNSS data, BP-IRI-2016 data, DORIS data, COSMIC data and Jason-2 / 3 data; L i Represents various observation vectors, B i represents the design matrix, and X represents the spherical harmonic coefficients to be estimated. i Represents the weight matrix of different types of observations; the same type of observation data determines the weight according to the satellite altitude angle, p i =sin2 (ele), where ele is the satellite altitude angle; represents the variance of each category of observations, calculated using the Helmert variance component estimate.

[0105] The execution process of step S6 includes the following steps:

[0106] (6.1) The parameters to be estimated can be obtained from the normal equation in step 5.3 The adjustment value is expressed as:

[0107]

[0108] The definitions of the parameters in the formula are consistent with those described in step 5.3.

[0109] (6.2) The Helmert variance component estimation method is used to determine the weights of the current five different sources and systems of observation data: P1, P2, P3, P4, and P5; the covariance of the parameters to be estimated for GNSS data, IRI data, DORIS data, Jason-2 / 3 data, and COSMIC data is determined, that is, the prior variance is expressed as: The diagonal elements of D represent the variance of the estimated spherical harmonic coefficients.

[0110] (6.3) Perform the first least squares adjustment, and the residual vector corresponding to the observed data is V i , there is a relationship Calculate the weighted residual sum of squares corresponding to each type of observation data P i V i , which is used for the subsequent variance component estimation.

[0111] (6.4) Since the Helmert variance component estimation formula is relatively complicated, the simplified formula of variance component estimation is often used in practice to calculate the posterior unit weight variance of various observations, expressed as:

[0112]

[0113] Where n i For degrees of freedom.

[0114] (6.5)Then the weight is determined according to the following formula:

[0115]

[0116] In the formula, i = 1, 2...5; c is an arbitrary constant, usually selected A value of .

[0117] (6.6) The adjustment-variance component estimation-weight determination cycle is repeated until the post-test unit weight variances of all types of observations are equal, that is,

[0118] (6.7) After determining the weight ratio between observations of different systems using the Helmert variance component estimate, the spherical harmonic coefficient X can be solved according to the estimated parameter formula in step 6.1. Then, the global ionospheric VTEC grid of the current epoch can be calculated in combination with the spherical harmonic coefficient model in step 5.1.

[0119] What is disclosed above is only a preferred embodiment of the present invention, and it certainly cannot be used to limit the scope of rights of the present invention. Ordinary technicians in this field can understand that all or part of the processes of the above embodiment and equivalent changes made according to the claims of the present invention still fall within the scope of the invention.

Claims

1. A global ionosphere inversion method based on BP neural network fusion of multi-source data, characterized in that: The following steps are involved: Step 1: Preprocess the multi-source and multi-system GNSS observation data, invert the oblique ionospheric delay of each system through the observation data, calculate the coordinates of the ionospheric piercing point of each system through the ephemeris data, and use the projection function to project to the vertical direction of the ionospheric piercing point to obtain VTEC; Step 2: Using the international reference ionosphere IRI-2016 model as the background model, the BP neural network is used to fuse the IRI-2016 model and the ionospheric VTEC obtained by ground-based inversion to obtain ΔVTEC; Step 3: Obtain ΔVTEC at the location to be determined by Kriging interpolation method, and use the IRI-2016 model at the location to be determined as the background model to obtain VTEC at the location to be determined. Combined with the VTEC obtained by ground foundation inversion, high-precision VTEC can be obtained. Step 4: Invert the VTEC data of ocean satellite track points through the Ku- and C-band dual-frequency data of ocean altimetry satellites, and invert the absolute VTEC data through the DORIS dual-frequency observation data and the VTEC data obtained by the COSMIC occultation system; Step 5: Use spherical harmonics combined with the VTEC data obtained in step 4 to establish the global ionosphere estimation observation equation and construct the normal equation in combination with the observation equation; Step 6: Estimate the weights of observation data from different sources and systems through the Helmert variance component estimation method, and solve the model parameters to be estimated, thereby obtaining the global ionospheric VTEC grid data.

2. The global ionosphere inversion method based on BP neural network fusion multi-source data as claimed in claim 1, characterized in that: The execution process of step 1 includes the following steps: Step 1.1: Preprocessing of ground-based multi-system GNSS observation data: cycle slip detection and repair of observation data by combining MW combination method and ionospheric residual method; setting sampling rate to 30s and cutoff elevation angle to 20°; Step 1.2: Preprocessing of Jason-2 / 3 data: Obtain ionospheric correction data of Jason-2 / 3 satellite altimeter through GDR product files, obtain TEC, remove data on ocean and ice surface, and then perform a 25s window smoothing. Step 1.3: The preprocessing method of DORIS system observation data is the same as that of ground-based multi-system GNSS observation data; Step 1.4: Obtain ground-based dual-frequency pseudorange and carrier phase observation data, without considering the influence of multipath effect and observation noise, subtract the pseudorange observation values ​​and carrier phase observation values ​​on the two frequencies to form an observation value without geometric distance combination, and use the carrier phase smoothed pseudorange algorithm to introduce the smoothing amount to solve STEC; Step 1.5: Select the ionospheric single layer model height as 450km, select the commonly used SLM projection function as the projection function, and project the ionospheric VTEC.

3. The global ionosphere inversion method based on BP neural network fusion multi-source data as claimed in claim 2, characterized in that: The execution process of step 2 is as follows: the BP neural network generates the predicted ΔVTEC by using the VTEC value provided by the IRI model and the VTEC value obtained by the foundation inversion measurement as the forward transmission signal, and then calculates the prediction error, and then adjusts the network weights by back propagation of the error, so that the prediction result is closer to the expected output ΔVTEC.

4. The global ionosphere inversion method based on BP neural network fusion multi-source data as claimed in claim 3 is characterized in that: In step 3, the Kriging interpolation method is used to first perform optimal and unbiased estimation on the values ​​of the variables to determine the weights, and then the electron content difference at the location to be determined is solved.

5. The global ionosphere inversion method based on BP neural network fusion multi-source data as claimed in claim 4 is characterized in that: The data used by the ocean altimetry satellite in step 4 are Jason-2 and Jason-3 ocean altimetry satellite data, and the data is obtained from the GDR product of CNSE.

6. The global ionosphere inversion method based on BP neural network fusion of multi-source data as claimed in claim 5, characterized in that: The DORIS data acquisition method in step 4 includes the following steps: Step 4.1: Invert the GIM VTEC at the DORIS ionospheric penetration point using ground-based GNSS observation data; Step 4.2: Project the GIM VTEC onto the signal propagation path using the SLM projection function to obtain the GIM STEC, and then calculate the average value of the difference between the GIM STEC and the relative DORIS STEC in each DORIS observation arc; Step 4.3: Add this average value to the DORIS STEC value to get the absolute DORIS STEC, and finally project the DORIS STEC using the SLM projection function to get the DORIS VTEC; Step 4.4: Use COSMIC occultation data to obtain VTEC, where the VTEC data source is extracted from the ionPrf product provided by UCAR, the Universities Consortium for Atmospheric Research in the United States.

7. The global ionosphere inversion method based on BP neural network fusion of multi-source data as claimed in claim 6, characterized in that: The observation equation establishment process in step 5 is to form the normal equation by superimposing the normal equations for the observation data from different sources. The variance of each type of observation data is calculated based on the Helmert variance component estimation.

8. The global ionosphere inversion method based on BP neural network fusion of multi-source data as claimed in claim 7, characterized in that: The process of determining the weights of observation data from different sources and systems using the Helmert variance component estimation method in step 6 includes the following steps: Step 6.1: Determine the prior variance of GNSS data, DORIS data, Jason-2 / 3 data and COSMIC data as D i =σ i 2 P i -1 ; Step 6.2: Perform the first adjustment to obtain the V corresponding to each type of observation data i T P i V i ; Step 6.3: Use the simplified formula for variance component estimation to calculate the posterior unit weight variance of each type of observation: Step 6.4: Determine the post-test variance of observation data from different sources based on the post-test unit weight variance Step 6.5: Repeat the adjustment, variance component estimation, and posterior variance calculation until the posterior unit weight variances of all types of observations are equal, that is,

Images