CN111273335A - Ionosphere tomography method based on vertical measurement data constraint - Google Patents

Abstract

The invention discloses an ionosphere tomography method based on vertical measurement data constraint, which comprises the following steps: step A, reconstructing a background ionosphere model F2 layer characteristic parameter based on vertical measurement data: and B, calculating the total inclined electron content of the ground-based GNSS ionosphere: step C, constructing an ionospheric tomography inversion matrix: the invention discloses an ionospheric tomography method based on vertical measurement data constraint, which adopts ionospheric vertical measurement data issued by a global radio ionospheric observation station and ground GNSS observation data issued by an International global navigation satellite system (IGS) as a data source of ionospheric tomography, updates an ionospheric tomography background model by using the vertical measurement data, corrects the background ionospheric model by using an improved multiplicative algebra iterative algorithm in combination with inclined total electron content data of the ground GNSS, and realizes regional and even global ionospheric tomography so as to acquire space-time variation characteristics of total electron content and electron density of the ionospheric.

Description

Ionosphere tomography method based on vertical measurement data constraint

Technical Field

The invention belongs to the field of space ionosphere environment monitoring, and particularly relates to an ionosphere tomography method based on vertical measurement data constraint in the field.

Background

The ionosphere, which is an important component of the space-day environment, has effects of refraction, reflection, scattering, absorption and the like on radio waves passing through the ionosphere, so that the performance of various radio information systems such as satellite navigation, communication, radar and the like is affected. The method utilizes various technical means to obtain the characteristic parameters of the ionized layer and research various phenomena in the ionized layer, thereby revealing the physical mechanism behind the ionized layer, and having important scientific, even political, military and economic significance. Among many characteristic parameters of the ionosphere, the ionosphere electron density is the most critical detection parameter, and the mastering of the time-space variation characteristics and rules thereof has important practical value for improving the performance of a radio information system, so that the accurate acquisition of the ionosphere electron density becomes one of the key directions of concern in the ionosphere research field of various countries in the world gradually, and the very important means for acquiring the ionosphere electron density is the ionosphere tomography technology.

With the deep ionosphere model research and the rapid improvement of the computer hardware performance, especially with the development of the satellite ionosphere detection technology, the ionosphere Tomography (Computerized ionosphere Tomography) technology has made great progress in the aspects of theoretical methods and applications in view of the current forecast requirements of different service guarantees on ionosphere parameters. The most typical technical result in the direction of ionospheric tomography is a multi-hand ionospheric tomography system MIDAS (Multi Instrument Data Analysis System) developed by Bath university in UK, the system introduces the latest international ionospheric tomography result, has very high operation stability and inversion accuracy, and can provide information of ionospheric TEC and electron density distribution in a large range of areas.

The major data sources for ionospheric tomography today are ground-based LEO beacon receivers and ground-based GNSS receivers. The ground-based LEO beacon receiver is mainly applied to two-dimensional ionosphere tomography, and the GNSS data is mainly applied to the field of three-dimensional ionosphere tomography. Three-dimensional ionospheric tomography is significantly different from two-dimensional ionospheric tomography in terms of measurement view angle, spatial coverage of observation data, angular velocity of CT scanning, error source of measurement TEC, number of unknowns to be solved, and the like, and three-dimensional ionospheric tomography is not a simple extension of two-dimensional ionospheric tomography in spatial dimension, and therefore three-dimensional ionospheric tomography needs to adopt a solution strategy different from that of two-dimensional ionospheric tomography. Particularly, due to the fact that the orbit height of the GNSS satellite is high, for a certain specific area, ray paths corresponding to most observation data are in a near-vertical direction compared with a CT imaging grid, and according to the ionospheric tomography theory, the vertical resolution which can be achieved by the ground-based GNSS imaging is low, that is, the ionospheric peak value accuracy obtained by inversion is poor.

Disclosure of Invention

The invention aims to provide an ionospheric tomography method based on vertical measurement data constraint.

The invention adopts the following technical scheme:

in an ionospheric tomography method based on vertical survey data constraints, the improvement comprising the steps of:

step A, reconstructing a background ionosphere model F2 layer characteristic parameter based on vertical measurement data:

step A1: downloading vertical measurement automatic interpretation data issued by a GIRO (global radio ionospheric observation station) according to the appointed date, and extracting the geographical latitude of the vertical measurement station

Longitude lambda, ionosphere F2 layer critical frequency

Peak height of ionosphere F2 layer

Wherein i is the position number of the vertical measuring station, t,

λ_iRespectively corresponding to the ith vertical measuring station at the time t, latitude and longitude;

step A2: selecting international reference ionosphere IRI as a background ionosphere model, and calculating the ionosphere F2 layer critical frequency of the IRI model at the same time and position

And peak height

Calculating the difference between all the observed values and the model values, wherein the calculation method comprises the following steps:

△Z_foF2＝foF2_ion-foF2_IRI

△Z_hmF2＝hmF2_ion-hmF2_IRI

wherein △ Z_foF2、△Z_hmF2Represents the difference between ionosphere foF2 and hmF2 measured by the plumb bob, respectively, relative to the IRI model calculations;

step A3: calculating the equivalent distance of the ionized layer by the following method:

in the formula of_iAnd

respectively, as longitude and latitude of a vertical measuring point i, in degrees, lambda_jAnd

respectively longitude and latitude of a vertical measuring point j, wherein F is a distance weighting factor, and can take a value of 3.5 for low-latitude and high-latitude areas and a value of 2.0 for medium-latitude areas;

step A4: for any time t, calculating the difference value between ionosphere foF2 and ionosphere hmF2 observed vertically at the time and the calculated value of the IRI model according to the step A2, and setting the difference value at the ith sample point as

And

calculating an arbitrary position of an imaging region

Difference of ionosphere foF2 and hmF2 with respect to IRI model

The calculation method is as follows:

wherein: w is a_foF2,i、w_hmF2,iWeighting coefficients representing ionosphere foF2 and hmF2, respectively, with N representing the total number of samples observed;

step A5: solving the weighting coefficient according to the kriging unbiased estimation and the optimal estimation theory, and jointly solving the following (N +1) equation sets, wherein the equation set expression method is as follows:

wherein: d_ijThe ionospheric equivalent distance between the ith sample point and the jth sample point is represented, and lambda is a Lagrange multiplier;

step A6: converting the joint solution equation into a linear equation, wherein a matrix solution expression of the linear equation is as follows:

Aw＝b

in the formula: a is a matrix of (N +1) × (N +1), w is a vector of (N +1) × 1, is an unknown coefficient to be solved, and comprises N weighting coefficients and a vector of lambda, and b is (N +1) × 1;

step A7: using least squares estimation, arbitrary positions can be calculated

The weighting coefficient is calculated as follows:

w＝(A^TA)^-1A^Tb

step A8: according to the calculation method in the step A4, the solved weighting coefficients are substituted into the formula again to reconstruct the increment value

And

foF for the region is obtained by adding the output of IRI model₂And hmF₂Distribution:

step A9: traversing all ionospheric grids in the region to be imaged, and respectively reconstructing ionospheric layers foF2 and hmF2 of all grids in the region according to the steps A4-A8;

and B, calculating the total inclined electron content of the ground-based GNSS ionosphere:

step B1: according to the same date as the vertical observation data, acquiring a GNSS observation file and a GNSS satellite ephemeris file issued by an IGS (International Global navigation System) organization, decompressing the data, extracting the longitude, the latitude and the height of an observation station, and recording the longitude, the latitude and the height as

Extracting pseudo-range between GNSS satellite and observation station

And carrier phase

Wherein i is the corresponding observation frequency band, i is 1,2, j is the corresponding satellite number, and r is the corresponding receiver number;

step B2: respectively calculating a non-geometric distance linear combination for all pseudo-ranges and carrier phase observed values, wherein the calculation method comprises the following steps:

wherein:

is a geometry-free linear combination of pseudoranges,

is a geometric distance-free linear combination of carrier phases;

step B3: in the same observation arc segment, smoothing filtering is carried out on the non-geometric distance linear combination P4 of the pseudo-range by utilizing the non-geometric distance linear combination L4 of the carrier phase, and the calculation method is expressed as follows:

wherein: t represents epoch time, N_arcRepresenting the total number of observation epochs in an overhead arc section;

step B4: calculating an elevation angle E and an azimuth angle A of the receiver relative to the GNSS satellite according to the coordinate of the receiver and the GNSS satellite coordinate obtained by satellite ephemeris resolving;

step B5: assuming that total concentration TEC of ionized layer electrons is concentrated in a thin layer at the height of 450km above the earth, calculating the geographic coordinates of puncture points of the observation station and the GNSS satellite which are connected on the height of the thin layer

The calculation formula is as follows:

wherein:

latitude and longitude coordinates of the GNSS observation station are shown, and A is an azimuth angle; psi₀For the geocentric angle between the observation station and the GNSS satellite, the calculation formula is as follows:

wherein: r_e6371.2km is the average radius of the earth; h is₀Is the station height; h is_I450km is the height of an ionosphere thin layer, and E is the elevation angle between a receiver and a satellite;

step B6: and calculating a mapping function value of the total inclined electron content and the total vertical electron content by the following calculation method:

M(E)＝[1-(R_ecosE/(R_e+h_I))²]^-1/2

step B7: a model of the vertical total electron content of the regional ionized layer along with the spatial distribution is constructed by utilizing a quadratic Taylor expansion series, and the calculation method is as follows:

wherein: n is_maxThe maximum order of the expansion is represented,

and

respectively the geographic latitude and the geographic longitude of the puncture point between the r-th receiver and the j-th satellite,

and λ_sSet as the center point geographical latitude and geographical longitude, E, respectively, of the tomographic imaging area_nmIs an unknown coefficient to be solved;

step B8: geometry-free linear combination using smoothed pseudoranges

Constructing a total electron content mapping solution equation, wherein the calculation method comprises the following steps:

wherein: f. of₁Is the first frequency of the satellite, f₂Being the second frequency of the satellite, DCB^rFor the hardware delay of the r-th receiver, DCB_jFor the hardware delay value of the jth satellite, in order to separate the hardware delays of the satellite and the receiver, the following four constraints are added:

wherein:

hardware delay value, N, representing the jth GPS satellite_GThe total number of GPS satellites;

denotes the hardware delay value, N, of the jth GLONASS satellite_RThe total number of GLONASS satellites;

hardware delay value N representing the jth Beidou satellite_CThe total number of the Beidou satellites;

hardware delay value, N, for the jth GALILEO satellite_EThe total number of GALILEO satellites;

step B9: and (3) selecting 1-2 hours as a time window, traversing all satellites and receivers in the time window, and establishing a linear matrix of the vertical total electron content according to the equation in the step B8 as follows:

Y＝AX

wherein: y is the observed quantity between all receivers and satellite

A is a coefficient matrix related to spherical harmonics, receivers and satellite hardware delays, and X is an unknown coefficient vector to be solved, and the specific expression mode is as follows:

step B10: solving a fitting equation by using a least square algorithm, and solving to obtain an unknown coefficient X, wherein the calculation method comprises the following steps:

X_LS＝(A^TA)^-1A^TY

wherein: x_LSIs a least squares solution comprising Taylor series expansion coefficients E_nmReceiver and satellite hardware delays;

step B11: and (3) substituting the resolved satellite and receiver hardware delay DCB into the Taylor expansion series again to resolve and obtain the ionized layer tilt TEC value, wherein the calculation method comprises the following steps:

wherein: f. of₁Is the first frequency of the satellite, f₂Is the second frequency of the satellite or satellites,

calculating to obtain the total electron content of the inclined ionized layer between the r receiver and the j satellite;

step C, constructing an ionospheric tomography inversion matrix:

step C1: defining the geographic coordinate range and grid size of ionospheric tomography region, and defining the latitude range of

Step size of latitude grid is

Longitude range of [ lambda ]_min,λ_max]The longitude grid step is d lambda; height range of [ h_min,h_max]The height grid step is dh, and the calculation method of the ith longitude, latitude and height grid point is as follows:

λ_grid,i＝λ_min+(i-1)dλ，i＝1,2...,N_λ，N_λ＝(λ_max-λ_min)/dλ

h_grid,i＝h_min+(k-1)dh，k＝1,2...,N_h，N_h＝(h_max-h_min)/dh

wherein: n is a radical of_λ，

N_hThe numbers of grid points respectively representing the longitudinal, latitudinal and elevational directions of the tomographic region;

step C2: GNSS receiver

And GNSS satellite longitude and latitude high coordinate

Converted into spatial rectangular XYZ coordinates, respectively denoted (X)_r,Y_r,Z_r) And (X)_s,Y_s,Z_s) The conversion expression is expressed as:

wherein

R_eRepresenting the radius of the earth, e²＝0.00669437999013；

Step C3: and establishing a ray equation between the receiver and the satellite according to the relative position between the receiver and the satellite, and calculating the expression as follows:

wherein (X, Y, Z) represents the rectangular coordinate of any point on the ray, K represents a proportionality constant, and any intersection point coordinate of the ray path can be determined as long as the value of K is determined;

step C4: according to the relative position between the receiver and the GNSS satellite, the longitude plane L, the latitude plane B and the altitude plane H of the intersection of the ray and the tomography area are determined, and the determination method is as follows:

wherein: the longitude plane L, the latitude plane B and the altitude plane H must satisfy

B∈[λ_min,λ_max]、 H∈[h_min,h_max]；

Step C5: calculating the intersection (X) of a ray with the longitudinal, latitudinal and elevational planes of a tomographic grid_i,Y_i,Z_i) Wherein i is the intersection number;

step C6: the coordinates of the cross point are expressed by rectangular coordinates (X)_i,Y_i,Z_i) Conversion to a geographical coordinate system

The calculation formula is as follows:

in the formula: e' is the second eccentricity of the reference ellipsoid;

a＝6378.137，b＝6356.752；

step C7: according to the geographic coordinates of the intersection

Calculating the midpoints between the intersections in order

The calculation method is as follows:

step C8: determining the index value of the lattice according to the midpoint

The determination method comprises the following steps:

IDλ＝count(λ_o≥λ_grid)

ID_h＝count(h_o≥h_grid)

wherein: lambda [ alpha ]_gridRepresenting longitude grid points;

representing all latitude grid points; h is_gridRepresents height grid points, count (-) represents statistical counts;

step C9: and calculating the grid number of the midpoint of the cross point, wherein the number calculation method in the chromatography inversion matrix comprises the following steps:

wherein: n is a radical of_λIndicates the number of grids divided in the longitudinal direction,

a grid number indicating a division in a latitudinal direction;

step C10: calculating the intercept of the ray intersection point in the grid through a distance formula between two points, and filling the intercept into the chromatography inversion matrix A, wherein the calculation method comprises the following steps:

wherein: p represents the p-th ray used for tomography; subscripts i and i +1 respectively represent a starting point and an end point of the rays in the same grid, and all the rays are traversed, so that an inversion matrix of ionospheric tomography can be obtained;

step D, time-varying three-dimensional ionosphere electron density inversion:

step D1: acquiring the total inclined electron content data of all GNSS receivers, and constructing a solving equation of ionosphere electron density inversion according to the geometric position relation between the receivers and satellites:

wherein: STEC denotes the ionospheric tilt total electron content, N, extracted by the GNSS receiver_eThe electron density value of the ionized layer is represented, and the upper and lower integral limits are the positions of the satellite and the receiver respectively;

step D2: partitioning inversion regions into

The grid adopts a series expansion algorithm and utilizes a group of space basis functions h_k(r) discretizing the electron density distribution by the following calculation method:

wherein: h is_k(r) is a basis function, a_kFor the weighting coefficients, the basis functions select the pixel basis functions, the calculation method is as follows:

step D3: selecting a particular time window, at which time a can be identified_kThe ionosphere electron density tomography equation is constructed without changing along with time, and the calculation method is as follows:

R_inrepresenting the intercept length of the ith ray within the nth grid point;

step D4: traversing ionized layer TEC data of all GNSS, and converting an ionized layer tomography equation into a problem of solving a linear equation set, wherein the calculation method comprises the following steps:

d＝Ax

wherein: d is a vector formed by ionized layer inclination TEC data, a matrix A is a chromatography inversion matrix given in the step C, and x is the electron density of each grid point to be solved;

step D5: obtaining the ionized layers foF2 and hmF2 of all grids in the region reconstructed in the step A, inputting the ionized layers foF2 and hmF2 into the IRI model again, calculating to obtain the electron density value of the ionized layer, and taking the electron density value as the initial value X of the background field of the tomography⁽⁰⁾The initial value satisfies the following condition:

wherein: n represents the total number of grids of the region to be tomography;

step D6: iterative solution is carried out on the tomography equation by adopting an improved multiplication algebra reconstruction algorithm, and the calculation formula is as follows:

wherein: the iteration sequence of the equation is generated by random numbers in the range of 1, N](ii) a m is the total number of the ionosphere inclined total electron content observed values, and i is the inclined TEC ray number; j is a grid number;

and

(ii) electron density of the jth grid for iterations k +1 and k, respectively; a is_ijIs the element in the inversion matrix A; | | a_iThe | | is the total length of the ith ray path; d_iIs the total electron content of the ith ray path; λ is pineA relaxation factor;

step D7: take the larger relaxation factor to iterate, λ_kThe value is 3.0, and the cut-off threshold value in the iterative process meets tol<10^-8And then stopping iteration, wherein the threshold value calculation method comprises the following steps:

wherein: tol represents an iteration cut-off threshold;

step D8: replacing X with the final X value obtained in the step D6⁽⁰⁾As an initial value of the iteration, a smaller relaxation factor lambda can be used_kAnd (5) iterating again for a second iteration when the electron density of the ionized layer is 0.9, and taking the resolved electron density of the ionized layer as a final ionized layer tomography data product.

Further, the step B4 specifically includes:

step B401: longitude and latitude high coordinates of receiver and GNSS satellite

Converting into space rectangular XYZ coordinates, and expressing the conversion expression as:

wherein

R_eWhich is the radius of the earth, is,

step B402: converting the space rectangular coordinate into a station center rectangular coordinate (N, E, H), wherein the conversion expression is as follows:

wherein: t is a rotation matrix, (X)₀,Y₀,Z₀) Is the rectangular coordinate of the receiver, (X)₁,Y₁,Z₁) For the GNSS satellite rectangular coordinate, the calculation method of the rotation matrix T is as follows:

wherein:

longitude and latitude coordinates of the GNSS receiver;

step B403: calculating the elevation angle E between the receiver and the satellite, wherein the calculation expression is as follows:

step B404: calculating the azimuth angle A between the receiver and the satellite, wherein the calculation expression is as follows:

further, the step C5 specifically includes:

step C501: calculating the intersection point of the ray and the longitude plane, wherein the tomography forms a curved surface along the longitude direction, namely a plane passing through the Z axis, and the included angle between the curved surface and the initial meridian plane is longitude L, and the equation of the curved surface is expressed as follows:

Y＝X tan(L)

step C502: and C501, solving the curved surface equation in the step C501 and the ray equation in the step C3 in a simultaneous manner to obtain a proportionality constant K corresponding to the intersection point of the ray and the longitude plane, wherein the calculation method is as follows:

step C503: the proportionality constant K is substituted into the ray equation of the step C3 again to obtain the coordinates of the intersection point of the ray and the longitude plane L, the coordinates are added into the intersection point sequence of the ray and the grid plane,

step C504: calculating the intersection point of the ray and the latitude plane, wherein the curved surface formed along the latitude direction is a conical surface with the vertex at the origin and the rotating axis as the z axis, the half angle of the cone of the conical surface is complementary with the latitude B, and the equation of the curved surface is expressed as follows:

(X²+Y²)·tan²W-Z²＝0

step C505: and C504, solving the curved surface equation in the step C504 and the ray equation in the step C3 in a simultaneous manner to obtain a proportionality constant K corresponding to the intersection point of the ray and the longitude plane, wherein the calculation method is as follows:

A＝(△X tanB)²+(△Y tanB)²-(△Z)²

B＝2×((△X·X_r·tan B)²+(△Y·Y_r·tan B)²-△Z·Z_r)

C＝(X_r·tan B)²+(Y_r·tan B)²-(Z_r)²

wherein: A. b, C are root coefficients of a quadratic equation;

step C506: substituting the proportionality constant K into the ray equation in the step C3 again to obtain the coordinates of the intersection point of the ray and the latitude plane B, judging whether the intersection point is between the starting point and the ending point, if so, adding the intersection point into the effective point sequence of the ray, otherwise, rejecting the point;

step C507: calculating the intersection point of the ray and the height surface H, wherein in the region to be chromatographically imaged, a curved surface formed by the height surface and the same height point from the ground is a spherical surface with the center of the sphere at the origin, and the equation of the curved surface is as follows:

(X²+Y²+Z²)-(R_e+H)²＝0

in the formula: r_eIs the average earth radius;

step C508: and C507, solving the curved surface equation in the step C507 and the ray equation in the step C3 in a simultaneous manner to obtain a proportionality constant K corresponding to the intersection point of the ray and the longitude plane, wherein the calculation method is as follows:

α＝(△X)²+(△Y)²+(△Z)²

β＝2×(△X·X_r+△Y·Y+△Z·Z_r)

γ＝(X_r)²+(Y_r)²+(Z_r)²-(R_e+H)²

wherein α, β and gamma are root coefficients of a quadratic equation of one unit;

step C509: substituting the proportionality constant K into the ray equation of the step C3 again to obtain the coordinates of the intersection point of the ray and the height surface H, judging whether the intersection point is between the starting point and the ending point, if so, adding the intersection point into the effective point sequence of the ray, otherwise, rejecting the point;

step C510, classifying all the intersection points into a sequence, and calculating the distance △ R from all the intersection points to the receiver, wherein the calculation formula is as follows:

step C511, the distances △ R from the effective point sequences of all the cross points to the receiver are sorted according to a method from small to large, and the coordinates of the corresponding points are respectively put into a coordinate sequence (X)_i,Y_i,Z_i) 1, 2.., n.

The invention has the beneficial effects that:

the invention discloses an Ionospheric tomography method based on vertical measurement data constraint, which adopts Ionospheric vertical measurement data issued by a Global Radio Ionospheric observation station (GIRO) and ground-based GNSS observation data issued by an International Global navigation satellite system (IGS) organization as a data source of Ionospheric tomography, updates an Ionospheric tomography background model by using the vertical measurement data, corrects the background Ionospheric model by using an improved multiplicative algebra iterative algorithm in combination with inclined total electron content data of the ground-based GNSS, and realizes regional and Global Ionospheric tomography so as to acquire the time-space variation characteristics of the total electron content and the electron density of the Ionospheric. Compared with the traditional ionosphere tomography algorithm, the method can effectively improve the vertical resolution of ionosphere tomography, thereby meeting the requirement of accurate monitoring of regional and even global high-resolution three-dimensional ionosphere electron density, and related data products can serve in multiple fields of space science research, ionosphere space weather forecast, ionosphere refraction error correction of a radio information system and the like.

Drawings

FIG. 1 is a schematic flow chart of the method disclosed in example 1 of the present invention;

fig. 2 is a schematic diagram of a reconstruction result of a critical frequency of a layer F2 of a background ionosphere model based on vertical measurement data in the method disclosed in embodiment 1 of the present invention;

FIG. 3 is a diagram illustrating the total electron content of the ionospheric tilt of the ground-based GNSS calculated by the method disclosed in embodiment 1 of the present invention;

FIG. 4 is a schematic diagram of a structure of an ionospheric tomography inversion matrix in the method disclosed in embodiment 1 of the present invention;

fig. 5 is a schematic diagram of a time-varying three-dimensional ionospheric electron density inversion result in the method disclosed in embodiment 1 of the present invention.

Detailed Description

In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention will be described in further detail below with reference to the accompanying drawings and examples. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.

In order to improve the vertical resolution of ionospheric three-dimensional CT imaging, an important method is to add other observation data with higher vertical resolution for imaging constraint, for example, to add ionospheric vertical measurement data to improve the vertical resolution of ionospheric tomography. With the construction of global GPS, GLONASS, Beidou, GALILEO and other navigation systems, measurement data based on the global GNSS system is increasingly abundant, and GNSS observation data of more than 300 stations can be acquired online at present; meanwhile, global foundation vertical measurement stations are continuously increased, and acquisition of ionosphere vertical measurement data is more convenient. For example, the global ionospheric vertical survey network, initiated and established at the university of Lowell, massachusetts, usa, has incorporated data for 46 altimeter stations worldwide, covering 27 countries and regions worldwide, including china.

The embodiment provides an ionospheric tomography method based on vertical measurement data constraint aiming at the defect of ionospheric tomography by singly using a ground GNSS, solves the key problems of reconstruction of characteristic parameters of a background ionospheric model F2 layer based on vertical measurement data, calculation of total electron content of inclination of the ionospheric GNSS ionosphere of a ground GNSS, construction of an ionospheric tomography inversion matrix, inversion of electron density of a time-varying three-dimensional ionosphere and the like, and provides a whole set of solution for high-precision and high-resolution ionospheric tomography measurement.

Specifically, embodiment 1, as shown in fig. 1, discloses an ionospheric tomography method based on vertical measurement data constraint, which includes the following steps:

step A, reconstructing a background ionosphere model F2 layer characteristic parameter based on vertical measurement data:

step A1: downloading automatic interpretation data of vertical survey (Ionsonde) issued by Global ionosphere observation station (Global Ionospheresis observer, GIRO) according to appointed date, and extracting geographic latitude of the vertical survey station

Longitude lambda, ionosphere F2 layer critical frequency

Peak height of ionosphere F2 layer

Wherein i is the vertical measurementStation position number, t,

λ_iRespectively corresponding to the ith vertical measuring station at the time t, latitude and longitude;

step A2: selecting an International Reference Ionosphere (IRI) as a background Ionosphere model, and calculating the Ionosphere F2 layer critical frequency of the IRI model at the same time and position

And peak height

Calculating the difference between all the observed values and the model values, wherein the calculation method comprises the following steps:

△Z_foF2＝foF2_ion-foF2_IRI

△Z_hmF2＝hmF2_ion-hmF2_IRI

wherein △ Z_foF2、△Z_hmF2Represents the difference between ionosphere foF2 and hmF2 measured by the plumb bob, respectively, relative to the IRI model calculations;

step A3: calculating the equivalent distance of the ionized layer by the following method:

in the formula of_iAnd

respectively, as the longitude and latitude (in degrees), λ, of the vertical point i_jAnd

respectively longitude and latitude of a vertical measuring point j, wherein F is a distance weighting factor, and can take a value of 3.5 for low-latitude and high-latitude areas and a value of 2.0 for medium-latitude areas;

step A4: for theAt any time t, calculating the difference value between ionosphere foF2 and ionosphere hmF2 observed vertically at the time and the calculated value of the IRI model according to the step A2, and setting the difference value at the ith sample point as

And

calculating an arbitrary position of an imaging region

Difference of ionosphere foF2 and hmF2 with respect to IRI model

The calculation method is as follows:

wherein: w is a_foF2,i、w_hmF2,iWeighting coefficients representing ionosphere foF2 and hmF2, respectively, with N representing the total number of samples observed;

step A5: solving the weighting coefficient according to the unbiased estimation and the optimal estimation theory of Kriging (Kriging), and jointly solving the following (N +1) equation sets, wherein the expression method of the equation sets is as follows:

wherein: d_ijThe ionospheric equivalent distance between the ith sample point and the jth sample point is represented, and lambda is a Lagrange multiplier;

step A6: converting the joint solution equation into a linear equation, wherein a matrix solution expression of the linear equation is as follows:

Aw＝b

in the formula: a is a matrix of (N +1) × (N +1), w is a vector of (N +1) × 1, is an unknown coefficient to be solved, and comprises N weighting coefficients and a vector of lambda, and b is (N +1) × 1;

step A7: using least squares estimation, arbitrary positions can be calculated

The weighting coefficient is calculated as follows:

w＝(A^TA)^-1A^Tb

step A8: according to the calculation method in the step A4, the solved weighting coefficients are substituted into the formula again to reconstruct the increment value

And

foF for the region is obtained by adding the output of IRI model₂And hmF₂Distribution:

step A9: traversing all ionospheric grids in the region to be imaged, and respectively reconstructing ionospheric layers foF2 and hmF2 of all grids in the region according to the steps A4-A8;

the reconstruction result of the critical frequency of the background ionosphere model F2 layer based on the vertical measurement data is shown in FIG. 2.

And B, calculating the total inclined electron content of the ground-based GNSS ionosphere:

step B1: according to the same date as the vertical observation data, acquiring a GNSS observation file and a GNSS satellite ephemeris file issued by an IGS (International Global navigation System) organization, decompressing the data, extracting the longitude, the latitude and the height of an observation station, and recording the longitude, the latitude and the height as

Extracting pseudo-range between GNSS satellite and observation station

And carrier phase

Wherein i is the corresponding observation frequency band, i is 1,2, j is the corresponding satellite number, and r is the corresponding receiver number;

step B2: for all pseudorange and carrier phase observations, a Geometry-Free linear Combination (Geometry-Free Combination) is calculated, respectively, as follows:

wherein:

is a geometry-free linear combination of pseudoranges,

is a geometric distance-free linear combination of carrier phases;

step B3: in the same observation arc segment, smoothing filtering is carried out on the non-geometric distance linear combination P4 of the pseudo-range by utilizing the non-geometric distance linear combination L4 of the carrier phase, and the calculation method is expressed as follows:

wherein: t represents epoch time, N_arcRepresenting the total number of observation epochs in an overhead arc section;

step B4: calculating an elevation angle E and an azimuth angle A of the receiver relative to the GNSS satellite according to the coordinate of the receiver and the GNSS satellite coordinate obtained by satellite ephemeris resolving;

the step B4 specifically includes:

step B401: longitude and latitude high coordinates of receiver and GNSS satellite

Converting into space rectangular XYZ coordinates, and expressing the conversion expression as:

wherein

R_eWhich is the radius of the earth, is,

step B402: converting the space rectangular coordinate into a station center rectangular coordinate (N, E, H), wherein the conversion expression is as follows:

wherein: t is a rotation matrix, (X)₀,Y₀,Z₀) Is the rectangular coordinate of the receiver, (X)₁,Y₁,Z₁) For the GNSS satellite rectangular coordinate, the calculation method of the rotation matrix T is as follows:

wherein:

longitude and latitude coordinates of the GNSS receiver;

step B403: calculating the elevation angle E between the receiver and the satellite, wherein the calculation expression is as follows:

step B404: calculating the azimuth angle A between the receiver and the satellite, wherein the calculation expression is as follows:

step B5: generally, the total concentration TEC of the ionized layer electrons is assumed to be concentrated in a thin layer at the height of 450km above the earth, and the geographic coordinates of puncture points of the observation station and the GNSS satellite which are connected on the height of the thin layer are calculated

The calculation formula is as follows:

wherein:

latitude and longitude coordinates of the GNSS observation station are shown, and A is an azimuth angle; psi₀For the geocentric angle between the observation station and the GNSS satellite, the calculation formula is as follows:

wherein: r_e6371.2km is the average radius of the earth; h is₀Is the station height; h is_I450km is the height of an ionosphere thin layer, and E is the elevation angle between a receiver and a satellite;

step B6: and calculating a mapping function value of the total inclined electron content and the total vertical electron content by the following calculation method:

M(E)＝[1-(R_ecosE/(R_e+h_I))²]^-1/2

step B7: a model of the vertical total electron content of the regional ionized layer along with the spatial distribution is constructed by utilizing a quadratic Taylor (Tayler) expansion series, and the calculation method is as follows:

wherein: n is_maxThe maximum order of the expansion is represented,

and

respectively the geographic latitude and the geographic longitude of the puncture point between the r-th receiver and the j-th satellite,

and λ_sSet as the center point geographical latitude and geographical longitude, E, respectively, of the tomographic imaging area_nmIs an unknown coefficient to be solved;

step B8: geometry-free linear combination using smoothed pseudoranges

Constructing a total electron content mapping solution equation, wherein the calculation method comprises the following steps:

wherein: f. of₁Is the first frequency of the satellite, f₂Being the second frequency of the satellite, DCB^rFor the hardware delay of the r-th receiver, DCB_jFor the hardware delay value of the jth satellite, in order to separate the hardware delays of the satellite and the receiver, the following four constraints are added:

wherein:

hardware delay value, N, representing the jth GPS satellite_GThe total number of GPS satellites;

denotes the hardware delay value, N, of the jth GLONASS satellite_RThe total number of GLONASS satellites;

hardware delay value N representing the jth Beidou satellite_CThe total number of the Beidou satellites;

hardware delay value, N, for the jth GALILEO satellite_EThe total number of GALILEO satellites;

step B9: and (3) selecting 1-2 hours as a time window, traversing all satellites and receivers in the time window, and establishing a linear matrix of the vertical total electron content according to the equation in the step B8 as follows:

Y＝AX

wherein: y is the observed quantity between all receivers and satellite

A is a coefficient matrix related to spherical harmonics, receivers and satellite hardware delays, and X is an unknown coefficient vector to be solved, and the specific expression mode is as follows:

step B10: solving a fitting equation by using a least square algorithm, and solving to obtain an unknown coefficient X, wherein the calculation method comprises the following steps:

X_LS＝(A^TA)^-1A^TY

wherein: x_LSIs a least squares solution comprising coefficients of expansion E in the Taylor (Taylor) series_nmReceiver and satellite hardware delays;

step B11: and (3) substituting the resolved satellite and receiver hardware delay DCB into the Taylor expansion series again to resolve and obtain the ionized layer tilt TEC value, wherein the calculation method comprises the following steps:

wherein: f. of₁Is the first frequency of the satellite, f₂Is the second frequency of the satellite or satellites,

calculating to obtain the total electron content of the inclined ionized layer between the r receiver and the j satellite;

the calculated total electron content of the ground-based GNSS ionosphere tilt is shown in FIG. 3

Step C, constructing an ionospheric tomography inversion matrix:

step C1: defining the geographic coordinate range and grid size of ionospheric tomography region, and defining the latitude range of

Step size of latitude grid is

Longitude range of [ lambda ]_min,λ_max]The longitude grid step is d lambda; height range of [ h_min,h_max]The height grid step is dh, and the calculation method of the ith longitude, latitude and height grid point is as follows:

λ_grid,i＝λ_min+(i-1)dλ，i＝1,2...,N_λ，N_λ＝(λ_max-λ_min)/dλ

h_grid,i＝h_min+(k-1)dh，k＝1,2...,N_h，N_h＝(h_max-h_min)/dh

wherein: n is a radical of_λ，

N_hThe numbers of grid points respectively representing the longitudinal, latitudinal and elevational directions of the tomographic region;

step C2: GNSS receiver

And GNSS satellite longitude and latitude high coordinate

Converted into spatial rectangular XYZ coordinates, respectively denoted (X)_r,Y_r,Z_r) And (X)_s,Y_s,Z_s) The conversion expression is expressed as:

wherein

R_eRepresenting the radius of the earth, e²＝0.00669437999013；

Step C3: and establishing a ray equation between the receiver and the satellite according to the relative position between the receiver and the satellite, and calculating the expression as follows:

wherein (X, Y, Z) represents the rectangular coordinate of any point on the ray, K represents a proportionality constant, and any intersection point coordinate of the ray path can be determined as long as the value of K is determined;

step C4: according to the relative position between the receiver and the GNSS satellite, the longitude plane L, the latitude plane B and the altitude plane H of the intersection of the ray and the tomography area are determined, and the determination method is as follows:

wherein: the longitude plane L, the latitude plane B and the altitude plane H must satisfy

B∈[λ_min,λ_max]、 H∈[h_min,h_max]；

Step C5: calculating the intersection (X) of a ray with the longitudinal, latitudinal and elevational planes of a tomographic grid_i,Y_i,Z_i) Wherein i is the intersection number;

the step C5 specifically includes:

step C501: calculating the intersection point of the ray and the longitude plane, wherein the tomography forms a curved surface along the longitude direction, namely a plane passing through the Z axis, and the included angle between the curved surface and the initial meridian plane is longitude L, and the equation of the curved surface is expressed as follows:

Y＝X tan(L)

step C502: and C501, solving the curved surface equation in the step C501 and the ray equation in the step C3 in a simultaneous manner to obtain a proportionality constant K corresponding to the intersection point of the ray and the longitude plane, wherein the calculation method is as follows:

step C503: the proportionality constant K is substituted into the ray equation of the step C3 again to obtain the coordinates of the intersection point of the ray and the longitude plane L, the coordinates are added into the intersection point sequence of the ray and the grid plane,

step C504: calculating the intersection point of the ray and the latitude plane, wherein the curved surface formed along the latitude direction is a conical surface with the vertex at the origin and the rotating axis as the z axis, the half angle of the cone of the conical surface is complementary with the latitude B, and the equation of the curved surface is expressed as follows:

(X²+Y²)·tan²W-Z²＝0

step C505: and C504, solving the curved surface equation in the step C504 and the ray equation in the step C3 in a simultaneous manner to obtain a proportionality constant K corresponding to the intersection point of the ray and the longitude plane, wherein the calculation method is as follows:

A＝(△X tan B)²+(△Y tan B)²-(△Z)²

B＝2×((△X·X_r·tan B)²+(△Y·Y_r·tan B)²-△Z·Z_r)

C＝(X_r·tan B)²+(Y_r·tan B)²-(Z_r)²

wherein: A. b, C are root coefficients of a quadratic equation;

step C506: substituting the proportionality constant K into the ray equation in the step C3 again to obtain the coordinates of the intersection point of the ray and the latitude plane B, judging whether the intersection point is between the starting point and the ending point, if so, adding the intersection point into the effective point sequence of the ray, otherwise, rejecting the point;

step C507: calculating the intersection point of the ray and the height surface H, wherein in the region to be chromatographically imaged, a curved surface formed by the height surface and the same height point from the ground is a spherical surface with the center of the sphere at the origin, and the equation of the curved surface is as follows:

(X²+Y²+Z²)-(R_e+H)²＝0

in the formula: r_eIs the average earth radius;

step C508: and C507, solving the curved surface equation in the step C507 and the ray equation in the step C3 in a simultaneous manner to obtain a proportionality constant K corresponding to the intersection point of the ray and the longitude plane, wherein the calculation method is as follows:

α＝(△X)²+(△Y)²+(△Z)²

β＝2×(△X·X_r+△Y·Y+△Z·Z_r)

γ＝(X_r)²+(Y_r)²+(Z_r)²-(R_e+H)²

wherein α, β and gamma are root coefficients of a quadratic equation of one unit;

step C509: substituting the proportionality constant K into the ray equation of the step C3 again to obtain the coordinates of the intersection point of the ray and the height surface H, judging whether the intersection point is between the starting point and the ending point, if so, adding the intersection point into the effective point sequence of the ray, otherwise, rejecting the point;

step C510, classifying all the intersection points into a sequence, and calculating the distance △ R from all the intersection points to the receiver, wherein the calculation formula is as follows:

step C511, the distances △ R from the effective point sequences of all the cross points to the receiver are sorted according to a method from small to large, and the coordinates of the corresponding points are respectively put into a coordinate sequence (X)_i,Y_i,Z_i) 1, 2.·, n;

step C6: the coordinates of the cross point are expressed by rectangular coordinates (X)_i,Y_i,Z_i) Conversion to a geographical coordinate system

The calculation formula is as follows:

in the formula: e' is the second eccentricity of the reference ellipsoid;

a＝6378.137，b＝6356.752；

step (ii) ofC7: according to the geographic coordinates of the intersection

Calculating the midpoints between the intersections in order

The calculation method is as follows:

step C8: determining the index value of the lattice according to the midpoint

The determination method comprises the following steps:

ID_λ＝count(λ_o≥λ_grid)

ID_h＝count(h_o≥h_grid)

wherein: lambda [ alpha ]_gridRepresenting longitude grid points;

representing all latitude grid points; h is_gridRepresents height grid points, count (-) represents statistical counts;

step C9: and calculating the grid number of the midpoint of the cross point, wherein the number calculation method in the chromatography inversion matrix comprises the following steps:

wherein: n is a radical of_λIndicates the number of grids divided in the longitudinal direction,

grid representing latitudinal divisionsCounting;

step C10: calculating the intercept of the ray intersection point in the grid through a distance formula between two points, and filling the intercept into the chromatography inversion matrix A, wherein the calculation method comprises the following steps:

wherein: p represents the p-th ray used for tomography; subscripts i and i +1 respectively represent a starting point and an end point of the rays in the same grid, and all the rays are traversed, so that an inversion matrix of ionospheric tomography can be obtained;

the composition structure of the ionospheric tomography inversion matrix is shown in fig. 4.

Step D, time-varying three-dimensional ionosphere electron density inversion:

step D1: acquiring the total inclined electron content data of all GNSS receivers, and constructing a solving equation of ionosphere electron density inversion according to the geometric position relation between the receivers and satellites:

wherein: STEC denotes the ionospheric tilt total electron content, N, extracted by the GNSS receiver_eThe electron density value of the ionized layer is represented, and the upper and lower integral limits are the positions of the satellite and the receiver respectively;

step D2: partitioning inversion regions into

The grid adopts a series expansion algorithm and utilizes a group of space basis functions h_k(r) discretizing the electron density distribution by the following calculation method:

wherein: h is_k(r) is a basis function, a_kThe basis functions are selected as pixel basis functions for the weight coefficients, and are calculated as follows：

Step D3: selecting a particular time window, at which time a can be identified_kThe ionosphere electron density tomography equation is constructed without changing along with time, and the calculation method is as follows:

R_inrepresenting the intercept length of the ith ray within the nth grid point;

step D4: traversing ionized layer TEC data of all GNSS, and converting an ionized layer tomography equation into a problem of solving a linear equation set, wherein the calculation method comprises the following steps:

d＝Ax

wherein: d is a vector formed by ionized layer inclination TEC data, a matrix A is a chromatography inversion matrix given in the step C, and x is the electron density of each grid point to be solved;

step D5: obtaining the ionized layers foF2 and hmF2 of all grids in the region reconstructed in the step A, inputting the ionized layers foF2 and hmF2 into the IRI model again, calculating to obtain the electron density value of the ionized layer, and taking the electron density value as the initial value X of the background field of the tomography⁽⁰⁾The initial value satisfies the following condition:

wherein: n represents the total number of grids of the region to be tomography;

step D6: iterative solution is carried out on the tomography equation by adopting an improved multiplication algebra reconstruction algorithm, and the calculation formula is as follows:

wherein: the iteration sequence of the equation is generated by random numbers, ranging from 1,N](ii) a m is the total number of the ionosphere inclined total electron content observed values, and i is the inclined TEC ray number; j is a grid number;

and

(ii) electron density of the jth grid for iterations k +1 and k, respectively; a is_ijIs the element in the inversion matrix A; | | a_iThe | | is the total length of the ith ray path; d_iIs the total electron content of the ith ray path; λ is a relaxation factor;

step D7: take the larger relaxation factor to iterate, λ_kThe value is generally 3.0, and the cut-off threshold value in the iterative process meets tol<10^-8And then stopping iteration, wherein the threshold value calculation method comprises the following steps:

wherein: tol represents an iteration cut-off threshold;

step D8: replacing X with the final X value obtained in the step D6⁽⁰⁾As an initial value of the iteration, a smaller relaxation factor lambda can be used_kAnd (5) iterating again for a second iteration when the electron density of the ionized layer is 0.9, and taking the resolved electron density of the ionized layer as a final ionized layer tomography data product.

The time-varying three-dimensional ionospheric electron density inversion results are shown in fig. 5.

Claims (3)

1. An ionospheric tomography method based on vertical measurement data constraint is characterized by comprising the following steps:

step A, reconstructing a background ionosphere model F2 layer characteristic parameter based on vertical measurement data:

step A1: downloading vertical measurement automatic interpretation data issued by a GIRO (global radio ionospheric observation station) according to the appointed date, and extracting the geographical latitude of the vertical measurement station

Longitude lambda, ionosphere F2 layer critical frequency

Peak height of ionosphere F2 layer

Wherein i is the position number of the vertical measuring station, t,

λ_iRespectively corresponding to the ith vertical measuring station at the time t, latitude and longitude;

step A2: selecting international reference ionosphere IRI as a background ionosphere model, and calculating the ionosphere F2 layer critical frequency of the IRI model at the same time and position

And peak height

Calculating the difference between all the observed values and the model values, wherein the calculation method comprises the following steps:

△Z_foF2＝foF2_ion-foF2_IRI

△Z_hmF2＝hmF2_ion-hmF2_IRI

wherein △ Z_foF2、△Z_hmF2Represents the difference between ionosphere foF2 and hmF2 measured by the plumb bob, respectively, relative to the IRI model calculations;

step A3: calculating the equivalent distance of the ionized layer by the following method:

in the formula of_iAnd

respectively, as longitude and latitude of a vertical measuring point i, in degrees, lambda_jAnd

respectively longitude and latitude of a vertical measuring point j, wherein F is a distance weighting factor, and can take a value of 3.5 for low-latitude and high-latitude areas and a value of 2.0 for medium-latitude areas;

step A4: for any time t, calculating the difference value between ionosphere foF2 and ionosphere hmF2 observed vertically at the time and the calculated value of the IRI model according to the step A2, and setting the difference value at the ith sample point as

And

calculating an arbitrary position of an imaging region

Difference of ionosphere foF2 and hmF2 with respect to IRI model

The calculation method is as follows:

wherein: w is a_foF2,i、w_hmF2,iWeighting coefficients representing ionosphere foF2 and hmF2, respectively, with N representing the total number of samples observed;

step A5: solving the weighting coefficient according to the kriging unbiased estimation and the optimal estimation theory, and jointly solving the following (N +1) equation sets, wherein the equation set expression method is as follows:

wherein: d_ijThe ionospheric equivalent distance between the ith sample point and the jth sample point is represented, and lambda is a Lagrange multiplier;

step A6: converting the joint solution equation into a linear equation, wherein a matrix solution expression of the linear equation is as follows:

Aw＝b

in the formula: a is a matrix of (N +1) × (N +1), w is a vector of (N +1) × 1, is an unknown coefficient to be solved, and comprises N weighting coefficients and a vector of lambda, and b is (N +1) × 1;

step A7: using least squares estimation, arbitrary positions can be calculated

The weighting coefficient is calculated as follows:

w＝(A^TA)^-1A^Tb

step A8: according to the calculation method in the step A4, the solved weighting coefficients are substituted into the formula again to reconstruct the increment value

And

foF for the region is obtained by adding the output of IRI model₂And hmF₂Distribution:

step A9: traversing all ionospheric grids in the region to be imaged, and respectively reconstructing ionospheric layers foF2 and hmF2 of all grids in the region according to the steps A4-A8;

and B, calculating the total inclined electron content of the ground-based GNSS ionosphere:

step B1: according to the same date as the vertical observation data, acquiring a GNSS observation file and a GNSS satellite ephemeris file issued by an IGS (International Global navigation System) organization, decompressing the data, extracting the longitude, the latitude and the height of an observation station, and recording the longitude, the latitude and the height as

Extracting pseudo-range between GNSS satellite and observation station

And carrier phase

Wherein i is the corresponding observation frequency band, i is 1,2, j is the corresponding satellite number, and r is the corresponding receiver number;

step B2: respectively calculating a non-geometric distance linear combination for all pseudo-ranges and carrier phase observed values, wherein the calculation method comprises the following steps:

wherein:

is a geometry-free linear combination of pseudoranges,

is a geometric distance-free linear combination of carrier phases;

step B3: in the same observation arc segment, smoothing filtering is carried out on the non-geometric distance linear combination P4 of the pseudo-range by utilizing the non-geometric distance linear combination L4 of the carrier phase, and the calculation method is expressed as follows:

wherein: t represents epoch time, N_arcRepresenting the total number of observation epochs in an overhead arc section;

step B4: calculating an elevation angle E and an azimuth angle A of the receiver relative to the GNSS satellite according to the coordinate of the receiver and the GNSS satellite coordinate obtained by satellite ephemeris resolving;

step B5: assuming that total concentration TEC of ionized layer electrons is concentrated in a thin layer at the height of 450km above the earth, calculating the geographic coordinates of puncture points of the observation station and the GNSS satellite which are connected on the height of the thin layer

The calculation formula is as follows:

wherein:

latitude and longitude coordinates of the GNSS observation station are shown, and A is an azimuth angle; psi₀For the geocentric angle between the observation station and the GNSS satellite, the calculation formula is as follows:

wherein: r_e6371.2km is the average radius of the earth; h is₀Is the station height; h is_I450km is the height of an ionosphere thin layer, and E is the elevation angle between a receiver and a satellite;

step B6: and calculating a mapping function value of the total inclined electron content and the total vertical electron content by the following calculation method:

M(E)＝[1-(R_ecosE/(R_e+h_I))²]^-1/2

step B7: a model of the vertical total electron content of the regional ionized layer along with the spatial distribution is constructed by utilizing a quadratic Taylor expansion series, and the calculation method is as follows:

wherein: n is_maxThe maximum order of the expansion is represented,

and

respectively the geographic latitude and the geographic longitude of the puncture point between the r-th receiver and the j-th satellite,

and λ_sSet as the center point geographical latitude and geographical longitude, E, respectively, of the tomographic imaging area_nmIs an unknown coefficient to be solved;

step B8: geometry-free linear combination using smoothed pseudoranges

Constructing a total electron content mapping solution equation, wherein the calculation method comprises the following steps:

wherein: f. of₁Is the first frequency of the satellite, f₂Being the second frequency of the satellite, DCB^rFor the hardware delay of the r-th receiver, DCB_jFor the hardware delay value of the jth satellite, in order to separate the hardware delays of the satellite and the receiver, the following four constraints are added:

wherein:

hardware delay value, N, representing the jth GPS satellite_GThe total number of GPS satellites;

denotes the hardware delay value, N, of the jth GLONASS satellite_RThe total number of GLONASS satellites;

hardware delay value N representing the jth Beidou satellite_CThe total number of the Beidou satellites;

hardware delay value, N, for the jth GALILEO satellite_EThe total number of GALILEO satellites;

step B9: and (3) selecting 1-2 hours as a time window, traversing all satellites and receivers in the time window, and establishing a linear matrix of the vertical total electron content according to the equation in the step B8 as follows:

Y＝AX

wherein: y is the observed quantity between all receivers and satellite

A is a coefficient matrix related to spherical harmonics, receivers and satellite hardware delays, and X is an unknown coefficient vector to be solved, and the specific expression mode is as follows:

step B10: solving a fitting equation by using a least square algorithm, and solving to obtain an unknown coefficient X, wherein the calculation method comprises the following steps:

X_LS＝(A^TA)^-1A^TY

wherein: x_LSIs a least squares solution comprising Taylor series expansion coefficients E_nmReceiver and satellite hardware delays;

step B11: and (3) substituting the resolved satellite and receiver hardware delay DCB into the Taylor expansion series again to resolve and obtain the ionized layer tilt TEC value, wherein the calculation method comprises the following steps:

wherein: f. of₁Is the first frequency of the satellite, f₂Is the second frequency of the satellite or satellites,

calculating to obtain the total electron content of the inclined ionized layer between the r receiver and the j satellite;

step C, constructing an ionospheric tomography inversion matrix:

step C1: defining the geographic coordinate range and grid size of ionospheric tomography region, and defining the latitude range of

Step size of latitude grid is

Longitude range of [ lambda ]_min,λ_max]The longitude grid step is d lambda; height range of [ h_min,h_max]The height grid step is dh, and the calculation method of the ith longitude, latitude and height grid point is as follows:

λ_grid,i＝λ_min+(i-1)dλ，i＝1,2...,N_λ，N_λ＝(λ_max-λ_min)/dλ

h_grid,i＝h_min+(k-1)dh，k＝1,2...,N_h，N_h＝(h_max-h_min)/dh

wherein: n is a radical of_λ，

N_hThe numbers of grid points respectively representing the longitudinal, latitudinal and elevational directions of the tomographic region;

step C2: GNSS receiver

And GNSS satellite longitude and latitude high coordinate

Converted into spatial rectangular XYZ coordinates, respectively denoted (X)_r,Y_r,Z_r) And (X)_s,Y_s,Z_s) The conversion expression is expressed as:

wherein

R_eRepresenting the radius of the earth, e²＝0.00669437999013；

Step C3: and establishing a ray equation between the receiver and the satellite according to the relative position between the receiver and the satellite, and calculating the expression as follows:

wherein (X, Y, Z) represents the rectangular coordinate of any point on the ray, K represents a proportionality constant, and any intersection point coordinate of the ray path can be determined as long as the value of K is determined;

step C4: according to the relative position between the receiver and the GNSS satellite, the longitude plane L, the latitude plane B and the altitude plane H of the intersection of the ray and the tomography area are determined, and the determination method is as follows:

wherein: the longitude plane L, the latitude plane B and the altitude plane H must satisfy

B∈[λ_min,λ_max]、H∈[h_min,h_max]；

Step C5: calculating the intersection (X) of a ray with the longitudinal, latitudinal and elevational planes of a tomographic grid_i,Y_i,Z_i) Wherein i is the intersection number;

step C6: the coordinates of the cross point are expressed by rectangular coordinates (X)_i,Y_i,Z_i) Conversion to a geographical coordinate system

The calculation formula is as follows:

in the formula: e' is the second eccentricity of the reference ellipsoid;

a＝6378.137，b＝6356.752；

step C7: according to the geographic coordinates of the intersection

Calculating the midpoints between the intersections in order

The calculation method is as follows:

step C8: determining the index value of the lattice according to the midpoint

The determination method comprises the following steps:

ID_λ＝count(λ_o≥λ_grid)

ID_h＝count(h_o≥h_grid)

wherein: lambda [ alpha ]_gridRepresenting longitude grid points;

representing all latitude grid points; h is_gridRepresents height grid points, count (-) represents statistical counts;

step C9: and calculating the grid number of the midpoint of the cross point, wherein the number calculation method in the chromatography inversion matrix comprises the following steps:

wherein: n is a radical of_λIndicates the number of grids divided in the longitudinal direction,

a grid number indicating a division in a latitudinal direction;

step C10: calculating the intercept of the ray intersection point in the grid through a distance formula between two points, and filling the intercept into the chromatography inversion matrix A, wherein the calculation method comprises the following steps:

wherein: p represents the p-th ray used for tomography; subscripts i and i +1 respectively represent a starting point and an end point of the rays in the same grid, and all the rays are traversed, so that an inversion matrix of ionospheric tomography can be obtained;

step D, time-varying three-dimensional ionosphere electron density inversion:

step D1: acquiring the total inclined electron content data of all GNSS receivers, and constructing a solving equation of ionosphere electron density inversion according to the geometric position relation between the receivers and satellites:

wherein: STEC denotes the ionospheric tilt total electron content, N, extracted by the GNSS receiver_eThe electron density value of the ionized layer is represented, and the upper and lower integral limits are the positions of the satellite and the receiver respectively;

step D2: partitioning inversion regions into

The grid adopts a series expansion algorithm and utilizes a group of space basis functions h_k(r) discretizing the electron density distribution by the following calculation method:

wherein: h is_k(r) is a basis function, a_kFor the weighting coefficients, the basis functions select the pixel basis functions, the calculation method is as follows:

step D3: selecting a particular time window, at which time a can be identified_kDoes not change with time, thereby constructing an ionosphere electron density tomography equation,the calculation method is as follows:

R_inrepresenting the intercept length of the ith ray within the nth grid point;

step D4: traversing ionized layer TEC data of all GNSS, and converting an ionized layer tomography equation into a problem of solving a linear equation set, wherein the calculation method comprises the following steps:

d＝Ax

wherein: d is a vector formed by ionized layer inclination TEC data, a matrix A is a chromatography inversion matrix given in the step C, and x is the electron density of each grid point to be solved;

step D5: obtaining the ionized layers foF2 and hmF2 of all grids in the region reconstructed in the step A, inputting the ionized layers foF2 and hmF2 into the IRI model again, calculating to obtain the electron density value of the ionized layer, and taking the electron density value as the initial value X of the background field of the tomography⁽⁰⁾The initial value satisfies the following condition:

wherein: n represents the total number of grids of the region to be tomography;

step D6: iterative solution is carried out on the tomography equation by adopting an improved multiplication algebra reconstruction algorithm, and the calculation formula is as follows:

wherein: the iteration sequence of the equation is generated by random numbers in the range of 1, N](ii) a m is the total number of the ionosphere inclined total electron content observed values, and i is the inclined TEC ray number; j is a grid number;

and

(ii) electron density of the jth grid for iterations k +1 and k, respectively; a is_ijIs the element in the inversion matrix A; | | a_iThe | | is the total length of the ith ray path; d_iIs the total electron content of the ith ray path; λ is a relaxation factor;

step D7: take the larger relaxation factor to iterate, λ_kThe value is 3.0, and the cut-off threshold value in the iterative process meets tol<10^-8And then stopping iteration, wherein the threshold value calculation method comprises the following steps:

wherein: tol represents an iteration cut-off threshold;

step D8: replacing X with the final X value obtained in the step D6⁽⁰⁾As an initial value of the iteration, a smaller relaxation factor lambda can be used_kAnd (5) iterating again for a second iteration when the electron density of the ionized layer is 0.9, and taking the resolved electron density of the ionized layer as a final ionized layer tomography data product.

2. The ionospheric tomography method based on vertical measurement data constraint according to claim 1, wherein the step B4 specifically comprises:

step B401: longitude and latitude high coordinates of receiver and GNSS satellite

Converting into space rectangular XYZ coordinates, and expressing the conversion expression as:

wherein

R_eWhich is the radius of the earth, is,

step B402: converting the space rectangular coordinate into a station center rectangular coordinate (N, E, H), wherein the conversion expression is as follows:

wherein: t is a rotation matrix, (X)₀,Y₀,Z₀) Is the rectangular coordinate of the receiver, (X)₁,Y₁,Z₁) For the GNSS satellite rectangular coordinate, the calculation method of the rotation matrix T is as follows:

wherein:

longitude and latitude coordinates of the GNSS receiver;

step B403: calculating the elevation angle E between the receiver and the satellite, wherein the calculation expression is as follows:

step B404: calculating the azimuth angle A between the receiver and the satellite, wherein the calculation expression is as follows:

3. the ionospheric tomography method based on vertical measurement data constraint according to claim 1, wherein the step C5 specifically comprises:

step C501: calculating the intersection point of the ray and the longitude plane, wherein the tomography forms a curved surface along the longitude direction, namely a plane passing through the Z axis, and the included angle between the curved surface and the initial meridian plane is longitude L, and the equation of the curved surface is expressed as follows:

Y＝Xtan(L)

step C502: and C501, solving the curved surface equation in the step C501 and the ray equation in the step C3 in a simultaneous manner to obtain a proportionality constant K corresponding to the intersection point of the ray and the longitude plane, wherein the calculation method is as follows:

step C503: the proportionality constant K is substituted into the ray equation of the step C3 again to obtain the coordinates of the intersection point of the ray and the longitude plane L, the coordinates are added into the intersection point sequence of the ray and the grid plane,

step C504: calculating the intersection point of the ray and the latitude plane, wherein the curved surface formed along the latitude direction is a conical surface with the vertex at the origin and the rotating axis as the z axis, the half angle of the cone of the conical surface is complementary with the latitude B, and the equation of the curved surface is expressed as follows:

(X²+Y²)·tan²W-Z²＝0

step C505: and C504, solving the curved surface equation in the step C504 and the ray equation in the step C3 in a simultaneous manner to obtain a proportionality constant K corresponding to the intersection point of the ray and the longitude plane, wherein the calculation method is as follows:

A＝(△XtanB)²+(△YtanB)²-(△Z)²

B＝2×((△X·X_r·tanB)²+(△Y·Y_r·tanB)²-△Z·Z_r)

C＝(X_r·tanB)²+(Y_r·tanB)²-(Z_r)²

wherein: A. b, C are root coefficients of a quadratic equation;

step C506: substituting the proportionality constant K into the ray equation in the step C3 again to obtain the coordinates of the intersection point of the ray and the latitude plane B, judging whether the intersection point is between the starting point and the ending point, if so, adding the intersection point into the effective point sequence of the ray, otherwise, rejecting the point;

step C507: calculating the intersection point of the ray and the height surface H, wherein in the region to be chromatographically imaged, a curved surface formed by the height surface and the same height point from the ground is a spherical surface with the center of the sphere at the origin, and the equation of the curved surface is as follows:

(X²+Y²+Z²)-(R_e+H)²＝0

in the formula: r_eIs the average earth radius;

step C508: and C507, solving the curved surface equation in the step C507 and the ray equation in the step C3 in a simultaneous manner to obtain a proportionality constant K corresponding to the intersection point of the ray and the longitude plane, wherein the calculation method is as follows:

α＝(△X)²+(△Y)²+(△Z)²

β＝2×(△X·X_r+△Y·Y+△Z·Z_r)

γ＝(X_r)²+(Y_r)²+(Z_r)²-(R_e+H)²

wherein α, β and gamma are root coefficients of a quadratic equation of one unit;

step C509: substituting the proportionality constant K into the ray equation of the step C3 again to obtain the coordinates of the intersection point of the ray and the height surface H, judging whether the intersection point is between the starting point and the ending point, if so, adding the intersection point into the effective point sequence of the ray, otherwise, rejecting the point;

step C510, classifying all the intersection points into a sequence, and calculating the distance △ R from all the intersection points to the receiver, wherein the calculation formula is as follows:

step C511, the distances △ R from the effective point sequences of all the cross points to the receiver are sorted according to a method from small to large, and the coordinates of the corresponding points are respectively put into a coordinate sequence (X)_i,Y_i,Z_i) 1, 2.., n.

Images