Swarm / GNSS ionosphere chromatography method based on mixed grid and weighted horizontal constraint

Abstract

The invention discloses a mixed grid and weighted horizontal constraint-based Swarm / GNSS ionosphere chromatography method. The method comprises the following steps of: firstly, establishing an effective observation ray database of an inversion region; then establishing an effective observation ray coefficient matrix; according to the observation ray coefficient matrix and an electron density initial value obtained by an IRI empirical model, calculating to obtain a TEC predicted value; obtaining inversion area electron density correction based on Swarm satellite observation data; and finally, correcting the electron density of the inversion region and adding vertical and horizontal constraints, and taking the processed electron density value as a final ionized layer chromatography result. According to the method, the problems of low inversion precision, serious initial value dependence and the like caused by neglecting the influence of resolution selection of grids in different regions, horizontal constraint weight distribution, lack of single data source information and the like can be effectively solved.

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 Art

[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 plays a decisive role in the performance of technical systems such as short-wave communication and radar detection. Therefore, studying 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 needs to be solved urgently in current ionospheric tomography.

[0003] Currently, there are mainly three types of ionospheric tomography methods. The first type of method 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 it cannot meet large-scale measurements. The second type of method 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 type 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, so as to obtain the ionospheric electron density over a large range 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 dependable than grids with fewer or even 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 object 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 factors such as 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 in layers, 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 appears in the inversion process. Through the electron density detector carried by the Swarm satellite, the electron density of the movement 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 grids, effectively overcoming the problem of dependence on the initial value and improving the smoothness of the reconstructed ionosphere.

[0005] To achieve the above object, the specific technical solution adopted by the present invention is 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 chronological order by station and by satellite to establish an effective observation ray database, which is stored in a 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 at a low resolution in three dimensions of longitude, latitude, and altitude to obtain a number of three-dimensional volume grids. It is assumed that the electron density within each three-dimensional volume 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) in the inversion area with each volume grid in the tomographic model grid space. Since an excessive or zero number of effective observation rays in a single volume grid will directly affect the subsequent iterative update of the volume grid by the observation rays, the resolution of the volume grid space is adaptively adjusted in layers. The volume grids at the same altitude are divided into the same volume grid layer. The volume grid layer with a small number of effective observation rays passing through remains at the original low resolution, while the volume grid layer with a larger number of passing rays is changed to a higher resolution, and the intercept values of the effective observation rays in the volume grids of this layer are recalculated. Then, code according to the division order of the volume grids, and construct the observation ray coefficient matrix using the intercept lengths of the effective observation rays in the volume grids. Finally, according to the inversion area, inversion time, and volume 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. By the intercept values of the effective observation rays passing through each volume grid (i.e., the observation ray coefficient matrix established in the second step) and the initial values of the electron density of each volume 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 volume grids passed through by the observation ray in a certain proportion. Updating all the volume 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 volume 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 using 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 selected. 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 that 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 within 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 Chapman function relationship between the electron density and height of the three-dimensional grids in each grid bundle. During the fitting process, the Chapman function parameters are adjusted to make the generated Chapman function formula better fit 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 the vertical direction based on its 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 corrected by 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 within 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, it can effectively solve the problems of 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, effectively improving the real-time performance, accuracy, and stability of ionospheric tomography, providing 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 lack of single data source observation information on the accuracy of ionospheric tomography. Effectively solve the problems such as low accuracy of ionospheric tomography results and strong dependence of inversion on initial values. It can not only meet the requirements of ionospheric tomography in terms of time and space resolution, but also meet the requirements of tomography in terms of accuracy, real-time performance, and stability, 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 by using existing Swarm satellites / GNSS navigation system satellites. Compared with the existing tomography methods based on ground sounders and radio occultation, it can not only solve the problems of insufficient horizontal resolution and inability to meet large-scale measurement requirements, but also solve the problems of poor spatial continuity restricted by satellite orbits. Compared with the existing tomography methods based on GNSS navigation satellites, it can not only solve the problems such as the influence of the selection of three-dimensional grid resolution in different regions on the accuracy of iterative inversion, but also solve the problems of 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 in 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 vertical constraint is adopted to weaken the influence of observation noise and improve the vertical accuracy of the reconstructed ionosphere. Non-equal weight horizontal constraints are adopted for each layer of 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 and updated by the observation rays and the problem of dependence on the initial value of inversion, and 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 inversion result optimization based on Chapman function vertical constraint 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 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 in 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 in 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 in the inversion area can be directly extracted from observation files. The connection line between the ground-based receiving station and the GNSS navigation satellite is called an effective observation ray, and 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 belongs to the common knowledge in the field of positioning, and the calculation process will not be listed in detail here. Calculate the actual value of the high-precision ionospheric TEC of each effective observation ray separately by time, measurement station and 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 the 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 the 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: Divide the inversion area into grids and establish an 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 a reduction in tomography accuracy. Therefore, a hierarchical adaptive three-dimensional grid resolution considering the distribution of effective observation rays is adopted here, as Figure 2 shown. The specific process is as follows:

[0026] First, divide the inversion area into low-resolution three-dimensional grids. 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 low-resolution three-dimensional grids in the longitude, latitude, and altitude directions respectively. The low-resolution three-dimensional grid specifically refers to a three-dimensional grid 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 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 , and 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 6,378,137 meters and the semi-minor axis b is 6,356,752.314 meters. The specific process of this step of the algorithm is as follows:

[0028] 1) Calculate the straight-line equation of the effective observation rays. 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 rays 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 is perpendicular to the latitude plane and both pass through the Z-axis of the above-mentioned geodetic coordinate system. Determine the angle between the longitude plane to be found 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 point in the space region and the geocenter of the ellipsoid of revolution around the Z-axis of the geodetic coordinate system. Determine the latitude of the latitude plane to be found 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 altitude plane to be found 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 in formula (3).

[0041] Simultaneously solve the equations of the longitude plane, latitude plane, and altitude plane, namely equations (2), (3), and (4), with the space straight line equation (1) formed by the effective observation ray, and obtain the intersection coordinates of each effective observation ray in the inversion region 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, and 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 encoding order of the three-dimensional grid. 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 encoding 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 penetrate the three-dimensional grids in each layer to the total number of three-dimensional grids in that layer is obtained, 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 penetrate relatively densely. To improve the situation where the number of rays penetrating the three-dimensional grid is too dense and exceeds 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 refined high-resolution three-dimensional grid layers, 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 penetrate all the three-dimensional grids in each three-dimensional grid layer is counted. The top 50% of the three-dimensional grid layers with a higher number of times the effective observation rays penetrate 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] Re-numbering based on the original three-dimensional grid numbers is adopted 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, the coefficient matrix of the remaining low-resolution three-dimensional grids and the coefficient matrix of the high-resolution three-dimensional grids are combined 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 in 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 iterative algorithm in this case selects MART based on the principle of maximum entropy, which effectively improves the rank deficiency problem in the parameter solution process and ensures that the electron density does not appear negative during the iteration. 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 terminates. 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 region 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 Swarm satellite motion path database. The original observation data extracted mainly includes the motion time of the three satellites A, B, and C of 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 region. Compare the positions of the Swarm satellites in the motion path database ( ) with the longitude, latitude, and altitude ranges of the inversion region, and delete the data located outside the inversion region 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 crossed 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. The longitude, latitude, and height of all low-resolution three-dimensional grid layers in the three-dimensional grid space are subdivided into half of the original, and the subdivided three-dimensional grids are renumbered. 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 mixed-resolution three-dimensional grid space is adjusted to a unified high-resolution three-dimensional grid space.

[0061] 2) Fit and correct the electron density with the Chapman function. 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. Taking 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, the Chapman function corresponding to each three-dimensional grid beam is fitted by 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 least squares to solve 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 to fit the Swarm measured values better, improving the accuracy and reliability of ionospheric tomography. There is a ready-made function module in the Matlab platform for the non-linear weighted least squares algorithm that can be called. Here we only give a simple description of the algorithm 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] In the formula, 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] In the formula, , 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] Next, 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 operations of linear algebra, λ is the damping factor, I is the identity matrix, which belongs to the common 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 change of the objective function and 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 of each stereo grid bundle fitting 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 mixed resolution stereo grid space. For the original high-resolution stereo grid, its electron density value remains unchanged. For the eight high-resolution stereo grids subdivided from 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 subdivided from the low-resolution stereo grid numbered j. Thus, the electron density value of the mixed resolution stereo grid with Chapman function vertical constraint is obtained .

[0077] Step 7: Apply 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 in 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, 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 referred to as 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, and smoothing processing is performed accordingly. 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 of this industry 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 It includes the following steps: Step 1: Establish an effective observation ray database for the inversion area. The effective observation rays are the connections between 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 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 data arrangement order; 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. Consider all three-dimensional grids at the same height in the inversion area 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 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 using the ratio of the predicted TEC value to the measured TEC value, and perform multiple rounds of iterative updates 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 in the inversion area based on Swarm satellite observation data; Step 5: Add vertical and horizontal constraints to the electron density correction in 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 a hybrid grid and weighted horizontal constraints 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 a hybrid grid and weighted horizontal constraints according to claim 2, characterized in that In the said step 2-1, the method for performing a 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 altitude 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 in ascending order from small to large in the order of first along the longitude, then along the latitude, and finally along the altitude. The position of each three-dimensional grid is represented by the longitude , latitude , altitude .

5. A Swarm / GNSS ionospheric tomography method based on hybrid grid and weighted horizontal constraint 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 for adaptively adjusting to high-resolution three-dimensional grids is as follows: All the three-dimensional grids at the same height within the inversion region are regarded as the same three-dimensional grid layer. The number of times 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 into three-dimensional grids, 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 high-resolution three-dimensional grid with a longitude span of 1 degree, a latitude span of 1 degree, and a height span of 25 km.

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, select the three-dimensional grids with two intersections, and the two intersections are denoted as P1 and P2 , and calculate the intercept ΔL of the ray within the three-dimensional grid 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 by 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, and γ 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 grid 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 7, characterized in that, In step 4-3, the electron density correction method is as follows: Obtain the three-dimensional grid where it is located according to the position of the Swarm satellite in the motion path database, and replace the inversely calculated electron density value obtained from the GNSS observation data originally corresponding to 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, the mean value of all its measured values is taken 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 region: 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 subdivided; Then, the three-dimensional grids with continuous heights at the same longitude and latitude positions are divided into the same grid bundle. The inversely calculated electron density values and the central height values of the three-dimensional grids 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 by 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, the mean value of their electron densities is taken 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 hybrid grid and weighted horizontal constraint according to claim 11, characterized in that, The specific formula of the Chapman function is as follows: ； Wherein, is the central height value of the three-dimensional grid, is the electron density of the three-dimensional grid at the central height , that is, , is the maximum value of the electron density of each three-dimensional grid 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 region: 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 referred to as 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. 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: ； Wherein, is the total number of neighborhood three-dimensional grids, is the number of the th three-dimensional grid in the neighborhood three-dimensional grids passing 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