Global ionospheric inversion method based on bp neural network fusion of multi-source data

Abstract

The present application relates to the technical field of satellite positioning and navigation enhancement, in particular to a global ionospheric inversion method based on BP neural network fusion of multi-source data, STEC is inverted through ground-based GNSS observation data, then a projection function is used to project it to the vertical direction of the ionospheric piercing point to obtain VTEC, the ground-based observation data is fused into the IRI-2016 ionospheric model by using the BP neural network and the Kriging interpolation method to solve the problem of uneven global distribution of GNSS tracking stations, so as to improve the inversion accuracy of the ionosphere, spherical harmonics is used to express the spatial distribution of VTEC, then the missing problem of ground-based GNSS observation data in the ocean area is improved by using ocean altimetry satellite observation data, DORIS observation data and COSMIC occultation data, a global ionospheric observation equation is established, a normal equation is constructed, and finally the global ionospheric grid data VTEC is generated by accurate weight calculation through Helmert variance estimation, so as to improve the global ionospheric observation accuracy.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of satellite positioning navigation enhancement, and particularly relates to a global ionospheric inversion method based on BP neural network fusion of multi-source data. BACKGROUND

[0002] The ionosphere is an important part of the near space, and the research on the ionosphere can help people understand the space system between the sun and the earth, and thus better serve the space activities of human beings. It is of great significance to study the physical characteristics of the ionosphere, and the requirement for improving the ionospheric inversion accuracy is getting higher and higher. It is particularly important to obtain high-precision global ionospheric VTEC. In practical applications, dual-frequency receiver users can use dual-frequency correction to eliminate the main error of the ionosphere, but for single-frequency receiver users, a high-precision ionospheric model is extremely important for ionospheric error correction.

[0003] At present, the global ionospheric TEC products of the ionospheric analysis center of each international GNSS service (International GNSS Service, IGS) still only use ground-based GNSS observation data. Since the spherical harmonic function used to establish the global ionospheric model requires that the global distribution of GNSS observation stations should be uniform to obtain ideal results, but the GNSS tracking stations are not uniformly distributed in the world, especially in the vast ocean area, there is a phenomenon of tracking station blanking, which leads to low model accuracy in the ocean area. Therefore, researchers use ionospheric empirical models (IRI model, NeQuick model, Klobuchar model, etc.) as virtual observations to fill in the blank observation area. The VTEC value calculated by the model is used as a virtual observation to constrain the blank area around the equator or the southern hemisphere, which can improve the problem of uneven distribution of GNSS tracking stations, but limited by the accuracy of the empirical model, the improved accuracy of the fused model is limited. To solve this problem, researchers have used a series of algorithms for smoothing and constraint, such as Kalman filter, sequential least squares, and piecewise linear, but the global VTEC accuracy still needs to be improved. On the other hand, researchers fuse multi-source data and ground-based GNSS observation data (ocean altimetry satellites, DORIS, COSMIC, etc.) to estimate global ionospheric VTEC, which provides the possibility of improving the accuracy of global ionospheric model in the ocean area to some extent. However, due to the different design principles, satellite altitudes and tracking methods of each satellite system, the ionospheric data accuracy of each system is different, and there are systematic biases between systems. It is necessary to accurately weight each type of data. If the classic adjustment method is used to weight each system observation data, it may cause inaccurate weighting, which will affect the accuracy and make it difficult to meet the accuracy requirement. SUMMARY

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

[0005] To achieve the above-mentioned purpose, the application provides a global ionospheric inversion method based on BP neural network fusion of multi-source data, comprising the following steps:

[0006] Step 1: Preprocessing of multi-source multi-system GNSS observation data, inverting the slant ionospheric delay of each system through observation data and calculating the ionospheric piercing point coordinates of each system through ephemeris data, projecting to the vertical direction of the ionospheric piercing point using the projection function, and obtaining the ionospheric VTEC;

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

[0008] Step 3: Obtaining the ΔVTEC at the to-be-solved point by Kriging interpolation method, and obtaining the VTEC at the to-be-solved point by taking the IRI-2016 model at the to-be-solved point as the background model, and combining the ground-based inversion VTEC to obtain high-precision VTEC;

[0009] Step 4: Inverting the VTEC data of the ocean satellite track point through the Ku and C band dual-frequency data of the ocean altimetry satellite and inverting the absolute VTEC data through the DORIS dual-frequency observation data and obtaining the VTEC data through the COSMIC occultation system;

[0010] Step 5: Using spherical harmonics combined with other obtained VTEC data to establish a global ionospheric estimation observation equation and constructing a normal equation combined with the observation equation;

[0011] Step 6: Estimating the weight of the observation data of different sources and different systems by Helmert variance component estimation method, and solving the model to be estimated parameters, and then obtaining the global ionospheric VTEC grid data.

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

[0013] Step 1.1: Ground-based multi-system GNSS observation data preprocessing, using the MW combination method and the ionospheric residual method to detect and repair the cycle slip of the observation data; setting the sampling rate to 30s and the cutoff elevation angle to 20°;

[0014] Step 1.2: The ionospheric correction data of Jason-2 / 3 satellite altimeter is obtained from the GDR product file, and then the TEC is obtained, the data on the sea and ice surface is removed, and a 25s window smoothing is performed;

[0015] Step 1.3: The DORIS system observation data is preprocessed in the same way as the ground-based multi-system GNSS observation data;

[0016] Step 1.4: Obtain the ground-based dual-frequency pseudorange and carrier phase observation data, without considering the influence of multipath effect and observation noise, and the pseudorange observation value and carrier phase observation value at two frequencies are subtracted to form the observation value of the non-geometric distance combination, and the carrier phase smoothing pseudorange algorithm is used to introduce the smoothing amount to solve the high-precision STEC;

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

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

[0019] Optionally, in step 3, the variable value is first optimally and unbiasedly estimated by the Kriging interpolation method to determine the weight value, and then the electron content difference is solved.

[0020] Optionally, in step 4, the data used by the ocean altimetry satellite is Jason-2 and Jason-3 ocean altimetry satellite data, and the data acquisition method is to obtain it from the GDR product of CNES.

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

[0022] Step 4.1: The GIM VTEC where the DORIS ionospheric piercing point is located is inversed from the ground-based GNSS observation data;

[0023] Step 4.2: Project the GIM VTEC to the signal propagation path by 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 segment;

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

[0025] Step 4.4: Obtain VTEC using COSMIC occultation data, wherein the VTEC data is obtained from the ionPrf product provided by the University Corporation for Atmospheric Research (UCAR).

[0026] Optionally, the process of establishing the observation equation in step 5 is specifically that the normal equation is formed by stacking the normal equations of different sources of observation data, and the variances of various observation data are calculated according to the Helmert variance component estimation.

[0027] Optionally, the process of determining the weights of observation data of different sources and different systems in step 6 by using the Helmert variance component estimation method 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 VTEC corresponding to each type of observation data i T P i V i ;

[0030] Step 6.3: Calculate the posterior unit weight variance of each type of observation value by using the simplified formula of variance component estimation:

[0031]

[0032] Step 6.4: Determine the posterior variance of observation data of different sources according to the posterior unit weight variance

[0033] Step 6.5: Perform adjustment, variance component estimation and posterior variance calculation in a loop until the posterior unit weight variances of various observation values are equal, that is,

[0034] The application provides a global ionosphere inversion method based on BP neural network fusion of multi-source data, which inverses oblique ionosphere TEC through ground-based GNSS observation data, then projects it to the vertical direction of the ionosphere piercing point to obtain vertical ionosphere TEC by using a projection function, and solves the problem of uneven global distribution of GNSS tracking stations by fusing ground-based observation data into the IRI-2016 ionosphere model by using BP neural network and Kriging interpolation method, so as to improve the inversion accuracy of the ionosphere, adopts spherical harmonic function to express the spatial distribution of VTEC, and then uses ocean altimetry satellite observation data, DORIS observation data and COSMIC occultation data to improve the problem of lack of ground-based GNSS observation data in the ocean area, establishes a global ionosphere observation equation, constructs a normal equation, and finally generates global ionosphere grid data VTEC by accurate weight calculation through Helmert variance estimation, so as to improve the global ionosphere observation accuracy. BRIEF DESCRIPTION OF DRAWINGS

[0035] In order to more clearly illustrate the technical solutions of the embodiments of the present application or the prior art, the following will briefly introduce the drawings needed to be used in the embodiments or prior art description. Obviously, the drawings in the following description are only some embodiments of the present application, and for those skilled in the art, other drawings can also be obtained without creative labor.

[0036] Figure 1 is a step flowchart of the global ionosphere inversion method based on BP neural network fusion of multi-source data of the present application.

[0037] Figure 2 is a specific scheme flowchart of the global ionosphere inversion method based on BP neural network fusion of multi-source data of the present application.

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

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

[0040] The embodiments of the present application will be described in detail below, and examples of the embodiments are shown in the drawings, wherein the same or similar reference signs represent the same or similar elements or elements having the same or similar functions throughout. The embodiments described below by referring to the drawings are exemplary and are intended to explain the present application, and cannot be understood as a limitation of the present application.

[0041] The English abbreviation terms appearing in the present application have the following Chinese and English meanings:

[0042] BP: backpropagation, backpropagation algorithm;

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

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

[0045] ΔVTEC, ionospheric electron content difference;

[0046] DORIS: Doppler Orbitography and Radio-positioning Integrated by Satellite, 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, the Constellation Observing System for Meteorology, Ionosphere and Climate;

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

[0049] CNES: Centre National d'Etudes Spatiales, Centre National d'Etudes Spatiales;

[0050] Please refer to Figure 1 The present application provides a global ionospheric inversion method based on BP neural network fusion of multi-source data, comprising the following steps:

[0051] S1: Preprocessing multi-source multi-system GNSS observation data, and inversing the ionospheric delay of each system by observation data and calculating the ionospheric piercing point coordinates of each system by ephemeris data, and projecting to the vertical direction of the ionospheric piercing point using the projection function to obtain the ionospheric VTEC;

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

[0053] S3: Obtain the ΔVTEC at the to-be-solved place by the Kriging interpolation method, and obtain the VTEC at the to-be-solved place based on the IRI-2016 model as a background model, and obtain the high-precision VTEC by combining the ground-based inverted VTEC;

[0054] S4: Invert the VTEC data of the ocean satellite track point from the Ku and C band dual-frequency data of the ocean altimetry satellite, and invert the absolute VTEC data from the DORIS dual-frequency observation data and obtain the VTEC data from the COSMIC occultation system;

[0055] S5: Establish a global ionospheric estimation observation equation by using the spherical harmonic function combined with other obtained VTEC data, and construct a method equation combined with the observation equation;

[0056] S6: Estimate the weight of the observation data of different sources and different systems by the Helmert variance component estimation method, solve the model to-be-estimated parameters, and then obtain the global ionospheric VTEC grid data.

[0057] The specific ionospheric inversion method flow is shown in Figure 2 The following will be further described in combination with specific embodiments and execution steps:

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

[0059] (1.1) Ground-based multi-system GNSS observation data preprocessing: The MW combination method and the ionospheric residual method are combined to detect and repair the cycle slip of the observation data; the sampling rate is set to 30s, and the cutoff elevation angle is set to 20°.

[0060] (1.2) Jason-2 / 3 data solving VTEC value processing: Obtain the ionospheric correction data of the Jason3 satellite altimeter through the GDR product file to obtain the VTEC, then remove the data on the ocean and ice surface, and then perform a 25s window smoothing.

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

[0062] (1.4) Use the carrier phase smoothing pseudorange method to obtain the total electron content STEC on the propagation path by using the pseudorange and carrier phase observation data, and the specific steps are as follows:

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

[0064]

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

[0066] Without considering the effects of multipath and observation noise, the difference between the pseudorange observations and the carrier phase observations at two frequencies yields observations without geometric distance combination, and the observation equation is as follows:

[0067]

[0068] In the formula, and These represent pseudorange and phase-without-geometric-distance combined observations, respectively.

[0069] A smoothing amount D is introduced using a carrier phase smoothing pseudorange algorithm. rs The expression for solving high-precision STEC is:

[0070]

[0071] Where, N TECp STEC, N is calculated from dual-frequency pseudorange observations. TECl The STEC is calculated from the dual-frequency carrier phase observations, where N is the number of effective samples. Therefore, a high-precision STEC can be calculated as: STEC = N TECl +D rs

[0072] It should be noted that the data preprocessing can select the pseudo-range and carrier phase observation data obtained by different satellite systems, adopt the combination of wide lane MW and ionosphere residual method to detect and repair the cycle slip of the observation data, set the sampling rate to 30s, and the cut-off elevation angle is 20°; and the preprocessing of Jason-2 / 3 ocean altimetry satellite data is to remove the data on the ocean and ice surface, and then perform 25s window smoothing on the VTEC value obtained by solving the ocean altimetry data; the preprocessing method of the dual-frequency observation data of the DORIS system is the same as that of the GNSS ground observation data; since the N TECp value obtained by the pseudo-range observation data is the absolute total electron content, but the pseudo-range has low precision and large noise, and the N TECl value obtained by the carrier phase observation data is the relative total electron content, although the carrier phase has high precision, but is affected by the integer ambiguity and can only reflect the relative change of the total electron content, in view of the above problems, the carrier phase smoothing pseudo-range method is adopted to solve the STEC, which combines the advantages of pseudo-range and carrier phase, and obtains the absolute change of high-precision total electron content.

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

[0074]

[0075] Wherein, z is the zenith distance at IPP, E is the satellite elevation angle, R is the earth radius, and H is the single-layer model height of 450km.

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

[0077] (2.1) The IRI model is an empirical model obtained according to a large amount of ionospheric detection data, the IRI-2016 model data is downloaded from the international reference ionosphere website, and the electron density in the range of 50-2000km above sea level can be calculated, and the electron density above 2000km to the satellite orbit height range can be obtained by extrapolation. The equation for calculating the electron density below 2000km from the IRI model is:

[0078]

[0079] In the formula, f(h) is the electron density corresponding to the height h, y i is the electron density corresponding to the height h i , and k is the coefficient to be solved. After solving the coefficient k by the least square method, the electron density above 2000km can be extrapolated, and finally the electron density in the range of 50km to the satellite orbit height is accumulated, so that the total electron content is obtained.

[0080] (2.2) Select BP neural network to fuse IRI model and VTEC obtained by ground-based inversion to obtain ΔVTEC, the core part of BP neural network is to continuously correct the weights of hidden layer and output layer in the continuous training process, the training process of neural network is: the VTEC value generated by the given IRI model is the input vector and the VTEC value inverted by the ground is the expected output, the error between the expected output and the actual output is calculated, 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 region to be solved is output. Fuse the error ΔVTEC with the IRI background model to obtain the VTEC value of the selected ionospheric inversion region, thereby solving the problem of low VTEC inversion accuracy caused by the blank of tracking stations in some areas.

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

[0082] (3.1) The method for obtaining ΔVTEC at the to-be-solved point by Kriging interpolation method 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 Jason-2 / 3 ocean altimetry satellite data, and the data acquisition method is to extract it from the GDR product of CNES (Centre National d’Etudes Spatiales). The VTEC below the satellite height is inverted by using Jason-2 / 3 altimetry data Jas2 / 3 , and the calculation formula is:

[0087]

[0088] Where, f Ku = 13.6 GHz, f C = 5.3 GHz, h Ku , h C are the distances from the satellite to the sea surface measured by Ku and C bands respectively, D Ku , D C are the sea state deviation corrections of Ku and C bands respectively, B Ku , B C are the instrument deviation corrections of Ku and C bands respectively.

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

[0090] (4.2) DORIS data acquisition is to obtain the ranging code and phase observation value in RINEX format through the DORIS receiver, and then to inverse STEC by using the DORIS satellite dual-frequency data, and its expression is as follows:

[0091]

[0092] In the formula, STEC is the relative STEC, f1=2036.25 MHz, f2=401.25 MHz, λ1 and λ2 are the respective corresponding wavelengths, φ1 and φ2 are the phase observation values, and ΔD=D1-D2 contains the integer ambiguity in the two frequency phase observation values and the dual-frequency phase observation value bias.

[0093] In order to obtain the absolute TEC, the GIM VTEC inverted by the ground-based data is needed to realize the conversion from the relative TEC to the absolute TEC, and the calculation flow chart is as shown in Figure 4 The specific steps are as follows: the GIM VTEC at the DORIS ionospheric piercing point is inverted by the ground-based GNSS observation data; the GIM VTEC is projected to the signal propagation path by the SLM projection function to obtain the GIM STEC, and then the average value of the difference between the GIM STEC and the relative DORIS STEC in each DORIS observation arc segment is calculated; then the average value is added to the DORIS STEC value to obtain the absolute DORIS STEC, and finally the DORIS STEC is projected by using the SLM projection function to obtain the DORIS VTEC. It should be noted that since the average value already contains the TEC part above the orbit height of the DORIS satellite, the finally obtained DORIS VTEC covers the entire ionospheric height. Specifically, the DORIS system inversion VTEC flow chart is as shown in Figure 4 .

[0094] (4.3) The VTEC is obtained by using the COSMIC occultation data, and 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) Global ionospheric model adopts spherical harmonic function model, in which the order of spherical harmonic function model adopted is 15x15 order, the time resolution is 2 hours, the latitude direction interval is 2.5°, and the longitude direction interval is 5°. The spherical harmonic function expression is as follows:

[0097]

[0098] In the formula, β is the geomagnetic latitude at IPP, n and m are the order of spherical harmonic function respectively, n max is the maximum order of expansion, are model parameters to be solved respectively, is the normalized Lagrange function of n order m times, MC(n, m) is the regularization function, and δ om represents the Kronecker function. s is the IPP longitude in the daily fixed system, and the calculation formula is: s = λ-λ0≈UT+λ-π, wherein 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] In the formula, VTEC(β, s) is expressed by spherical harmonic function; VTEC BP-IRI , VTEC Jas2 / 3 , VTEC DORIS , and VTEC COSMIC respectively represent BP-IRI2016 VTEC, Jason-2 / 3 altimetric satellite VTEC, DORIS system VTEC, and COSMIC occultation system VTEC; DCB r , DCB s respectively represent the hardware delay of the receiver and the hardware delay of the satellite; and bias represents the systematic bias of various types of observation data systems.

[0102] (5.3) The normal equation is formed in the form of superposition of normal equations for different source observation data, which is represented as:

[0103]

[0104] In the formula, m represents m types of observation values, which include GNSS data, BP-IRI-2016 data, DORIS data, COSMIC data and Jason-2 / 3 data; L i represents various types of observation vectors, B i represents the design matrix, and X represents the spherical harmonic coefficients to be estimated. P i represents the weight matrix of different types of observation values; the weight of the same type of observation data is determined according to the satellite elevation angle, and p i= sin 2 (ele), where ele is the satellite elevation angle; denotes the variance of each type of observation data, which is calculated according to the Helmert variance component estimation.

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

[0106] (6.1) The adjustment value of the to-be-estimated parameter obtained from the normal equation of step 5.3 is expressed as:

[0107]

[0108] The definitions of the parameters 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 different systems of observation data: P1, P2, P3, P4, and P5; and to determine the covariance of the to-be-estimated parameters of GNSS data, IRI data, DORIS data, Jason-2 / 3 data, and COSMIC data, i.e., the prior variance, which is expressed as: The diagonal elements of D represent the variances of the estimated spherical harmonic coefficients.

[0110] (6.3) First-order least squares adjustment is performed, and the residual vector corresponding to the observation data is V i , which satisfies the relationship The weighted residual sum of squares V i corresponding to each type of observation data is calculated. T i V i , which is used for subsequent variance component estimation.

[0111] (6.4) Since the Helmert variance component estimation formula is relatively complex, the simplified formula for variance component estimation is often used in practice to calculate the posterior unit weight variance of each type of observation value, which is expressed as:

[0112]

[0113] In the formula, n i is the degree of freedom.

[0114] (6.5) The weights are then determined according to the following formula:

[0115]

[0116] In the formula, i = 1, 2,..., 5; c is an arbitrary constant, and is generally selected as a certain value of .

[0117] (6.6) Repeat the adjustment-variance component estimation-weighting until the posterior unit weight variances of each type of observations are equal, i.e.

[0118] (6.7) After the weight ratio between different system observations is determined by the Helmert variance component estimation, the spherical harmonic coefficients X can be solved according to the formula of the parameters to be estimated in step 6.1, and then the global ionospheric VTEC grid at the current epoch can be calculated by combining the spherical harmonic coefficient model in step 5.1.

[0119] The above only discloses a preferred embodiment of the present application, and of course cannot limit the scope of the present application, and those skilled in the art can understand that all or part of the processes of the above-mentioned embodiments are implemented, and equivalent changes made according to the claims of the present application still belong to the scope covered by the present application.

Claims

1. A global ionospheric inversion method based on BP neural network fusion of multi-source data, characterized in that, Comprise the following steps: Step 1: preprocessing of multi-source multi-system GNSS observation data, through the inversion of observation data of each system oblique ionospheric delay and through ephemeris data to calculate the ionospheric piercing point coordinates of each system, projection function is projected to the vertical direction of ionospheric piercing point, and VTEC is obtained; The execution process of step 1 comprises the following steps: Step 1.1: ground-based multi-system GNSS observation data preprocessing, the combination of MW method and ionospheric residual method is used to detect and repair the observation data; The sampling rate is set to 30 seconds, and the cut-off elevation angle is 20°; Step 1.2: Jason-2 / 3 data preprocessing of ocean altimetry satellite, obtaining ionospheric correction data of Jason-2 / 3 satellite altimeter through GDR product file, removing data on ocean and ice surface after obtaining TEC, and then performing 25 seconds 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: obtaining ground-based dual-frequency pseudorange and carrier phase observation data, without considering the influence of multipath effect and observation noise, the pseudorange observation value and carrier phase observation value at two frequencies are subtracted to form the observation value of non-geometric distance combination, and the carrier phase smoothing pseudorange algorithm is used to introduce smoothing quantity to solve STEC; Step 1.5: selecting ionospheric single-layer model height as 450km, and selecting commonly used SLM projection function for projection to obtain ionospheric VTEC; Step 2: taking the international reference ionosphere IRI-2016 model as the background model, using BP neural network to fuse the IRI-2016 model and the ionospheric VTEC obtained by ground-based inversion to obtain ΔVTEC; Step 3: obtaining the ΔVTEC to be solved by Kriging interpolation method, and obtaining the VTEC to be solved by taking the IRI-2016 model to be solved as the background model, and combining the ground-based inversion VTEC to obtain high-precision VTEC; Step 4: through the Ku and C band dual-frequency data of ocean altimetry satellite, the VTEC data of ocean satellite track point is inversed, and through the DORIS dual-frequency observation data, the absolute VTEC data and COSMIC occultation system are obtained VTEC data; Step 5: using spherical harmonic function to establish global ionospheric estimation observation equation combined with the VTEC data obtained in step 4, and constructing the equation of law combined with the observation equation; Step 6: estimating the weight of observation data of different sources and different systems by Helmert variance component estimation method, and solving the model to be estimated parameter, and then obtaining the global ionospheric VTEC grid data; The process of determining the weight of observation data of different sources and different systems by Helmert variance component estimation method in step 6 comprises 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: The first adjustment is made to obtain the V i T P i V i ; Step 6.3: the simplified formula of variance component estimation is used to calculate the posteriori unit weight variance of each type of observation value: Step 6.4: Determine the post-filter variance of different source observations based on the post-filter unit weight variance Step 6.5: The adjustment, variance component estimation, and post-adjustment variance calculation are iterated until the post-adjustment unit weight variances for each type of observation are equal, i.e.

2. The global ionospheric inversion method based on BP neural network fusion of multiple sources according to claim 1, characterized in that, The execution process of step 2 is specifically that the BP neural network generates a predicted AVTEC by taking the VTEC value provided by the IRI model and the VTEC value obtained by ground-based inversion measurement as a forward transmission signal, then calculates a prediction error, and adjusts the network weight through back propagation of the error, so that the prediction result is closer to the expected output AVTEC.

3. The global ionospheric inversion method based on BP neural network fusion of multi-source data according to claim 2, characterized in that, In step 3, the Kriging interpolation method is used to determine the weight value by first performing optimal and unbiased estimation on the variable value, and then solving the electron content difference value.

4. The global ionospheric inversion method based on BP neural network fusion of multi-source data according to claim 3, characterized in that, In step 4, the data used by the ocean altimetry satellite is the Jason-2 and Jason-3 ocean altimetry satellite data, and the data acquisition method is to obtain it from the GDR product of CNES.

5. The global ionospheric inversion method based on BP neural network fusion of multi-source data according to claim 4, characterized in that, The DORIS data acquisition method in step 4 includes the following steps: Step 4.1: The GIM VTEC at the DORIS ionospheric piercing point is inverted from the ground-based GNSS observation data; Step 4.2: The GIM VTEC is projected to the signal propagation path by the SLM projection function to obtain the GIM STEC, and then the average value of the difference between the GIM STEC and the relative DORIS STEC in each DORIS observation arc segment is calculated; Step 4.3: The average value is added to the DORIS STEC value to obtain the absolute DORIS STEC, and finally the DORIS STEC is projected by using the SLM projection function to obtain the DORIS VTEC; Step 4.4: The VTEC is obtained by using the COSMIC occultation data, wherein the VTEC data is extracted from the ionPrf product provided by the University Corporation for Atmospheric Research (UCAR) in the United States.

6. The global ionospheric inversion method based on BP neural network fusion of multi-source data according to claim 5, characterized in that, The observation equation establishment process in step 5 is specifically that the normal equation is formed by using the normal equation superposition method for different source observation data, and the variance of each type of observation data is calculated according to the Helmert variance component estimation.

Images