A method for extracting three-dimensional ionospheric disturbances

Abstract

The application discloses a three-dimensional ionospheric disturbance extraction method, which comprises the following steps: acquiring double-frequency observation data of a research area from a GNSS receiving station; then, calculating STEC from the GNSS double-frequency signal and converting the STEC into dSTEC to highlight the ionospheric disturbance signal; secondly, constructing an observation matrix containing prior ionospheric model information and utilizing an improved three-dimensional tomography algorithm to iteratively correct and gradually approach the real electron density distribution; thirdly, in the iteration process, adopting inverse distance weighted interpolation to correct the pixels not being passed through, so as to ensure the overall smoothness; finally, inputting the dSTEC data into the improved three-dimensional tomography model, obtaining the three-dimensional ionospheric electron density distribution through iterative solution, and comparing the three-dimensional ionospheric electron density distribution with the actual observation data to verify the accuracy of the method. The method can effectively extract the three-dimensional structure of the ionospheric disturbance, improve the reconstruction accuracy and resolution, and provide an effective tool for studying the response of the ionosphere to events such as earthquakes.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of ionospheric electron density three-dimensional monitoring, and particularly relates to a three-dimensional ionospheric disturbance extraction method, which is suitable for high-precision three-dimensional reconstruction of ionospheric disturbances. BACKGROUND

[0002] The ionosphere is an atmospheric layer with an altitude of 60 to 1000 kilometers above the earth, and is an important part of the earth's sphere. The dynamic changes of its electron density distribution not only reflect the coupling process of the magnetosphere-ionosphere-thermosphere system, but also have a direct and significant impact on satellite communication, navigation and positioning, and space weather research. Therefore, high-precision and high-resolution reconstruction of the ionospheric electron density distribution is crucial for deepening the study of ionospheric physical mechanisms, ensuring the stability of radio communication, improving the accuracy of satellite navigation, and optimizing the ability of space weather monitoring.

[0003] The emergence of the Global Navigation Satellite System (GNSS) has provided a powerful tool for ionospheric monitoring. GNSS provides Slant Total Electron Content (STEC) observation data, making it possible to use 3DCIT technology based on GNSS for three-dimensional computerized ionospheric tomography. 3DCIT technology can effectively utilize GNSS data to reconstruct the three-dimensional electron density distribution in the ionosphere; however, existing ionospheric monitoring methods have many limitations: (1) When dealing with ionospheric disturbances caused by events such as earthquakes, they are often affected by background wind, magnetic field and nonlinear phenomena, resulting in unstable and low-precision reconstruction results. (2) When dealing with ionospheric disturbances, it is usually assumed that the ionosphere remains constant over a period of time, which limits its ability to monitor rapidly changing ionospheric disturbances. (3) There are deficiencies in vertical resolution and horizontal coverage, making it difficult to accurately describe the three-dimensional structure of ionospheric disturbances. (4) When dealing with pixels that are not traversed by GNSS rays, the Laplace operator is often used for smoothing, which may cause over-smoothing and mask the true characteristics of ionospheric disturbances.

[0004] During earthquakes, Rayleigh waves at the Earth's surface, sound waves generated by ruptures, and gravity waves triggered by tsunamis all induce ionospheric electron density disturbances, known as seismo-traveling ionospheric disturbances (STIDs). The propagation characteristics of STIDs are influenced by various factors, including background winds, magnetic fields, and nonlinear effects, making it difficult for numerical models to accurately simulate their propagation process. While existing research has revealed the horizontal velocity and vertical distribution of STIDs to some extent, detailed three-dimensional analyses of their propagation and attenuation characteristics at different altitudes are still lacking. Summary of the Invention

[0005] The purpose of this invention is to provide a three-dimensional ionospheric disturbance extraction method to solve many problems existing in the current ionospheric monitoring technology mentioned in the background art.

[0006] To achieve the above objectives, the present invention provides a method for extracting three-dimensional ionospheric perturbations, comprising the following steps:

[0007] S1. Obtain dual-frequency observation data of the study area from the GNSS receiving station;

[0008] S2. Calculate the total electron content (STEC) of the ionosphere using the carrier phase smoothed pseudorange value of the GNSS dual-frequency signal; subtract the moving average from the STEC time series to obtain the detrended total electron content (STEC), i.e., the dSTEC measurement value.

[0009] S3. The ionosphere is divided into pixels. When dividing the pixels, a multi-resolution grid is used in the vertical direction, and the division resolution in the peak height range of the F2 layer is higher than the division resolution in other height ranges. The electron density value is initialized using a priori ionosphere model, and the ionosphere is iteratively corrected using synchronous algebraic reconstruction technology. Then, for pixels that have not been penetrated by GNSS rays, inverse distance weighted interpolation is used for correction. Finally, an improved three-dimensional tomographic model is constructed.

[0010] S4. Input the dSTEC measurement value obtained in step S2 into the improved three-dimensional tomography model constructed in step S3, and obtain the three-dimensional ionospheric perturbation through iterative calculation.

[0011] Furthermore, in step S2, the formula for calculating the total electron content (STEC) of the strabismus is:

[0012]

[0013] In the formula, and These are the two carrier signal frequencies of the GPS satellites; and For and the carrier phase smoothed pseudo-range value of the frequency signal; and respectively the differential code bias of the GNSS receiver and the GPS satellite.

[0014] Further, in step S2, the calculation formula of dSTEC is:

[0015]

[0016] wherein m is the number of STEC measurements of the GPS satellite; n is the total number of pixels of the reconstruction area; is the observation matrix; represents the electron density of STIDs; ε is the STEC measurement noise error.

[0017] Further, in each iteration of step S3, the electron density of each pixel is gradually corrected according to the intercept of the ray path and the current electron density estimation value; the iteration expression is:

[0018]

[0019]

[0020]

[0021] wherein, is the ray-corrected electron density value of the jth pixel after the k+1th iteration, P is the number of pixels through which the ith ray passes, and 0≤P≤n; λ is a relaxation parameter, 0≤λ≤1; is the total ionospheric electron content of the ith ray; is the intercept of the ith ray in the ath pixel, and 0≤a≤n; is the electron density value of the ath pixel after the kth iteration; Δ is the correction amount of the ith ray path; W is the weight of the jth pixel in the TEC correction distribution of the ith ray; represents the total number of pixels through which the ray passes on a specific ray path, i.e. the number of pixels in the calculation of the correction amount.

[0022] Further, in step S2, the correction formula for interpolating the pixels not through the ray is:

[0023]

[0024] wherein, represents the electron density value of the pixel to be interpolated, represents the electron density value of the pixel that has been corrected; denotes the distance between the pixel to be interpolated and the already corrected pixel; denotes the number of already corrected pixels involved in the interpolation.

[0025] Compared with the prior art, the present application has the following beneficial effects:

[0026] 1. The three-dimensional ionospheric disturbance extraction method of the present application uses the detrended slant total electron content (dSTEC) instead of the traditional STEC, and combines the improved three-dimensional tomography technology to effectively remove the interference of background changes and accurately capture the ionospheric disturbance signals caused by events such as earthquakes. At the same time, the dSTEC has been processed by detrending and other processes on the basis of the STEC, and has higher precision than the original STEC, which can reach 0.5TECU, which is 1-2 orders of magnitude higher than the traditional STEC, and is sufficient to capture small amplitude ionospheric disturbances caused by earthquakes; the high precision and high resolution reconstruction effect of the present application provides strong technical support for in-depth research on the response mechanism of the ionosphere to earthquakes.

[0027] 2. The three-dimensional ionospheric disturbance extraction method of the present application uses the improved constrained synchronous iterative reconstruction technology algorithm, which introduces the product of the in-pixel intercept and the electron density to calculate the electron density correction amount, so that the correction amount at different altitudes is more reasonable, and uses inverse distance weighted interpolation to correct the pixels not traversed by GNSS rays, which can effectively avoid over-smoothing and effectively preserve the true characteristics of the ionospheric disturbance. The multi-resolution grid design is adopted in the vertical direction, which especially improves the reconstruction accuracy near the F2 layer peak height (HmF2), and solves the problem of insufficient accuracy in the key area of the traditional method.

[0028] In summary, the present application not only improves the accuracy and reliability of ionospheric disturbance extraction, but also expands its application range in space weather research and ionospheric monitoring.

[0029] In addition to the purposes, features and advantages described above, the present application has other purposes, features and advantages. The present application will be further described in detail below with reference to the accompanying drawings. BRIEF DESCRIPTION OF DRAWINGS

[0030] The accompanying drawings are used to provide further understanding of the embodiments of the present application, and constitute a part of the specification, and are used together with the following specific embodiments to explain the embodiments of the present application, but do not constitute a limitation on the embodiments of the present application. In the drawings:

[0031] Figure 1 is a flow chart of an embodiment of the three-dimensional ionospheric disturbance extraction method provided by the present application;

[0032] Figure 2 is a ray path correction assignment diagram of a certain latitude plane in the improved three-dimensional tomography technology involved in the present application. DETAILED DESCRIPTION

[0033] The present application will be described in detail below with reference to the embodiments shown in the drawings, but it should be noted that these embodiments are not limiting of the present application, and equivalent transformations or substitutions of function, method, or structure made by those of ordinary skill in the art based on these embodiments are within the scope of the present application.

[0034] Reference should be made to Figure 1 The present embodiment provides a three-dimensional ionospheric disturbance extraction method, which uses detrended slant total electron content (dSTEC) data instead of traditional slant total electron content (STEC) data, and combines an improved three-dimensional tomography method to highlight ionospheric disturbance signals by calculating STEC from GNSS dual-frequency signals and converting it to dSTEC, constructing an observation matrix containing prior ionospheric model information, and using an improved three-dimensional tomography algorithm for iterative correction to gradually approach the real electron density distribution. In the iteration process, the inverse distance weighted interpolation is used to correct the pixels not passed through to ensure the overall smoothness. Finally, the dSTEC data is input into the improved three-dimensional tomography model, and the three-dimensional ionospheric electron density distribution is obtained by iterative solution, and the accuracy of the method is verified by comparison with the actual observation data. This method can effectively extract the three-dimensional structure of ionospheric disturbances, improve the reconstruction accuracy and resolution, and provide an effective tool for studying the response of the ionosphere to events such as earthquakes. The method includes the following specific steps:

[0035] 1. Obtain dual-frequency observation data of the study area from GNSS receiving stations;

[0036] 2. When the GNSS signal passes through the ionosphere, the free electrons in the ionosphere affect the propagation speed of the signal, causing the signal to delay. This delay is inversely proportional to the square of the signal frequency. By measuring the pseudorange and carrier phase of the GNSS signal and using dual-frequency observation data, the effect of ionospheric delay can be eliminated to estimate the electron content in the ionosphere. Specifically, STEC reflects the electron density integral when the satellite signal passes through the ionosphere, and its calculation formula is:

[0037]

[0038] In formula (1), and are the frequencies of the two carrier signals of the GPS satellite, = 1575.42 MHz, = 1227.60 MHz; and are and the carrier phase smoothed pseudorange values of the frequency signals; and Differential Code Bias (DCB) of GNSS receiver and GPS satellite, respectively.

[0039] It is noted that DCB is one of the important error sources for estimating precise TEC using GNSS data. The International GNSS Service Analysis Center has routinely provided DCB estimates for GNSS satellites and IGS ground receivers, but not for regional and local network receivers. Moreover, it is not accurate to assume that the DCB values of GNSS satellites or receivers remain constant within 1 day or 1 month. In the present embodiment, M_DCB software is used to estimate the DCB of GNSS satellites and receivers with time intervals from hours to days; the M_DCB software is developed in Matlab (version: 2010b). The input of M_DCB is GPS RINEX observation files and precise ephemeris, which here default to contain P1 and P2 observations. The output is the DCB estimates of satellites and receivers.

[0040] 3. Next, the moving average is subtracted from the STEC time series to obtain the detrended slant total electron content STEC, i.e. dSTEC measurements, to remove background variations (e.g. diurnal variations and extreme ultraviolet flux variations) and highlight ionospheric disturbance signals. Specifically:

[0041]

[0042] where t is the moving average window; preferably, t is 10 min.

[0043] In actual calculation, the reconstruction region is divided into a plurality of small pixels, and it is assumed that the electron density in one pixel is the same. Therefore, dSTEC can be expressed as:

[0044]

[0045] where m is the number of STEC measurements of the GPS satellite, and n is the total number of pixels of the reconstruction region; is an observation matrix constructed based on the number of STEC measurements m and the total number of pixels n of the reconstruction region, where the element represents the intercept of the ith ray in the jth pixel; represents the electron density of the STIDs, and ε is the STEC measurement noise error.

[0046] 4、Pixel division of ionosphere, in order to improve the reconstruction accuracy near the height of HmF2 (peak height of F2 layer), multi-resolution grid is used in the vertical direction when dividing pixels, the resolution is 10 km in the range of peak height of F2 layer (200~420 km), and the resolution is 50 km in other heights.

[0047] 5、Using prior ionospheric model to initialize electron density value, iterative correction of SART (Synchronous Algebraic Reconstruction Technique) is carried out; in each iteration, the electron density of each pixel is gradually corrected according to the intercept of ray path and the current electron density estimate value; the iteration expression is as follows:

[0048]

[0049]

[0050]

[0051] In the formula, is the ray-corrected electron density value of the jth pixel after the k+1th iteration, is the number of pixels crossed by the ith ray, and 0≤P≤n; λ is a relaxation parameter, 0≤λ≤1; in the embodiment of the present application, λ is set to 0.2. is the total ionospheric electron content of the ith ray; is the intercept of the ith ray in the ath pixel, and 0≤a≤n; is the electron density value of the ath pixel after the kth iteration; Δ is the correction amount of the ith ray path; W is the weight of the jth pixel in the TEC correction distribution of the ith ray; represents the total number of pixels crossed by the ray on a certain ray path, that is, the number of pixels in the calculation of the correction amount. In some cases, and may be equal, that is, all the pixels crossed by the ray are used to calculate the correction amount Δ. But in other cases, may be greater than , that is, more pixels may be considered in the calculation of the correction amount Δ, which may include some pixels not directly crossed by the ray but indirectly affected.

[0052] As shown in Figure 2 , a ray path correction distribution diagram of a longitude-latitude plane is shown. If A i,j =A i,k, according to the traditional method of distributing the correction according to the ray intercept, the correction of the jth pixel and the kth pixel is the same. In fact, the correction should be different because the electron density and its variation in pixels j and k are very different. Therefore, this correction assignment process is unreasonable. In the present application, the denominator part of the weight W is changed from the traditional to , that is, the traditional correction according to the ray intercept is changed to the multiplication according to the intercept and the electron density, which can effectively avoid the over-correction problem of the low-height area (small electron density).

[0053] In the embodiment of the present application, when the maximum value of the difference between the current iteration and the last iteration result is less than 0.03 ×10 11 el / m 3 or the iteration number is greater than 20 times, the iteration will stop.

[0054] 6. For the pixels not penetrated by the GNSS ray, inverse distance weighted interpolation (IDW) is used for correction. This method assumes that the values of adjacent pixels are more similar, thereby ensuring the smoothness of the reconstructed area. The correction formula is as follows:

[0055]

[0056] wherein, represents the electron density value of the pixel to be interpolated, represents the electron density value of the pixel that has been corrected; represents the distance between the pixel to be interpolated and the pixel that has been corrected; represents the number of pixels that have been corrected and participate in the interpolation. The electron density value of the ray-corrected pixel is not affected by the inverse distance weighting, and only the pixels not penetrated by the ray are corrected.

[0057] 7. The obtained dSTEC measurement value is input into the constructed improved three-dimensional tomographic model, and the three-dimensional ionospheric disturbance is obtained through iteration.

[0058] The above only describes the preferred embodiments of the present application and is not used to limit the present application. For those skilled in the art, the present application can have various modifications and changes. Any modification, equivalent replacement, improvement, etc. made within the spirit and principles of the present application shall be included in the protection scope of the present application.

Claims

1. A method for extracting three-dimensional ionospheric perturbations, characterized in that, Includes the following steps: S1. Obtain dual-frequency observation data of the study area from the GNSS receiving station; S2. Calculate the outward total electron content (STEC) of the ionosphere using the carrier phase smoothing pseudorange value of the GNSS dual-frequency signal. Subtracting the moving average from the STEC time series yields the detrended slanted total electron content STEC, also known as dSTEC measurement. S3. The ionosphere is divided into pixels. A multi-resolution grid is used in the vertical direction for pixel division, with higher resolution within the peak height range of layer F2 than in other height ranges. The electron density value is initialized using a priori ionospheric model, and synchronous algebraic reconstruction technology is used for iterative correction of the ionosphere. Pixels not penetrated by GNSS rays are then corrected using inverse distance weighted interpolation. Finally, an improved three-dimensional tomographic model is constructed. Wherein: In each iteration, the electron density of each pixel is progressively adjusted based on the ray path intercept and the current electron density estimate; the iterative expression is: In the formula, is the ray-corrected electron density value of the j-th pixel after the (k+1)-th iteration; P is the number of pixels traversed by the i-th ray, and 0≤P≤n; λ is the relaxation parameter, 0≤λ≤1; It is the total ionospheric electron content of the i-th ray; It is the intercept of the i-th ray in the a-th pixel, and 0≤a≤n; Δ is the electron density value of the a-th pixel after the k-th iteration; Δ is the correction amount for the i-th ray path; W is the weight of the j-th pixel in the TEC correction allocation for the i-th ray. This represents the total number of pixels that are traversed by the ray along the ray path, i.e., the number of pixels used in calculating the correction amount; The correction formula for interpolating pixels that were not passed through by the ray is: In the formula, This represents the electron density value of the pixel to be interpolated. This indicates the electron density value of the pixel that has been corrected; This represents the distance between the pixel to be interpolated and the pixel that has been corrected; This indicates the number of pixels that have been corrected in the interpolation process; S4. Input the dSTEC measurement value obtained in step S2 into the improved three-dimensional tomography model constructed in step S3, and obtain the three-dimensional ionospheric perturbation through iterative calculation.

2. The three-dimensional ionospheric perturbation extraction method according to claim 1, characterized in that, In step S2, the formula for calculating the total electron content (STEC) of the strabismus is: ； In the formula, and These are the two carrier signal frequencies of the GPS satellites; and for and The carrier phase smoothing pseudorange value of the frequency signal; and These are the differential code offsets for GNSS receivers and GPS satellites, respectively.

3. The three-dimensional ionospheric perturbation extraction method according to claim 1, characterized in that, In step S2, the formula for calculating dSTEC is: In the formula, m is the number of STEC measurements by GPS satellites; n is the total number of pixels in the reconstructed area; The observation matrix; ε represents the electron density of STIDs; ε is the STEC measurement noise error.

Images