CN110244365A - A Multi-resolution Ionospheric Chromatography Method - Google Patents

Abstract

The invention provides a multi-resolution ionospheric tomography method, including gridding the inversion area and setting a time window function; preprocessing the GNSS observation data and calculating the TEC required for the first inversion; creating experience Orthogonal function; use the empirical orthogonal function to construct the mapping matrix, and regularize the matrix; use the processed matrix and the calculated TEC to construct the normal equation, find the optimal solution of the equation, and get an inversion result; set Set the value of level, create an empirical orthogonal function again, construct a mapping matrix, and regularize the matrix; use the results of the first inversion and the matrix after the second processing to construct a normal equation, find the optimal solution of the equation, and obtain The result of the second inversion. The invention performs secondary reconstruction on the ionospheric structure, greatly improving the accuracy of the ionospheric tomography; by setting the level, the ionospheric tomography results with different and smaller resolutions can be realized.

Description

A Multi-resolution Ionospheric Chromatography Method

technical field

The invention belongs to the technical field of ionospheric detection, in particular to a multi-resolution ionospheric chromatography method.

Background technique

As an all-weather, large-scale ionospheric detection technology, ionospheric tomography has many advantages such as low cost, simple operation, and wide detection range. However, due to the limitation of ground-based GNSS data coverage and the uneven distribution of ground stations, the tomographic inversion data source is lacking and unevenly distributed, which limits the reconstruction accuracy of the ionospheric structure, especially for the reconstruction of small and medium scale structures. severe constraints.

In recent years, with the increase of land-based GNSS ground stations and the construction of multi-constellation GNSS in some areas, it has provided a good opportunity for the development of ionospheric imaging algorithms. Zhao Haishan et al. proposed a new adaptive ionospheric tomography algorithm for the problem of low iteration accuracy of the traditional multiplicative algebraic reconstruction algorithm (MART). The comprehensive influence of the electron density can be reasonably assigned to the iterative difference; on the other hand, an adaptive relaxation factor related to the electron density value is proposed to effectively overcome the influence of the propagation noise on the electron density inversion. Aiming at the ionospheric electron density reconstruction problem, Tang Jun et al. proposed an ionospheric tomography algorithm that comprehensively uses the total variation minimization and multiplicative algebraic reconstruction algorithm. The algorithm imposes total variation constraints on the parameters of the inversion model to improve the stability of the inversion process and the accuracy of the results. Mingfeng et al. proposed an ionospheric tomography method with an additional scale factor, taking the ratio of the measured STEC to the STEC calculated by the current electron density in the iteration as the basic scale, based on the spatial correlation of the electron density and the distance from the grid center to the observed ray , the corresponding scale is calculated and adjusted for the grids that no signal passes through, so that it gradually approaches the true value and improves the inversion accuracy. The above methods have improved the accuracy of ionospheric tomographic inversion to a certain extent, but these tomographic methods are based on one inversion to reconstruct the ionospheric structure, resulting in the inversion accuracy of ionospheric tomography cannot reach a more ideal state .

Contents of the invention

In view of this, the present invention aims to propose a multi-resolution ionospheric tomography method to realize multi-resolution ionospheric structure tomography.

In order to achieve the above object, technical solution of the present invention is achieved in that way:

A multi-resolution ionospheric chromatography method, comprising the following steps:

(1) The ionosphere in the inversion area is gridded, and each grid is divided into three-dimensional volume pixels. The electron density integral on the signal path passing through the grid is TEC, and the inversion time function is set;

(2) Preprocessing the GNSS observation data, calculating the required TEC of step (1);

(3) create an empirical orthogonal function for the signal path passing through the grid in step (1);

(4) build compensation mapping matrix M to the empirical orthogonal function that step (3) obtains, and do regularization processing to matrix;

(5) Construct a normal equation for the matrix after the regularization process obtained in step (4) and the TEC calculated in step (2), and estimate the solution of the equation to obtain the optimal solution of the equation, which is the result of an inversion;

(6) Set the value of level, and perform steps (3) and (4) again;

(7) According to the primary inversion result obtained in step (5) as input data and the processed matrix obtained in step (6), the normal equation is constructed again, the solution of the equation is estimated, and the optimal solution of the secondary inversion is obtained.

Further, the step (1) specifically includes: using the setgrid function for the first inversion to set the grid to a rougher 4°×4°, and set the time interval to 30 minutes.

Further, the step (2) specifically includes: performing cycle slip and integer ambiguity detection on the observation data, and then calculating the TEC on the signal path passing through the ionospheric grid by using the carrier phase observation method.

Further, the step (3) specifically includes: using the Chapman function or the International Reference Model (IRI) to create an Empirical Orthogonal Function (EOF) from a group of radial distributions of ionospheric electron density.

Further, the step (4) specifically includes: in the inversion process, the electron density is assumed to be a linear combination of a given set of three-dimensional EOF, that is, the mapping matrix uses a set of EOF to map the radial change of the electron density, The longitude and latitude map the horizontal direction of the electron density, and the inverted time window maps the time axis of the electron density, and then perform secondary smoothing on the electron density obtained by using the mapping matrix in the horizontal direction, radial direction, and time respectively.

Further, the step (5) specifically includes: using the minimum residual method to introduce the minimum norm residual solution into the normal equation to obtain the optimal solution or the quadratic programming method to obtain the optimal solution of the equation.

Further, the step (6) specifically includes: setting level=[2,1], level represents the spatial and temporal resolution, and the second inversion is based on the first inversion with the spatial and temporal resolution divided by 2 , the ionospheric grid is 1°×1° in the second inversion, and the time is 15 minutes.

Further, the step (7) specifically includes: interpolating the primary inversion result into the inversion grid and time function set by level, then constructing a normal equation, and using the minimum residual method or quadratic programming method to solve The optimal solution of the normal equation this time.

Another object of the present invention is to propose a multi-resolution ionospheric chromatography device, and the specific scheme is realized in this way:

A multi-resolution ionospheric chromatography device comprising:

It is used to grid the ionosphere in the inversion area, divide each grid into three-dimensional volume pixels, and the integral of electron density on the signal path passing through the grid is TEC, and set the ionosphere of the inversion time function Grid module;

A TEC calculation module for preprocessing the GNSS observation data and calculating the TEC required by the ionospheric gridding module;

an empirical orthogonal function creation module for creating empirical orthogonal functions for signal paths through the grid in the ionospheric gridding module;

An empirical orthogonal function correction module for constructing a compensation mapping matrix M for the empirical orthogonal function obtained by the empirical orthogonal function creation module, and performing regularization processing on the matrix;

It is used to construct the normal equation from the regularized matrix obtained by the empirical orthogonal function correction module and the TEC calculated by the TEC calculation module, and estimate the solution of the equation to obtain the optimal solution of the equation, that is, the first inversion result of one inversion result module;

A data processing module for setting the value of level, and again using the empirical orthogonal function creation module and empirical orthogonal function correction module to perform empirical orthogonal function creation and empirical orthogonal function correction operations;

The secondary inversion is used to construct the normal equation again based on the primary inversion result obtained by the primary inversion result module as input data and the processed matrix obtained by the data processing module, estimate the solution of the equation, and obtain the optimal solution of the secondary inversion Results module.

Compared with the prior art, a kind of multi-resolution ionospheric chromatography method and system of the present invention has the following advantages:

(1) The present invention carries out secondary reconstruction to ionospheric structure, has improved the precision of ionospheric tomography to a large extent;

(2) The present invention can realize ionospheric tomography results with different and smaller resolutions by setting the level.

Description of drawings

The drawings constituting a part of the present invention are used to provide a further understanding of the present invention, and the schematic embodiments of the present invention and their descriptions are used to explain the present invention, and do not constitute an improper limitation of the present invention. In the attached picture:

Fig. 1 is the flowchart of the multi-resolution ionospheric tomography technique method provided by the present invention;

Figure 2 is a schematic diagram of the ionospheric grid structure;

Figure 3(a) is a comparison of no-multi, multi and ionos at 4°×4° on March 16, 2015 at MHJ45 site;

Figure 3(b) is a comparison of MHJ45 no-multi, multi and ionos at 1°×1° on March 16, 2015 at the MHJ45 site;

Figure 4(a) is a comparison of no-multi, multi and ionos at 4°×4° on March 16, 2015 at IF843 site;

Figure 4(b) is a comparison of no-multi, multi and ionos at the IF843 site on March 16, 2015 at 1°×1°.

Detailed ways

It should be noted that, in the case of no conflict, the embodiments of the present invention and the features in the embodiments can be combined with each other.

The present invention will be described in detail below with reference to the accompanying drawings and examples.

As shown in Figure 1, the multi-resolution ionospheric chromatography method provided by the present invention carries out the following steps in order:

(1) The ionosphere in the inversion area is gridded, and each grid is divided into three-dimensional volume pixels. The electron density integral on the signal path passing through the grid is TEC, and the inversion time function is set;

Use the setgrid function to grid the ionosphere, then set the transmission path between the satellite and the receiver as h, and the intercepts in the vertical direction when the ray passes through each voxel are recorded as H ₁ , H ₂ ,...,H _i . ..assuming that the electron density in each grid is a fixed value x ₁ , x ₂ ,..., _xi ... then the TEC value in the unit grid is H _i x _i , then the TEC value on a transmission path It can be considered as the sum of the product of the intercept of the ray passing through the grid and the electron density in the grid, that is, TEC=H ₁ x ₁ +H ₂ x ₂ +…+H _i x _i +…; then all the rays The upper electron content is Hx=z, and x is the desired ionospheric electron density. And the inversion time period of one day is set to a time interval of 30 minutes, so as to perform the first inversion of 4°×4°.

(2) Preprocessing the GNSS observation data, calculating the required TEC of step (1);

Specific steps are as follows:

A RINEX file for a GPS receiver contains the following observables.

where P1 and P2 are the _pseudoranges obtained from the refined code, and _L1 and _L2 represent the carrier phase _of the signal. P ₀ is the ionospheric free pseudorange, n and λ represent the integer ambiguity and carrier wavelength, respectively. ε represents receiver and satellite hardware error components.

Comparing the difference between the observables of phase and code distance, two equations for TEC can be obtained

In summary, the TEC value can be calculated by the noise term in Equation 2.5, or by the bias term of the integer ambiguity associated with it in Equation 2.6. When the satellite is within the visible range (except during the cycle slip period, such as a large jump or sudden change), the ambiguity of the entire cycle is constant, and the offset when a large cycle slip occurs is equal to the weighted average of the difference between formula 2.5 and formula 2.6 , and the weight is related to the signal-to-noise ratio. Compared with using pseudorange to calculate TEC value, the accuracy of TEC calculated by carrier phase is higher, so the carrier phase method is usually used to calculate the TEC on the signal path passing through the ionospheric grid.

(3) Solving the incomplete problem of the path passing through the grid in step (1), creating an empirical orthogonal function;

(a) The Chapman function creates an empirical orthogonal function (EOF):

The general equation of the electron density distribution Chapman function is:

Among them, χ is the solar zenith angle; N is the electron density; k is a constant with value, when representing the Chapmanβ layer, k=1, representing the Chapmanα layer, z is given by:

where h and h _m stand for height and peak density height, respectively; H is elevation. We usually use Chapman's β function to construct an orthogonal function, which generally requires the following parameters: the peak height within the range of 260km to 320km during the day, the peak height within the range of 260km to 360km at night, and the elevation of 60km to 90km. In addition, the ionospheric grid within a certain range of latitude and longitude, that is, the grid set by the setgrid function and level, needs to be input to obtain the three-dimensional space empirical orthogonal function set required by the electron density inversion reconstruction algorithm.

(b) The IRI model creates an Empirical Orthogonal Function (EOF):

The International Reference Ionosphere (IRI) is a standard specification for plasma parameters in the Earth's ionosphere. The IRI model synthesizes the observation data of many ionospheric sounding instruments such as ionospheric sounding sonar and scattered radar on a global scale, and describes the monthly average change of ionospheric density and temperature in the range of 50km to 1500km. To create empirical orthogonal functions with IRI model coefficients, the following parameters are generally required: the number of required empirical orthogonal functions, the date, and the specified ionospheric spatial grid. Specific steps are as follows:

1. Obtain the ionospheric profile time series at the center of the survey area and at the same time as the observation data through the IRI model;

2. Find the covariance matrix of this time series;

3. Decompose the eigenvalue of this covariance matrix to obtain the corresponding EOF.

(4) Construct the compensation mapping matrix M for the empirical orthogonal function obtained in step (2), and perform regularization processing on the matrix.

During the inversion, the electron density is assumed to be a linear combination of empirical orthogonal functions (EOF), i.e. the mapping matrix uses the empirical orthogonal function (EOF) calculated in step (3) to map the radial variation of electron density, longitude and The latitude maps the horizontal direction of the electron density, and the inverted time window maps the time axis of the electron density, wherein the latitude-longitude interval and the time window are respectively set according to step (1) to form a mapping matrix M, that is, the equation Hx=z becomes (HM)X = z, and then introduce the regularization matrix R to perform quadratic smoothing on the solution of the equation in the horizontal, radial and time directions.

(5) Construct a normal equation for the processed matrix obtained in step (4) and the TEC calculated in step (2), and estimate the solution of the equation to obtain the optimal solution of the equation, which is the result of an inversion.

Specific steps:

(a) Construct the normal equation H ^T H+cR ^T R=H ^T z, where c is a user-defined regularization constant

(b) Calculate the solution of the equation using the minimum residual method

The minimum residual method (minres) tries to introduce the minimum norm residual solution x into the linear equation system, for example: Ax=b, in the minimum residual method, the n-order coefficient matrix A must be symmetric but not necessarily positive definite, and the right side of the equation The rows or columns of vector b must be of length n.

(c) Calculate the solution of the equation using the quadratic programming method (quadprog)

The quadratic programming method (quadprog) solves the problem in the form:

Among them, HX≤b

In this way, the vector X that minimizes the above problem can be obtained. In order to make this problem have a finite minimum solution, it is necessary to ensure that H is a positive definite matrix. The Quadprog algorithm requires two or more inputs, such as upper and lower bounds, equality and inequality constraints, etc., to impose constraints on the solution of the inversion problem. In order to make the solution of the inverse problem in the transformation basis non-negative, the quadprog algorithm imposes constraints on the inverse problem through the inequality MX≥0, which ensures that the vector of EOF coefficients is a set of positive numbers. Then, the final solution is based on the original basis, and the finally obtained electron density profile distribution is: x=MX

Choose any one of these two methods to find the optimal solution of the equation, which is the result of an inversion.

(6) Set the value of level, and perform steps (3) and (4) again;

Set level=[2,1], that is, the ionospheric grid in the inversion area becomes 4°×4°/2^2, the inversion time interval is 30min/2^1, and the empirical orthogonal function is created and the mapping matrix is constructed , and the matrix regularization process is the same as steps (3) and (4).

(7) According to the primary inversion result obtained in step (5) as input data and the processed regular matrix obtained in step (6), the normal equation is constructed again, the solution of the equation is estimated, and the optimal solution of the secondary inversion is obtained.

Use the method of step (5) again to construct the normal equation, and then solve the optimal solution of the equation according to the minimum residual method or the quadratic programming method to obtain the result of the quadratic inversion.

The effect of the multi-resolution chromatography method provided by the present invention can be further illustrated by the following experimental results, and the results calculated by the method of the present invention and the altimeter data are compared and analyzed.

Experimental data description: The GPS observation data is provided by UNAVCO, the precise ephemeris data is downloaded from the IGS official website, and the altimeter data is from the Digital Ionogram Database official website. One of the altimeters is located in the Idaho National Lab in the west of the United States (code: IF843, location: northern latitude 43.81°, west longitude 112.68°), its local time LT = universal time UT-7; another Millstone Hill (code: MHJ45, location: north latitude 42.6°, west longitude 71.5°) in the eastern United States, its local time LT = UT-5.

Table one has provided the inventive method and the single resolution method and the altimeter data comparison result, from the statistical result, the inventive method calculation result is higher precision with respect to the altimeter data; In addition the layer of the inventive method Higher resolution imaging.

Table 1 Comparison of multi, no-multi chromatography results and the NmF2 error of the altimeter

The above descriptions are only preferred embodiments of the present invention, and are not intended to limit the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present invention shall be included in the scope of the present invention. within the scope of protection.

Claims (9)

1. a kind of Ionospheric Tomography method of multiresolution, it is characterised in that: the following steps are included:

(1) by the ionosphere grid in inverting region, each grid is divided into three-dimensional volume pixel, across the signal of grid Electronics density integral is TEC on path, and the function of time of inverting is arranged；

(2) data prediction is observed to GNSS, TEC needed for calculating step (1)；

(3) Empirical Orthogonal Function is created to the signal path for passing through grid in step (1)；

(4) compensation map matrix M is constructed to the Empirical Orthogonal Function that step (3) obtains, and Regularization is done to matrix；

(5) TEC that matrix and step (2) after obtaining Regularization to step (4) are calculated constructs normal equation, and estimates Non trivial solution is counted, equation optimal solution, i.e. an inversion result are obtained；

(6) value of level is set, and carries out the operation of step (3), (4) again；

(7) inversion result obtained according to step (5) as input data and step (6) matrix that obtains that treated again Normal equation is constructed, non trivial solution is estimated, obtains the optimal solution of quadratic inversion.

2. a kind of Ionospheric Tomography method of multiresolution according to claim 1, it is characterised in that: the step (1) Specifically include: first time inverting sets more rough 4 ° × 4 ° for grid with setgrid function, and sets time interval to 30min。

3. a kind of Ionospheric Tomography method of multiresolution according to claim 1, it is characterised in that: the step (2) It specifically includes: cycle slip being carried out to observation data and integer ambiguity detects, is then calculated through electricity using carrier phase observation method TEC on the signal path of absciss layer grid.

4. a kind of Ionospheric Tomography method of multiresolution according to claim 1, it is characterised in that: the step (3) It specifically includes: being created using Chapman function or International Reference model (IRI) from the radial distribution of one group of ionospheric electron density Build Empirical Orthogonal Function (EOF).

5. a kind of Ionospheric Tomography method of multiresolution according to claim 4, it is characterised in that: the step (4) Specifically include: in refutation process, electron density is assumed the linear combination of one group of given three-dimensional EOF, i.e. mapping matrix Map the radial variations of electron density using one group of EOF, the horizontal direction of longitude and latitude mapping electron density, inverting when Between window mapping electron density time shaft, then respectively in the horizontal direction, radial direction, on the time to the electricity found out using mapping matrix Sub- density carries out secondary smoothing processing.

6. a kind of Ionospheric Tomography method of multiresolution according to claim 1, it is characterised in that: the step (5) It specifically includes: minimum norm residual error solution being introduced into normal equation using minimum residual method and seeks optimal solution or quadratic programming Seek equation optimal solution.

7. a kind of Ionospheric Tomography method of multiresolution according to claim 1, it is characterised in that: the step (6) Specifically include: setting level=[2,1], level represent room and time resolution ratio, and quadratic inversion is in first time inverting base Room and time resolution ratio is reduced again with 2 power on plinth, i.e., ionosphere grid is 1 ° × 1 ° in second of inverting, the time For 15min.

8. a kind of Ionospheric Tomography method of multiresolution according to claim 1, it is characterised in that: the step (7) It specifically includes: will be inserted into an inversion result in inverting grid and the function of time by level setting, then construct regular side Journey solves the optimal solution of this normal equation using minimum residual method or quadratic programming.

9. according to a kind of Ionospheric Tomography device of multiresolution of the claims, it is characterised in that: include:

For each grid to be divided into three-dimensional volume pixel, across the letter of grid for the ionosphere grid in inverting region Electronics density integral is TEC on number path, and the ionospheric grid module of the function of time of inverting is arranged；

For observing data prediction to GNSS, the TEC computing module of TEC needed for calculating ionospheric grid module；

For the Empirical Orthogonal Function to the signal path creation Empirical Orthogonal Function for passing through grid in ionospheric grid module Creation module；

Empirical Orthogonal Function for obtaining to Empirical Orthogonal Function creation module constructs compensation map matrix M, and does to matrix The Empirical Orthogonal Function correction module of Regularization；

For obtaining the matrix after Regularization to Empirical Orthogonal Function correction module and TEC computing module is calculated TEC constructs normal equation, and estimates non trivial solution, obtains equation optimal solution, i.e., an inversion result mould of one time inversion result Block；

For setting the value of level, and progress Empirical Orthogonal Function creation module, Empirical Orthogonal Function correction module are utilized again Carry out the data processing module of Empirical Orthogonal Function creation, Empirical Orthogonal Function amendment operation；

An inversion result for being obtained according to an inversion result module is obtained as input data and data processing module Treated, and matrix again pulls up normal equation, estimates non trivial solution, obtains the quadratic inversion result of the optimal solution of quadratic inversion Module.