Hierarchical constraint ionosphere chromatography method based on hybrid interpolation and adaptive relaxation factor

Abstract

The invention discloses a hierarchical constraint ionosphere chromatography method based on hybrid interpolation and an adaptive relaxation factor. The method comprises the following steps: firstly, determining an inversion area, setting a virtual observation station in the inversion area, generating virtual STEC data of the virtual observation station through a hybrid interpolation method, and integrating actual STEC data and the virtual STEC data to obtain observation STEC data; calculating the intercept of each ray in the voxel according to the actual observation value in the inversion area and the coordinates of the virtual observation station, and constructing a coefficient matrix; obtaining an electron density initial value through an IRI-2016 model, and then iterating the electron density value by using the coefficient matrix and observation STEC data based on a multiplication algebraic reconstruction algorithm; adding constraints in the horizontal direction and applying constraints in the vertical direction; and taking the constrained electron density as an electron density initial value of a new round of iteration until the root-mean-square error is smaller than a preset threshold value.

Description

Technical Field

[0001] The present invention relates to the field of remote sensing inversion of navigation satellites, and particularly to a hierarchical constrained ionospheric tomography method based on hybrid interpolation and adaptive relaxation factors. Background Technique

[0002] The ionosphere is located at a distance of 60 km to 2000 km from the earth's surface. As an important part of the solar-terrestrial space environment, it has an important impact on human production and life. A large number of free electrons existing in the ionosphere form the geomagnetic space. On the one hand, it can reduce the effects of solar ultraviolet radiation and cosmic high-energy particles, protecting humans and other organisms on the earth from being harmed and enabling them to live normally; on the other hand, the ionosphere will have effects such as reflection, scattering, absorption, and refraction on electromagnetic waves; finally, due to the influence of solar activities and the geomagnetic field, the ionosphere will show irregular perturbation phenomena and produce drastic changes, which may lead to the instability or even damage of aerospace, communication, navigation, and many ground technology systems, resulting in spacecraft malfunctions or damage, wireless communication interruptions, satellite navigation service failures, etc. Therefore, studying the fine structure and spatio-temporal variation characteristics of the ionosphere, establishing an accurate and reliable ionosphere model, and realizing ionospheric refraction error correction are of great significance for improving the positioning accuracy of navigation satellite systems, ensuring the stability of wireless communication, and guaranteeing the safety of aerospace activities.

[0003] Ionospheric tomography imaging technology combines computer tomography (CT) technology with satellite radio signals and uses the total electron content on the ray path to obtain the spatio-temporal distribution of ionospheric electron density. The ionosphere reconstructed by three-dimensional ionospheric tomography technology can simultaneously reflect the horizontal and vertical structure information of the ionosphere. However, due to objective conditions such as a small number of ground receivers and uneven distribution, as well as limited observation angles, three-dimensional ionospheric tomography based on GNSS (Global Navigation Satellite System) is a typical ill-posed problem.

[0004] Currently, ionospheric tomography algorithms are mainly divided into two categories. The first category is non-iterative algorithms mainly represented by the singular value decomposition algorithm, the Kalman filtering method, the regularization method, etc. Non-iterative algorithms do not rely on any initial values during inversion, but the solution methods are relatively complex. Moreover, due to the existence of GNSS observation noise and discrete errors, the obtained inversion results are usually only approximate solutions. In the case where the distribution of observation stations is uneven or the number of observation rays is insufficient, that is, when the rank deficiency of the coefficient matrix of the tomography equation is relatively serious, the obtained solutions will deviate significantly from the true values, resulting in the distortion of the reconstructed ionosphere.

[0005] The second category is iterative reconstruction algorithms mainly represented by the additive algebraic reconstruction algorithm and the multiplicative algebraic reconstruction algorithm. When using the iterative reconstruction algorithm for ionospheric inversion, the ionospheric model is used to provide the initial parameter values. The residual or numerical ratio between the iteratively reconstructed slant total electron content (STEC) and the slant total electron content calculated from the observed data is used to circularly correct the parameter to be estimated, so as to obtain the reconstruction result of the three-dimensional distribution of the ionospheric electron density. However, the existing algorithms have a strong dependence on the initial values. Moreover, due to the insufficient total amount of existing observed data, the electron density reconstructed by the iterative reconstruction algorithm for the voxels without observed rays still retains the initial values, resulting in a reduction in the accuracy of the final three-dimensional ionospheric inversion result. At the same time, the existing algorithms do not consider the influence of the relaxation factor increasing with the number of iterations during the iterative reconstruction process, but take the relaxation factor as a constant, thus leading to a reduction in the algorithm efficiency and limited ionospheric reconstruction accuracy. In addition, the existing iterative reconstruction algorithms often use a single method for vertical constraint, ignoring the significant differences in the electron density distribution characteristics at different heights of the ionosphere, resulting in the inversion result of the electron density distribution in the vertical direction not conforming to the actual physical characteristics, and reducing the reliability and applicability of the inversion result. Summary of the Invention

[0006] The purpose of the present invention is to overcome the deficiencies of the prior art and propose a hierarchical constrained ionospheric tomography method based on hybrid interpolation and adaptive relaxation factor. By studying the effects of observed data completion, adaptive adjustment of the relaxation factor, and different vertical constraints on the three-dimensional inversion accuracy of the ionospheric electron density value, it effectively solves problems such as the strong dependence of traditional iterative reconstruction algorithms on initial values and sensitivity to noise. It can not only meet the requirements of the spatial resolution of ionospheric tomography, but also meet the requirements of tomography accuracy, stability, etc., providing strong support for ionospheric tomography inversion.

[0007] To achieve the above purpose, the specific technical solution adopted by the present invention is as follows:

[0008] A hierarchical constrained ionospheric tomography method based on hybrid interpolation and adaptive relaxation factor, comprising the following steps:

[0009] Step 1: Generate virtual observation station data based on hybrid interpolation method. Determine the size of the inversion area and the voxel resolution in all three directions. According to this resolution, the ionospheric space in the inversion area is divided into regular discrete grids. Each grid is called a voxel, and the voxels are numbered from large to small in the order of longitude, latitude, and finally altitude. First, use the actual observation stations in the inversion area to receive GNSS observation data, mainly including raw observation values such as pseudorange and carrier phase, and solve the high-precision actual STEC data of each station. Secondly, a series of virtual observation stations are selected in the inversion area with a longitude and latitude resolution of 0.5°×0.5°. With each virtual observation station as the center, several actual observation stations within a certain cutoff range are selected as interpolation stations. Finally, hybrid interpolation is used to generate virtual STEC data based on the number of selected interpolation stations. If there are more than six observation stations, the nearest six are selected as interpolation stations, and the multi-quadratic surface fitting method is used to fit and generate the data of the virtual observation stations. If there are more than three but less than six observation stations, it is necessary to divide the circle with the cutoff range as the radius into sectors (interval 2π / 3), select the actual observation station closest to the virtual observation station in each sector as the interpolation station, and ensure that the virtual observation station is in the triangle formed by the three interpolation stations to optimize the geometric distribution. With these three observation stations as interpolation stations, triangular linear interpolation is used to generate virtual STEC data. If there are two observation stations, the sectors are also divided to ensure that the two selected stations are distributed in different azimuth directions of the virtual station to ensure the stability and accuracy of the interpolation. With these two observation stations as interpolation stations, inverse distance weighted interpolation is used to generate virtual STEC data. Virtual observation stations with less than two observation stations within the cutoff range are discarded. Because the actual STEC data are directly calculated from the original observation data, and the virtual STEC data are generated by interpolation from the actual STEC data, the actual STEC data and virtual STEC data are integrated and called observed STEC data.

[0010] Step 2: Solve the coefficient matrix. First, obtain the instantaneous coordinates of each satellite through the GPS (Global Positioning System) satellite precise ephemeris file, and calculate the spatial position of each observation ray and the corresponding straight line equation by combining the coordinate file of the ground observation station in the inversion area and the coordinates of the virtual observation station. Solve the coordinates of the intersection of the observation ray and each voxel, and determine whether the intersection is within the inversion area, and discard the intersection outside the inversion area. Secondly, according to the total number of intersections between each ray and the height plane of all voxels, screen out the effective rays that pass through the bottom and top of the inversion area at the same time during the inversion time, and eliminate the rays whose number of intersections is less than the number of voxels in the height direction to ensure the correctness of the tomography results. Finally, calculate the intercept of each effective ray in the voxel it passes through based on the intersection of the ray and the voxel, and store the intercept value of each ray in each voxel in the coefficient matrix in sequence according to the voxel coding order.

[0011] Step 3: Iteratively correct the electron density value based on the MART algorithm with an adaptive relaxation factor. First, use the open-source code of the IRI-2016 model in Fortran language and solar magnetic index data. After compilation and conversion, input the inversion time, longitude, latitude, and altitude range of the inversion area to obtain the electron density distribution corresponding to the inversion area as the initial value of the iterative algorithm. Second, substitute the initial electron density value, the coefficient matrix solved in the second step, and the observed STEC data into the tomography equation, and use the MART iterative formula to invert the electron density value. Based on the ratio of the reconstructed STEC to the observed STEC calculated from the electron density value obtained in each iteration, the electron density value is corrected accordingly until the reconstruction result converges. Finally, adjust the relaxation factor in the iterative formula, and adaptively perform the next iteration based on the residual value between the reconstructed STEC and the observed STEC in each iteration, so as to effectively accelerate the convergence speed of the inversion and avoid overshoot when the residual is small, ensuring the stability of the inversion.

[0012] Step 4: Apply horizontal smooth constraint and vertical stratification constraint. First, perform mean smoothing constraint on the output value of each iteration of the MART algorithm in Step 3 in the horizontal direction. For the voxels on the same height plane in the horizontal direction, calculate the mean value of the electron density values of all voxels within the neighborhood range, and use this mean value as the corrected electron density value of the central voxel, so as to control the difference between the electron density values of the central voxel and adjacent voxels and ensure the smoothness of the reconstructed ionosphere in the horizontal direction. Second, based on the difference in the distribution characteristics of electron density values at different heights, adopt a stratification constraint strategy to perform vertical constraint on the iterative output value applied with horizontal constraint. In the altitude range of 100 - 200 km, the change of electron density is relatively gentle, and mean smoothing constraint is performed on the iterative electron density value. The same as the mean smoothing constraint method in the horizontal direction, calculate the mean value of the electron density of all voxels within the neighborhood as the corrected value of the central voxel. In the altitude range of 200 - 400 km, vertical constraint is performed based on the Chapman function. Use the least squares method and fit the parameters of the Chapman function according to the electron density values and altitude values of the voxel bundles at different heights at the same horizontal position. Substitute different altitude values into the Chapman function with the fitted parameters to obtain the constrained electron density value. In the altitude range above 400 km, vertical constraint is performed based on the empirical orthogonal function. Based on the empirical data set of the vertical electron density profile (EDP) constructed by the IRI-2016 model, use the singular value decomposition technique (SVD) to extract the empirical orthogonal basis functions (EOF). Select the first three EOFs and reconstruct the vertical electron density distribution through the linear combination of EOF coefficients to achieve vertical direction constraint on the inverted electron density value in the area above 400 km.

[0013] Finally, the electron density values that have completed the horizontal and vertical constraints are used as the initial values for the next iteration, and the next round of iteration is carried out. The process is repeated until the reconstruction result converges, and then the final electron density value is output.

[0014] After being processed by the above three-dimensional ionospheric tomography algorithm, it can effectively solve the problems of low inversion accuracy, poor stability, and sensitivity to noise caused by the existing methods ignoring the sparse satellite observation data, fixed relaxation factor, and differences in the electron density distribution characteristics at different altitudes, effectively improve the accuracy, stability, and efficiency of three-dimensional ionospheric tomography, and provide guarantee for the delay correction of satellite positioning.

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

[0016] By using the method of the present invention, the three-dimensional electron distribution of the ionosphere can be tomographically reconstructed with higher precision using existing GNSS satellite observation stations. Compared with existing non-iterative algorithms, it can not only solve problems such as relatively complex solution methods, and the inversion results are usually approximate solutions due to observation noise and discrete errors, but also solve problems such as when the distribution of observation stations is uneven or the number of rays is insufficient, the obtained solution will deviate seriously from the true value, resulting in distorted reconstruction of the ionosphere. Compared with existing iterative algorithms based on MART, it can not only solve problems such as strong dependence on initial values due to insufficient total amount of existing observation data, but also solve problems such as the relaxation factor being taken as a constant during the iterative process, resulting in reduced algorithm efficiency and limited ionosphere reconstruction accuracy. At the same time, applying vertical constraints using a layering strategy can ensure the accuracy and stability of the reconstruction of the electron density distribution in the vertical direction, thus meeting the requirements of high-precision and high-stability three-dimensional ionospheric tomography. By screening effective virtual observation stations based on the number of actual observation stations within the cut-off range centered on the virtual observation station, spatial interpolation errors are reduced on the premise of ensuring the accuracy of virtual data. For virtual observation stations with different numbers of actual observation stations within the cut-off range, three different interpolation methods, namely multiquadric surface fitting, triangular linear interpolation, and inverse distance weighted interpolation, are used to generate virtual observation data, reducing the number of voxels without ray penetration and improving the quality of three-dimensional ionospheric tomography. Using the International Reference Ionosphere IRI-2016 ionosphere reference model as the initial value of the iteration, iterative reconstruction is carried out through the multiplicative algebraic reconstruction algorithm to ensure the non-negativity of the solution of the electron density value. For the relaxation factor during the iterative reconstruction process, it is adjusted according to the residual between the reconstructed STEC and the observed STEC after each iteration, accelerating the convergence of the algorithm and enhancing the stability of three-dimensional ionospheric tomography. Introducing a mean smoothing constraint in the horizontal direction improves the accuracy of ionospheric iterative reconstruction. By applying vertical constraints to different heights of the ionosphere using a layering strategy, the accuracy and reliability of three-dimensional ionospheric tomography inversion are significantly improved. Taking the electron density value after completing the constraints as the initial value of the next iteration and repeating until the reconstruction result converges, ensuring the high precision, stability, and reliability of three-dimensional ionospheric tomography. Description of the Drawings

[0017] Figure 1 Flowchart of generating virtual observation station data based on a hybrid interpolation method.

[0018] Figure 2 Flowchart of solving the coefficient matrix based on the coordinates of actual and virtual observation stations.

[0019] Figure 3 Flowchart of iteratively correcting the electron density based on the MART algorithm with an adaptive relaxation factor.

[0020] Figure 4 Flowchart of applying horizontal smooth constraints and vertical layering constraints based on the MART algorithm. Detailed Implementation Manner

[0021] The present invention will be described in detail below in conjunction with 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 and features in the embodiments of the present invention can be combined with each other.

[0022] A hierarchical constrained ionospheric tomography method based on hybrid interpolation and adaptive relaxation factor includes the following steps:

[0023] Step 1: Determine the range of the inversion area and the voxel resolution in the longitude, latitude, and height directions. Divide the ionospheric space of the area to be inverted into regular discrete grids, and each grid is a voxel (specifically referring to each regular discrete grid in the discretized inversion area space), and number the voxels according to their spatial positions. Assume that the ionospheric area is divided into α, β, and γ voxels in the longitude, latitude, and height directions respectively, then the entire area is divided into α×β×γ voxels, and then number the voxels in the order of longitude first, then latitude, and finally height.

[0024] Step 2: Generate virtual STEC data by using hybrid interpolation. This process can be divided into three steps. The first step is to calculate the STEC data of the actual observation stations, the second step is to select the available actual observation stations and determine their corresponding interpolation stations, and the third step is to generate the data of the virtual observation stations by using hybrid interpolation according to the number of interpolation stations.

[0025] Specifically, as Figure 1 shown, the method for generating virtual STEC data by using hybrid interpolation is as follows:

[0026] First, receive GNSS observation data through the actual observation stations in the inversion area, mainly including original observations such as pseudorange, carrier phase, satellite and receiver clock errors, etc., and calculate the actual STEC data of each station through these data. This process belongs to the common knowledge in the field of ionospheric tomography, and the calculation process will not be listed in detail here.

[0027] Secondly, a series of virtual observation stations distributed in a grid are taken within the inversion area, with a distribution density of 0.5°×0.5° in longitude and latitude. Actual observation stations are selected as interpolation stations within a circle centered on each virtual observation station with a cut-off distance as the radius. The cut-off distance is the maximum distance between a virtual observation station and available interpolation stations. If the location of an actual observation station exceeds the cut-off distance, there will be a large interpolation error when it is used to interpolate and generate virtual observation data. Therefore, the determination of the cut-off distance is a compromise between the number of available interpolation stations and the interpolation error. In this embodiment, the cut-off distance is set to 150 km. When there are more than six actual observation stations within the cut-off range, the six closest to the center are selected as interpolation stations. When there are 3 - 5 actual observation stations within the cut-off range, three sector areas are drawn at intervals of 2π / 3 in the circle, and then one actual observation station closest to the virtual observation station in each sector is selected as an interpolation station, while ensuring that the virtual observation station is within the triangle formed by the connection lines of the three interpolation stations. When there are two actual observation stations within the cut-off range, four sector areas are divided at intervals of π / 4 in the circle to ensure that the two observation stations are in non-adjacent sectors; otherwise, the virtual observation station is discarded. When the number of actual observation stations within the cut-off range is less than two, the virtual observation station is discarded.

[0028] Finally, different interpolation methods are used for virtual observation stations with different numbers of interpolation stations. The hybrid interpolation process is divided into three steps:

[0029] If the number of actual observation stations within the cut-off range ≥ 6, the six actual observation stations closest to the virtual observation station are taken, and multiquadric fitting is used to generate virtual STEC data;

[0030] Specifically, in this embodiment, multiquadric fitting is used to generate virtual STEC data for virtual observation stations with six interpolation stations. Multiquadric fitting is usually used to interpolate data in a multi-dimensional space, and its theoretical basis is that any smooth mathematical surface can always be approximated with arbitrary accuracy by the sum of a series of regular mathematical surfaces. The specific process is as follows: First, the weights of the interpolation stations are determined from the actual STEC data and the distance between two stations, and the calculation formula is as follows:

[0031] (1)

[0032] (2)

[0033] In the formula, STEC p represents the actual STEC data of the p-th interpolation station, c q is the weight of each quadratic surface, and the sum of the weights is assumed to be 0. v is the number of interpolation stations. In this embodiment, the MSF method is used to generate virtual STEC data for six interpolation stations, that is, v = 6. g pqThe distance parameter between two interpolation stations can be expressed as a hyperboloid, a paraboloid, or a conical surface. Since it is difficult to determine the optimal constant term a of the hyperboloid and there are certain problems in the calculation of the paraboloid, a conical surface is usually selected for calculation.

[0034] Then, the virtual STEC data of the virtual observation station is calculated by solving the weights and the distances between the virtual observation station and each interpolation station. The calculation formula is as follows:

[0035] (3)

[0036] In the formula, STEC E is the virtual STEC data of the calculated virtual observation station, g 0q is the distance parameter between interpolation station q and the virtual observation station, and c q and v are the same as those in formula (1).

[0037] If the number of actual observation stations is 3 - 5, after dividing the sectors, select the three stations closest to the virtual observation station, and ensure that the virtual station is inside the triangle formed by the three stations' connection lines. Triangular linear interpolation is used to generate the virtual STEC data;

[0038] Specifically, in this embodiment, triangular linear interpolation is used to generate the virtual STEC data of the virtual observation station with three interpolation stations. Triangular linear interpolation can generate the virtual data of a virtual observation station located in the triangle formed by the connection lines of three interpolation stations. Using the coordinates of interpolation stations A, B, C, and virtual observation station E, the interpolation can be calculated as follows:

[0039] (4)

[0040] In the formula, D is the intersection point of line AE and line BC, AD represents the distance between interpolation station A and point D, AE represents the distance between interpolation station A and virtual station E, BC represents the distance between interpolation station B and interpolation station C, BD represents the distance between interpolation station B and point D, DC represents the distance between point D and interpolation station C, ED represents the distance between virtual station E and point D, and STEC E represents the virtual STEC data of the required virtual station, and STEC A , STEC B , STEC C are the actual STEC data of interpolation stations A, B, and C respectively.

[0041] If the number of actual observation stations is 2, inverse distance weighted interpolation is used to generate the virtual STEC data after dividing the sectors;

[0042] Specifically, in this embodiment, inverse distance weighted interpolation is used to generate the virtual STEC data of the virtual observation station with two interpolation stations. The calculation formula is as follows:

[0043] (5)

[0044] Wherein, STEC E represents the virtual STEC data of the virtual station to be obtained, and STEC r represents the actual STEC data of the r-th interpolation station, d r is the distance from the virtual observation point to the r-th interpolation station, w is the power parameter, and in this embodiment, w = 2 is taken. v is the same as in formula (1) and both represent the number of interpolation stations, so here v = 2 is taken.

[0045] Virtual observation stations with the number of actual observation stations < 2 are discarded.

[0046] After generating the virtual STEC data by interpolating the actual STEC data, the actual STEC data and the virtual STEC data are integrated, and then collectively referred to as the observed STEC data.

[0047] It should be noted that generating the virtual STEC data actually generates the observation rays between the virtual observation stations and the satellites. Coupled with the actual observation rays, the total number of rays increases, and the number of voxels penetrated by the rays in the inversion area increases accordingly. In this way, more values can be solved in the coefficient matrix, and the reconstruction accuracy of the same area will be higher than before. The "virtual STEC data" is already included in the observed STEC data, so what is mentioned in the subsequent steps is the observed STEC data

[0048] Step 3: Calculate the intercept value of each effective observation ray passing through each voxel and construct the coefficient matrix. It should be noted that the ray refers to the connection line between the observation station and the satellite.

[0049] This process can be divided into three steps. The first step is to find the straight-line equation corresponding to each ray, the second step is to find the intersection points of the straight line and the six faces of each voxel, and the third step is to calculate the intercept length of the straight line within each voxel, so as to construct the coefficient matrix.

[0050] Specifically, as Figure 2 shown, the method for constructing the coefficient matrix is as follows:

[0051] First, find the straight-line equation corresponding to each ray. Obtain the instantaneous coordinates of the satellite through the precise ephemeris file, and at the same time read the coordinates of all actual observation stations within the inversion area and the coordinates of the virtual observation stations selected in the second step. In the space rectangular coordinate system (WGS84 coordinate system), let the ground station coordinates be (X t , Y t , Z t ), and the satellite coordinates be (X s , Y s , Z s), then the ray equation can be expressed as follows:

[0052] (6)

[0053] Where X, Y, Z, and k are the variables in the equation.

[0054] Secondly, solve the coordinates of the intersection of the ray and each voxel. Each voxel is composed of six planes, including two longitude planes, two latitude planes and two altitude planes. The equations corresponding to the intersection of a straight line with different curved surfaces are different. The equations corresponding to the three curved surfaces are given below. The coordinates of the intersection are obtained by calculating the intersection of the straight line and the curved surface.

[0055] Longitude plane: refers to the plane passing through the Z axis of the coordinate system. The angle between this plane and the zero meridian plane is longitude J. The equation is as follows:

[0056] (7)

[0057] Where X and Y are the variables in the equation.

[0058] Latitude plane: refers to the conical surface with the center of mass of the earth as the vertex, and the angle between its generatrix and the Z axis is the latitude W. The equation is as follows:

[0059] (8)

[0060] Where X, Y, and Z are the variables in the equation.

[0061] Height plane: refers to the surface parallel to the earth's spherical surface. When the height is H, its equation is expressed as:

[0062] (9)

[0063] In the formula, R represents the radius of the earth, which is 6372 kilometers, and X, Y, and Z are the variables of the equation.

[0064] Combine equation (6) with equation (7), equation (8), and equation (9) to find the intersection coordinates (X i , Y i , Z i ). Determine whether the intersection point is in the inversion area. If not, discard the intersection point. After obtaining the coordinates of all intersection points, filter out valid rays by calculating the total number of intersections between each ray and all voxel height planes. When the number of intersection points is equal to the number of voxels in the height direction, it means that the ray passes through the entire inversion area, and the ray can be used for tomography; when the number of intersection points is less than the number of voxels in the height direction, it means that the ray does not pass through the entire inversion area and cannot be used as tomography data. The ray should be removed from the inversion data, otherwise it will lead to incorrect tomography results.

[0065] Finally, calculate the intercept and construct the coefficient matrix based on it. After obtaining the intersection coordinates of the straight line and each voxel, filter out the two intersection points P1 and P2 that fall on the surface of each voxel. The coordinates of P1 are (X1, Y1, Z1) and the coordinates of P2 are (X2, Y2, Z2). Calculate the intercept ∆L of the ray within each voxel according to the distance formula between two points. The calculation formula is as follows:

[0066] (10)

[0067] Store the intercept values of each ray within each voxel in the coefficient matrix in the order of voxel encoding. The intercept of the voxel not penetrated by the ray is 0. Each row of the coefficient matrix represents a ray, and each column corresponds to a voxel, that is, the arrangement order of the elements in each row corresponds to the encoding order of the voxels. Obtain the coefficient matrix in the inversion equation through the above processing. In the subsequent steps, it will be represented by A m×n where m is the total number of observed rays, including actual observed rays and virtual observed rays, and n is the total number of voxels.

[0068] Step 4: Obtain the initial value of electron density through the IRI-2016 model. Call the open-source Fortran code of the IRI-2016 ionosphere model, input the inversion time, the longitude and latitude of the inversion area, and the height range to obtain the electron density values of a series of coordinate points. At the same time, set the grid spacing of the output data in the longitude and latitude directions and the height direction, so that the output electron density values exactly correspond to the central point coordinates of each voxel, and use this as the average initial value of electron density within each voxel. In the subsequent steps, it is called .

[0069] Step 5: Iteratively invert the electron density based on the multiplicative algebraic reconstruction algorithm.

[0070] Specifically, as Figure 3 shown, the multiplicative algebraic reconstruction algorithm gradually improves the estimated value of each voxel in an iterative manner based on the initial value of electron density obtained in Step 4. During the correction, use the ratio of the reconstructed STEC data calculated from the electron density value calculated in each iteration to the observed STEC data (including actual STEC data and virtual STEC data), and correct the ionospheric electron density (IED) one by one through multiplication. The specific iterative correction formula is as follows:

[0071] (11)

[0072] In the formula, k represents the number of iterations, is the electron density iteration value of the th voxel at the th iteration, is the electron density value of the voxel before the k-th iteration. i represents the i-th observation ray, and a is the intercept coefficient matrix A i of the m×n k-th row, that is, the intercept value of the i-th ray passing through each voxel. m is the total number of observation rays, including actual and virtual observation rays, and n is the total number of voxels. a ij is the intercept distance of the i-th observation ray passing through the j-th voxel. λ (k-1) is the relaxation factor of the k-th iteration, where 0 < λ < 1.

[0073] Specifically, the method for obtaining the value of λ is as follows:

[0074] It can be understood that the larger the value of the relaxation factor, the higher the iteration efficiency, that is, the faster the algorithm converges, but the greater the influence of noise; the smaller the relaxation factor, the smoother the iteration process, with less noise interference, but the number of iterations increases and the operation time becomes longer. In traditional methods, the relaxation factor usually takes a fixed value, making it difficult to balance between fast convergence in the early stage and stable optimization in the later stage. Therefore, an adaptive relaxation factor based on the residual is adopted here. The specific process is as follows:

[0075] First, after obtaining the electron density values of all voxels after the k-th iteration through step 5, the reconstructed STEC data based on the electron density iteration values of this iteration can be calculated. The reconstructed STEC expression corresponding to each ray is as follows:

[0076] (12)

[0077] In the formula, represents the reconstructed STEC data of the i-th observation ray after the k-th iteration. i, j, k, a ij , are the same as those in formula (11), and n is the total number of voxels. After obtaining the reconstructed STEC data, the residual (i.e., the difference between the reconstructed STEC data and the observed STEC data) can be calculated according to the following formula:

[0078] (13)

[0079] In the formula, STEC i represents the STEC value of the i-th observation ray.

[0080] Next, the relaxation factor λ for the next iteration can be generated using the relaxation adjustment formula based on the residual adaption. It should be noted that λ (k) ,(k) It corresponds to the (k + 1)-th iteration, and the specific calculation formula is as follows:

[0081] (14)

[0082] In the formula, λ0 is the initial value, and in this embodiment, λ0 = 1 is taken; μ is a positive parameter that controls the exponential decay rate, and in this embodiment, μ = 0.12 is taken; λ min is the lower limit of the relaxation factor value, and in this embodiment, λ min = 0.15, is the initial residual vector, is the residual vector at the k-th iteration. Because when → 0, → ∞, which causes λ (k) → 0. This is reasonable in theory, but in practical applications, this may cause λ (k) to decay too fast, affecting the later convergence. Therefore, a lower limit of the value is added to the formula to avoid this situation. Through the above processing, the relaxation factor based on the residual adaption can be obtained and used in the next iteration formula.

[0083] Step 6: Add horizontal constraints to the electron density values using the smoothing constraint method. The smoothing constraint corrects the electron density value of the central voxel by using the mean value of the electron density values of all voxels within the neighborhood range during the iteration process to ensure the smoothness of the reconstructed ionosphere in the horizontal direction. This method is essentially similar to the mean smoothing idea in digital image processing methods.

[0084] Specifically, as Figure 4 shown, the specific process of the horizontal mean smoothing constraint algorithm is as follows:

[0085] First, the voxels in the inversion area are divided into different horizontal layers according to different heights, and voxel smoothing is performed on each horizontal layer after each iteration. For the electron density value of any voxel in a certain horizontal layer, the smoothing method is as follows:

[0086] (15)

[0087] In the formula, and respectively represent the electron density values of the j-th voxel before and after smoothing after the k-th iteration correction; is the electron density value of the l-th voxel within the neighborhood range of the j-th voxel before smoothing; k represents the number of iterations, n is the total number of voxels, and I is the total number of voxels within the neighborhood range.

[0088] In this embodiment, according to the different positions of the voxels, the value of I is divided into the following three cases:

[0089] 1) When the voxel is inside a horizontal layer, the electron density values of itself and the four voxels in front, behind, left, and right are averaged. At this time, I = 5 in Equation (15).

[0090] 2) When the voxel is at the edge of a horizontal layer, the electron density values of itself and the three adjacent voxels around it are averaged. At this time, I = 4 in Equation (15).

[0091] 3) When the voxel is at the corner, the electron density values of itself and the two adjacent voxels around it are averaged. At this time, I = 3 in Equation (15).

[0092] Through the above processing, a horizontal voxel smoothing constraint is performed at each iteration, which can effectively alleviate the error brought by MART, and at the same time provides a more stable and reliable input value for applying the vertical constraint in the next step.

[0093] Step 7: Apply vertical constraints to the electron density values using a layering strategy based on the differences in the distribution characteristics of electron density values at different heights. As Figure 4 shown, the algorithm mainly consists of three parts:

[0094] 1) Use the smoothing constraint method to add vertical constraints to the voxels in the height range of 100 - 200 km. The voxels in the height range of 100 - 200 km are divided into different voxel layers according to different longitudes, and smoothing constraints are performed on each layer separately. For the electron density value of any voxel in a certain longitude layer, the specific operation of the smoothing constraint is the same as that in Step 7.

[0095] 2) Use the Chapman function to add vertical constraints to the voxels in the height range of 200 - 400 km. The Chapman function is a function model used to represent the vertical profile of the ionospheric electron density, and it gives an empirical model of the vertical distribution of the ionosphere. The overall vertical distribution of the ionosphere can be expressed as:

[0096] (16)

[0097] In the formula, x(h) is the electron density value at height h, h is the height at which the electron density is currently calculated, N m F2 is the peak density of the F2 layer of the ionosphere, which can be abbreviated as N m , h m F2 is the peak height of the F2 layer of the ionosphere, h m That is, h m The abbreviation of F2, the F2 layer is the area above 200 km in height in the ionosphere, H1 is the scale height at the bottom of the ionosphere, and H2 is the scale height at the top of the ionosphere.

[0098] According to the spatial positions of voxels, voxels located at consecutive heights at the same horizontal position are divided into a voxel string. Repeating this step divides all voxels in the 200 - 400 km range in the entire inversion area into a series of voxel strings. Taking the MART reconstructed electron density values with horizontal mean smoothing constraints and the central height values of voxels corresponding to all voxels in each voxel string as input values, the N of the Chapman function is obtained through least - squares fitting. m , h m , H1, H2 four parameter values, and then substituting the parameter values obtained by fitting back into the Chapman function formula to reconstruct the vertical distribution function of this voxel number. Finally, by inputting the height values of each voxel in this voxel bundle, the corrected electron density value can be obtained, thus achieving vertical constraint on the output value of the MART algorithm using the Chapman function.

[0099] Among them, there is a ready - made function module for the nonlinear least - squares fitting algorithm in the Matlab platform that can be called. Here, only a simple algorithm description is given as follows:

[0100] For a nonlinear fitting problem, the goal is to obtain the optimal model parameters by minimizing the sum of the squares of the errors between the observed data and the model - predicted data. Assume the model is in the form of a nonlinear function:

[0101] (17)

[0102] In the formula, y is the observed value, f(x,θ) is a nonlinear function of the independent variable x and the parameter θ to be estimated, and ϵ is the error term representing the observation error. Through the fitting algorithm, we hope to find the parameter θ such that the model - predicted value is as close as possible to the observed data. The goal of nonlinear least - squares fitting is to minimize the following objective function:

[0103] (18)

[0104] In the formula, y δ is the δ - th observed value, x δ is the δ - th independent variable value, f(x δ ,θ) is the predicted value of the model, G represents the total amount of observed values, and θ is the parameter to be estimated. The specific implementation of the algorithm can be divided into the following three processes:

[0105] First, according to the parameter estimate value θ k , calculate the error between the predicted value of the model f(x,θ) and the actual observed value. For each data point x δ , the prediction error is:

[0106] (19)

[0107] In the formula, ϵ pred,δdenotes the prediction error for the δ-th data point, y obs,δ denotes the observed data of the δ-th data point, f(x δ ,θ k ) denotes the predicted value calculated for the δ-th data point with the parameter estimate value θ k .

[0108] Secondly, linearize the objective function through Taylor expansion to transform the non-linear problem into a linear problem. Calculate the Jacobian matrix J, which represents the partial derivative of the model function f(x,θ) with respect to the parameter θ. The expression is as follows:

[0109] (20)

[0110] Use the Levenberg-Marquardt or Gauss-Newton method to update the parameter θ, and adjust the parameter by calculating the gain matrix K. The calculation formula for the gain matrix K is as follows:

[0111] (21)

[0112] In the formula, J T represents the transpose of the Jacobian matrix J. Then, use the gain matrix K to update the parameter estimate. The update formula is:

[0113] (22)

[0114] In the formula, θ k is the current parameter estimate value, θ k+1 is the updated parameter estimate value, K is the gain matrix, and ϵ pred is the prediction error for all data points.

[0115] Finally, perform convergence judgment. If the change in ϵ pred is very small or the objective function S(θ) reaches a predetermined threshold, it means that the model has converged. At this time, the final θ k+1 is the optimal fitting parameter, denoted as θ*, terminate the algorithm and output the optimal fitting parameter θ*. Through the above process, the optimal non-linear model parameters can be obtained by using the non-linear least squares fitting algorithm.

[0116] 3) Add vertical constraints to voxels in the range above 400 km based on Empirical Orthogonal Function. The Empirical Orthogonal Function (EOF), also known as Principal Component Analysis, is a powerful method for extracting the characteristics of a data set. It can eliminate redundant information in the original data, extract the main information in the data set, and complete data analysis and dimensionality reduction. When a specific data set is given, the Singular Value Decomposition technique (SVD) can be used to decompose the data into several orthogonal basis functions, and most of the information of the data set can be approximately expressed based on the linear combination of these orthogonal basis functions. Specifically as follows:

[0117] (23)

[0118] In the formula, X P×Q represents the original data set, P is the number of data samples, Q is the number of variables, 、 and represent the left singular matrix, the singular value matrix, and the right singular matrix respectively. Through the SVD operation, EOF is stored in the U P×P matrix, and each column of this matrix represents an EOF; Σ P×Q is a diagonal matrix, where the non-zero elements are singular values and are arranged in size. The column corresponding to the largest singular value in the U P×P matrix is EOF-1, the next is EOF-2, and so on. The contribution rates of different EOFs can be calculated as:

[0119] (24)

[0120] In the formula, C W is the contribution rate of the Wth EOF, Σ W,W is the Wth singular value, is the sum of all singular values, M represents the total number of singular values, and Q is the same as in formula (23).

[0121] The specific process of adding vertical constraints based on EOF is as follows:

[0122] First, extract the vertical pixel bundle in the height range above 400 km in the inversion area. To obtain appropriate vertical direction constraints, in this embodiment, an empirical data set is established based on the vertical electron density profile (EDP) provided by the IRI-2016 model, and EOF is obtained from the empirical data set using the Singular Value Decomposition technique, so that it can be used to express the EDP to be solved. At this time, P in formula (23) is the number of voxels in a vertical direction, and Q represents the number of empirical electron density EDPs obtained based on the IRI-2016 model. Therefore, based on the obtained EOF, the EDP to be solved can be expressed as:

[0123] (25)

[0124] In the formula, E represents the EDP to be solved, and e u represents the u-th EOF, and φ u represents the projection of E on the u-th EOF mode e u . N represents the total number of EOFs selected to participate in the constraint. According to the contribution rates of different EOFs and their constraint information, in this embodiment, the first three EOFs are selected for vertical constraint, that is, N = 3. Through the above process, the EDP of each vertical pixel bundle can be obtained, that is, the vertical constraint on the reconstructed electron density value is completed.

[0125] So far, the vertical constraint has been applied to the ionospheric inversion regions at each altitude, making the final tomographic result closer to the true value in the vertical direction.

[0126] Step 8: Replace the IED provided by the IRI-2016 model with the electron density values applied with horizontal and vertical constraints as the initial value of the new round of iteration, and substitute it into the MART algorithm formula. Repeat steps 5 to 7 until the root mean square error between the reconstructed STEC data of each ray of the iterative inversion result and the observed STEC data is less than the tolerance σ, and then stop the iteration.

[0127] Step 9: Save and output the inversion result, and calculate the root mean square error between the inversion result and the electron density value data measured by the ionosonde, so as to judge the accuracy of the inversion result.

[0128] 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 hierarchical constrained ionospheric tomography method based on hybrid interpolation and adaptive relaxation factor, characterized in that The method includes the following steps: Step 1: First, determine the inversion area, set virtual observation stations within the inversion area, determine the number of actual observation stations according to the cut-off range of the virtual observation stations within the inversion area and use them as interpolation stations, and then generate virtual STEC data of the virtual observation stations through a hybrid interpolation method, and integrate the actual STEC data and the virtual STEC data to obtain the observed STEC data; Step 2: Calculate the intercept of each ray within the voxel according to the coordinates of the actual observation stations and virtual observation stations within the inversion area, and then construct a coefficient matrix, where each row of the coefficient matrix represents a ray and each column corresponds to a voxel; Step 3: Obtain the initial electron density value through the IRI-2016 model, where the initial electron density value corresponds to the central point coordinates of each voxel, and then iteratively update the electron density value based on the multiplicative algebraic reconstruction algorithm using the coefficient matrix and the observed STEC data; Step 4: Apply horizontal constraints to the iterated electron density value using a smoothing constraint method, and apply vertical constraints to the electron density value according to the stratification strategy of the distribution characteristics of the electron density values at different heights; Step 5: Use the constrained electron density as the initial electron density value for the new round of iteration, and repeat steps (3)-(4) until the root mean square error is less than the preset threshold.

2. The method according to claim 1, wherein In step 1, the method for generating virtual STEC data is as follows: First, receive GNSS observation data through the actual observation stations within the inversion area and calculate the actual STEC data of each station; Then, take a series of virtual observation stations distributed in a grid within the inversion area, with a distribution density of 0.5°×0.5° in longitude and latitude. Select the actual observation stations as interpolation stations within the circle with each virtual observation station as the center and the cut-off distance as the radius, and then determine the number of interpolation stations; Finally, use different interpolation methods for virtual observation stations with different numbers of interpolation stations.

3. The method according to claim 1, wherein The specific hybrid interpolation method is as follows: If the number of actual observation stations within the cut-off range ≥ 6, take the 6 actual observation stations closest to the virtual observation station and use multi-quadric surface fitting to generate virtual STEC data; If the number of actual observation stations is 3-5, divide the sector and select the three stations closest to the virtual observation station, and ensure that the virtual station is inside the triangle formed by the three stations, and use triangular linear interpolation to generate virtual STEC data; If the number of actual observation stations is 2, use inverse distance weighted interpolation to generate virtual STEC data after dividing the sector; Virtual observation stations with the number of actual observation stations < 2 are discarded.

4. The method according to claim 1, wherein In step 2, the method for constructing the coefficient matrix is as follows: First, read the coordinates of all actual observation stations and the selected virtual observation stations within the inversion area simultaneously. In the space rectangular coordinate system, let the ground station coordinates be ( , , ), and the satellite coordinates be ( , , ). Then, construct the equation of the ray: ； In the formula, , , and are the variables of the equation; Then, calculate the intersection coordinates of the ray and each voxel; Finally, screen out the two intersection points P1 and P2 that fall on the surface of each voxel, calculate the intercept of the ray within each voxel according to the coordinates of the two intersection points, and then store the intercept value of each ray within each voxel in the coefficient matrix in the order of voxel coding.

5. The method according to claim 1, wherein In step 3, the open-source Fortran code of the IRI-2016 ionosphere model is called. By inputting the inversion time, the longitude and latitude of the inversion area, and the height range, a series of electron density values of coordinate points are obtained. At the same time, the grid spacing in the longitude and latitude directions and the height direction of the output data is determined, so that the output electron density values exactly correspond to the central point coordinates of each voxel, and these are used as the initial electron density values within each voxel.

6. The method according to claim 1, characterized in that The multiplicative algebraic reconstruction algorithm has the following expression: ； In the formula, represents the number of iterations, is the iterative value of the electron density of the th voxel at the th iteration, is the electron density value of the th voxel before the th iteration, represents the th observation ray, m×n is the th row of the intercept coefficient matrix A i.e., the intercept value of the th ray passing through each voxel. m is the total number of observation rays, including actual and virtual observation rays, and n is the total number of voxels. is the intercept distance of the (k-1) th observation ray passing through the th voxel. λ (k-1) is the relaxation factor at the kth iteration, where 0 < λ < 1.

7. The method according to claim 6, characterized in that The relaxation factor is adaptively adjusted by the residual. The specific method is as follows: ； In the formula, is the initial value, is a positive parameter for controlling the exponential decay rate, is the lower limit of the relaxation factor value, is the initial residual vector, is the residual vector at the k-th iteration.

8. The method according to claim 1, characterized in that, The smoothing constraint method is as follows: The voxels within the inversion area are divided into different horizontal layers according to different heights. Voxel smoothing is performed on each horizontal layer after each iteration. For the electron density value of any voxel in a certain horizontal layer, the smoothing method is as follows: ； In the formula, and respectively represent the electron density values of the th voxel before and after smoothing after the kth iteration correction; is the electron density value of the th voxel within the neighborhood range of the th voxel before smoothing; represents the number of iterations, is the total number of voxels, is the total number of voxels within the neighborhood range.

9. The method according to claim 8, wherein The smoothing constraint method is used to add vertical constraints to the voxels in the height range of 100 - 200 km.

10. The method according to claim 1, characterized in that, In step 4, the Chapman function is used to add vertical constraints to the voxels in the height range of 200 - 400 km: ； wherein, is the electron density value at the altitude , is the altitude at which the electron density is currently calculated, is the peak density of the ionospheric layer, is the peak altitude of the ionospheric layer, layer is the region in the ionosphere with an altitude range above 200 km, is the scale height at the bottom of the ionosphere, is the scale height at the top of the ionosphere, where the , , and are obtained by least squares fitting.

11. The method according to claim 1, wherein In step 4, the empirical orthogonal function is used to add vertical constraints to the voxels above 400 km: Based on the vertical electron density profile data set of the IRI-2016 model, the empirical orthogonal function is extracted by the singular value decomposition method; The first three empirical orthogonal functions are selected for vertical constraints.

12. The method according to claim 1, characterized in that In step 5, the reconstructed STEC data is calculated according to the electron density values after iteration. Iteration stops until the root mean square error between the reconstructed STEC data of each ray in the iteration inversion result and the observed STEC data is less than the limit difference σ.

Images