Functional basis ionosphere chromatography modeling method

Abstract

The invention discloses a functional basis ionosphere chromatography modeling method, which comprises the following steps of: acquiring longitude and latitude grid information of a research area, and calculating the electron content of a plasma layer at a preset height of a central point of each grid by utilizing an IRI-plas model. And fitting the data by adopting a third-order polynomial, and correcting the ionized layer delay observation value according to the data. Next, a non-ray penetration grid is identified based on GNSS data, a virtual ray is added, and an observation value thereof is calculated. And solving a model coefficient in combination with the corrected ionosphere delay and the virtual ray observation value, and constructing an initial ionosphere chromatography model based on the spherical harmonic function. And finally, calculating an ionosphere residual value, and performing residual correction on the initial model to obtain a corrected ionosphere chromatography model. The technical scheme of the invention has the advantages of high modeling precision and low calculation amount, and the corrected electron density value is closer to actual measurement data.

Description

Technical Field

[0001] The invention belongs to the technical field of GNSS atmosphere inversion, and in particular relates to a function-based ionospheric tomography modeling method. Background Art

[0002] The ionosphere is located in the range of 60km to 2000km above the ground. It is one of the important components of the atmosphere and a hot research topic in the field of satellite navigation. Under the action of solar radiation and cosmic rays, neutral atmospheric molecules in the atmosphere are ionized to produce free electrons that are constantly moving. The area where a large number of free electrons exist constitutes the ionosphere. The ionosphere is closely related to human production and life. On the one hand, the ionosphere absorbs ultraviolet radiation and cosmic high-energy ions, avoiding direct radiation from the sun and harmful rays to the earth, and also provides a propagation medium for long-distance radio transmission; on the other hand, electromagnetic waves propagating through the ionosphere will produce refraction and other effects, resulting in changes in signal speed and propagation path. When drastic changes occur, it may even cause communication interruption between satellites and ground equipment, greatly affecting the development of navigation, communication, aerospace and other fields. Therefore, studying the ionosphere is not only of great application value in ensuring the normal progress of aerospace activities, improving satellite navigation positioning accuracy and radio long-distance communication performance, but also of great scientific significance for studying the changing characteristics of the space atmospheric environment including the ionosphere.

[0003] Traditional ionospheric monitoring methods, such as GNSS satellites, occultation, altimeters, incoherent scattering radars, etc., can only obtain one-dimensional or two-dimensional ionospheric parameters, making it difficult to study the three-dimensional spatiotemporal changes of the ionosphere. Ionospheric tomography technology overcomes the limitations of the two-dimensional ionospheric model and provides the possibility for studying the three-dimensional ionospheric structure and disturbances. However, due to the influence of the geometric distribution of GNSS satellites and the density of ground observation stations, ionospheric tomography inversion faces ill-posed problems caused by insufficient observation data, uneven ray distribution, and limited observation angle.

[0004] However, the existing ionospheric tomography inversion technology still has the following shortcomings: 1) In order to solve the ill-posed problem, empirical models such as IRI2020 and NeQuick are often used as prior values ​​in areas where data is scarce. However, the accuracy of these empirical models is limited, which affects the accuracy of the inverted electron density; 2) In function-based ionospheric tomography, it is usually assumed that the plasmaspheric electron content is a fixed constant, such as 3TECU. However, the plasmaspheric electron content is not constant, and it varies with latitude, longitude, time and solar activity. Summary of the invention

[0005] The purpose of the present invention is to provide a function-based ionospheric tomography modeling method to solve the problems existing in the above-mentioned prior art.

[0006] To achieve the above object, the present invention provides a function-based ionospheric tomography modeling method, comprising:

[0007] Obtain the latitude and longitude grid division information of the study area, calculate the plasmasphere electron content at the preset height of each grid center point based on the IRI-plas model; fit the calculated plasmasphere electron content based on the third-order polynomial function, perform plasmasphere correction on the ionospheric delay observation value based on the fitted data, and obtain the corrected ionospheric delay observation value;

[0008] Based on GNSS observation data, grids without ray penetration are identified, virtual rays are added to the grids without ray penetration, and virtual ray observation values ​​are calculated;

[0009] Solving the model coefficients based on the corrected ionospheric delay observations and the virtual ray observations;

[0010] An initial ionospheric tomography model based on spherical harmonics is constructed based on the model coefficients, an ionospheric residual value between a corrected ionospheric delay observation value and a predicted value of the initial ionospheric tomography model is calculated, and a residual correction is performed on the initial ionospheric tomography model based on the ionospheric residual value to obtain a corrected ionospheric tomography model.

[0011] Optionally, the process of obtaining the latitude and longitude grid division information specifically includes:

[0012] The study area is divided horizontally to obtain a number of grids.

[0013] Optionally, the plasma layer correction is performed on the ionospheric delay observation value based on the fitted data, and the specific calculation formula is:

[0014]

[0015] In the formula, is the corrected ionospheric delay observation, STEC is the original ionospheric delay observation, VPEC is the plasmaspheric electron content, and Z′ is the zenith distance at the puncture point.

[0016] Optionally, the calculation process of the virtual ray observation value specifically includes:

[0017] The electron content of the plasma layer at the center point of the grid corresponding to the grid after adding the virtual ray is calculated at the preset height to obtain the virtual ray observation value.

[0018] Optionally, the model coefficients are solved based on the corrected ionospheric delay observation value and the virtual ray observation value, and the specific calculation formula is:

[0019]

[0020] In the formula, is the corrected ionospheric delay observation, VTEC is the added virtual ray observation, A1 and A2 are the matrices composed of the intercepts of the rays in the pixels, B is the coefficient matrix related to the horizontal and vertical direction functions, x is the column vector of the coefficients to be determined, and e is the observation noise.

[0021] Optionally, the calculation process of the ionospheric residual value specifically includes:

[0022]

[0023] In the formula, is the observed ionospheric delay, STEC' is the inverted ionospheric delay, STEC res is the ionospheric residual value.

[0024] Optionally, performing residual correction on the initial ionospheric tomography model based on the ionospheric residual value specifically includes:

[0025] Based on the ionospheric residual value, residual correction is performed on the ionospheric delay and electron density corresponding to each layer of each ray in the initial ionospheric tomography model to obtain a corrected ionospheric tomography model.

[0026] Optionally, the calculation formula for performing residual correction specifically includes:

[0027]

[0028] In the formula, STEC i is the ionospheric delay of the ray in layer i, STEC cor(i) is the ionospheric delay correction value of the ray in the i-th layer, D(i) is the intercept of the ray in the grid in the i-th layer, ELEC cor(i) is the electron density correction value of a certain ray in layer i, and is the final value of the ionospheric delay and the final value of the electron density of the ray in the i-th layer, ELEC i is the electron density value of the ray in the i-th layer.

[0029] The technical effects of the present invention are:

[0030] After the present invention uses a cubic polynomial to remove the plasma layer electron content, the ionospheric tomography result has higher accuracy than the direct modeling without eliminating the plasma layer; the initial model of the present invention uses a function-based ionospheric tomography model to significantly reduce the amount of calculation; the present invention adjusts the electron density in each voxel based on the STEC residual, and the corrected electron density value is closer to the actual measurement data. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0032] The drawings constituting a part of the present application are used to provide a further understanding of the present application. The illustrative embodiments and descriptions of the present application are used to explain the present application and do not constitute an improper limitation on the present application. In the drawings:

[0033] Figure 1 is a schematic diagram of an ionospheric tomography model in an embodiment of the present invention;

[0034] Figure 2 is the electron content of the plasma layer above 2000 km in the study area at a certain moment in the embodiment of the present invention;

[0035] Figure 3 Schematic diagram of plasma layer correction in an embodiment of the present invention;

[0036] Figure 4 is a ray distribution diagram observed at a certain moment in an embodiment of the present invention;

[0037] Figure 5 This is a diagram of the ionospheric tomography internal compliance accuracy result in an embodiment of the present invention;

[0038] Figure 6 This is a diagram of the ionospheric tomography external coincidence accuracy result in an embodiment of the present invention;

[0039] Figure 7 It is a schematic diagram of the comparison between the tomographic inversion result in the embodiment of the present invention and the digital altimeter hmF2, NmF2;

[0040] Figure 8 Schematic diagram of the comparison results of the electron density of tomographic inversion and the occultation data on August 22, 2021 in an embodiment of the present invention;

[0041] Fig. 9 Schematic diagram of the comparison results of the electron density of tomographic inversion and the occultation data on August 27, 2021 in an embodiment of the present invention;

[0042] Fig.10 This is a modeling flow chart in an embodiment of the present invention. DETAILED DESCRIPTION

[0043] Various exemplary embodiments of the present invention will now be described in detail. This detailed description should not be considered as limiting the present invention, but should be understood as a more detailed description of certain aspects, features, and embodiments of the present invention.

[0044] It should be understood that the terms described in the present invention are only for describing special embodiments and are not intended to limit the present invention. In addition, for the numerical range in the present invention, it should be understood that each intermediate value between the upper and lower limits of the scope is also specifically disclosed. Each smaller range between the intermediate value in any stated value or stated range and any other stated value or intermediate value in the described range is also included in the present invention. The upper and lower limits of these smaller ranges can be independently included or excluded in the scope.

[0045] Unless otherwise indicated, all technical and scientific terms used herein have the same meanings as those of ordinary skill in the art generally understood by the present invention. Although the present invention describes only preferred methods, any method similar or equivalent to that described herein may also be used in the implementation or testing of the present invention. All documents mentioned in this specification are incorporated by reference to disclose and describe the methods related to the documents. In the event of a conflict with any incorporated document, the content of this specification shall prevail.

[0046] It will be apparent to those skilled in the art that various modifications and variations may be made to the specific embodiments of the present invention description without departing from the scope or spirit of the present invention. Other embodiments derived from the present invention description will be apparent to those skilled in the art. The present application description and examples are exemplary only.

[0047] The words “include,” “including,” “have,” “contain,” etc. used in this article are open-ended terms, meaning including but not limited to.

[0048] It should be noted that, in the absence of conflict, the embodiments and features in the embodiments of the present application can be combined with each other. The present application will be described in detail below with reference to the accompanying drawings and in combination with the embodiments.

[0049] Embodiment 1

[0050] like Figure 1 - Fig.10As shown, in this embodiment, a function-based ionospheric tomography modeling method is provided, including: obtaining the latitude and longitude grid division information of the study area, calculating the plasma spheroid electron content at the center point of each grid at a preset height based on the IRI-plas model; fitting the calculated plasma spheroid electron content based on a third-order polynomial function, performing plasma spheroid correction on the ionospheric delay observation value based on the fitted data, and obtaining a corrected ionospheric delay observation value; identifying grids without ray penetration based on GNSS observation data, adding virtual rays to the grids without ray penetration, and calculating virtual ray observation values; solving model coefficients based on the corrected ionospheric delay observation values ​​and the virtual ray observation values; constructing an initial ionospheric tomography model based on spherical harmonic functions based on the model coefficients, calculating the ionospheric residual value between the corrected ionospheric delay observation value and the initial ionospheric tomography model prediction value, and performing residual correction on the initial ionospheric tomography model based on the ionospheric residual value to obtain a corrected ionospheric tomography model.

[0051] This embodiment, through in-depth research on high-precision function-based ionospheric tomography modeling, can not only effectively reduce the impact of ionospheric errors on the navigation system and improve the accuracy and reliability of positioning, but also help provide more accurate ionospheric information for multimodal navigation systems, and promote the integration of satellite navigation and other navigation technologies, such as inertial navigation, visual navigation, etc., to further improve the robustness and global applicability of the navigation system. The ionosphere may undergo abnormal changes due to factors such as solar activity, which may have a serious impact on the navigation system. High-precision ionospheric modeling research helps to identify and respond to these abnormal events in a timely manner, and also provides valuable data support for areas such as early warning of catastrophic geological events and space weather forecasting.

[0052] STEC is the integral of the ionospheric electron density from the receiver position to the satellite position and is expressed as follows:

[0053]

[0054] Among them, r, s are the positions of the receiver and satellite respectively, Ne is the electron density, λ, φ, h are the longitude, latitude and altitude respectively.

[0055] In this embodiment, the distribution of electron density in the horizontal direction is represented by spherical harmonic functions, and the vertical direction is represented by three EOFs. Then the electron density at any point above the surface can be expressed as:

[0056]

[0057] In ionospheric tomography, the study area is usually divided into pixels of equal size, assuming that the electron density in each pixel remains constant. The STEC value on each ray can be considered to be equal to the sum of the product of the intercept of the ray in the corresponding pixel it passes through and the electron density in the pixel, such as Figure 1 Assume that the total number of pixels in the ionosphere region to be inverted after gridding is , there are a total of observation rays passing through the region to be inverted, and the intercept of the ray in each pixel is a i (i=1,2,3,…n), the electron density of each pixel is x i (i=1,2,3,…n), its TEC value is STEC i , which can be expressed as:

[0058]

[0059] Among them, e i is the observation noise.

[0060] Written in matrix form:

[0061] STEC m×1 =A m×n X n×1 +e m×1 #(4)

[0062] Among them, A is the matrix composed of the intercepts of each ray in the pixel, and X is the electron density matrix. The electron density is expressed as a mathematical function in the horizontal and vertical directions, and the above formula can be written as:

[0063] STEC m×1 =A m×n B n×u x u×1 +e m×1 #(5)

[0064] Among them, B is the coefficient matrix related to the horizontal and vertical direction functions, m is the number of rays, n is the number of pixels in the tomography area, u is the number of coefficients to be determined, and x is the column vector of coefficients to be determined.

[0065] (1) Plasmasphere electron content correction: Under normal circumstances, the plasmasphere electron content is considered to be a fixed constant 3TECU, and this constant is divided by the sine of the altitude angle as the plasma electron content correction on the oblique path. However, the plasmasphere electron content changes with the latitude, longitude, time and solar activity. Therefore, this embodiment divides the study area into 10*13 grids horizontally (10 in the latitude direction and 13 in the longitude direction), and uses the IRI-plas 2020 model of assimilation CODE-GIM to calculate the plasmasphere electron content in the zenith direction above 2000km at the center of each grid, recorded as VPEC. IRI-plas extends the International Reference Ionosphere (IRI) to the plasmasphere, and the altitude can reach the GPS orbit altitude of 20,200 kilometers. IRI-plas can input external information such as the critical frequency of the F2 layer (foF2), the maximum ionization altitude (hmF2) and the total electron content (TEC) to update the ionospheric state. Studies have shown that after inputting external TEC, the TEC and foF2 of IRI-plas will be more accurate. The daily VPEC data are fitted with a third-order polynomial function, as shown in the figure.

[0066] After obtaining the VPEC at a certain moment, the plasmasphere correction can be performed on the ionospheric delay observation value:

[0067]

[0068] in, is the corrected ionospheric delay observation, STEC is the original ionospheric delay observation, VPEC is the plasma spheric electron content, and Z′ is the zenith distance at the puncture point. Figure 3 shown.

[0069] This embodiment introduces the IRI-plas model into ionospheric tomography for the first time, and analyzes the spatiotemporal characteristics of electron content above 2000 kilometers based on the IRI-plas model. Plasma layer correction is performed before ionospheric tomography modeling to eliminate the contribution of the plasma layer to the total electron content and improve the modeling accuracy.

[0070] (2) Assimilation of GIM’s IRI-plas virtual ray enhancement: Due to the limitations of the location distribution of GNSS observation stations and satellite observation angles, GNSS rays often cannot penetrate all grids in the area to be inverted. At the same time, the electron density distribution in the ionosphere is constantly changing, and ionospheric tomography is often based on the assumption that the ionosphere remains unchanged over a period of time. Therefore, the GNSS observation time interval for ionospheric tomography cannot be too long. This results in the electron density of many grids that are not penetrated by rays in ionospheric tomography not being effectively corrected, thus affecting the modeling accuracy. Figure 4Shown is the distribution of observed rays for 30 minutes at UT14:00 on August 20, 2021.

[0071] Depend on Figure 4 It can be seen that even after accumulating 30 minutes of observation data, there are still a large number of grids that have no ray penetration. The reason for the grids not having ray penetration is that there are a lack of observation stations in the areas near these grids. Therefore, the distribution of observation stations is a major factor affecting the accuracy of ionospheric tomography modeling. If some virtual rays are added in areas where there are no observation stations on the ground, the grids without ray penetration can be covered.

[0072] The vertical total electron content (VTEC) in the region of 80-2000 km in the absence of observational data is calculated by integrating the IRI-plas 2020 model with GIM, and compared with the corrected ionospheric delay observations. If we solve them together, formula (5) can be expressed as:

[0073]

[0074] in, is the corrected ionospheric delay observation, VTEC is the added virtual ray observation, A1 and A2 are matrices of the intercepts of the rays in the pixels, B is the coefficient matrix related to the horizontal and vertical functions, and x is the column vector of the coefficients to be calculated. The model coefficient x is solved using weighted least squares, where The weight ratio of the observed value and the VTEC observed value is 2:1, and e is the observation noise.

[0075] By setting up virtual observatories in areas lacking ground observatories and adding a series of virtual rays, on the one hand, the voxels that were originally not penetrated by rays can be effectively corrected, and on the other hand, the number of voxels penetrated by cross-rays can be increased, alleviating the rank deficiency and ill-conditioned problems of the tomographic equation. This method can effectively invert the three-dimensional ionospheric electron density, and the inversion accuracy is significantly improved in areas lacking rays.

[0076] (3) Correction of residual errors in the tomographic model: Since the spherical harmonic function is a relatively smooth function model, it can only describe relatively gentle changes in ionospheric delay. If the ionospheric delay changes dramatically or the amplitude is large, such as during a magnetic storm, the value of the oblique ionospheric delay STEC may be several times that of a calm period. At this time, it is difficult to accurately describe the change trend of STEC using only a spherical harmonic function. Therefore, an initial function base analysis model is first established, and then the model is corrected. The correction method is as follows:

[0077] is the observed value of ionospheric delay, STEC' is the inverted ionospheric delay, STEC res is the ionospheric residual value, see formula (7).

[0078]

[0079] According to formulas (8), (9), (10), (11), the ionospheric delay and electron density of each layer of each ray are corrected in turn:

[0080]

[0081] Among them, STEC i is the ionospheric delay of a ray in layer i, STEC cor(i) is the ionospheric delay correction value of a ray in the i-th layer, D(i) is the intercept of a ray in the grid of the i-th layer, ELEC cor(i) is the electron density correction value of a certain ray in layer i, and is the final value of ionospheric delay and electron density of a certain ray in the i-th layer.

[0082] After the cubic polynomial is used to remove the electron content of the plasma layer in this embodiment, the ionospheric tomography result is more accurate than the direct modeling without eliminating the plasma layer; the initial model in this embodiment adopts the function-based ionospheric tomography model to significantly reduce the amount of calculation; the electron density in each voxel is adjusted based on the STEC residual in this embodiment, and the corrected electron density value is closer to the actual measurement data.

[0083] A function-based ionospheric tomography modeling method with additional plasma correction and virtual ray enhancement proposed in this embodiment is used to achieve regional high-precision ionospheric modeling. The results show that the new method can significantly improve the accuracy of ionospheric regional modeling. Comparison of the modeling results with the ionospheric altimeter and COSMIC2 occultation data shows that this method can effectively reconstruct the three-dimensional electron density and improve the inversion accuracy in areas with limited ray penetration.

[0084] (1) Internal and external accuracy of the ionospheric tomography model: The root mean square error of 48 inversion periods in one day is taken as the overall internal accuracy evaluation index for one day. The internal accuracy results are as follows: Figure 5 The black columns in the figure represent the internal accuracy of the ionospheric tomography model established using only spherical harmonics in the horizontal direction and empirical orthogonal decomposition functions in the vertical direction. The red columns represent the internal accuracy after the ionospheric residual correction is performed on the basis of the previous step. If only the first step is used, the average internal accuracy for 30 days is 3.38 TECU. If the second step of ionospheric residual correction is performed on the basis of the first step, the average internal accuracy for 30 days is 1.34 TECU, and the accuracy of the tomography model is significantly improved.

[0085] In order to objectively evaluate the overall accuracy of the inverted ionospheric tomography model, the inversion results are compared with the observation data that does not participate in ionospheric tomography to evaluate its external accuracy. The root mean square error of all time periods in a day is counted and the average value is taken as the evaluation index of the overall external accuracy of the day. This embodiment counts the average value of the external accuracy for 30 days from August 1 to 30, 2021, as shown in Figure 6 Shown

[0086] Figure 6 The blue column in the middle represents the external accuracy of the ionospheric tomography model established using only the horizontal spherical harmonics and the vertical empirical orthogonal decomposition function. The green column represents the external accuracy after the ionospheric residual correction is performed on the basis of the previous step. If only the first step is used, the average external accuracy for 30 days is 2.83TECU. If the second step of ionospheric residual correction is performed on the basis of the first step, the average external accuracy for 30 days is 1.55TECU, and the external accuracy within 30 days fluctuates slightly, all below 2TECU.

[0087] (2) Digital altimeter verification: Since the peak density (NmF2) and peak height (hmF2) of the ionospheric F2 layer are important parameters reflecting the ionospheric electron density profile, if the inverted electron density NmF2 and hmF2 are close to the true value, it can be proved that the inversion result is accurate enough and the algorithm is relatively reliable. Therefore, it can be used as an important reference indicator for evaluating the accuracy of ionospheric tomography inversion results. The absolute errors between NmF2 and hmF2 obtained by ionospheric tomography inversion and the data of the digital altimeter at the same time are counted, and the daily average root mean square error RMSE of the absolute errors in all time periods within a day is calculated. The results are as follows: Figure 7 shown.

[0088] Depend on Figure 7 It can be seen that after using the second step of ionospheric residual correction, the root mean square error of the peak density NmF2 from August 1 to 30, 2021 is significantly reduced, the peak density is closer to the true value, and the change in the peak height hmF2 is not obvious. This proves that the model accuracy is improved after using residual correction.

[0089] (3) COSMIC2 occultation data verification: The COSMIC series of atmospheric detection satellites were jointly developed by the United States and Taiwan, China. They are equipped with GNSS receivers, ionospheric photometers and other instruments, which can provide parameters such as atmospheric humidity, refractive index and electron density profiles. They are mainly used for meteorological environment monitoring, ionospheric detection and weather forecasting. IonPrfs is a level-2 atmospheric inversion product of COSMIC 1 / 2, which mainly provides electron density data obtained by inverting the occultation event using the Abel transformation method. COSMIC2 occultation data has the characteristics of large data volume, high accuracy, all-weather, and not affected by the distribution of observation sites. The actual measured data can better reflect the small-scale daily changes in the ionosphere. Therefore, the electron density inverted by ionospheric tomography can be compared with the occultation data to verify the inversion results.

[0090] On August 22, 2021, geomagnetic activity was relatively calm. There were 4721 occultation events on that day. After removing the inaccurate electron density files, there were still 4506 occultation events, of which 38 occurred over the study area. The electron density profile of each event was compared with the electron density inverted at the same time and location. The difference between the inverted electron density and the occultation electron density at each altitude from 80 to 560 km was compared at an interval of 20 km, and the root mean square error of all events at the same altitude was calculated. The results are as follows: Figure 8 As shown in Figure 2, at altitudes above 200 km, the RMS error after the second-step residual correction is reduced, and the results above 300 km are more obvious, indicating that the electron density inverted using the two-step method is more accurate.

[0091] On August 27, 2021, the geomagnetic activity was relatively active. There were 4893 occultation events on this day. After removing the inaccurate electron density files, there were still 4650 occultation events, of which 38 occurred over the study area. The root mean square error of this day was calculated according to the previous method, and the results are as follows Fig. 9 The root mean square error after the second step residual correction is reduced, and the results above 300 km are more obvious, indicating that whether the ionosphere is relatively calm or relatively active, the electron density inverted using this method is closer to the true value, and the model has good robustness.

[0092] The above is only a preferred specific implementation of the present application, but the protection scope of the present application is not limited thereto. Any changes or substitutions that can be easily thought of by a person skilled in the art within the technical scope disclosed in the present application should be included in the protection scope of the present application. Therefore, the protection scope of the present application should be based on the protection scope of the claims.

Claims

1. A function-based ionospheric tomography modeling method, characterized in that: include: Obtain the latitude and longitude grid division information of the study area, calculate the plasmasphere electron content at the preset height of each grid center point based on the IRI-plas model; fit the calculated plasmasphere electron content based on the third-order polynomial function, perform plasmasphere correction on the ionospheric delay observation value based on the fitted data, and obtain the corrected ionospheric delay observation value; Based on GNSS observation data, grids without ray penetration are identified, virtual rays are added to the grids without ray penetration, and virtual ray observation values ​​are calculated; Solving the model coefficients based on the corrected ionospheric delay observations and the virtual ray observations; An initial ionospheric tomography model based on spherical harmonics is constructed based on the model coefficients, an ionospheric residual value between a corrected ionospheric delay observation value and a predicted value of the initial ionospheric tomography model is calculated, and a residual correction is performed on the initial ionospheric tomography model based on the ionospheric residual value to obtain a corrected ionospheric tomography model.

2. A function-based ionospheric tomography modeling method according to claim 1, characterized in that: The process of obtaining the latitude and longitude grid division information specifically includes: The study area is divided horizontally to obtain a number of grids.

3. A function-based ionospheric tomography modeling method according to claim 1, characterized in that: The plasma layer correction is performed on the ionospheric delay observation value based on the fitted data, and the specific calculation formula is: In the formula, is the corrected ionospheric delay observation, STEC is the original ionospheric delay observation, VPEC is the plasmaspheric electron content, and Z' is the zenith distance at the puncture point.

4. The function-based ionospheric tomography modeling method according to claim 1, characterized in that: The calculation process of the virtual ray observation value specifically includes: The electron content of the plasma layer at the center point of the grid corresponding to the grid after adding the virtual ray is calculated at the preset height to obtain the virtual ray observation value.

5. The function-based ionospheric tomography modeling method according to claim 1, characterized in that: The model coefficients are solved based on the corrected ionospheric delay observation value and the virtual ray observation value, and the specific calculation formula is: In the formula, is the corrected ionospheric delay observation, VTEC is the added virtual ray observation, A1 and A2 are the matrices composed of the intercepts of the rays in the pixels, B is the coefficient matrix related to the horizontal and vertical direction functions, x is the column vector of the coefficients to be determined, and e is the observation noise.

6. A function-based ionospheric tomography modeling method according to claim 1, characterized in that: The calculation process of the ionospheric residual value specifically includes: In the formula, is the observed ionospheric delay, STEC' is the inverted ionospheric delay, STEC res is the ionospheric residual value.

7. The function-based ionospheric tomography modeling method according to claim 1, characterized in that: The performing residual correction on the initial ionospheric tomography model based on the ionospheric residual value specifically includes: Based on the ionospheric residual value, residual correction is performed on the ionospheric delay and electron density corresponding to each layer of each ray in the initial ionospheric tomography model to obtain a corrected ionospheric tomography model.

8. A function-based ionospheric tomography modeling method according to claim 7, characterized in that: The calculation formula for performing residual correction specifically includes: In the formula, STEC i is the ionospheric delay of the ray in layer i, STEC cor(i) is the ionospheric delay correction value of the ray in the i-th layer, D(i) is the intercept of the ray in the grid in the i-th layer, ELEC cor(i) is the electron density correction value of the ray in layer i, and is the final value of the ionospheric delay and the final value of the electron density of the ray in the i-th layer, ELEC i is the electron density value of the ray in the i-th layer.

Images