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

Ionosphere tomography method based on vertical measurement data constraint Download PDF

Info

Publication number
CN111273335B
CN111273335B CN201911330615.XA CN201911330615A CN111273335B CN 111273335 B CN111273335 B CN 111273335B CN 201911330615 A CN201911330615 A CN 201911330615A CN 111273335 B CN111273335 B CN 111273335B
Authority
CN
China
Prior art keywords
satellite
ray
equation
follows
longitude
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201911330615.XA
Other languages
Chinese (zh)
Other versions
CN111273335A (en
Inventor
欧明
韩超
陈亮
冯健
熊雯
於晓
刘钝
陈龙江
宋方雷
陈丽
甄卫民
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Institute of Radio Wave Propagation CETC 22 Research Institute
Original Assignee
China Institute of Radio Wave Propagation CETC 22 Research Institute
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China Institute of Radio Wave Propagation CETC 22 Research Institute filed Critical China Institute of Radio Wave Propagation CETC 22 Research Institute
Priority to CN201911330615.XA priority Critical patent/CN111273335B/en
Publication of CN111273335A publication Critical patent/CN111273335A/en
Application granted granted Critical
Publication of CN111273335B publication Critical patent/CN111273335B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01T—MEASUREMENT OF NUCLEAR OR X-RADIATION
    • G01T1/00—Measuring X-radiation, gamma radiation, corpuscular radiation, or cosmic radiation
    • G01T1/29—Measurement performed on radiation beams, e.g. position or section of the beam; Measurement of spatial distribution of radiation
    • G01T1/2914—Measurement of spatial distribution of radiation
    • G01T1/2985—In depth localisation, e.g. using positron emitters; Tomographic imaging (longitudinal and transverse section imaging; apparatus for radiation diagnosis sequentially in different planes, steroscopic radiation diagnosis)
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/35—Constructional details or hardware or software details of the signal processing chain
    • G01S19/37—Hardware or software details of the signal processing chain
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01T—MEASUREMENT OF NUCLEAR OR X-RADIATION
    • G01T1/00—Measuring X-radiation, gamma radiation, corpuscular radiation, or cosmic radiation
    • G01T1/29—Measurement performed on radiation beams, e.g. position or section of the beam; Measurement of spatial distribution of radiation
    • G01T1/2914—Measurement of spatial distribution of radiation
    • G01T1/2992—Radioisotope data or image processing not related to a particular imaging system; Off-line processing of pictures, e.g. rescanners

Landscapes

  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Physics & Mathematics (AREA)
  • Molecular Biology (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Signal Processing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

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
Figure GDA0002477926470000021
Longitude lambda, ionosphere F2 layer critical frequency
Figure GDA0002477926470000022
Peak height of ionosphere F2 layer
Figure GDA0002477926470000023
Wherein i is the position number of the vertical measuring station, t,
Figure GDA0002477926470000024
λ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
Figure GDA0002477926470000025
And peak height
Figure GDA0002477926470000026
Calculating the difference between all the observed values and the model values, wherein the calculation method comprises the following steps:
△ZfoF2＝foF2ion-foF2IRI
△ZhmF2＝hmF2ion-hmF2IRI
wherein, Delta ZfoF2、△ZhmF2Represents 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:
Figure GDA0002477926470000027
in the formula ofiAnd
Figure GDA0002477926470000028
respectively, as longitude and latitude of a vertical measuring point i, in degrees, lambdajAnd
Figure GDA0002477926470000029
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
Figure GDA00024779264700000210
And
Figure GDA00024779264700000211
calculating an arbitrary position of an imaging region
Figure GDA00024779264700000212
Difference of ionosphere foF2 and hmF2 with respect to IRI model
Figure GDA00024779264700000213
The calculation method is as follows:
Figure GDA0002477926470000031
Figure GDA0002477926470000032
wherein: w is afoF2,i、whmF2,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:
Figure GDA0002477926470000033
wherein: dijThe 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
Figure GDA0002477926470000034
The weighting coefficient is calculated as follows:
w＝(ATA)-1ATb
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
Figure GDA0002477926470000035
And
Figure GDA0002477926470000036
foF for the region is obtained by adding the output of IRI model2And hmF2Distribution:
Figure GDA0002477926470000037
Figure GDA0002477926470000038
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
Figure GDA0002477926470000041
Extracting pseudo-range between GNSS satellite and observation station
Figure GDA0002477926470000042
And carrier phase
Figure GDA0002477926470000043
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:
Figure GDA0002477926470000044
Figure GDA0002477926470000045
wherein:
Figure GDA0002477926470000046
is a geometry-free linear combination of pseudoranges,
Figure GDA0002477926470000047
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:
Figure GDA0002477926470000048
wherein: t represents epoch time, NarcRepresenting 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 total ionospheric electron concentrationThe content TEC is concentrated in a thin layer with the height of 450km above the earth, and the geographic coordinates of the puncture points of the connecting lines of the observation station and the GNSS satellite on the height of the thin layer are calculated
Figure GDA0002477926470000049
The calculation formula is as follows:
Figure GDA00024779264700000410
wherein:
Figure GDA00024779264700000411
latitude and longitude coordinates of the GNSS observation station are shown, and A is an azimuth angle; psi0For the geocentric angle between the observation station and the GNSS satellite, the calculation formula is as follows:
Figure GDA00024779264700000412
wherein: re6371.2km is the average radius of the earth; h is0Is the station height; h isI450km 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-(RecosE/(Re+hI))2]-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:
Figure GDA0002477926470000051
wherein: n ismaxThe maximum order of the expansion is represented,
Figure GDA0002477926470000052
and
Figure GDA0002477926470000053
respectively the geographic latitude and the geographic longitude of the puncture point between the r-th receiver and the j-th satellite,
Figure GDA0002477926470000054
and λsSet as the center point geographical latitude and geographical longitude, E, respectively, of the tomographic imaging areanmIs an unknown coefficient to be solved;
step B8: geometry-free linear combination using smoothed pseudoranges
Figure GDA0002477926470000055
Constructing a total electron content mapping solution equation, wherein the calculation method comprises the following steps:
Figure GDA0002477926470000056
wherein: f. of1Is the first frequency of the satellite, f2Being the second frequency of the satellite, DCBrFor the hardware delay of the r-th receiver, DCBjFor 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:
Figure GDA0002477926470000057
wherein:
Figure GDA0002477926470000058
hardware delay value, N, representing the jth GPS satelliteGThe total number of GPS satellites;
Figure GDA0002477926470000059
denotes the hardware delay value, N, of the jth GLONASS satelliteRThe total number of GLONASS satellites;
Figure GDA00024779264700000510
hardware delay value N representing the jth Beidou satelliteCThe total number of the Beidou satellites;
Figure GDA0002477926470000061
hardware delay value, N, for the jth GALILEO satelliteEThe 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
Figure GDA0002477926470000062
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:
Figure GDA0002477926470000063
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:
XLS＝(ATA)-1ATY
wherein: xLSIs a least squares solution comprising Taylor series expansion coefficients EnmReceiver 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:
Figure GDA0002477926470000064
wherein: f. of1Is the first frequency of the satellite, f2Is the second frequency of the satellite or satellites,
Figure GDA0002477926470000069
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
Figure GDA0002477926470000065
Step size of latitude grid is
Figure GDA0002477926470000066
Longitude range of [ lambda ]min,λmax]The longitude grid step is d lambda; height range of [ hmin,hmax]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λ
Figure GDA0002477926470000067
hgrid,i＝hmin+(k-1)dh，k＝1,2...,Nh，Nh＝(hmax-hmin)/dh
wherein: n is a radical ofλ，
Figure GDA0002477926470000068
NhThe numbers of grid points respectively representing the longitudinal, latitudinal and elevational directions of the tomographic region;
step C2: GNSS receiver
Figure GDA0002477926470000071
And GNSS satellite longitude and latitudeHigh coordinate
Figure GDA0002477926470000072
Converted into spatial rectangular XYZ coordinates, respectively denoted (X)r,Yr,Zr) And (X)s,Ys,Zs) The conversion expression is expressed as:
Figure GDA0002477926470000073
wherein
Figure GDA0002477926470000074
ReRepresenting the radius of the earth, e2＝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:
Figure GDA0002477926470000075
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:
Figure GDA0002477926470000076
wherein: the longitude plane L, the latitude plane B and the altitude plane H must satisfy
Figure GDA0002477926470000077
B∈[λmin,λmax]、H∈[hmin,hmax]；
Step C5: computed ray and tomography gridCross point (X) of longitude plane, latitude plane and altitude plane ofi,Yi,Zi) Wherein i is the intersection number;
step C6: the coordinates of the cross point are expressed by rectangular coordinates (X)i,Yi,Zi) Conversion to a geographical coordinate system
Figure GDA0002477926470000078
The calculation formula is as follows:
Figure GDA0002477926470000081
in the formula: e' is the second eccentricity of the reference ellipsoid;
Figure GDA0002477926470000082
Figure GDA0002477926470000083
a＝6378.137，b＝6356.752；
step C7: according to the geographic coordinates of the intersection
Figure GDA0002477926470000084
Calculating the midpoints between the intersections in order
Figure GDA0002477926470000085
The calculation method is as follows:
Figure GDA0002477926470000086
step C8: determining the index value of the lattice according to the midpoint
Figure GDA0002477926470000087
The determination method comprises the following steps:
IDλ＝count(λo≥λgrid)
Figure GDA0002477926470000088
IDh＝count(ho≥hgrid)
wherein: lambda [ alpha ]gridRepresenting longitude grid points;
Figure GDA0002477926470000089
representing all latitude grid points; h isgridRepresents 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:
Figure GDA00024779264700000810
wherein: n is a radical ofλIndicates the number of grids divided in the longitudinal direction,
Figure GDA00024779264700000811
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:
Figure GDA0002477926470000091
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:
Figure GDA0002477926470000092
wherein: STEC denotes the ionospheric tilt total electron content, N, extracted by the GNSS receivereThe 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
Figure GDA0002477926470000093
The grid adopts a series expansion algorithm and utilizes a group of space basis functions hk(r) discretizing the electron density distribution by the following calculation method:
Figure GDA0002477926470000094
wherein: h isk(r) is a basis function, akFor the weighting coefficients, the basis functions select the pixel basis functions, the calculation method is as follows:
Figure GDA0002477926470000095
step D3: selecting a particular time window, at which time a can be identifiedkThe ionosphere electron density tomography equation is constructed without changing along with time, and the calculation method is as follows:
Figure GDA0002477926470000096
Rinrepresenting 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(0)The initial value satisfies the following condition:
Figure GDA0002477926470000101
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:
Figure GDA0002477926470000102
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;
Figure GDA0002477926470000103
and
Figure GDA0002477926470000104
(ii) electron density of the jth grid for iterations k +1 and k, respectively; a isijIs the element in the inversion matrix A; | | aiThe | | is the total length of the ith ray path; diIs 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:
Figure GDA0002477926470000105
wherein: tol represents an iteration cut-off threshold;
step D8: replacing X with the final X value obtained in the step D6(0)As an initial value of the iteration, a smaller relaxation factor lambda can be usedkAnd (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
Figure GDA0002477926470000106
Converting into space rectangular XYZ coordinates, and expressing the conversion expression as:
Figure GDA0002477926470000111
wherein
Figure GDA0002477926470000112
ReWhich is the radius of the earth, is,
Figure GDA0002477926470000113
step B402: converting the space rectangular coordinate into a station center rectangular coordinate (N, E, H), wherein the conversion expression is as follows:
Figure GDA0002477926470000114
wherein: t is a rotation matrix, (X)0,Y0,Z0) Is the rectangular coordinate of the receiver, (X)1,Y1,Z1) For the GNSS satellite rectangular coordinate, the calculation method of the rotation matrix T is as follows:
Figure GDA0002477926470000115
wherein:
Figure GDA0002477926470000116
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:
Figure GDA0002477926470000117
step B404: calculating the azimuth angle A between the receiver and the satellite, wherein the calculation expression is as follows:
Figure GDA0002477926470000118
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:
Figure GDA0002477926470000121
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,
Figure GDA0002477926470000122
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:
(X2+Y2)·tan2W-Z2＝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:
Figure GDA0002477926470000123
A＝(△X tanB)2+(△Y tanB)2-(△Z)2
B＝2×((△X·Xr·tan B)2+(△Y·Yr·tan B)2-△Z·Zr)
C＝(Xr·tan B)2+(Yr·tan B)2-(Zr)2
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:
(X2+Y2+Z2)-(Re+H)2＝0
in the formula: reIs 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:
Figure GDA0002477926470000131
α＝(△X)2+(△Y)2+(△Z)2
β＝2×(△X·Xr+△Y·Y+△Z·Zr)
γ＝(Xr)2+(Yr)2+(Zr)2-(Re+H)2
wherein: the alpha, beta and gamma are root coefficients of a quadratic equation;
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 cross points into a sequence, and calculating the distances delta R from all the cross points to the receiver, wherein the calculation formula is as follows:
Figure GDA0002477926470000132
step C511: sorting the distances delta R from the effective point sequences of all cross points to the receiver according to a method from small to large, and respectively putting the coordinates of corresponding points into a coordinate sequence (X)i,Yi,Zi) 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 Ionospheral Radio observer, GIRO) according to appointed date, and extracting geographic latitude of the vertical survey station
Figure GDA0002477926470000151
Longitude lambda, ionosphere F2 layer critical frequency
Figure GDA0002477926470000152
Peak height of ionosphere F2 layer
Figure GDA0002477926470000153
Wherein i is the position number of the vertical measuring station, t,
Figure GDA0002477926470000154
λ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
Figure GDA0002477926470000155
And peak height
Figure GDA0002477926470000156
Calculating the difference between all the observed values and the model values, wherein the calculation method comprises the following steps:
△ZfoF2＝foF2ion-foF2IRI
△ZhmF2＝hmF2ion-hmF2IRI
wherein, Delta ZfoF2、△ZhmF2Represents 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:
Figure GDA0002477926470000157
in the formula ofiAnd
Figure GDA0002477926470000158
respectively, as the longitude and latitude (in degrees), λ, of the vertical point ijAnd
Figure GDA0002477926470000159
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 between ionosphere foF2 and ionosphere hmF2 observed vertically at that time relative to the IRI model calculation value according to the step A2, and settingDefining the difference value at the ith sample point as
Figure GDA00024779264700001510
And
Figure GDA00024779264700001511
calculating an arbitrary position of an imaging region
Figure GDA00024779264700001512
Difference of ionosphere foF2 and hmF2 with respect to IRI model
Figure GDA00024779264700001513
The calculation method is as follows:
Figure GDA0002477926470000161
Figure GDA0002477926470000162
wherein: w is afoF2,i、whmF2,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:
Figure GDA0002477926470000163
wherein: dijThe 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
Figure GDA0002477926470000164
The weighting coefficient is calculated as follows:
w＝(ATA)-1ATb
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
Figure GDA0002477926470000165
And
Figure GDA0002477926470000166
foF for the region is obtained by adding the output of IRI model2And hmF2Distribution:
Figure GDA0002477926470000167
Figure GDA0002477926470000168
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
Figure GDA0002477926470000171
Extracting pseudo-range between GNSS satellite and observation station
Figure GDA0002477926470000172
And carrier phase
Figure GDA0002477926470000173
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:
Figure GDA0002477926470000174
Figure GDA0002477926470000175
wherein:
Figure GDA0002477926470000176
is a geometry-free linear combination of pseudoranges,
Figure GDA0002477926470000177
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:
Figure GDA0002477926470000178
wherein: t represents epoch time, NarcRepresenting 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
Figure GDA0002477926470000179
Converting into space rectangular XYZ coordinates, and expressing the conversion expression as:
Figure GDA00024779264700001710
wherein
Figure GDA0002477926470000181
ReWhich is the radius of the earth, is,
Figure GDA0002477926470000182
step B402: converting the space rectangular coordinate into a station center rectangular coordinate (N, E, H), wherein the conversion expression is as follows:
Figure GDA0002477926470000183
wherein: t is a rotation matrix, (X)0,Y0,Z0) Is the rectangular coordinate of the receiver, (X)1,Y1,Z1) For the GNSS satellite rectangular coordinate, the calculation method of the rotation matrix T is as follows:
Figure GDA0002477926470000184
wherein:
Figure GDA0002477926470000185
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:
Figure GDA0002477926470000186
step B404: calculating the azimuth angle A between the receiver and the satellite, wherein the calculation expression is as follows:
Figure GDA0002477926470000187
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
Figure GDA0002477926470000188
The calculation formula is as follows:
Figure GDA0002477926470000189
wherein:
Figure GDA00024779264700001810
latitude and longitude coordinates of the GNSS observation station are shown, and A is an azimuth angle; psi0For the geocentric angle between the observation station and the GNSS satellite, the calculation formula is as follows:
Figure GDA0002477926470000191
wherein: re6371.2km is the average radius of the earth; h is0Is the station height; h isI450km 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-(RecosE/(Re+hI))2]-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:
Figure GDA0002477926470000192
wherein: n ismaxThe maximum order of the expansion is represented,
Figure GDA0002477926470000193
and
Figure GDA0002477926470000194
respectively the geographic latitude and the geographic longitude of the puncture point between the r-th receiver and the j-th satellite,
Figure GDA0002477926470000195
and λsSet as the center point geographical latitude and geographical longitude, E, respectively, of the tomographic imaging areanmIs an unknown coefficient to be solved;
step B8: geometry-free linear combination using smoothed pseudoranges
Figure GDA0002477926470000196
Constructing a total electron content mapping solution equation, wherein the calculation method comprises the following steps:
Figure GDA0002477926470000197
wherein: f. of1Is the first frequency of the satellite, f2Being the second frequency of the satellite, DCBrFor the hardware delay of the r-th receiver, DCBjFor 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:
Figure GDA0002477926470000198
wherein:
Figure GDA0002477926470000201
hardware delay value, N, representing the jth GPS satelliteGThe total number of GPS satellites;
Figure GDA0002477926470000202
denotes the hardware delay value, N, of the jth GLONASS satelliteRThe total number of GLONASS satellites;
Figure GDA0002477926470000203
hardware delay value N representing the jth Beidou satelliteCThe total number of the Beidou satellites;
Figure GDA0002477926470000204
hardware delay value, N, for the jth GALILEO satelliteEThe 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
Figure GDA0002477926470000205
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:
Figure GDA0002477926470000206
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:
XLS＝(ATA)-1ATY
wherein: xLSIs a least squares solution comprising coefficients of expansion E in the Taylor (Taylor) seriesnmReceiver 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:
Figure GDA0002477926470000207
wherein: f. of1Is the first frequency of the satellite, f2Is the second frequency of the satellite or satellites,
Figure GDA0002477926470000208
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
Figure GDA0002477926470000209
Step size of latitude grid is
Figure GDA00024779264700002010
Longitude range of [ lambda ]min,λmax]The longitude grid step is d lambda; height range of [ hmin,hmax]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λ
Figure GDA0002477926470000211
hgrid,i＝hmin+(k-1)dh，k＝1,2...,Nh，Nh＝(hmax-hmin)/dh
wherein: n is a radical ofλ，
Figure GDA0002477926470000212
NhThe numbers of grid points respectively representing the longitudinal, latitudinal and elevational directions of the tomographic region;
step C2: GNSS receiver
Figure GDA0002477926470000213
And GNSS satellite longitude and latitude high coordinate
Figure GDA0002477926470000214
Converted into spatial rectangular XYZ coordinates, respectively denoted (X)r,Yr,Zr) And (X)s,Ys,Zs) The conversion expression is expressed as:
Figure GDA0002477926470000215
wherein
Figure GDA0002477926470000216
ReRepresenting the radius of the earth, e2＝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:
Figure GDA0002477926470000217
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:
Figure GDA0002477926470000218
wherein: the longitude plane L, the latitude plane B and the altitude plane H must satisfy
Figure GDA0002477926470000219
B∈[λmin,λmax]、H∈[hmin,hmax]；
Step C5: calculating the intersection (X) of a ray with the longitudinal, latitudinal and elevational planes of a tomographic gridi,Yi,Zi) 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:
Figure GDA0002477926470000221
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,
Figure GDA0002477926470000222
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:
(X2+Y2)·tan2W-Z2＝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:
Figure GDA0002477926470000223
A＝(△X tan B)2+(△Y tan B)2-(△Z)2
B＝2×((△X·Xr·tan B)2+(△Y·Yr·tan B)2-△Z·Zr)
C＝(Xr·tan B)2+(Yr·tan B)2-(Zr)2
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:
(X2+Y2+Z2)-(Re+H)2＝0
in the formula: reIs 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:
Figure GDA0002477926470000231
α＝(△X)2+(△Y)2+(△Z)2
β＝2×(△X·Xr+△Y·Y+△Z·Zr)
γ＝(Xr)2+(Yr)2+(Zr)2-(Re+H)2
wherein: the alpha, beta and gamma are root coefficients of a quadratic equation;
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 cross points into a sequence, and calculating the distances delta R from all the cross points to the receiver, wherein the calculation formula is as follows:
Figure GDA0002477926470000232
step C511: sorting the distances delta R from the effective point sequences of all cross points to the receiver according to a method from small to large, and respectively putting the coordinates of corresponding points into a coordinate sequence (X)i,Yi,Zi) 1, 2.·, n;
step C6: the coordinates of the cross point are expressed by rectangular coordinates (X)i,Yi,Zi) Conversion to a geographical coordinate system
Figure GDA0002477926470000233
The calculation formula is as follows:
Figure GDA0002477926470000241
in the formula:e' is the second eccentricity of the reference ellipsoid;
Figure GDA0002477926470000242
Figure GDA0002477926470000243
a＝6378.137，b＝6356.752；
step C7: according to the geographic coordinates of the intersection
Figure GDA0002477926470000244
Calculating the midpoints between the intersections in order
Figure GDA0002477926470000245
The calculation method is as follows:
Figure GDA0002477926470000246
step C8: determining the index value of the lattice according to the midpoint
Figure GDA0002477926470000247
The determination method comprises the following steps:
IDλ＝count(λo≥λgrid)
Figure GDA0002477926470000248
IDh＝count(ho≥hgrid)
wherein: lambda [ alpha ]gridRepresenting longitude grid points;
Figure GDA0002477926470000249
representing all latitude grid points; h isgridRepresents 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:
Figure GDA00024779264700002410
wherein: n is a radical ofλIndicates the number of grids divided in the longitudinal direction,
Figure GDA00024779264700002411
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:
Figure GDA0002477926470000251
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:
Figure GDA0002477926470000252
wherein: STEC denotes the ionospheric tilt total electron content, N, extracted by the GNSS receivereThe 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
Figure GDA0002477926470000253
Personal netLattice, using a series expansion algorithm, using a set of spatial basis functions hk(r) discretizing the electron density distribution by the following calculation method:
Figure GDA0002477926470000254
wherein: h isk(r) is a basis function, akFor the weighting coefficients, the basis functions select the pixel basis functions, the calculation method is as follows:
Figure GDA0002477926470000255
step D3: selecting a particular time window, at which time a can be identifiedkThe ionosphere electron density tomography equation is constructed without changing along with time, and the calculation method is as follows:
Figure GDA0002477926470000256
Rinrepresenting 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(0)The initial value satisfies the following condition:
Figure GDA0002477926470000261
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:
Figure GDA0002477926470000262
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;
Figure GDA0002477926470000263
and
Figure GDA0002477926470000264
(ii) electron density of the jth grid for iterations k +1 and k, respectively; a isijIs the element in the inversion matrix A; | | aiThe | | is the total length of the ith ray path; diIs 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:
Figure GDA0002477926470000265
wherein: tol represents an iteration cut-off threshold;
step D8: replacing X with the final X value obtained in the step D6(0)As an initial value of the iteration, a smaller relaxation factor lambda can be usedkAnd (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 characteristic parameters of a background ionosphere F2 layer 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
Figure FDA0003176611340000011
Longitude lambda, ionosphere F2 layer critical frequency
Figure FDA0003176611340000012
Peak height of ionosphere F2 layer
Figure FDA0003176611340000013
Wherein i is the position number of the vertical measuring station, t,
Figure FDA0003176611340000014
λ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
Figure FDA0003176611340000015
And peak height
Figure FDA0003176611340000016
Calculating the difference between all the observed values and the model values, wherein the calculation method comprises the following steps:
ΔZfoF2＝foF2ion-foF2IRI
ΔZhmF2＝hmF2ion-hmF2IRI
wherein, Delta ZfoF2、ΔZhmF2Represents 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:
Figure FDA0003176611340000017
in the formula ofiAnd
Figure FDA0003176611340000018
respectively, as longitude and latitude of a vertical measuring point i, in degrees, lambdajAnd
Figure FDA0003176611340000019
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
Figure FDA00031766113400000110
And
Figure FDA00031766113400000111
calculating an arbitrary position of an imaging region
Figure FDA00031766113400000112
Difference of ionosphere foF2 and hmF2 with respect to IRI model
Figure FDA00031766113400000113
The calculation method is as follows:
Figure FDA00031766113400000114
Figure FDA0003176611340000021
wherein: w is afoF2,i、whmF2,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:
Figure FDA0003176611340000022
wherein: dijThe 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
Figure FDA0003176611340000027
The weighting coefficient is calculated as follows:
w＝(ATA)-1ATb
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
Figure FDA0003176611340000023
And
Figure FDA0003176611340000024
foF for the region is obtained by adding the output of IRI model2And hmF2Distribution:
Figure FDA0003176611340000025
Figure FDA0003176611340000026
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
Figure FDA0003176611340000031
Extracting pseudo-range between GNSS satellite and observation station
Figure FDA0003176611340000032
And carrier phase
Figure FDA0003176611340000033
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:
Figure FDA0003176611340000034
Figure FDA0003176611340000035
wherein:
Figure FDA0003176611340000036
is a geometry-free linear combination of pseudoranges,
Figure FDA0003176611340000037
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:
Figure FDA0003176611340000038
wherein: t represents epoch time, NarcRepresenting 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
Figure FDA0003176611340000039
The calculation formula is as follows:
Figure FDA00031766113400000310
wherein:
Figure FDA00031766113400000311
latitude and longitude coordinates of the GNSS observation station are shown, and A is an azimuth angle; psi0For the geocentric angle between the observation station and the GNSS satellite, the calculation formula is as follows:
Figure FDA00031766113400000312
wherein: re6371.2km is the average radius of the earth; h is0Is the station height; h isI450km 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-(RecosE/(Re+hI))2]-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:
Figure FDA0003176611340000041
wherein: n ismaxThe maximum order of the expansion is represented,
Figure FDA0003176611340000042
and
Figure FDA0003176611340000043
respectively the geographic latitude and the geographic longitude of the puncture point between the r-th receiver and the j-th satellite,
Figure FDA0003176611340000044
and λsSet as the center point geographical latitude and geographical longitude, E, respectively, of the tomographic imaging areanmAs unknown to be solvedA coefficient;
step B8: geometry-free linear combination using smoothed pseudoranges
Figure FDA0003176611340000045
Constructing a total electron content mapping solution equation, wherein the calculation method comprises the following steps:
Figure FDA0003176611340000046
wherein: f. of1Is the first frequency of the satellite, f2Being the second frequency of the satellite, DCBrFor the hardware delay of the r-th receiver, DCBjFor 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:
Figure FDA0003176611340000047
wherein:
Figure FDA0003176611340000048
hardware delay value, N, representing the jth GPS satelliteGThe total number of GPS satellites;
Figure FDA0003176611340000049
denotes the hardware delay value, N, of the jth GLONASS satelliteRThe total number of GLONASS satellites;
Figure FDA00031766113400000410
hardware delay value N representing the jth Beidou satelliteCThe total number of the Beidou satellites;
Figure FDA00031766113400000411
hardware delay value, N, for the jth GALILEO satelliteEThe 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
Figure FDA0003176611340000051
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:
Figure FDA0003176611340000052
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:
XLS＝(ATA)-1ATY
wherein: xLSIs a least squares solution comprising Taylor series expansion coefficients EnmReceiver 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:
Figure FDA0003176611340000053
wherein: f. of1Is the first frequency of the satellite, f2Is the second frequency of the satellite or satellites,
Figure FDA0003176611340000054
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
Figure FDA0003176611340000055
Step size of latitude grid is
Figure FDA0003176611340000056
Longitude range of [ lambda ]min,λmax]The longitude grid step is d lambda; height range of [ hmin,hmax]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λ
Figure FDA0003176611340000057
hgrid,i＝hmin+(k-1)dh，k＝1,2...,Nh，Nh＝(hmax-hmin)/dh
wherein: n is a radical ofλ，
Figure FDA0003176611340000058
NhThe numbers of grid points respectively representing the longitudinal, latitudinal and elevational directions of the tomographic region;
step C2: GNSS receiver
Figure FDA0003176611340000059
And GNSS satellite longitude and latitude high coordinate
Figure FDA00031766113400000510
Converted into spatial rectangular XYZ coordinates, respectively denoted (X)r,Yr,Zr) And (X)s,Ys,Zs) The conversion expression is expressed as:
Figure FDA0003176611340000061
wherein
Figure FDA0003176611340000062
ReRepresenting the radius of the earth, e2＝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:
Figure FDA0003176611340000063
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:
Figure FDA0003176611340000064
wherein: the longitude plane L, the latitude plane B and the altitude plane H must satisfy
Figure FDA0003176611340000065
B∈[λmin,λmax]、H∈[hmin,hmax]；
Step C5: calculating the intersection (X) of a ray with the longitudinal, latitudinal and elevational planes of a tomographic gridi,Yi,Zi) Wherein i is the intersection number;
step C6: the coordinates of the cross point are expressed by rectangular coordinates (X)i,Yi,Zi) Conversion to a geographical coordinate system
Figure FDA0003176611340000066
The calculation formula is as follows:
Figure FDA0003176611340000067
in the formula: e' is the second eccentricity of the reference ellipsoid;
Figure FDA0003176611340000071
Figure FDA0003176611340000072
a＝6378.137，b＝6356.752；
step C7: according to the geographic coordinates of the intersection
Figure FDA0003176611340000073
Calculating the midpoints between the intersections in order
Figure FDA0003176611340000074
The calculation method is as follows:
Figure FDA0003176611340000075
step C8: determining the index value of the lattice according to the midpoint
Figure FDA0003176611340000076
The determination method comprises the following steps:
IDλ＝count(λo≥λgrid)
Figure FDA0003176611340000077
IDh＝count(ho≥hgrid)
wherein: lambda [ alpha ]gridRepresenting longitude grid points;
Figure FDA0003176611340000078
representing all latitude grid points; h isgridRepresents 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:
Figure FDA0003176611340000079
wherein: n is a radical ofλIndicates the number of grids divided in the longitudinal direction,
Figure FDA00031766113400000710
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:
Figure FDA00031766113400000711
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:
Figure FDA0003176611340000081
wherein: STEC denotes the ionospheric tilt total electron content, N, extracted by the GNSS receivereThe 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
Figure FDA0003176611340000082
The grid adopts a series expansion algorithm and utilizes a group of space basis functions hk(r) discretizing the electron density distribution by the following calculation method:
Figure FDA0003176611340000083
wherein: h isk(r) is a basis function, akFor the weighting coefficients, the basis functions select the pixel basis functions, the calculation method is as follows:
Figure FDA0003176611340000084
step D3: selecting a particular time window, at which time a can be identifiedkThe ionosphere electron density tomography equation is constructed without changing along with time, and the calculation method is as follows:
Figure FDA0003176611340000085
ΔRinrepresenting 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(0)The initial value satisfies the following condition:
Figure FDA0003176611340000086
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:
Figure FDA0003176611340000091
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;
Figure FDA0003176611340000092
and
Figure FDA0003176611340000093
(ii) electron density of the jth grid for iterations k +1 and k, respectively; a isijIs the element in the inversion matrix A; | | aiThe | | is the total length of the ith ray path; diIs 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:
Figure FDA0003176611340000094
wherein: tol represents an iteration cut-off threshold;
step D8: replacing X with the final X value obtained in the step D6(0)As an initial value of the iteration, a smaller relaxation factor lambda can be usedkAnd (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
Figure FDA0003176611340000095
Converting into space rectangular XYZ coordinates, and expressing the conversion expression as:
Figure FDA0003176611340000096
wherein
Figure FDA0003176611340000097
ReWhich is the radius of the earth, is,
Figure FDA0003176611340000098
step B402: converting the space rectangular coordinate into a station center rectangular coordinate (N, E, U), wherein the conversion expression is as follows:
Figure FDA0003176611340000101
wherein: t is a rotation matrix, (X)0,Y0,Z0) Is the rectangular coordinate of the receiver, (X)1,Y1,Z1) For the GNSS satellite rectangular coordinate, the calculation method of the rotation matrix T is as follows:
Figure FDA0003176611340000102
wherein:
Figure FDA0003176611340000103
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:
Figure FDA0003176611340000104
step B404: calculating the azimuth angle A between the receiver and the satellite, wherein the calculation expression is as follows:
Figure FDA0003176611340000105
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: the curved surface equation in the step C501 and the ray equation in the step C3 are solved simultaneously to obtain a proportionality constant K corresponding to the intersection point of the ray and the longitude planeLThe calculation method is as follows:
Figure FDA0003176611340000111
step C503: will be constant of proportionality KLAnd re-substituting into the ray equation of the step C3 to obtain the coordinates of the intersection point of the ray and the longitude plane L, adding the coordinates into the intersection point sequence of the ray and the grid plane,
Figure FDA0003176611340000112
step C504: calculating the intersection point of the ray and the latitude surface, 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 W, and the equation of the curved surface is expressed as follows:
(X2+Y2)·tan2W-Z2＝0
step C505: c504 surface equation and C3 ray equation are solved simultaneously to obtain proportional constant K corresponding to the intersection point of the ray and the longitude planeWThe calculation method is as follows:
Figure FDA0003176611340000113
A＝(ΔX tan W)2+(ΔY tan W)2-(ΔZ)2
B＝2×((ΔX·Xr·tan W)2+(ΔY·Yr·tan W)2-ΔZ·Zr)
C＝(Xr·tan W)2+(Yr·tan W)2-(Zr)2
wherein: A. b, C are root coefficients of a quadratic equation;
step C506: will be constant of proportionality KWSubstituting the obtained coordinates of the intersection point of the ray and the latitude plane B into the ray equation of the step C3 again, 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:
(X2+Y2+Z2)-(Re+H)2＝0
in the formula: reIs the average earth radius;
step C508: c507, solving the curved surface equation in the step C507 and the ray equation in the step C3 simultaneously to obtain a proportionality constant K corresponding to the intersection point of the ray and the longitude surfaceHThe calculation method is as follows:
Figure FDA0003176611340000121
α＝(ΔX)2+(ΔY)2+(ΔZ)2
β＝2×(ΔX·Xr+ΔY·Y+ΔZ·Zr)
γ＝(Xr)2+(Yr)2+(Zr)2-(Re+H)2
wherein: the alpha, beta and gamma are root coefficients of a quadratic equation;
step C509: will be constant of proportionality KHSubstituting the obtained coordinates of the intersection point of the ray and the height plane H into the ray equation of the step C3 again, 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 distances delta R from all the intersection points to the receiver, wherein the calculation formula is as follows:
Figure FDA0003176611340000122
step C511: sequencing the distances delta R from the effective point sequences of all the cross points to the receiver according to a method from small to large, and respectively placing the coordinates of corresponding pointsInto a coordinate sequence (X)i,Yi,Zi) 1, 2.., n.
CN201911330615.XA 2019-12-20 2019-12-20 Ionosphere tomography method based on vertical measurement data constraint Active CN111273335B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911330615.XA CN111273335B (en) 2019-12-20 2019-12-20 Ionosphere tomography method based on vertical measurement data constraint

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911330615.XA CN111273335B (en) 2019-12-20 2019-12-20 Ionosphere tomography method based on vertical measurement data constraint

Publications (2)

Publication Number Publication Date
CN111273335A CN111273335A (en) 2020-06-12
CN111273335B true CN111273335B (en) 2021-09-17

Family

ID=70996934

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911330615.XA Active CN111273335B (en) 2019-12-20 2019-12-20 Ionosphere tomography method based on vertical measurement data constraint

Country Status (1)

Country Link
CN (1) CN111273335B (en)

Families Citing this family (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111505702B (en) * 2020-06-15 2023-08-11 华东交通大学 Ionosphere chromatography method based on vertical boundary truncated rays
CN111781594B (en) * 2020-06-23 2022-03-11 中国空间技术研究院 Ionosphere tomography method based on satellite-borne ice radar
CN111830582B (en) * 2020-07-07 2021-12-24 中国矿业大学 Average value constraint-based chromatographic inversion method
CN111985117B (en) * 2020-09-01 2022-03-22 华东交通大学 ACMART method suitable for GNSS ionosphere chromatography
CN112272067B (en) * 2020-10-15 2022-04-08 天津大学 Short wave broadcast frequency scheduling method based on multi-source data processing
CN112649899B (en) * 2020-11-19 2023-01-24 中国电波传播研究所（中国电子科技集团公司第二十二研究所） Global ionosphere data assimilation and forecasting method
CN112526617A (en) * 2020-11-19 2021-03-19 中国电波传播研究所（中国电子科技集团公司第二十二研究所） Ionospheric tomography system observation data simulation method based on multi-source satellite signals
CN112782727B (en) * 2020-11-20 2022-03-25 中国电波传播研究所（中国电子科技集团公司第二十二研究所） Station distribution design method based on observation weak area compensation
CN112882066A (en) * 2021-01-12 2021-06-01 中国人民解放军63770部队 Polar region ionized layer large-scale plasma structure image acquisition method
CN113093224B (en) * 2021-02-18 2022-08-02 北京航空航天大学 Edge-enhanced ionosphere chromatography method
CN113219405B (en) * 2021-04-15 2022-08-12 南京理工大学 Indoor dynamic multi-target passive positioning and quantity estimation method
CN113534213B (en) * 2021-07-26 2022-06-10 中国电子科技集团公司第五十四研究所 High-precision modeling and correcting method for atmospheric phase inconsistency in kilometer-level region
CN114384556B (en) * 2021-12-31 2024-05-14 中国电波传播研究所（中国电子科技集团公司第二十二研究所） Regional high-resolution ionosphere TEC map reconstruction method
CN117008154B (en) * 2023-08-03 2024-03-26 北京航空航天大学 Rapid ionosphere chromatography method based on relaxation factor reverse time decay function
CN118036665B (en) * 2024-04-12 2024-08-23 天津云遥宇航科技有限公司 Multi-intelligent meteorological model fusion method based on GNSS occultation data correction

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5428358A (en) * 1994-05-03 1995-06-27 The United States Of America As Represented By The Secretary Of The Navy Apparatus and method for ionospheric mapping
CN104536019A (en) * 2014-12-12 2015-04-22 中国电子科技集团公司第二十二研究所 GNSS ionized layer delay correction method based on ionized layer spatial correlation
CN105954764A (en) * 2016-04-27 2016-09-21 东南大学 GNSS computerized ionospheric tomography projection matrix acquisition method based on ellipsoid
CN108491616A (en) * 2018-03-19 2018-09-04 东南大学 A kind of vertical total electron content modeling method in ionosphere based on ellipsoid harmonic function theory
CN109613565A (en) * 2019-01-14 2019-04-12 中国人民解放军战略支援部队信息工程大学 Ionospheric Tomography method and system based on more constellation GNSS

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5428358A (en) * 1994-05-03 1995-06-27 The United States Of America As Represented By The Secretary Of The Navy Apparatus and method for ionospheric mapping
CN104536019A (en) * 2014-12-12 2015-04-22 中国电子科技集团公司第二十二研究所 GNSS ionized layer delay correction method based on ionized layer spatial correlation
CN105954764A (en) * 2016-04-27 2016-09-21 东南大学 GNSS computerized ionospheric tomography projection matrix acquisition method based on ellipsoid
CN108491616A (en) * 2018-03-19 2018-09-04 东南大学 A kind of vertical total electron content modeling method in ionosphere based on ellipsoid harmonic function theory
CN109613565A (en) * 2019-01-14 2019-04-12 中国人民解放军战略支援部队信息工程大学 Ionospheric Tomography method and system based on more constellation GNSS

Also Published As

Publication number Publication date
CN111273335A (en) 2020-06-12

Similar Documents

Publication Publication Date Title
CN111273335B (en) Ionosphere tomography method based on vertical measurement data constraint
Chen et al. Voxel-optimized regional water vapor tomography and comparison with radiosonde and numerical weather model
Rama Rao et al. On the validity of the ionospheric pierce point (IPP) altitude of 350 km in the Indian equatorial and low-latitude sector
Yao et al. Global ionospheric modeling based on multi-GNSS, satellite altimetry, and Formosat-3/COSMIC data
CN114690208B (en) Ionosphere three-dimensional electron density sparse chromatography method and device thereof
CN109613565A (en) Ionospheric Tomography method and system based on more constellation GNSS
CN111125609A (en) Ionized layer three-dimensional electron density reconstruction method based on double-exponential driving
Razin et al. Regional ionosphere modeling using spherical cap harmonics and empirical orthogonal functions over Iran
Seemala Estimation of ionospheric total electron content (TEC) from GNSS observations
Jin et al. 3-D ionospheric tomography from dense GNSS observations based on an improved two-step iterative algorithm
Li et al. Real‐time sensing of precipitable water vapor from beidou observations: hong kong and CMONOC networks
CN112526617A (en) Ionospheric tomography system observation data simulation method based on multi-source satellite signals
Haji-Aghajany et al. The effect of function-based and voxel-based tropospheric tomography techniques on the GNSS positioning accuracy
Ren et al. Electron density reconstruction by ionospheric tomography from the combination of GNSS and upcoming LEO constellations
Liu et al. Evaluation of HY-2A satellite-borne water vapor radiometer with shipborne GPS and GLONASS observations over the Indian Ocean
CN114384564B (en) Ionosphere tomography method based on multi-source data driving
Zhang et al. AN IMPROVED TROPOSPHERIC TOMOGRAPHY METHOD BASED ON THE DYNAMIC NODE PARAMETERIZED ALGORITHM.
CN113093224B (en) Edge-enhanced ionosphere chromatography method
Yang et al. GNSS water vapor tomography based on Kalman filter with optimized noise covariance
Yan et al. Correction of atmospheric delay error of airborne and spaceborne GNSS-R sea surface altimetry
Bahadur et al. Real-time single-frequency multi-GNSS positioning with ultra-rapid products
Lu et al. Virtual reference station-based computerized ionospheric tomography
Cheng et al. Global monitoring of geomagnetic storm-induced ionosphere anomalies using 3-D ionospheric modeling with multi-GNSS and COSMIC measurements
Chen et al. Global ionosphere modeling based on GNSS, satellite altimetry, radio occultation, and DORIS data considering ionospheric variation
Tang et al. Adaptive regularization method for 3-D GNSS ionospheric tomography based on the U-curve

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant