A method and device for calculating the electron density of the ionosphere in a region, and a computer device

Abstract

The application discloses a regional ionospheric electron density calculation method and device and computer equipment. The method comprises the following steps: according to all absolute STEC values, all station coordinate information and the coordinate information of the GNSS satellite, the signal ray intercept length through each grid in the to-be-tomographic region is calculated, and a tomographic observation equation is established; according to the position of each grid point, an air-based occultation ionospheric electron density empirical model is used to generate the prior initial value of the electron density of each grid point at the tomographic moment; and the prior initial value of the grid point electron density is iteratively corrected according to the tomographic observation equation, so that the regional ionospheric electron density distribution is obtained. According to the application, the constrained ambiguity float solution is used as the ambiguity initial value of the next epoch, so that the influence of the decimal part of the ambiguity parameter is eliminated, and the initial value of the air-based ionospheric empirical model is greatly accelerated, and the convergence speed is greatly accelerated.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of GNSS atmosphere monitoring, and particularly relates to a regional ionospheric electron density calculation method and device and computer equipment. BACKGROUND

[0002] When GNSS satellite signals pass through the ionosphere, the propagation speed and the propagation path of the signals will change. The degree of change in the propagation speed mainly depends on the signal frequency and the electron density in the ionosphere. The slight bending of the propagation path has little effect on the ranging result and is generally ignored; however, the delay of the signal propagation time will seriously affect the positioning accuracy. However, because of the dispersion characteristics of the ionosphere, the time delay of the satellite signal is related to the electron content in the ionosphere and the frequency of the satellite signal, and therefore, the ionospheric delay can be detected by using multi-frequency GNSS observation data. With the development of space-based and ground-based GNSS, sufficient multi-source observation data has been accumulated, and how to process and fuse the multi-source observation data has become a major problem in current research. In addition, in the traditional precise point positioning (PPP) positioning solution, the ambiguity parameter is a real number solution due to the influence of the phase fractional bias. This limits the convergence speed of PPP and makes it unable to be widely used in real-time positioning services. How to eliminate the influence of the fractional part of the ambiguity parameter and restore its integer characteristics is also a difficult problem in current research. SUMMARY

[0003] The embodiments of the present application provide a regional ionospheric electron density calculation method and device and computer equipment, an empirical model of the space-based ionosphere is constructed as a background field as an initial value to ensure the convergence speed of PPP, and a tomographic algorithm is used to realize joint space / ground-based GNSS detection of the ionospheric electron density.

[0004] To achieve the above object, a first aspect of the embodiments of the present application provides a regional ionospheric electron density calculation method, comprising:

[0005] performing grid division on the to-be-tomographed region to obtain the positions of each grid point;

[0006] calculating the absolute STEC values between each ground station and the GNSS satellite at the tomo-graphing time by using non-difference non-combined PPP according to the GNSS observation data of each ground station and the terminal operation data, and obtaining the coordinate information of each station and the coordinate information of the GNSS satellite;

[0007] calculating the intercept length of the signal ray passing through each grid in the to-be-tomographed region according to all the absolute STEC values, the coordinate information of all the stations and the coordinate information of the GNSS satellite, and establishing a tomographic observation equation;

[0008] According to the positions of the grid points, an empirical model of ionospheric electron density based on occultation is used to generate initial values of electron density at the grid points at the tomography time;

[0009] The initial values of electron density at the grid points are iteratively corrected according to the tomography observation equation to obtain the distribution of ionospheric electron density in the region.

[0010] In a possible implementation of the first aspect, the grid division of the region to be tomographed specifically includes:

[0011] The region to be tomographed is divided into grids in three dimensions of latitude, longitude and height, and each grid point is assigned a number.

[0012] In a possible implementation of the first aspect, the absolute STEC values between each ground station and GNSS satellite at the tomography time are calculated by using non-difference non-combination PPP according to the GNSS observation data of each ground station and the terminal operation data, and specifically includes:

[0013] The observation equation of non-difference non-combination PPP and MW combined observation values are obtained according to the GNSS observation data of each ground station and the terminal operation data.

[0014] The MW combined observation values are processed by multi-epoch smoothing to obtain wide-lane ambiguities and perform FCB correction, and the inter-satellite single-difference wide-lane ambiguity is constructed.

[0015] The inter-satellite single-difference wide-lane ambiguity is rounded and fixed to obtain fixed inter-satellite single-difference wide-lane ambiguity.

[0016] The inter-satellite single-difference narrow-lane ambiguity is obtained by using the fixed inter-satellite single-difference wide-lane ambiguity and I F floating ambiguity, and narrow-lane FCB correction is performed.

[0017] The inter-satellite single-difference narrow-lane ambiguity is processed by using the LAMBDA algorithm to obtain fixed inter-satellite single-difference narrow-lane ambiguity.

[0018] The fixed inter-satellite single-difference wide-lane ambiguity and the fixed inter-satellite single-difference narrow-lane ambiguity are linearly combined to obtain ionosphere-free combined ambiguity.

[0019] The ionospheric TEC observation value is obtained in combination with the ionosphere-free combined ambiguity and the observation equation of non-difference non-combination PPP.

[0020] In a possible implementation of the first aspect, after the inter-satellite single-difference wide-lane ambiguity is rounded and fixed to obtain the fixed inter-satellite single-difference wide-lane ambiguity, the method further includes:

[0021] The correctness of the fixed inter-star single-difference wide-lane ambiguity is verified by using a probability judgment function, and the correct fixed inter-star single-difference wide-lane ambiguity is stored.

[0022] In a possible implementation manner of the first aspect, after the inter-star single-difference narrow-lane ambiguity is processed by using the LAMBDA algorithm to obtain the fixed inter-star single-difference narrow-lane ambiguity, the method specifically comprises the following steps.

[0023] The correctness of the fixed inter-star single-difference narrow-lane ambiguity is verified by using a ratio value, and the correct fixed inter-star single-difference narrow-lane ambiguity is stored.

[0024] In a possible implementation manner of the first aspect, the signal ray intercept length through each grid in the tomographic region to be tomographed is calculated according to all absolute STEC values, coordinate information of all stations and coordinate information of the GNSS satellite, and a tomographic observation equation is established, which specifically comprises the following steps.

[0025] A plurality of intersection points of the signal ray through the height plane, the longitude plane and the latitude plane are respectively obtained;

[0026] The plurality of intersection points are sorted in ascending order, and the distance between adjacent two points is obtained as the signal ray intercept through the corresponding grid from low to high;

[0027] According to the grid point number, the signal ray intercept through the corresponding grid is assigned to the intercept matrix by the grid;

[0028] The tomographic observation equation is established according to the intercept matrix.

[0029] In a possible implementation manner of the first aspect, the establishment process of the air-based occultation ionospheric electron density empirical model specifically comprises the following steps.

[0030] Key parameters are extracted from GNSS observation data of the air-based station based on least squares;

[0031] The ionospheric key parameter matrix is reorganized according to the key parameters after the space-time grid is divided;

[0032] The ionospheric key parameter matrix is subjected to space-time empirical orthogonal decomposition to obtain an empirical orthogonal basis function and a time coefficient;

[0033] The time coefficient corresponding to each order mode is subjected to Fourier series fitting to establish an ionospheric key parameter empirical model;

[0034] The ionospheric key parameter empirical model value is input into a Vary-Chapman model to calculate a first electron density;

[0035] The second electron density is calculated by the form of exponential function extrapolation of the COSMIC series of low-orbit satellite orbit height above the part of the electron density, and the second electron density is obtained.

[0036] The first electron density and the second electron density are added to obtain the electron density distribution of the space-based occultation ionospheric electron density empirical model.

[0037] In a possible implementation manner of the first aspect, the iterative correction of the grid point electron density prior initial value according to the tomographic observation equation specifically includes:

[0038] According to the tomographic observation equation, the iterative calculation is continuously performed until the standard deviation of the residual between the STEC of the iterative calculation and the measured STEC is less than a preset threshold.

[0039] The second aspect of the embodiment of the application provides a regional ionospheric electron density calculation device, which comprises:

[0040] The division module is configured to divide the tomographic region into grids to obtain the positions of the grid points.

[0041] The calculation module is configured to calculate the absolute STEC values between each ground station and a GNSS satellite at a tomographic time according to the GNSS observation data of each ground station and the end running data, and obtain the coordinate information of each station and the coordinate information of the GNSS satellite.

[0042] The tomographic module is configured to calculate the intercept length of the signal ray passing through each grid in the tomographic region according to all the absolute STEC values, the coordinate information of all the stations and the coordinate information of the GNSS satellite, and establish a tomographic observation equation.

[0043] The initial module is configured to generate the grid point electron density prior initial value at the tomographic time according to the positions of the grid points by using the space-based occultation ionospheric electron density empirical model.

[0044] The iteration module is configured to perform iterative correction on the grid point electron density prior initial value according to the tomographic observation equation to obtain the regional ionospheric electron density distribution.

[0045] The second aspect of the embodiment of the application provides a computer device comprising a processor and a memory, wherein the memory is configured to store a computer program, and the computer program is configured to be executed by the processor to implement the above-mentioned regional ionospheric electron density calculation method.

[0046] Compared with the prior art, the regional ionospheric electron density calculation method, device and computer equipment provided by the embodiment of the present application can eliminate the influence of the decimal part of the ambiguity parameter, restore the integer characteristic, improve the extraction accuracy of the ionospheric observation value, adopt the strategy of obtaining better PPP solving results based on the satellite attitude quaternion and phase bias products provided by IGS / MGEX to solve the processing problem of satellite yaw attitude, and adopt the fix-and-hold mode as the fixing method of the integer ambiguity in PPP, that is, the ambiguity of the ionosphere-free (IF) combination is uniformly decomposed into wide-lane and narrow-lane ambiguities for fixing in sequence, the IF combination ambiguity in the undifferenced and uncombined PPP model is obtained from the original ambiguity of the dual-frequency, after the fixed IF combination ambiguity is obtained, the ambiguity is used as a virtual observation value and a virtual observation equation is constructed to strongly constrain the Kalman filtering state, then the floating point solution of the constrained ambiguity is used as the initial value of the ambiguity of the next epoch, so that the influence of the decimal part of the ambiguity parameter is eliminated, and the convergence speed is greatly accelerated as the initial value of the empirical ionospheric model: the IonPrf product provided by COSMIC is used, the Vary-Chapman and top-layer exponential electron density model are used, the empirical orthogonal decomposition method and Fourier series are used to construct a regional ionospheric empirical model considering the longitude, latitude, local time, annual day and solar activity index. Then the constructed space-based ionospheric empirical model is used as a background field, and the tomographic algorithm is used to realize the joint space-based / ground-based GNSS detection of ionospheric electron density.

[0047] In summary, the embodiment of the present application uses high-precision phase observation values, accurately estimates the ambiguity value, and fixes the ambiguity parameter, so as to improve the accuracy of extracting the ionospheric total electron content observation value by using the undifferenced and uncombined precise point positioning technology. The satellite attitude quaternion and phase bias products provided by IGS / MGEX are used to obtain better PPP solving results and accelerate the fixing speed of the ambiguity. BRIEF DESCRIPTION OF DRAWINGS

[0048] Figure 1 is a flowchart of the regional ionospheric electron density calculation method provided by the embodiment of the present application;

[0049] Figure 2 is a flowchart of the process of fixing the ambiguity by using the FCB product in the embodiment of the present application;

[0050] Figure 3 is a flowchart of the process of establishing the regional space-based ionospheric empirical model by using the occultation data in the embodiment of the present application. DETAILED DESCRIPTION

[0051] With reference to the drawings of the embodiments of the present application, the technical solutions in the embodiments of the present application will be described clearly and completely. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments of the present application. Based on the embodiments of the present application, all other embodiments obtained by a person of ordinary skill in the art without creative work fall within the protection scope of the present application.

[0052] Please refer to Figure 1 An embodiment of the present application provides a method for calculating regional ionospheric electron density, comprising:

[0053] S10, performing grid division on the to-be-tomographic region to obtain positions of grid points.

[0054] S11, according to GNSS observation data of each ground station and each end running data, using non-difference non-combination PPP to calculate absolute STEC values between each ground station and GNSS satellites at a tomographic moment, and obtaining coordinate information of each station and coordinate information of GNSS satellites.

[0055] S12, according to all absolute STEC values, coordinate information of all stations and coordinate information of the GNSS satellites, calculating an intercept length of a signal ray passing through each grid in the to-be-tomographic region, and establishing a tomographic observation equation.

[0056] S13, according to the positions of the grid points, using an air-based occultation ionospheric electron density empirical model to generate prior initial values of electron density of each grid point at the tomographic moment.

[0057] S14, performing iterative correction on the prior initial values of electron density of the grid points according to the tomographic observation equation to obtain regional ionospheric electron density distribution.

[0058] The specific content of the embodiment is that total electron content of the ionosphere is obtained based on non-difference non-combination precise point positioning technology, satellite end pseudorange code bias is corrected by using DCB (Differential Code Bias) products provided by the Chinese Academy of Sciences on the basis of GNSS original observation equation, and the differential code bias of the receiver end is taken as an unknown parameter, and Kalman filtering is used for estimation, so that absolute path electron content is obtained.

[0059] However, the ambiguity parameter in the non-combination PPP method is a real number solution, and it usually takes a long convergence time to achieve sufficient accuracy, cannot fully utilize the observation data, and cannot use the integer constraint condition of ambiguity, thereby limiting the extraction accuracy of ionospheric TEC. If high-precision phase observations are directly used and the ambiguity value is accurately estimated and fixed, the extraction accuracy of ionospheric total electron content (TEC) observations will be greatly improved. In order to eliminate the influence of the decimal part of the ambiguity parameter and restore its integer characteristics, while improving the extraction accuracy of ionospheric observations, the strategy of obtaining better PPP solution results based on satellite attitude quaternions and phase bias products provided by IGS / MGEX is adopted to solve the problem of satellite yaw attitude processing, and the fix-and-hold mode is used as the fixing method of the integer ambiguity in PPP, that is, the ambiguity of the ionosphere-free (IF) combination is uniformly decomposed into wide-lane and narrow-lane ambiguities for fixing in turn. The IF combination ambiguity is obtained from the original dual-frequency ambiguity, and after the fixed IF combination ambiguity is obtained, it is used as a virtual observation value and a virtual observation equation is constructed to strongly constrain the Kalman filter state. The mean error of the virtual observation value is set to an empirical value of 0.003 m, and then the floating-point solution of the constrained ambiguity is used as the initial value of the ambiguity at the next epoch. Using the IonPrf product provided by COSMIC, based on the Vary-Chapman and top-layer exponential electron density model, an empirical model of regional ionosphere considering latitude, longitude, local time, annual day and solar activity index is constructed by using the empirical orthogonal decomposition method and Fourier series. Then, the constructed space-based ionospheric empirical model is used as the background field, and the method of joint space / ground-based GNSS probing of ionospheric electron density is realized by tomographic algorithm.

[0060] Illustratively, the grid division of the tomographic region to be tomographed specifically includes:

[0061] The to-be-tomographed region is grid divided in three dimensions of latitude, longitude and height, and each grid point is assigned a number.

[0062] Illustratively, the absolute STEC values between each ground station and GNSS satellite at the tomographic time are calculated by using non-difference non-combination PPP according to the GNSS observation data of each ground station and the terminal running data, specifically including:

[0063] The observation equation of non-difference non-combination PPP and the MW combination observation value are obtained according to the GNSS observation data of each ground station and the terminal running data.

[0064] Performing multi-epoch smoothing on the MW combination observations, obtaining wide-lane ambiguities and performing FCB correction, constructing inter-satellite single-difference wide-lane ambiguities;

[0065] Performing rounding fixing on the inter-satellite single-difference wide-lane ambiguities, obtaining fixed inter-satellite single-difference wide-lane ambiguities;

[0066] Obtaining inter-satellite single-difference narrow-lane ambiguities by using the fixed inter-satellite single-difference wide-lane ambiguities and IF float ambiguities, and performing narrow-lane FCB correction;

[0067] Performing processing on the inter-satellite single-difference narrow-lane ambiguities by using LAMBDA algorithm, obtaining fixed inter-satellite single-difference narrow-lane ambiguities;

[0068] Performing linear combination on the fixed inter-satellite single-difference wide-lane ambiguities and the fixed inter-satellite single-difference narrow-lane ambiguities, obtaining ionosphere-free combination ambiguities;

[0069] Combining the ionosphere-free combination ambiguities and the observation equation of undifferenced non-combination PPP, obtaining corresponding ionospheric TEC observation values.

[0070] Firstly, the value of TEC is solved by using the undifferenced non-combination PPP method, and since the precise clock difference provided by each IGS analysis center will have satellite-end uncorrected pseudo-range hardware delay bias, the observation equation of undifferenced non-combination PPP is expressed as follows:

[0071]

[0072] Wherein

[0073]

[0074] In the above formula (1) (2), and respectively represent pseudo-range and carrier phase observation values (m); the subscript i represents frequency index; is the geometric distance from the station to the satellite (m). Since the precise satellite clock difference product provided by IGS is generally estimated based on double-frequency IF observation values, there is a linear combination term of satellite-end double-frequency uncorrected pseudo-range hardware delay (Uncalibrated Code Delay, UCD) in the precise satellite clock difference product, that is, In order to unify, the receiver clock difference, ionospheric delay and ambiguity parameters are adjusted, and are respectively the reconstructed receiver clock difference, ionospheric delay on the slant path and ambiguity parameters; is the frequency-dependent ionospheric scaling factor; d r,i and are respectively uncorrected pseudo-range hardware delays of the receiver and the satellite (m); br,i and respectively represent the uncalibrated phase hardware delay (UPD) of the receiver and satellite end (m); λ i is the wavelength at the corresponding frequency (m); and respectively represent the unmodeled errors of the pseudorange and carrier phase observations (m).

[0075] In order to improve the extraction accuracy of undifferenced and uncombined PPP and the rapid fixing of integer ambiguity, the FCB products of GPS wide-lane and narrow-lane released by SGG (School of Geodesy and Geomatics, Wuhan University) since 2015 are used. FCB (Fractional Cycle Bias) is the fractional part of the uncalibrated phase hardware delay UPD, and accurate and reliable FCB estimation requires high-precision parameter estimation. First, the original dual-frequency ambiguity is converted into the linear combination of wide-lane and narrow-lane ambiguity for parameter estimation, as shown in equation (3):

[0076]

[0077] When the user uses the FCB product and its corresponding precise orbit clock error product to fix the PPP ambiguity, first, the MW combined observation (the pseudorange observation should be corrected by the satellite end DCB to keep consistent with the FCB product end) of equation (4) is processed by multi-epoch smoothing to reduce the influence of observation noise and multipath error, and the wide-lane ambiguity containing the satellite end and receiver end hardware delay is obtained, as shown in equation (4):

[0078]

[0079] Then, the satellite with the highest elevation angle is taken as the reference star to construct the inter-satellite single-difference wide-lane ambiguity, and the wide-lane FCB product is corrected; due to the longer wavelength of the wide-lane, the rounding method can be directly used for fixing, and the probability judgment function is used to check whether the fixing is correct; after the wide-lane is fixed, the inter-satellite single-difference narrow-lane ambiguity is obtained according to equation (5) using the fixed inter-satellite single-difference wide-lane ambiguity and IF floating ambiguity, and the narrow-lane FCB correction is performed, and then the LAMBDA algorithm is used to obtain the fixed inter-satellite single-difference narrow-lane ambiguity;

[0080]

[0081]

[0082] Finally, the more accurate ionosphere-free combined ambiguity is obtained by reducing the equation (6) to get the fixed solution. It is worth noting that the fixed narrow-lane ambiguity should be added to the FCB and substituted into the equation during the reduction process.

[0083]

[0084] The overall process is shown in Figure 2 .

[0085] Exemplarily, after the fixed inter-satellite single-difference wide-lane ambiguity is obtained by rounding the inter-satellite single-difference wide-lane ambiguity, the method further comprises:

[0086] The correctness of the fixed inter-satellite single-difference wide-lane ambiguity is verified by using a probability judgment function, and the correct fixed inter-satellite single-difference wide-lane ambiguity is stored.

[0087] Exemplarily, after the fixed inter-satellite single-difference narrow-lane ambiguity is obtained by processing the inter-satellite single-difference narrow-lane ambiguity using the LAMBDA algorithm, the method specifically comprises:

[0088] The correctness of the fixed inter-satellite single-difference narrow-lane ambiguity is verified by using a ratio value, and the correct fixed inter-satellite single-difference narrow-lane ambiguity is stored.

[0089] The observation equation for wide-lane and narrow-lane FCB estimation of different satellite navigation systems can be uniformly written as:

[0090]

[0091] Considering the linear correlation of the FCB at the receiver end and the satellite end, the sum of the wide-lane or narrow-lane FCB of all satellites of the same system is generally set to zero to solve the rank defect problem of equation (7).

[0092] In the use of FCB products and its corresponding precise track clock difference products to fix PPP ambiguity, firstly, the MW combination observation value (pseudo-range observation value should be corrected by satellite end DCB, consistent with FCB product end) of formula (4) is processed by multi-epoch smoothing to reduce the influence of observation noise and multi-path error, and wide-lane ambiguity containing satellite end and receiver end hardware delay is obtained; then the satellite with the highest elevation angle is taken as the reference star to construct inter-satellite single-difference wide-lane ambiguity, and the wide-lane FCB product is corrected; since the wide-lane wavelength is long, the rounding method can be directly used for fixing, and whether the fixing is correct is verified according to the probability judgment function; after the wide-lane is fixed, the inter-satellite single-difference narrow-lane ambiguity is obtained according to formula (5) by using the fixed inter-satellite single-difference wide-lane ambiguity and IF floating ambiguity, and the narrow-lane FCB is corrected, and the LAMBDA algorithm is used to obtain the fixed inter-satellite single-difference narrow-lane ambiguity; finally, more accurate ionosphere-free combination ambiguity is obtained by formula (3), and the fixed solution is obtained. In the reduction process, the fixed narrow-lane ambiguity is added to the FCB and substituted into the formula. The DCB of the satellite end is corrected, and the DCB of the receiver end is estimated as a parameter in the ambiguity parameter in the non-combination PPP model, so that the reconstructed IF ambiguity is consistent with the hardware delay in the IF ambiguity in the IF-PPP model. The fixed wide-lane and narrow-lane ambiguities are substituted into the ionosphere-free L4 combination observation value of formula (7), and the corresponding ionospheric TEC observation value is obtained, and then the ionospheric model can be established by using the observation value.

[0093]

[0094] Exemplarily, the intercept length of the signal ray passing through each grid in the to-be-tomographic region is calculated according to all absolute STEC values, coordinate information of all stations and coordinate information of the GNSS satellite, and a tomographic observation equation is established, and specifically, the method comprises the following steps:

[0095] A plurality of intersection points of the signal ray passing through the height plane, the longitude plane and the latitude plane are respectively obtained;

[0096] The plurality of intersection points are sorted in ascending order, and the distance between adjacent two points is taken as the intercept of the signal ray passing through the corresponding grid from low to high;

[0097] According to the grid point number, the intercept of the signal ray passing through the corresponding grid is assigned to the intercept matrix by the grid;

[0098] The tomographic observation equation is established according to the intercept matrix.

[0099] Exemplarily, the establishment process of the air-based occultation ionospheric electron density empirical model is specifically as follows:

[0100] Extracting key parameters from GNSS observation data of space-based observation station based on least square method;

[0101] Dividing space-time grid, reorganizing ionosphere key parameter matrix according to the key parameters;

[0102] Performing space-time empirical orthogonal decomposition on the ionosphere key parameter matrix to obtain empirical orthogonal basis function and time coefficient;

[0103] Performing Fourier series fitting on the time coefficient corresponding to each order mode to establish an empirical model of ionosphere key parameters;

[0104] Inputting the empirical model value of ionosphere key parameters into Vary-Chapman model to calculate the first electron density;

[0105] Calculating the electron density above the orbit height of COSMIC series low-orbit satellite by exponential function extrapolation to obtain the second electron density;

[0106] Adding the first electron density and the second electron density to obtain the electron density distribution of the space-based occultation ionosphere electron density empirical model.

[0107] Referring to Figure 3 , the process of establishing a regional space-based ionosphere empirical model using occultation data is described next. First, the COSMIC provided electron density profile data is preprocessed, and the top electron density gradient Slope, density average deviation MD and noise factor delta of the profile within the 420-490km altitude range are used as data rejection indicators. Among them, the density average deviation is defined as follows:

[0108]

[0109] In the above formula, N represents the total number of electron densities in the original profile, i represents the i-th electron density, represents the sliding average value of the electron density, and the sliding window size is set to 5. When MD>1, it represents that the electron density profile data is obviously disturbed and should be rejected. The screening indicator Slope is used to remove some profile data without obvious F2 layer peak electron density, which is defined as follows:

[0110]

[0111] Ionosphere F layer and top layer only need peak density NmF2, peak height HmF2, height at peak height H mand the slope a1 and a2 above and below the peak height can recover the reconstructed ionospheric electron density distribution. The key steps of constructing the regional space-based ionospheric empirical model using the OnPrf product data provided by COSMIC-1 / 2 are as follows:

[0112] 1) Data preprocessing, key parameters are extracted based on least squares;

[0113] 2) Spatial grid division, reorganize the ionospheric key parameter matrix. First, take 10, 20, 30, …, 360 as the center, and divide the original data into 36 data sets with a data window size of 20 days. Then, in each year accumulated day data set, further division is made. Given that the ionosphere changes little in a small range in the longitude direction, the data window is centered at 80°E, 88°E, 96°E, …, 136°E, and 16° is used for data windowing. Each year accumulated day data set is divided into 8 longitude data sets. Then, according to latitude, local time, and F10.7a, the longitude data set is further divided, with the division scale being 1° (10°N-50°N), 1 hr (1-24 h), and 20 sfu (70-150 sfu), respectively. Thus, 41x24x5 data are divided in each longitude data set. Considering that COSMIC-1 / 2 cannot guarantee data coverage at the above-mentioned latitude, local time, and F10.7a grid points, a spherical harmonic function is used with local time, latitude, and F10.7a as parameters to least square fit each longitude data set and then interpolate the data at the above-mentioned latitude, local time, and F10.7a grid points. The interpolation function expression is as follows:

[0114]

[0115] wherein,

[0116] The above formula f (lat, lt, F10.7) can represent NmF2, HmF2, H m and a1 and a2 any one key parameter, indicates the Legendre coefficient, N indicates the maximum order, and after multiple tests, a 6-order function is finally used. As known from the previous chapter, the solar activity index F10.7a has a linear relationship with the ionosphere size, so the coefficients of the spherical harmonic function model are expressed as a first-order linear relationship with F10.7a. The fitting coefficients h mn , h ′ mn , j mn and j ′ mn can be solved by least square estimation, and grid data interpolation is realized.

[0117] 3) Perform an empirical orthogonal decomposition of the ionospheric key parameter matrix in space and time to obtain the empirical orthogonal basis functions and time coefficients. Based on step 2), the 5-dimensional data structure is reorganized into a 2-dimensional dataset, where the "row data" of the matrix is ​​arranged in the grid size of latitude, local time, and longitude, and the "column data" is arranged in the grid size of annual accumulation day and F10.7a. Then, the above 2-dimensional space and time data are empirically decomposed according to the following form:

[0118]

[0119] In the above formula, u represents the maximum mode used in the empirical orthogonal decomposition. After the above decomposition, the empirical orthogonal basis function E of each order mode can be obtained. k The corresponding time coefficient A k .

[0120] 4) Perform Fourier series fitting on the time coefficients corresponding to each mode order and establish empirical models of key ionospheric parameters. In order to take into account the annual and diurnal variation characteristics of the ionosphere, a 6th-order Fourier series fitting of the time coefficient A is used. k , establish the time coefficient A of each mode related to the annual cumulative day, local time and solar activity index F10.7a k The empirical model of the ionosphere is then used to reconstruct the key parameters of the ionosphere at this moment. First, we establish A k The formula for the 6th-order Fourier model with annual variation and annual cumulative days as parameters is as follows:

[0121]

[0122] Then q m and b m It is again a function of local time lt and F10.7, reflecting its diurnal variation characteristics, as shown below:

[0123]

[0124]

[0125] 5) The empirical model values ​​of key ionospheric parameters are input into the Vary-Chapman model for calculation. The electron density above the orbital altitude of the COSMIC series low-orbit satellites is calculated by exponential function extrapolation. The electron density distribution of the space-based ionospheric empirical model can be obtained by adding the two.

[0126]

[0127] In the above formula, is the electron density at the COSMIC low-orbit satellite orbital altitude, h orbH is the height of the orbit P H is the height of the orbit

[0128] Exemplarily, the iteration correction of the grid point electron density prior initial value according to the tomographic observation equation specifically includes:

[0129] According to the tomographic observation equation, the iteration calculation is continuously performed until the standard deviation of the residual between the iteration calculated STEC and the measured STEC is less than a preset threshold.

[0130] The matrix form expression of the tomographic ionospheric electron density observation equation based on the pixel method is:

[0131] y = Ax + e (16)

[0132] Each row of the coefficient matrix A represents each ray, and each column represents the intercept size of the ray passing through the grid. The rest of the letters represent the same meaning as above. Therefore, the essence of establishing the tomographic observation equation is to determine the coefficient A matrix composed of the intercepts of the rays passing through each grid.

[0133] According to the spatial coordinates of the station (X C ,Y C ,Z C ) and the satellite (X S ,Y S ,Z S ), the spatial straight line equation of the signal ray in the WGS-84 spatial coordinate system can be determined:

[0134]

[0135] In a small range, the height surface can be regarded as an ellipsoid surface parallel to the earth reference ellipsoid. When the height surface is at a height of H from the ground, its ellipsoid surface equation can be expressed as follows:

[0136]

[0137] Where a and b are the long semi-axis and short semi-axis of the earth reference ellipsoid, respectively. By combining (17) and (18), the intersection point with the height surface can be obtained. In actual solving, the ray straight line equation can be combined with the lowest and highest height surface equations, respectively, to solve the two intersection points (P l and P h ) with the lowest and highest height surfaces. Then, it is determined whether the two intersection point coordinates are within the tomographic area. As long as any one of the intersection point coordinates is not within the tomographic area, it means that the ray does not completely pass through the tomographic area, and the ray should be discarded. When both intersection points are within the tomographic area, the ray straight line equation is combined with each layer height surface from low to high to solve the intersection points P l ,2,…,h .

[0138] According to the intersection P of the lowest and highest height planes l and P h The coordinate range can further determine the latitude and longitude planes through which the ray passes. First, the spatial rectangular coordinates of P l and P h are converted into geodetic latitude and longitude coordinates, respectively represented as For the latitude plane, it can be represented as a conical surface obtained by rotating the line connecting a spatial point and the center of the earth around the Z axis, and when the latitude is , its corresponding plane equation is in the following form:

[0139]

[0140] For the longitude plane, it can be represented as a plane passing through the Z axis, and when the longitude is λ, its corresponding plane equation is:

[0141] tanλ·X-Y=0 (20)

[0142] Therefore, in the latitude direction, as long as the latitude plane whose latitude is between and is found according to the corresponding latitude size of each latitude plane, and its corresponding latitude plane equation (20) is combined with the ray straight line equation (17), the intersection of the ray through each latitude plane can be obtained. After obtaining each intersection point of the height plane, the longitude plane and the latitude plane, then sort these intersection points from small to large according to the height size. After sorting, the distance between adjacent points is calculated, and the distance of each point is the intercept of the ray through the corresponding grid from low to high. Finally, according to the self-defined grid numbering rule, the intercept of each ray through the grid is assigned to the coefficient matrix A in turn, and for the grid which the ray does not pass through, it is assigned as 0.

[0143] There is a key problem in the pixel-based tomography, which is that the tomography equation matrix may be rank-deficient. This is mainly related to the lack of GNSS ground stations, uneven spatial distribution and the small number of observable satellites. The above factors result in that there is no ray passing through some grid in the tomography area, which leads to the non-uniqueness of the calculated grid electron density. To solve the ill-posed problem of tomography equation, the iterative algorithm is generally used. At the beginning of the iterative calculation, MART needs to assign an initial value to each ionospheric grid in the tomography area. The initial value is generally obtained by the priori empirical model of electron density or other detection methods. The iterative process is performed on each observation equation. Assuming that there are n observations in total, i.e. there are n GNSS rays in the tomography time, when n steps of iteration are completed, it is called one round of iteration. The correction of each iteration is based on the ratio of the STEC calculated by the electron density of the kth iteration result to the measured STEC, and then it is distributed to each grid electron density, so that it gradually converges. The correction formula of the kth iteration is:

[0144]

[0145] In the above formula, represents the (k+1)th iteration value of the jth grid, k represents the vector composed of the electron density of each grid in the kth step, represents the transpose of the row vector of the i-th row of the coefficient matrix A, and represents the inner product symbol, k represents the relaxation factor of the kth step. In actual iterative solution, the relaxation factor of each step is generally regarded as the same value, and its size is between 0 and 1. Through continuous iterative calculation, when the standard deviation of the residual between the STEC of the electron density iterative calculation and the measured STEC is less than a certain value, the iteration is terminated.

[0146] Compared with the prior art, the regional ionospheric electron density calculation method provided by the embodiment of the application eliminates the influence of the decimal part of the ambiguity parameter, restores the integer characteristic, improves the extraction accuracy of the ionospheric observation value, adopts the strategy of obtaining better PPP solving results based on the satellite attitude quaternion and phase bias products provided by IGSMGEX to solve the processing problem of satellite yaw attitude, and adopts the fix-and-hold mode as the fixing method of the integer ambiguity in PPP, that is, the ambiguity of the ionosphere-free (IF) combination is uniformly decomposed into wide-lane and narrow-lane ambiguities for fixing in sequence, the IF combination ambiguity in the undifferenced and uncombined PPP model is obtained from the original ambiguity of the dual-frequency, after the fixed IF combination ambiguity is obtained, it is used as a virtual observation value and a virtual observation equation is constructed to strongly constrain the Kalman filtering state, then the floating point solution of the constrained ambiguity is used as the initial value of the ambiguity of the next epoch, so that the influence of the decimal part of the ambiguity parameter is eliminated, and the convergence speed is greatly accelerated as the initial value of the empirical ionospheric model; the IonPrf product provided by COSMIC is used, based on the Vary-Chapman and top layer exponential electron density model, an empirical orthogonal decomposition method and a Fourier series are used to construct a regional ionospheric empirical model considering the longitude, latitude, local time, annual day and solar activity index. Then, the constructed space-based ionospheric empirical model is used as a background field, and the tomographic algorithm is used to realize the joint space / ground-based GNSS detection of the ionospheric electron density.

[0147] An embodiment of the present application provides a regional ionospheric electron density calculation device, which comprises a division module, a calculation module, a tomographic module, an initial module and an iteration module.

[0148] The division module is used for grid division on a to-be-tomographed region to obtain the positions of grid points.

[0149] The calculation module is used for calculating the absolute STEC values between each ground station and a GNSS satellite at a tomographic time by using undifferenced and uncombined PPP according to the GNSS observation data of each ground station and the terminal operation data, and obtaining the coordinate information of each station and the coordinate information of the GNSS satellite.

[0150] The tomographic module is used for calculating the intercept length of a signal ray passing through each grid in the to-be-tomographed region according to all the absolute STEC values, the coordinate information of all the stations and the coordinate information of the GNSS satellite, and establishing a tomographic observation equation.

[0151] The initial module is used for generating the prior initial value of the electron density of each grid point at the tomographic time by using a space-based occultation ionospheric electron density empirical model according to the positions of the grid points.

[0152] An iteration module is configured to iteratively correct the prior initial value of the electron density at the grid points according to the tomographic observation equation, so as to obtain the regional ionospheric electron density distribution.

[0153] Those skilled in the art can clearly understand that, for the convenience and brevity of description, the specific working process of the above-described device can refer to the corresponding process in the foregoing method embodiments, which will not be repeated here.

[0154] An embodiment of the present application provides a computer device, including a processor and a memory, the memory is used to store a computer program, the computer program is executed by the processor to realize the above-mentioned regional ionospheric electron density calculation method.

[0155] The computer device can be a smart phone, a tablet computer, a desktop computer, a cloud server and the like. The computer device can include but is not limited to a processor and a memory. Those skilled in the art can understand that the computer device can include an input / output device, a network access device and the like.

[0156] The processor can be a central processing unit (CPU), and can also be other general-purpose processors, digital signal processors (DSP), application specific integrated circuits (ASIC), field-programmable gate arrays (FPGA) or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The general-purpose processor can be a microprocessor or any conventional processor.

[0157] The memory can be an internal storage unit of the computer device in some embodiments, for example, a hard disk or a memory of the computer device. The memory can also be an external storage device of the computer device in other embodiments, for example, a plug-in hard disk, a smart media card (SMC), a secure digital (SD) card, a flash card and the like. Further, the memory can include both the internal storage unit and the external storage device of the computer device. The memory is used to store an operating system, an application program, a boot loader, data and other programs, for example, program codes of the computer program, etc. The memory can also be used to temporarily store data that has been output or will be output.

[0158] The above is the preferred embodiment of the present application, it should be noted that for those skilled in the art, without departing from the principles of the present application, can also make several improvements and refinements, these improvements and refinements are also considered the scope of protection of the present application.

Claims

1. A method for calculating regional ionospheric electron density, characterized in that: include: Divide the tomography area into grids and obtain the position of each grid point; Based on the GNSS observation data of each ground station and the operation data of each terminal, the absolute STEC value between each ground station and the GNSS satellite at the tomography time is calculated using non-difference non-combined PPP, and the coordinate information of each station and the coordinate information of the GNSS satellite are obtained; Calculating the intercept length of the signal ray passing through each grid in the to-be-to-tomography area based on all absolute STEC values, coordinate information of all measuring stations, and coordinate information of the GNSS satellite, and establishing a tomographic observation equation; According to the positions of the grid points, the empirical model of electron density of the space-based occultation ionosphere is used to generate a priori initial values ​​of the electron density of each grid point at the tomography moment, specifically: extracting key parameters from the GNSS observation data of the space-based station based on least squares; dividing the space-time grid, and reorganizing the key parameter matrices of the ionosphere according to the key parameters; performing space-time empirical orthogonal decomposition on the key parameter matrices of the ionosphere to obtain empirical orthogonal basis functions and time coefficients; performing Fourier series fitting on the time coefficients corresponding to each order mode to establish an empirical model of key parameters of the ionosphere; inputting the empirical model values ​​of the key parameters of the ionosphere into the Vary-Chapman model for calculation to obtain a first electron density; calculating the electron density of the part above the orbital altitude of the COSMIC series low-orbit satellites by exponential function extrapolation to obtain a second electron density; adding the first electron density and the second electron density to obtain the electron density distribution of the empirical model of electron density of the space-based occultation ionosphere; The grid point electron density priori initial value is iteratively corrected according to the tomographic observation equation to obtain the regional ionospheric electron density distribution.

2. The method for calculating regional ionospheric electron density according to claim 1, wherein: The grid division of the to-be-chromatographed region specifically includes: The area to be layered is divided into grids in three dimensions: latitude, longitude and altitude, and each grid point is assigned a number.

3. The method for calculating regional ionospheric electron density according to claim 1, wherein: The absolute STEC value between each ground station and the GNSS satellite at the tomographic time is calculated using non-differenced non-combined PPP based on the GNSS observation data of each ground station and the operation data of each terminal, specifically including: Based on the GNSS observation data of each ground station and the operation data of each terminal, the observation equation of non-difference non-combined PPP and the MW combined observation value are obtained; Multi-epoch smoothing is performed on the MW combined observations to obtain wide-lane ambiguities, which are then corrected with FCB to construct inter-satellite single-difference wide-lane ambiguities. rounding and fixing the inter-satellite single-difference wide-lane ambiguity to obtain a fixed inter-satellite single-difference wide-lane ambiguity; Obtaining inter-satellite single-difference narrow-lane ambiguity using the fixed inter-satellite single-difference wide-lane ambiguity and the IF floating-point ambiguity, and performing narrow-lane FCB correction; Processing the inter-satellite single-difference narrow-lane ambiguity using a LAMBDA algorithm to obtain a fixed inter-satellite single-difference narrow-lane ambiguity; performing a linear combination of the fixed inter-satellite single-difference wide-lane ambiguity and the fixed inter-satellite single-difference narrow-lane ambiguity to obtain an ionospheric-free combined ambiguity; The corresponding ionospheric TEC observation values ​​are obtained by combining the observation equations of the ionospheric-free combined ambiguity and the non-differenced non-combined PPP.

4. The method for calculating regional ionospheric electron density according to claim 3, wherein: After rounding and fixing the inter-satellite single-difference wide-lane ambiguity to obtain a fixed inter-satellite single-difference wide-lane ambiguity, the method further includes: The correctness of the fixed inter-satellite single-difference wide-lane ambiguity is checked by using a probability judgment function, and the correct fixed inter-satellite single-difference wide-lane ambiguity is stored.

5. The method for calculating regional ionospheric electron density according to claim 3, wherein: After the inter-satellite single-difference narrow-lane ambiguity is processed using the LAMBDA algorithm to obtain a fixed inter-satellite single-difference narrow-lane ambiguity, the method further includes: The correctness of the fixed inter-satellite single-difference narrow-lane ambiguity is checked using the Ratio value, and the correct fixed inter-satellite single-difference narrow-lane ambiguity is stored.

6. The method for calculating regional ionospheric electron density according to claim 1, wherein: The method of calculating the intercept length of the signal ray passing through each grid in the to-be-to-tomography area based on all absolute STEC values, the coordinate information of all measuring stations, and the coordinate information of the GNSS satellite, and establishing a tomographic observation equation specifically includes: Calculate multiple intersection points where the signal ray passes through the altitude plane, longitude plane and latitude plane respectively; Sort the multiple intersection points in ascending order and calculate the distance between two adjacent points as the intercept of the signal ray passing through the corresponding grid from low to high; According to the grid point number, the grid assigns the intercepts of the signal rays passing through the corresponding grid to the intercept matrix one by one; A tomographic observation equation is established according to the intercept matrix.

7. The method for calculating regional ionospheric electron density according to claim 1, wherein: The iterative correction of the a priori initial value of the grid point electron density according to the tomographic observation equation specifically includes: According to the tomographic observation equation, continuous iterative calculation is performed until the standard deviation of the residual between the iteratively calculated STEC and the measured STEC is less than a preset threshold.

8. A device for calculating regional ionospheric electron density, characterized in that: include: The division module is used to divide the tomography area into grids and obtain the position of each grid point; The calculation module is used to calculate the absolute STEC value between each ground station and the GNSS satellite at the tomography time using non-difference non-combined PPP based on the GNSS observation data of each ground station and the operation data of each terminal, and obtain the coordinate information of each station and the coordinate information of the GNSS satellite; a tomography module, configured to calculate the intercept length of a signal ray passing through each grid in the tomography area based on all absolute STEC values, coordinate information of all measuring stations, and coordinate information of the GNSS satellite, and establish a tomography observation equation; An initial module is used to generate a priori initial values ​​of the electron density of each grid point at the tomography moment based on the position of each grid point and using the space-based occultation ionospheric electron density empirical model, specifically: extracting key parameters from the space-based GNSS observation data based on least squares; dividing the space-time grid and reorganizing the ionospheric key parameter matrices according to the key parameters; performing space-time empirical orthogonal decomposition on the ionospheric key parameter matrices to obtain empirical orthogonal basis functions and time coefficients; performing Fourier series fitting on the time coefficients corresponding to each order mode to establish an empirical model of the ionospheric key parameters; inputting the empirical model values ​​of the ionospheric key parameters into the Vary-Chapman model to calculate and obtain a first electron density; calculating the electron density of the part above the orbital altitude of the COSMIC series low-orbit satellite by exponential function extrapolation to obtain a second electron density; adding the first electron density and the second electron density to obtain the electron density distribution of the space-based occultation ionospheric electron density empirical model; The iterative module is used to iteratively correct the prior initial value of the grid point electron density according to the tomographic observation equation to obtain the regional ionospheric electron density distribution.

9. A computer device, characterized in that: The method comprises a processor and a memory, wherein the memory is used to store a computer program, and when the computer program is executed by the processor, the method for calculating the regional ionospheric electron density according to any one of claims 1 to 7 is implemented.

Images