A Swarm / GNSS Ionospheric Tomography Method Based on Hybrid Grid and Weighted Horizontal Constraint

Abstract

The present invention discloses a Swarm / GNSS ionospheric tomography method based on hybrid grids and weighted horizontal constraints. First, an effective observation ray database for the inversion area is established; then an effective observation ray coefficient matrix is established; furthermore, the TEC predicted value is calculated according to the observation ray coefficient matrix and the initial electron density value obtained from the IRI empirical model; the electron density correction in the inversion area is obtained based on the Swarm satellite observation data; finally, the electron density correction in the inversion area is additionally vertically and horizontally constrained, and the processed electron density value is used as the final ionospheric tomography result. This method can effectively overcome the problems of low inversion accuracy and severe dependence on the initial value caused by neglecting the selection of grid resolution in different regions, the distribution of horizontal constraint weights, and the lack of information from a single data source.

Description

Technical Field

[0001] The present invention relates to the field of remote sensing inversion of navigation satellites, and specifically to a Swarm / GNSS ionospheric tomography method based on a hybrid grid and weighted horizontal constraint. Background Technique

[0002] The ionosphere is a high-altitude region in the Earth's atmosphere, about 60 kilometers to 2,000 kilometers above the Earth's surface. It is formed by the ionization of atmospheric molecules under the action of solar radiation and cosmic rays, and mainly contains free electrons, positive ions, and negative ions, with significant spatio-temporal variation characteristics. The ionosphere is not only an important part of the Earth's space environment but also a key medium for radio signal propagation. Its electron density distribution directly affects the propagation path and accuracy of satellite navigation signals, and at the same time has a decisive effect on the performance of technical systems such as short-wave communication and radar detection. Therefore, researching and constructing a high-precision ionospheric delay correction model can significantly reduce satellite navigation positioning errors, optimize radio communication link design, and reduce the risk of signal propagation interruption. This is an important problem that urgently needs to be solved in current ionospheric tomography.

[0003] Currently, there are mainly three categories of ionospheric tomography methods. The first category is the tomography method using ground sounders. By using an ionosonde or a side-looking sounder to emit electromagnetic waves and receive ionospheric reflection signals to obtain parameters such as electron density, although high precision can be obtained, the horizontal resolution is insufficient and cannot meet large-scale measurements. The second category is the tomography method of radio occultation. When an occultation event occurs, the signal transmitted by a GPS satellite to a LEO satellite is received after refraction. By analyzing information such as the phase delay of the received signal, ionospheric parameters can be obtained. Its advantages are global coverage and high vertical resolution, but the spatial distribution of occultation events is restricted by satellite orbits, and a single detection only covers a specific path, with poor spatial continuity. The third category is ionospheric tomography based on GNSS (Global Navigation Satellite System) navigation satellites. This method extracts the TEC (Total Electron Content) of the ionosphere using the observation data of GNSS satellites to obtain the ionospheric electron density over a large area in real time. However, the existing methods do not consider the non-uniform distribution of effective observation rays in the inversion area, and often use a unified grid resolution, resulting in an excessive or insufficient number of effective observation rays passing through individual grids, reducing the tomography accuracy. At the same time, the existing methods do not consider that grids with more rays passing through are more reliable and dependent than grids with fewer or no rays passing through, and use a mean smoothing horizontal constraint with a unified weight, resulting in a reduction in tomography accuracy. In addition, the existing methods often only use a single GNSS data source, ignoring the impact of insufficient observation information caused by the limited number and uneven distribution of ground stations, resulting in ill-posed inversion and low result accuracy. Summary of the Invention

[0004] The purpose of the present invention is to overcome the deficiencies of the prior art and propose a Swarm / GNSS ionospheric tomography method based on hybrid grids and weighted horizontal constraints. This method can effectively overcome problems such as low inversion accuracy and severe initial value dependence caused by neglecting the selection of grid resolution in different regions, the distribution of horizontal constraint weights, and the lack of information from a single data source. By adaptively adjusting the grid resolution layer by layer, the stability and reliability of ionospheric inversion are improved. Using the International Reference Ionosphere (IRI) as the initial value of iteration, the convergence speed and accuracy of inversion are increased. The multiplicative algebraic reconstruction technique (MART) is used to correct the grid electron density values, solving the rank deficiency problem that occurs during the inversion process. Through the electron density detector carried by the Swarm satellite, the electron density of the moving path is accurately corrected, reducing the error of the reconstructed ionosphere. The Chapman function vertical constraint is adopted to weaken the influence of observation noise and further improve the vertical accuracy of the reconstructed ionosphere. Finally, non-equal weight horizontal constraints are applied to each layer of the grid, effectively overcoming the problem of dependence on the initial value and improving the smoothness of the reconstructed ionosphere.

[0005] To achieve the above purpose, the specific technical solutions adopted by the present invention are as follows:

[0006] A Swarm / GNSS ionospheric tomography method based on hybrid grids and weighted horizontal constraints, comprising the following steps:

[0007] The first step: Establish an effective observation ray database. Use a geodetic receiver to collect the GNSS raw observation data of the measurement day. Extract the positions of GNSS navigation satellites in the inversion area on the measurement day from the precise ephemeris file, and extract the positions, times, pseudorange, carrier phase and other raw observation values of ground-based receiving stations in the inversion area from the observation file. And calculate the high-precision ionospheric TEC value of the satellite-receiving station connection line at the corresponding epoch through the pseudorange and carrier phase, as the measured ionospheric TEC value in the subsequent solution. Finally, arrange the GNSS navigation satellite positions, ground-based receiving station positions, and ionospheric TEC values of the satellite-receiving station connection line in time sequence by measuring station and satellite to establish an effective observation ray database, which is stored in the local computer in text form.

[0008] Step 2: Establish the observation ray coefficient matrix and simulate the initial value. First, the inversion area is gridded with low resolution in three dimensions of longitude, latitude, and altitude to obtain a number of three-dimensional solid grids. It is assumed that the electron density within each three-dimensional solid grid is constant in a short period. The longitude and latitude are divided into degrees, and the altitude is divided into kilometers. Second, calculate the intercept values of the effective observation rays (i.e., the effective observation ray database established in the first step) within the inversion area and each solid grid in the tomographic model grid space. Since an excessive or zero number of effective observation rays in a single solid grid will directly affect the subsequent iterative update of the solid grid by the observation rays, the resolution of the solid grid space is adaptively adjusted in layers. The solid grids at the same altitude are divided into the same solid grid layer. For the solid grid layer with a small number of effective observation rays passing through, the original low resolution is still retained, and for the solid grid layer with a larger number of passing rays, a higher resolution is adopted, and the intercept values of the effective observation rays within the solid grids of this layer are recalculated. Then, code according to the division order of the solid grids, and use the intercept lengths of the effective observation rays within the solid grids to construct the observation ray coefficient matrix. Finally, according to the inversion area, inversion time, and solid grid scale, extract the ionospheric electron density on the measurement day from the IRI empirical model as the initial estimate of ionospheric tomography.

[0009] Step 3: The iterative algorithm completes the reconstruction of the electron density in the inversion area. The iterative algorithm uses the multiplicative algebraic reconstruction method based on the principle of maximum entropy. Through the intercept values of the effective observation rays passing through each solid grid (i.e., the observation ray coefficient matrix established in the second step) and the initial values of the electron density of each solid grid (obtained in the second step), the predicted value of the ionospheric TEC on the propagation path of the observation ray can be obtained. Divide it by the measured value of the ionospheric TEC extracted from the GNSS observation data (obtained in the first step), and distribute the ratio to the electron density values of the solid grids passed through by the observation ray in a certain proportion. Updating all the solid grids passed through by the effective observation rays in the inversion area once is called one iteration. And so on, perform multiple rounds of iterative updates on the electron density in the solid grid space of the inversion area, gradually improving the initial estimate of the ionosphere in the inversion area until the iteration terminates when the convergence condition is reached, and obtain the reconstructed electron density value in the inversion area.

[0010] Step 4: Correction of the inverted regional electron density based on Swarm satellite observation data. The original observation data for the measurement day are obtained by the Swarm satellite receivers, mainly including the movement time, movement position of the three satellites A, B, and C of Swarm, and the in-situ electron density values measured by the Langmuir probes carried by them. According to the above information, the movement trajectories of the three Swarm satellites A, B, and C are analyzed, and the running paths within the three-dimensional grid space of the inversion region are screened out. The electron density of the three-dimensional grid through which its path passes is corrected to the high-precision electron density measured by the Langmuir probe at this moment, replacing the inversion value based on GNSS observation data (obtained in Step 3). Through the Swarm / GNSS multi-source data fusion, the estimation deviation problem caused by insufficient observation information of a single data source is improved, effectively enhancing the accuracy and reliability of the reconstructed ionosphere.

[0011] Step 5: Adding vertical and horizontal constraints to the inversion results. In view of the influence of observation noise, the Chapman function vertical constraint is adopted for the three-dimensional grid space in the inversion region. The three-dimensional grids at consecutive heights at the same latitude and longitude position are divided into grid bundles, and the least squares method is used to fit the relationship between the electron density and height of the three-dimensional grids in each grid bundle with the Chapman function. During the fitting process, the parameters of the Chapman function are adjusted to make the generated Chapman function formula fit better with the electron density values measured by the Swarm satellites (obtained in Step 4), improving the accuracy of the fitting. Each grid bundle then reconstructs the electron density in its vertical direction based on the fitted Chapman function to achieve precise vertical constraint. In view of the problem of the initial value dependence of the inversion results, a non-equal weight horizontal constraint considering the distribution of effective observation rays is adopted. The electron density value of the central three-dimensional grid is constrained and corrected using the electron density values of the three-dimensional grids within the neighborhood range, thereby controlling the difference in electron density between the central three-dimensional grid and the adjacent three-dimensional grids. At the same time, since the three-dimensional grids through which more observation rays pass are more credible and reliable than those through which fewer or even no rays pass, the constraint weight value for the central three-dimensional grid is divided according to the number of effective observation rays passing through the neighborhood three-dimensional grids. The horizontal constraint is performed on each three-dimensional grid in the inversion region according to the above method, and the processed electron density value is used as the final ionospheric tomography result.

[0012] After being processed by the above ionospheric tomography algorithm, the problems such as low inversion accuracy and serious initial value dependence caused by the existing methods ignoring the selection of three-dimensional grid resolution in different regions, the distribution of horizontal constraint weight values, and the lack of information from a single data source can be effectively solved, effectively improving the real-time performance, accuracy, and stability of ionospheric tomography, providing a guarantee for the reconstruction of the ionosphere.

[0013] Study the influence of the selection of three-dimensional grid resolution in different regions, the weight distribution of horizontal constraints, and the deficiency of single data source observation information on the accuracy of ionospheric tomography, and effectively solve problems such as low accuracy of ionospheric tomography results and strong dependence of inversion on initial values. This can not only meet the requirements of time and space resolution of ionospheric tomography, but also meet the requirements of accuracy, real-time performance, and stability of tomography, providing strong support for ionospheric tomography.

[0014] The present invention has the following characteristics and beneficial effects:

[0015] By using the method of the present invention, the ionosphere can be tomographed using existing Swarm satellites / GNSS navigation system satellites. Compared with existing tomography methods based on ground sounders and radio occultation, it can solve problems such as insufficient horizontal resolution and inability to meet large-scale measurement requirements, and can also solve problems such as poor spatial continuity due to satellite orbit limitations. Compared with existing tomography methods based on GNSS navigation satellites, it can solve problems such as the influence of the selection of three-dimensional grid resolution in different regions on the accuracy of iterative inversion, and can also solve problems such as low inversion accuracy and strong dependence on initial values caused by ignoring the credibility of three-dimensional grids and lack of single data source information during horizontal constraint. By adaptively adjusting the three-dimensional grid resolution layer by layer, the effective observation rays are evenly distributed in the three-dimensional grid, ensuring that each three-dimensional grid is fully iteratively corrected, and improving the stability and reliability of ionospheric inversion. By using the IRI model as the initial value of iteration, the convergence speed and accuracy of inversion are improved. The multiplicative algebraic reconstruction method is used to correct the electron density values of the three-dimensional grid, improving the rank deficiency problem that appears during the inversion process. Through the Langmuir electron density probe carried by the Swarm satellite, the electron density of the three-dimensional grid passed by it in the inversion area is accurately corrected, improving the accuracy of the reconstructed ionosphere. The Chapman function is used for vertical constraint, weakening the influence of observation noise and improving the vertical accuracy of the reconstructed ionosphere. Non-equal weight horizontal constraints are used for each layer of the three-dimensional grid, and the weight value is determined by the number of effective observation rays passing through the three-dimensional grid, so that more credible three-dimensional grids are given higher weight values, overcoming the negative constraints generated by the three-dimensional grids that are not passed through by the observation rays and updated, and the problem of dependence on the initial value of inversion, further ensuring the accuracy, real-time performance, and stability of ionospheric tomography. Description of the Drawings

[0016] Figure 1 Flowchart for establishing an effective observation ray database based on GNSS observation data

[0017] Figure 2 Flowchart for establishing an observation ray coefficient matrix based on a hierarchical adaptive three-dimensional grid

[0018] Figure 3 Flowchart for iterative inversion of electron density based on the multiplicative algebraic reconstruction method

[0019] Figure 4 Flow chart of electron density correction based on Swarm satellite observation data

[0020] Figure 5 Flow chart of optimizing the inversion result based on the vertical constraint of Chapman function and non - equal - weight horizontal constraint Specific implementation manner

[0021] The present invention will be described in detail below with reference to specific embodiments. The following embodiments will help those skilled in the art to further understand the present invention, but do not limit the present invention in any form. It should be noted that, without conflict, the embodiments in the present invention and the features in the embodiments can be combined with each other.

[0022] A Swarm / GNSS ionospheric tomography method based on a hybrid grid and weighted horizontal constraint, comprising the following steps:

[0023] Step 1: Extract the original observation data of GNSS navigation satellites and establish an effective observation ray database.

[0024] Specifically, as Figure 1 shown, the extracted original observation data mainly includes the longitude, latitude, and altitude positions of GNSS navigation satellites within the inversion area on the measurement day, the longitude, latitude, and altitude positions of ground - based receiving stations, time, pseudorange, and carrier phase. The positions of GNSS navigation satellites within the inversion area on the measurement day can be directly extracted from precise ephemeris files, and the positions, time, pseudorange, and carrier phase of ground - based receiving stations within the inversion area can be directly extracted from observation files. The connection line between a ground - based receiving station and a GNSS navigation satellite is called an effective observation ray. The high - precision ionospheric TEC of the effective observation ray corresponding to the epoch can be calculated through the carrier phase and pseudorange. This process is common knowledge in the positioning field, and the calculation process will not be listed in detail here. Calculate the actual values of the high - precision ionospheric TEC of each effective observation ray separately by time, by measurement station, and by satellite, arrange the data of each effective observation ray in chronological order, and number the effective observation rays in the arranged order. The data of the effective observation ray includes the position and name of the receiving station corresponding to this effective observation ray, the position and satellite number of the corresponding satellite, and the actual value of the ionospheric TEC. Establish a database and store it in a local computer in text form. The high - precision ionospheric TEC value extracted here is called the measured value of ionospheric TEC. For the convenience of subsequent description, it is represented by in the subsequent steps. The longitude, latitude, and altitude of the ground - based receiving station are respectively represented as 、 、 respectively, and the longitude, latitude, and altitude of the GNSS navigation satellite are respectively represented as 、 、 。

[0025] Step 2: Grid the inversion area and establish the observation ray coefficient matrix. The traditional method uses a uniform three-dimensional grid resolution, without considering that the number of effective observation rays passing through individual three-dimensional grids exceeds the optimal inversion amount or is zero, resulting in reduced tomography accuracy. Therefore, a hierarchical adaptive three-dimensional grid resolution considering the distribution of effective observation rays is adopted here, as shown in Figure 2 The specific process is as follows:

[0026] First, divide the inversion area into three-dimensional grids with low resolution. By using low resolution, the number of three-dimensional grids through which the observation rays pass is reduced to zero, and the problem of initial value dependence caused by the failure of the three-dimensional grids to be updated by the rays is improved. The inversion area is divided into d, f, and r three-dimensional grids with low resolution in the longitude, latitude, and altitude directions respectively. The three-dimensional grids with low resolution specifically refer to: three-dimensional grids with a longitude span of 2 degrees, a latitude span of 2 degrees, and an altitude span of 50 kilometers; that is, the entire inversion area is divided into d×f×r three-dimensional grids, where the longitude and latitude units are divided into degrees and the altitude unit is divided into kilometers. Then, number the three-dimensional grids in ascending order from small to large in the order of longitude first, then latitude, and finally altitude. The position of each three-dimensional grid is represented by the longitude , latitude , altitude of its center.

[0027] Secondly, solve the observation ray coefficient matrix, that is, find the intercept value of each ray in the three-dimensional grid and store it in the matrix corresponding to the position of the three-dimensional grid. During the solution process, in order to improve the inversion accuracy of the tomography model, the Earth is regarded as a regular ellipsoid under ideal conditions, and the solution is carried out in the Earth ellipsoid coordinate system, where the semi-major axis a of the ellipsoid is 6378137 meters and the semi-minor axis b is 6356752.314 meters. The specific process of this step of the algorithm is as follows:

[0028] 1) Calculate the linear equation of the effective observation ray. In the Earth ellipsoid coordinate system, according to the spatial coordinates ( ) (obtained in step 1) of the ground-based receiving station and the spatial coordinates ( ) (obtained in step 1) of the GNSS navigation satellite, the equation of the satellite effective observation ray can be obtained by the principle of solving the spatial straight line equation. The formula is as follows:

[0029] (1)

[0030] In the formula, X, Y, and Z are the coordinate values of any point on the effective observation ray, and c is the proportional parameter.

[0031] 2) Calculate the intersection coordinates of the effective observation ray with each three-dimensional grid. Each three-dimensional grid is composed of six faces, including two longitude planes, two latitude planes, and two altitude planes, and the equations corresponding to different surfaces are different.

[0032] The longitude plane and the latitude plane are perpendicular and both pass through the Z-axis of the above-mentioned geodetic coordinate system. Determine the angle between the required longitude plane and the prime meridian plane in the geodetic coordinate system as A, then the equation expression of the longitude plane is:

[0033] (2)

[0034] In the formula, X and Y are the coordinate values of any point on the longitude plane.

[0035] The latitude plane is the surface obtained by rotating the line connecting a certain point in the space region and the geodetic centroid around the Z-axis of the geodetic coordinate system. Determine the latitude of the required latitude plane as W, then the equation expression of the latitude plane is:

[0036] (3)

[0037] In the formula, X, Y, and Z are the coordinate values of any point on the latitude plane, a = 6378137 meters, which is the semi-major axis of the ellipsoid in the geodetic coordinate system, and b = 6356752.314 meters, which is the semi-minor axis of the ellipsoid in the geodetic coordinate system.

[0038] The altitude planes are all parallel to the Earth's surface. Determine the distance between the required altitude plane and the Earth's surface as H, then the equation expression of the altitude plane is:

[0039] (4)

[0040] In the formula, X, Y, and Z are the coordinate values of any point on the altitude plane, and the meanings of a and b are the same as those in formula (3).

[0041] Simultaneously solve the equations of the longitude plane, latitude plane, and altitude plane (2), (3), and (4) with the space straight line equation (1) formed by the effective observation ray to obtain the intersection coordinates of each effective observation ray in the inversion area with the longitude plane, latitude plane, and altitude plane of each three-dimensional grid .

[0042] 3) Calculate the intercept value of the effective observation ray. After obtaining the intersection coordinates of the effective observation ray with each three-dimensional grid, screen out the three-dimensional grids with two intersections. The two intersections are represented as P1 and P2 . According to the distance formula between two points, calculate the intercept ΔL of the ray in this three-dimensional grid. The specific formula is as follows:

[0043] (5)

[0044] 4) Construct the observation ray coefficient matrix. The intercept values of each effective observation ray within each three-dimensional grid are stored in the coefficient matrix in sequence according to the coding order of the three-dimensional grid, where the intercept of the three-dimensional grid not penetrated by the ray is 0. Each row of the observation ray coefficient matrix represents an observation ray, each column corresponds to a three-dimensional grid, and the arrangement order of the elements in each row corresponds to the coding order of the three-dimensional grid.

[0045] After solving the observation ray coefficient matrix for the low-resolution three-dimensional grid by the above method, within the inversion region, the three-dimensional grids at the same altitude are divided into three-dimensional grid layers. The proportion of the number of times the effective observation rays pass through the three-dimensional grids in each layer to the total number of three-dimensional grids in that layer is obtained by calculating the number of intercepts of the effective observation rays in the three-dimensional grids of each layer, which is called the ray penetration rate. The ray penetration rates of each three-dimensional grid layer are sorted from large to small. The three-dimensional grid layers in the top 50% are those where the effective observation rays pass through relatively densely. To improve the situation where the rays passing through the three-dimensional grid are too dense and exceed the optimal inversion amount, the three-dimensional grid layer is further subdivided, that is, the longitude, latitude, and altitude of each three-dimensional grid in the three-dimensional grid layer are scaled to half of the original, and changed to a higher-resolution three-dimensional grid. For the three-dimensional grids within these three-dimensional grid layers refined to high resolution, all the three-dimensional grids at the same altitude within the inversion region are regarded as the same three-dimensional grid layer. The number of times the effective observation rays pass through all the three-dimensional grids in each three-dimensional grid layer is counted, and the three-dimensional grid layers in the top 50% with a higher number of times the effective observation rays pass through are further divided, that is, the longitude, latitude, and altitude of each three-dimensional grid in this three-dimensional grid layer are all divided into half of the original, and adjusted to a three-dimensional grid with a longitude span of 1 degree, a latitude span of 1 degree, and a height span of 25 km, which is a high-resolution three-dimensional grid.

[0046] Adopt re-numbering based on the original three-dimensional grid numbers to ensure the uniqueness of each three-dimensional grid identifier. For example, the three-dimensional grid with the original number 3 is divided into eight smaller three-dimensional grids after refinement, and its numbers are represented as (3,1), (3,2) to (3,8), and so on, to obtain the re-numbering of all high-resolution three-dimensional grids. Recalculate the observation ray coefficient matrix of the high-resolution three-dimensional grid, that is, repeat steps 2), 3), and 4) of the above process. Finally, merge the coefficient matrix of the remaining low-resolution three-dimensional grids and the coefficient matrix of the high-resolution three-dimensional grids to obtain the final observation ray coefficient matrix.

[0047] Step 3: Use the IRI model to simulate the initial electron density of the three-dimensional grid. The IRI model provides important parameters such as the electron density of the non-polar ionosphere within the altitude range of 50 - 2000 km under quiet geomagnetic field conditions. By inputting the inversion region, inversion time, and three-dimensional grid scale into the IRI online model, the ionospheric electron density of each three-dimensional grid within the inversion region on the measurement day can be generated, which can be expressed as: , as the initial estimate for tomographic solution.

[0048] Step 4: Invert the electron density in the grid space through MART iteration.

[0049] Specifically, as Figure 3 shown, the MART based on the principle of maximum entropy is selected as the iterative algorithm for this case, which effectively improves the rank deficiency problem in the parameter solution process and ensures that the electron density will not be negative during the iteration process. MART is based on the initial value of the electron density in the three-dimensional grid space within the inversion region and gradually improves the estimated value of each three-dimensional grid in an iterative manner. Correcting the electron density of all three-dimensional grids passed through by an effective observation ray is called one iteration. During the correction, the intercept value of the observation ray within the three-dimensional grid (obtained in Step 2) is multiplied by the electron density value calculated in the previous iteration to obtain the predicted value of the ionospheric TEC along the propagation path of the observation ray. Using the ratio of this predicted value to the measured TEC value, the electron density values of all three-dimensional grids passed through by this observation ray are updated accordingly. Updating all the three-dimensional grids passed through by all the rays in the inversion region once is called one round of iteration. By analogy, multiple rounds of iteration are performed on the electron density in the three-dimensional grid space of the inversion region to gradually improve the initial estimate of the ionosphere to be reconstructed until the iteration terminates when the convergence condition is reached, and the reconstructed electron density value is obtained. The correction formula of the specific MART algorithm in the k-th iteration is as follows:

[0050] (6)

[0051] In the formula, is the electron density value of the j-th three-dimensional grid after the k-th iteration, where is the initial estimate of the electron density of the j-th three-dimensional grid, that is, = (obtained in Step 3), is the measured value of the ionospheric TEC of the i-th effective observation ray (obtained in Step 1), is the i-th row vector of the observation ray coefficient matrix (obtained in Step 2), corresponding to the intercept of the i-th observation ray within the three-dimensional grid, is the element in the i-th row and j-th column of the observation ray coefficient matrix, corresponding to the intercept of the i-th observation ray within the j-th three-dimensional grid. γ is the relaxation factor of the iteration, 0 < γ < 1, which is used to control the convergence speed. In this case, γ = 0.2. When the root mean square error between the predicted TEC value and the measured TEC value of each observation ray is less than the given tolerance σ, the iteration is terminated. In this case, σ is set to 0.01 TECU. TECU is a commonly used unit in the field, 1 TECU = 1×10 16 electrons per square meter. Through the above process, the electron density value obtained by iterative inversion using GNSS observation data is expressed as: 。

[0052] Step 5: Correction of the electron density in the inversion area based on Swarm satellites.

[0053] It should be noted that the Swarm constellation consists of three identical satellites, namely A, B, and C, all equipped with Langmuir probes, which can measure the high-precision electron density values at the positions of the satellites at corresponding times. By using the high-precision measurement values of the Swarm satellites to correct the electron density values in the three-dimensional grid space after inversion, the problem of reduced accuracy caused by insufficient observation data in the traditional method is improved. As Figure 4 shown, the specific process is divided into three steps:

[0054] 1) Extract the original observation data of the Swarm satellites and establish a database of the Swarm satellite motion paths. The original observation data extracted mainly includes the motion time of the three satellites A, B, and C of the Swarm on the measurement day, the corresponding longitude, latitude, and altitude positions, and the in-situ electron density values measured by the Langmuir probes they carry. All the above information can be directly obtained from the observation files. Arrange the electron density values measured by the Swarm satellites and the corresponding positions according to time and satellite name, and establish a motion path database, which is stored in the local computer in text form. For the convenience of subsequent description, the electron density value measured by the Swarm satellite is expressed as: , and the longitude, latitude, and altitude of the Swarm satellite are respectively expressed as 、 、 。

[0055] 2) Screen out the motion paths of the Swarm satellites within the inversion area. Compare the positions of the Swarm satellites in the motion path database ( ) with the longitude, latitude, and altitude ranges of the inversion area, and delete the data located outside the inversion area to obtain the screened Swarm satellite motion path database.

[0056] 3) Obtain the three-dimensional grid where the Swarm satellite is located according to the position of the Swarm satellite in the motion path database ( ), and replace the electron density inversion value (obtained in step 4) originally obtained from GNSS observation data in this three-dimensional grid with the high-precision electron density measured value of the Swarm satellite. If the same three-dimensional grid is passed through by multiple Swarm satellites, take the mean of all the measured values as the electron density correction value of this three-dimensional grid. Through the above process, the electron density value in the three-dimensional grid space corrected by the high-precision measured value of the Swarm satellite is obtained, which is expressed as: 。

[0057] Step 6: Add the Chapman function vertical constraint to the three-dimensional grid space within the inversion region. Due to the limitation of observational information, the vertical accuracy of the reconstructed ionosphere is generally low. Therefore, in this case, the Chapman function, which can better describe the characteristics of the ionospheric vertical profile, is used to add the vertical constraint to the electron density of the three-dimensional grid space to improve the accuracy of the reconstructed ionosphere. As Figure 5 shown, the specific process is divided into three steps:

[0058] 1) Unify the three-dimensional grid resolution. The three-dimensional grid space of ionospheric tomography in this case is a mixed combination of low-resolution and high-resolution three-dimensional grids. Since the Chapman function vertical constraint requires dividing the three-dimensional grid beam with a unified resolution three-dimensional grid, the low-resolution three-dimensional grids need to be subdivided before the constraint. Subdivide the longitude, latitude, and height of all low-resolution three-dimensional grid layers in the three-dimensional grid space into half of the original, and re-number the subdivided three-dimensional grids. The specific process is as described in Step 2. The electron density of the subdivided three-dimensional grid remains the same as that before subdivision, that is, one low-resolution three-dimensional grid is subdivided into eight high-resolution three-dimensional grids with the same electron density as the original. The specific formula is as follows:

[0059] (7)

[0060] In the formula, is the electron density of the m-th sub-three-dimensional grid obtained by subdividing the low-resolution three-dimensional grid numbered j. By subdividing the low-resolution three-dimensional grid, the original three-dimensional grid space with mixed resolution is adjusted to a unified high-resolution three-dimensional grid space.

[0061] 2) Fit the Chapman function and correct the electron density. According to the spatial position of the three-dimensional grid, the three-dimensional grids with continuous heights at the same longitude and latitude positions are divided into a three-dimensional grid beam. Then, the electron density values and corresponding height values of the three-dimensional grids in this three-dimensional grid beam need to satisfy the Chapman function distribution. Take the inverted electron density values and the three-dimensional grid center height values corresponding to all three-dimensional grids in each three-dimensional grid beam as inputs, and fit the Chapman function corresponding to each three-dimensional grid beam through the non-linear weighted least squares algorithm. The specific formula of the Chapman function is as follows:

[0062] (8)

[0063] In the formula, h is the center height value of the three-dimensional grid, that is, h = (obtained in Step 2), is the electron density of the three-dimensional grid at the center height h, that is, (obtained in Step 5), h mis the maximum value of the electron density of each three-dimensional grid in the three-dimensional grid bundle. The F2 layer is the region with the highest electron density in the ionosphere, N m F2 is the peak electron density value of the F2 layer, h m F2 is the peak height value of the F2 layer, H1 represents the range height of the bottom ionosphere, and H2 represents the range height of the top ionosphere. The above four parameters can be solved by the non-linear weighted least squares method. Since some three-dimensional grids in the three-dimensional grid space are corrected by the measured values of the Swarm satellites in step 5, when performing the least squares solution of the Chapman parameters for the three-dimensional grid bundle containing such three-dimensional grids, a weight value of 1 is assigned to the electron density value obtained from its GNSS observation data, and a weight value of 5 is assigned to the measured electron density value corrected by the Swarm satellites. By assigning a high weight value to the high-precision measured values of the Swarm satellites, the Chapman function fitted by the parameters is made more consistent with the Swarm measured values, improving the accuracy and reliability of ionospheric tomography. There is a ready-made function module for the non-linear weighted least squares algorithm in the Matlab platform. Here, we only give a simple algorithm description as follows:

[0064] The core idea of the non-linear weighted least squares algorithm is to find the optimal model parameters by minimizing the weighted squared error. The expression of its objective function is:

[0065] (9)

[0066] where n is the number of data points, is the weight of each data point, is the i-th observed value, is the predicted value calculated according to the model parameters θ and the input data The goal is to find a set of optimal parameters θ * , such that the objective function is minimized.

[0067] To achieve this goal, first, calculate the gradient of the objective function with respect to the parameters. The specific formula is as follows:

[0068] (10)

[0069] where , is the error of the i-th observed data point, is the partial derivative of the model predicted value with respect to the parameter , usually called the Jacobian matrix. The meanings of other parameters are the same as those in formula (9).

[0070] Then, use the Levenberg-Marquardt algorithm to update the parameters of the model. The specific formula is as follows:

[0071] (11)

[0072] In the formula, is the current parameter estimate, is the updated parameter estimate, J is the Jacobian matrix, T is the transpose operation, which belongs to the common sense operations of linear algebra, λ is the damping factor, I is the identity matrix, which belongs to the common sense matrix of linear algebra, and the meaning of e is the same as that in formula (10). After each update, it is judged whether to converge by checking the changes of the objective function and the parameters. After meeting the convergence conditions, the algorithm terminates and outputs the optimal parameter θ * .

[0073] Through the above non-linear weighted least squares algorithm, the four Chapman function parameters fitted for each stereo grid bundle can be accurately obtained, and the parameter values obtained by fitting are substituted back into the Chapman function formula to reconstruct the vertical distribution function of each stereo grid bundle. Thus, each stereo grid bundle inputs the central height value of each of its stereo grids , and the corrected electron density value of the Chapman function can be obtained

[0074] 3) Restore the unified stereo grid resolution and calculate the electron density value in the stereo grid space with mixed resolution. For the original high-resolution stereo grid, its electron density value remains unchanged. For the eight high-resolution stereo grids obtained by subdividing the same low-resolution stereo grid, the average value of their electron densities is taken as the final electron density value of the low-resolution stereo grid. The specific formula is as follows:

[0075] (12)

[0076] In the formula, is the electron density of the m-th sub-stereo grid obtained by subdividing the low-resolution stereo grid numbered j. Thus, the electron density value of the stereo grid with mixed resolution vertically constrained by the Chapman function is obtained .

[0077] Step 7: Add non-equal-weight horizontal constraints to the stereo grid space in the inversion area. The traditional method is to use the mean smoothing constraint with uniform weight for the stereo grids of each height layer. However, the stereo grids with fewer or even no rays passing through may have significantly lower or higher values because they are updated less by the observed rays, which is very likely to lead to negative constraint results contrary to the constraint intention. Therefore, non-equal-weight horizontal constraints considering the distribution of effective observed rays are adopted here

[0078] Such as Figure 5As shown in the figure, the specific process is as follows: The three-dimensional grid space is stratified according to different heights, and each three-dimensional grid is corrected layer by layer one by one. The three-dimensional grid to be corrected and its adjacent three-dimensional grids are collectively called the neighborhood three-dimensional grids, where the adjacent three-dimensional grids are all the three-dimensional grids in the same height layer that have a common edge or a common vertex with the three-dimensional grid to be corrected. During correction, a weight factor is set for the three-dimensional grid according to the number of effective observation rays passing through it, and the weighted average of the electron densities of the neighborhood three-dimensional grids is used as the corrected electron density value of the three-dimensional grid to be corrected, so as to perform smoothing processing. The specific formula is as follows:

[0079] (13)

[0080] In the formula, q is the total number of neighborhood three-dimensional grids, is the number of observation rays passing through the i-th three-dimensional grid in the neighborhood three-dimensional grids, is the weight factor of the i-th three-dimensional grid in the neighborhood three-dimensional grids, is the electron density of the i-th three-dimensional grid in its neighborhood three-dimensional grids when correcting the j-th three-dimensional grid. Thus, the electron density after adding the additional non-equal weight horizontal constraint is obtained as the final electron density inversion value.

[0081] Step 8: Save and output the inversion result.

[0082] The above shows and describes the basic principles, main features and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited by the above embodiments. The above embodiments and the descriptions in the specification are only preferred examples of the present invention and are not used to limit the present invention. Without departing from the spirit and scope of the present invention, the present invention will have various changes and improvements, and these changes and improvements all fall within the scope of the present invention claimed. The scope of protection claimed by the present invention is defined by the appended claims and their equivalents.

Claims

1. A Swarm / GNSS ionospheric tomography method based on hybrid grid and weighted horizontal constraint, characterized in that The method includes the following steps: Step 1: Establish an effective observation ray database for the inversion area. The effective observation rays are the lines connecting the ground-based receiving stations and GNSS navigation satellites. The data of the effective observation rays include the position and name of the receiving station corresponding to the effective observation ray, the position and satellite number of the corresponding satellite, and the actual value of the ionospheric TEC. Each row in the effective observation ray database represents the data of an effective observation ray, and each effective observation ray is encoded according to the order of data arrangement; Step 2: Establish an observation ray coefficient matrix; Step 2-1: Divide the inversion area into low-resolution three-dimensional grids; Step 2-2: Calculate the intercept values of the effective observation rays in the inversion area with all three-dimensional grids. All three-dimensional grids at the same height in the inversion area are regarded as the same three-dimensional grid layer. According to the number of times the effective observation rays pass through each three-dimensional grid layer, all the three-dimensional grids in this three-dimensional grid layer are adaptively adjusted to high-resolution three-dimensional grids; Step 2-3: Recalculate the intercept values of the effective observation rays with all high-resolution three-dimensional grids; Step 2-4: Arrange the three-dimensional grids according to the encoding and intercept length of the effective observation rays passing through them, and use the IRI empirical model as the initial value of the electron density of the three-dimensional grids. Finally, output the observation ray coefficient matrix. Each row of the observation ray coefficient matrix represents an observation ray, each column corresponds to a three-dimensional grid, and the arrangement order of the elements in each row corresponds to the encoding order of the three-dimensional grids; Step 3: Calculate the predicted TEC value based on the observation ray coefficient matrix and the initial electron density value obtained from the IRI empirical model. Use the multiplicative algebraic reconstruction method to correct the electron density of each three-dimensional grid by the ratio of the predicted TEC value to the measured TEC value, and perform multiple rounds of iterative update on the electron density of the three-dimensional grids until the convergence condition is reached to obtain the reconstructed electron density value of the inversion area; Step 4: Obtain the electron density correction of the inversion area based on the Swarm satellite observation data; Step 5: Add vertical and horizontal constraints to the electron density correction of the inversion area, and use the processed electron density value as the final ionospheric tomography result.

2. The Swarm / GNSS ionospheric tomography method based on hybrid grid and weighted horizontal constraint according to claim 1, wherein In step 1, a geodetic receiver is used to collect the GNSS raw observation data on the measurement day, and the positions of GNSS navigation satellites in the inversion area on the measurement day are extracted from the precise ephemeris file. The GNSS raw observation data includes the geodetic coordinates and height of GNSS navigation satellites in the inversion area on the measurement day, the geodetic coordinates and height of the ground-based receiving station, time, pseudorange, and carrier phase.

3. A Swarm / GNSS ionospheric tomography method based on a hybrid grid and weighted horizontal constraints according to claim 2, characterized in that, The actual TEC value of the effective observation ray is calculated through the carrier phase and pseudorange. The TEC values of the effective observation rays and the positions of the corresponding receiving stations and satellites are arranged according to time, receiving station name, and satellite number, and the effective observation rays are numbered in order to obtain the effective observation ray database.

4. A Swarm / GNSS ionospheric tomography method based on hybrid grids and weighted horizontal constraints according to claim 2, characterized in that, In the said step 2-1, the method for performing three-dimensional grid division on the inversion area is as follows: the inversion area is divided into d, f, and r three-dimensional grids in the longitude, latitude, and height directions respectively, that is, the entire inversion area is divided into d×f×r three-dimensional grids, and then the three-dimensional grids are numbered from small to large in the order of first along the longitude, then along the latitude, and finally along the height. The position of each three-dimensional grid is represented by the longitude , latitude , and height .

5. A Swarm / GNSS ionospheric tomography method based on a hybrid grid and weighted horizontal constraints according to claim 1, characterized in that, The low-resolution three-dimensional grid refers to a three-dimensional grid with a longitude span of 2 degrees, a latitude span of 2 degrees, and a height span of 50 kilometers.

6. A Swarm / GNSS ionospheric tomography method based on a hybrid grid and weighted horizontal constraints according to claim 5, characterized in that, In step 2-2, the method of adaptively adjusting to high-resolution three-dimensional grids is as follows: All the three-dimensional grids at the same height within the inversion area are regarded as the same three-dimensional grid layer. The number of times that all the three-dimensional grids in each three-dimensional grid layer are penetrated by the effective observation rays is counted. The first 50% of the three-dimensional grid layers with a higher number of times of being penetrated by the effective observation rays are re-divided, that is, the longitude, latitude, and height of each three-dimensional grid in this three-dimensional grid layer are all divided into half of the original, and adjusted to a three-dimensional grid with a longitude span of 1 degree, a latitude span of 1 degree, and a height span of 25 km, namely a high-resolution three-dimensional grid.

7. A Swarm / GNSS ionospheric tomography method based on a hybrid grid and weighted horizontal constraints according to claim 4, characterized in that, The calculation method of the intercept value is as follows: First, after obtaining the intersection coordinates of the effective observation rays and each three-dimensional grid, the three-dimensional grids with two intersections are screened out. The two intersections are denoted as P1 and P2 , and the intercept ΔL of the ray within the three-dimensional grid is calculated according to the distance formula between two points. The specific formula is as follows: 。 8. A Swarm / GNSS ionospheric tomography method based on a hybrid grid and weighted horizontal constraints according to claim 1, characterized in that, In step 3, the electron density in the three-dimensional grid space is inversely calculated through MART iteration. The correction formula of the MART algorithm in the k-th iteration is as follows: ； In the formula, is the electron density value after the k-th iteration of the j-th three-dimensional grid, where is the initial electron density value of the j-th three-dimensional grid, is the measured value of the ionospheric TEC of the i-th effective observation ray, is the i-th row vector of the observation ray coefficient matrix, corresponding to the intercept of the i-th observation ray within the three-dimensional grid, is the element in the i-th row and j-th column of the observation ray coefficient matrix, corresponding to the intercept of the i-th observation ray within the j-th three-dimensional grid. γ is the relaxation factor of the iteration, 0 < γ < 1, which is used to control the convergence speed.

9. A Swarm / GNSS ionospheric tomography method based on a hybrid grid and weighted horizontal constraints according to claim 1, characterized in that In step 4, the following sub-steps are included: Step 4-1: Collect the original observation data of the Swarm satellite on the measurement day, so as to obtain the motion trajectory of the Swarm satellite, and screen out the motion path of the Swarm satellite in the three-dimensional grid space; Step 4-2: Locate the three-dimensional grids penetrated by the motion path of the Swarm satellite; Step 4-3: Correct the electron density of the three-dimensional grids penetrated by the Swarm satellite to the measured high-precision electron density, and output the corrected electron density value.

10. A Swarm / GNSS ionospheric tomography method based on hybrid grid and weighted horizontal constraint according to claim 9, characterized in that, In step 4-3, the electron density correction method is as follows: According to the position of the Swarm satellite in the motion path database, obtain the three-dimensional grid where it is located, and replace the electron density inversion value obtained from the original corresponding GNSS observation data of the three-dimensional grid with the measured high-precision electron density value measured by the Swarm satellite. If the same three-dimensional grid is penetrated by multiple Swarm satellites, take the average value of all its measured values as the electron density correction value of this three-dimensional grid.

11. A Swarm / GNSS ionospheric tomography method based on a hybrid grid and weighted horizontal constraints according to claim 1, characterized in that In step 5, a Chapman function vertical constraint is added to the three-dimensional grid space within the inversion area: First, the original mixed-resolution three-dimensional grid space is adjusted to a unified high-resolution three-dimensional grid space by subdividing the low-resolution three-dimensional grids. The expression is as follows: ； wherein, is the electron density of the m-th sub-stereoscopic grid into which the low-resolution stereoscopic grid numbered j is divided; Then, the three-dimensional grids with continuous heights at the same longitude and latitude positions are divided into the same grid bundle. The electron density inversion values and the three-dimensional grid center height values corresponding to all the three-dimensional grids in each grid bundle are used as inputs, and the Chapman function corresponding to each grid bundle is fitted through the non-linear weighted least squares algorithm to obtain the electron density value after the Chapman function correction; Finally, for the original high-resolution three-dimensional grids, their electron density values remain unchanged. For the eight high-resolution three-dimensional grids subdivided from the same low-resolution three-dimensional grid, take the average value of their electron densities as the final electron density value of the low-resolution three-dimensional grid. The specific formula is as follows: ； In the formula, is the electron density of the m-th sub-stereo grid obtained by subdividing the low-resolution stereo grid numbered j. Thus, the electron density value of the mixed-resolution stereo grid vertically constrained by the Chapman function is obtained. .

12. A Swarm / GNSS ionospheric tomography method based on a hybrid grid and weighted horizontal constraints according to claim 11, characterized in that The specific formula of the Chapman function is as follows: ； In the formula, is the central height value of the three-dimensional grid, is the electron density of the three-dimensional grid at the central height at , that is , is the maximum value of the electron densities of the three-dimensional grids in the three-dimensional grid beam, layer is the region with the highest electron density in the ionosphere, is the peak electron density value of the layer, is the peak height value of the layer, represents the range height of the bottom ionosphere, represents the range height of the top ionosphere.

13. A Swarm / GNSS ionospheric tomography method based on a hybrid grid and weighted horizontal constraints according to claim 1, characterized in that, In step 5, a non-equal weight horizontal constraint is added to the three-dimensional grid space within the inversion area: First, the three-dimensional grid space is stratified according to different heights, and each three-dimensional grid is corrected layer by layer. The three-dimensional grid to be corrected and its adjacent three-dimensional grids are collectively called the neighborhood three-dimensional grids. When correcting, a weight factor is set for the three-dimensional grid according to the number of effective observation rays passing through it, and the weighted average of the electron densities of the neighborhood three-dimensional grids is used as the corrected electron density value of the three-dimensional grid to be corrected, and smoothing processing is carried out accordingly. The specific formula is as follows: ; In the formula, is the total number of neighborhood three-dimensional grids, is the number of the th three-dimensional grid in the neighborhood three-dimensional grids that passes through the observation ray, is the weight factor of the th three-dimensional grid in the neighborhood three-dimensional grids, is the electron density of the th three-dimensional grid in its neighborhood three-dimensional grids when correcting the th three-dimensional grid, and the electron density after obtaining the additional non-equal weight level constraint is used as the final electron density inversion value.

Images