Large-scale marching type ionosphere disturbance three-dimensional monitoring method and system based on GNSS (Global Navigation Satellite System)

Abstract

The invention discloses a GNSS-based large-scale marching type ionosphere disturbance three-dimensional monitoring method and system. The method comprises the following steps: calculating the total electron content of an oblique ionosphere between a satellite and a receiver by adopting non-difference non-combination precise point positioning; based on Savitzky-Golay smoothing and Butterworth band-pass filtering, the total electron content of the detrending oblique ionized layer is obtained; based on spatial and temporal distribution formed by multiple stations, multiple satellites and multiple elevation puncture points, a three-dimensional ionosphere electron density disturbance model is constructed, and ionosphere disturbance three-dimensional propagation characteristics are solved. According to the method, the total electron content of the detrended oblique ionized layer is obtained based on Savitzky-Golay smoothing and Butterworth band-pass filtering, trend removal and disturbance reservation can be better balanced, and the signal-to-noise ratio and reliability of the extracted large-scale marching type ionized layer disturbance signal are improved; meanwhile, the constructed three-dimensional ionosphere electron density disturbance model not only can reveal vertical coupling characteristics of ionosphere disturbance, but also can stably extract key propagation elements and improve monitoring accuracy.

Description

Technical Field

[0001] This invention relates to the field of ionospheric disturbance monitoring, and in particular to a large-scale traveling three-dimensional monitoring method and system for ionospheric disturbance based on GNSS. Background Technology

[0002] The ionosphere is a crucial component of Earth's space environment, and its spatiotemporal variations in electron density directly impact critical applications such as satellite navigation, satellite communication, and shortwave positioning. Triggered by space weather processes like solar activity and geomagnetic storms, as well as natural and man-made events such as earthquakes, volcanic eruptions, and rocket launches, traveling disturbances with significant spatiotemporal coherence frequently occur in the ionosphere. These disturbances have significant engineering impacts at the regional scale, particularly increasing GNSS ranging errors, causing carrier phase discontinuities, and difficulties in integer cycle calculation, thereby reducing the stability and reliability of precise positioning and timing services. Therefore, accurate monitoring of traveling disturbances in the ionosphere is essential for ensuring the continuous availability of navigation and communication systems.

[0003] Large-scale traveling ionospheric disturbances (LSTIDs) typically exhibit quasi-periodic fluctuations in electron density in the F-region, with horizontal scales ranging from hundreds to thousands of kilometers and characteristic periods from ten to hundreds of minutes, and show distinct propagation directions and velocities. Their propagation process exhibits significant three-dimensional coupling characteristics: in addition to in-plane phase advancement, it is accompanied by variations in vertical energy distribution and peak height. Detection of LSTIDs using only two-dimensional projection planes is insufficient to reveal their vertical coupling and propagation mechanisms, limiting a deeper understanding of the disturbance's nature and potential impacts, and failing to provide sufficient constraint information for navigation enhancement and anti-interference strategies.

[0004] Existing monitoring methods mostly rely on two-dimensional TEC / dTEC fields retrieved from ground-based GNSS networks or time-delay cross-correlation analysis from a small number of stations. While these methods can identify the occurrence and horizontal propagation direction of disturbances to some extent, they are insufficient for characterizing the vertical structure of disturbances and providing a complete quantitative representation of three-dimensional propagation parameters. Therefore, it is necessary to develop a three-dimensional monitoring method that can robustly extract key propagation elements of LSTID and provide a structured description under general ground-based GNSS observation conditions. This would enhance the physical understanding of ionospheric disturbance processes and provide support for related applied research. Summary of the Invention

[0005] The purpose of this invention is to overcome the above-mentioned defects and problems in the prior art and provide a large-scale traveling three-dimensional monitoring method and system for ionospheric disturbances based on GNSS. This method can not only reveal the vertical coupling characteristics of ionospheric disturbances, but also robustly extract key propagation elements and improve the accuracy of monitoring.

[0006] To achieve the above objectives, the technical solution of the present invention is:

[0007] In a first aspect, the present invention provides a method for large-scale, traveling, three-dimensional monitoring of ionospheric disturbances based on GNSS, comprising:

[0008] The total electron content of the oblique ionosphere between the satellite and the receiver is calculated using non-differential and non-combined precise single-point positioning.

[0009] The detrended oblique total electron content of the ionosphere is obtained based on Savitzky-Golay smoothing and Butterworth bandpass filtering;

[0010] Based on the spatiotemporal distribution of multiple stations, multiple satellites, and multiple elevation angle puncture points, a three-dimensional ionospheric electron density perturbation model is constructed, and the three-dimensional propagation characteristics of ionospheric perturbation are obtained by solving the model.

[0011] Preferably, the construction of a three-dimensional ionospheric electron density perturbation model based on the spatiotemporal distribution of multiple stations, multiple satellites, and multiple elevation angle puncture points includes:

[0012] Spatial gridding was performed on ionospheric puncture points from different stations, satellites, and elevation angles at the same time. The puncture points were mapped to corresponding voxels according to their geographical latitude, longitude, and reference altitude. By integrating the geographical distribution of multiple stations, the orbital distribution of multiple satellites, and the multi-directional and multi-path constraints formed by multi-elevation angle observations, a set of linear equations between the total electron content of oblique perturbation and the electron density perturbation of each voxel was established. Spatial smoothing constraints and physical rationality constraints were applied to solve the equations, resulting in a three-dimensional ionospheric electron density perturbation model.

[0013] Preferably, the three-dimensional ionospheric electron density perturbation model is as follows:

[0014] ;

[0015] ;

[0016] ;

[0017] In the formula, For the first After the nth iteration Electron density perturbation value of individual elements; For the first The number of rays in a voxel; It is a relaxation factor; For the first The total electron content of the de-orbiting ionosphere on the observed X-ray; For the first The observed ray is at the Intercept values ​​in each grid; The first one provided for the IRI2020 model Background electron density values ​​of individual elements; To observe the number of rays; This represents the number of grid cells in the tomographic region. For the first The observed ray is at the Intercept values ​​in each grid; For the first Electron density perturbation values ​​in each grid; To observe the total electron content of the de-trending oblique ionosphere along the ray path; This represents the ray path from the GNSS satellite to the ground-based GNSS receiver. For height The electron density perturbation value at that location.

[0018] Preferably, the method for calculating the total electron content of the de-trending oblique ionosphere is as follows:

[0019] The trend term was obtained using the ionospheric STEC sequence from the Savitzky-Golay smoothed ground-based GNSS receiver to each GNSS satellite:

[0020] ;

[0021] In the formula, For in position The total electron content of the ionosphere estimated by the smooth trend at the location; The center index for the current trend to be calculated. ; The sequence length; For the relative offset index within the sliding window, from arrive , Half the window length, Indicates the sample on the left side of the center. Indicates the sample on the right side of the center; These are the Savitzky-Golay convolution coefficients;

[0022] The result obtained through a single detrending process is:

[0023] ;

[0024] In the formula, The position obtained using the Savitzky-Golay method The total electron content of the ionosphere at the de-trending angle; For in position The total electron content of the original oblique ionosphere at that location;

[0025] The following frequency band filtering was performed using the Butterworth bandpass filtering method:

[0026] ;

[0027] In the formula, The sampling frequency; The sampling time interval; The lower cutoff frequency; This is the upper cutoff frequency; The upper limit of the target period; The lower limit of the target period; for The normalized cutoff frequency; for The normalized cutoff frequency;

[0028] Employing forward and backward zero-phase filtering Obtain the total electron content of the detrended oblique ionosphere:

[0029] ;

[0030] In the formula, This represents the order of the Butterworth bandpass filter.

[0031] Preferably, the method for obtaining the three-dimensional propagation features includes:

[0032] On a three-dimensional ionospheric electron density perturbation model, analysis was performed using a sliding time window of fixed duration, and the main period of each sliding time window was determined.

[0033] Perform a two-dimensional fast Fourier transform on the horizontal slice of each sliding time window to determine the propagation direction and horizontal wavelength of the ionospheric disturbance;

[0034] Two adjacent profiles are selected along the propagation direction on the horizontal plane, and the cross-correlation function is calculated based on the disturbance time series of the two adjacent profiles. The time lag corresponding to the maximum correlation coefficient is obtained, and the horizontal propagation speed of the disturbance is obtained by dividing the distance between the two adjacent profiles by the time lag.

[0035] Preferably, the step of performing a two-dimensional fast Fourier transform on the horizontal slices of each sliding time window to determine the propagation direction and horizontal wavelength of the ionospheric disturbance includes:

[0036] Perform a two-dimensional fast Fourier transform on each horizontal slice of the sliding time window:

[0037] ;

[0038] In the formula, The transformed horizontal wavenumber spectrum; and The horizontal wavenumber; For each horizontal slice of the time window; and The horizontal coordinate; The imaginary unit;

[0039] The direction angle of the wavenumber vector where the main spectral peak is located is taken as the propagation direction of the ionospheric disturbance, and the distance from the main spectral peak to the origin is converted into the horizontal wavelength.

[0040] Preferably, the step of calculating the total electron content of the oblique ionosphere between the satellite and the receiver using non-differential, non-combination precise single-point positioning includes:

[0041] Collect ground-based multi-system, multi-frequency GNSS observation data and simultaneously acquire precise ephemeris, precise clock error, and differential code deviation products;

[0042] Construct the GNSS primitive pseudorange and carrier phase observation equations;

[0043] IGS was introduced to release precise ephemeris and clock bias products. IGS uses an ionospheric desaturation combination to estimate satellite clock bias, and the final clock bias product absorbs the satellite pseudorange hardware delay.

[0044] Based on clock difference products and the GNSS original pseudorange and carrier phase observation equations, a non-differential and non-combined precise single-point positioning model is obtained. Then, the GNSS original pseudorange and carrier phase observation equations are solved using the least squares method to determine the oblique ionospheric delay, i.e., the total electron content of the oblique ionospheric.

[0045] Preferably, the GNSS primitive pseudorange and carrier phase observation equations are as follows:

[0046] ;

[0047] ;

[0048] In the formula, superscript Indicates satellite; subscript Indicates receiver; subscript Indicates the L1 or L2 frequency band of the carrier signal; These are pseudorange observations; These are carrier phase observations; This is the coefficient matrix of the unknowns after linearization; For receiver position parameters; The geometric distance between the satellite and the receiver; The speed of light; and These are the clock differences between the receiver and the satellite, respectively. This is a tropospheric slant delay; This is the ionospheric slack delay; and These are the pseudorange hardware delays at the receiver and satellite ends, respectively. The carrier phase wavelength; For carrier phase integer ambiguity; and These are the phase hardware delays at the receiver and satellite ends, respectively. and These represent the sum of multipath effects, observation noise, and other unmodeled errors in pseudorange and carrier observations, respectively.

[0049] Secondly, the present invention provides a large-scale traveling ionospheric disturbance three-dimensional monitoring system based on GNSS, the system being used to implement the aforementioned large-scale traveling ionospheric disturbance three-dimensional monitoring method based on GNSS, the system comprising:

[0050] The STEC acquisition module for the ionosphere is used to calculate the total electron content of the oblique ionosphere between the satellite and the receiver using non-differential and non-combined precise single-point positioning.

[0051] The detrended ionosphere STEC acquisition module is used to acquire the total electron content of the detrended oblique ionosphere based on Savitzky-Golay smoothing and Butterworth bandpass filtering;

[0052] The three-dimensional propagation feature acquisition module is used to construct a three-dimensional ionospheric electron density perturbation model based on the spatiotemporal distribution of multiple stations, multiple satellites, and multiple elevation angle puncture points, and solve for the three-dimensional propagation features of ionospheric perturbation.

[0053] Thirdly, the present invention provides a large-scale traveling three-dimensional monitoring device for ionospheric disturbances based on GNSS, including a memory and a processor;

[0054] The memory is used to store computer program code and to transmit the computer program code to the processor;

[0055] The processor is configured to execute, according to instructions in the computer program code, the large-scale traveling ionospheric disturbance three-dimensional monitoring method based on GNSS as described above.

[0056] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0057] This invention discloses a large-scale, traveling ionospheric disturbance three-dimensional monitoring method and system based on GNSS. It employs non-differential, non-combined precise single-point positioning to calculate the total oblique ionospheric electron content between the satellite and receiver, improving the extraction accuracy and stability of the total oblique ionospheric electron content. Based on Savitzky-Golay smoothing and Butterworth bandpass filtering, it obtains the detrended oblique ionospheric electron content, better balancing trend removal and disturbance preservation, improving the signal-to-noise ratio and reliability of the extracted large-scale traveling ionospheric disturbance signal, and providing high-quality input data for subsequent three-dimensional inversion. Based on the spatiotemporal distribution of multiple stations, multiple satellites, and multiple elevation angle penetration points, a three-dimensional ionospheric electron density disturbance model is constructed, and the three-dimensional propagation characteristics of the ionospheric disturbance are obtained by solving the model. The constructed three-dimensional ionospheric electron density disturbance model not only inherits the advantages of multi-source geometric constraints but also significantly improves the horizontal and vertical resolution of ionospheric disturbance three-dimensional monitoring. Therefore, this invention can not only reveal the vertical coupling characteristics of ionospheric disturbances but also robustly extract key propagation elements, improving monitoring accuracy. Attached Figure Description

[0058] Figure 1 This is a flowchart of the large-scale traveling three-dimensional monitoring method for ionospheric disturbances based on GNSS proposed in this embodiment of the invention.

[0059] Figure 2 This is a schematic diagram of two-dimensional ionospheric TEC disturbance in a certain area during a geomagnetic storm, as described in an embodiment of the present invention.

[0060] Figure 3 This is a schematic diagram of the three-dimensional electron density perturbation distribution in a certain area during a geomagnetic storm, as described in an embodiment of the present invention.

[0061] Figure 4 This is a schematic diagram of the propagation characteristics of three-dimensional electron density disturbance in a certain region during a geomagnetic storm, as described in an embodiment of the present invention.

[0062] Figure 5 This is a structural block diagram of the large-scale traveling ionospheric disturbance three-dimensional monitoring system based on GNSS proposed in the embodiments of the present invention.

[0063] Figure 6 This is a structural block diagram of the large-scale traveling ionospheric disturbance three-dimensional monitoring device based on GNSS proposed in the embodiments of the present invention. Detailed Implementation

[0064] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0065] See Figure 1 This invention provides a large-scale, traveling, three-dimensional monitoring method for ionospheric disturbances based on GNSS, comprising:

[0066] S1. The total electron content of the oblique ionosphere between the satellite and the receiver is calculated using non-differential and non-combined precise single-point positioning.

[0067] S2. Obtain the detrended oblique total electron content of the ionosphere based on Savitzky-Golay smoothing and Butterworth bandpass filtering;

[0068] S3. Based on the spatiotemporal distribution of multiple stations, multiple satellites, and multiple elevation angle puncture points, a three-dimensional ionospheric electron density perturbation model is constructed, and the three-dimensional propagation characteristics of ionospheric perturbation are obtained by solving the model.

[0069] Before step S1, ground-based multi-system multi-frequency GNSS observation data (RINEX format) is collected with a sampling interval of less than or equal to 30 seconds. The GFZRNX data processing software is used to perform data format standardization, cycle slip detection, and quality check analysis on the raw observation data, and precise ephemeris, precise clock error, and differential code deviation products are acquired simultaneously.

[0070] This invention employs non-differential, non-combined precise point positioning (PPP) to calculate the oblique total ionospheric electron content (STEC) between the satellite and receiver; it constructs detrended STEC (dSTEC) based on Savitzky-Golay (SG) smoothing and Butterworth bandpass filtering; and it combines the spatiotemporal distribution of multiple stations, multiple satellites, and multiple elevation angle puncture points (IPPs) to construct a three-dimensional ionospheric electron density perturbation model, which can effectively identify and quantify the LSTID propagation characteristics caused by space weather, natural disasters, and human events (such as geomagnetic storms, strong earthquakes / volcanoes, and rocket launches). Compared with existing monitoring methods that are mainly based on two-dimensional dTEC monitoring or single-station geostationary satellite (GEO) fitting, this invention achieves quantitative inversion from "two-dimensional projection" to "true three-dimensional propagation," effectively compensating for the shortcomings of existing two-dimensional methods in vertical coupling and propagation characterization.

[0071] Furthermore, non-differential, non-combined precise single-point positioning is used to calculate the total electron content of the oblique ionosphere between the satellite and the receiver, including:

[0072] Collect ground-based multi-system, multi-frequency GNSS observation data and simultaneously acquire precise ephemeris, precise clock error, and differential code deviation products;

[0073] Construct the GNSS primitive pseudorange and carrier phase observation equations;

[0074] IGS was introduced to release precise ephemeris and clock bias products. IGS uses an ionospheric desaturation combination to estimate satellite clock bias, and the final clock bias product absorbs the satellite pseudorange hardware delay.

[0075] Based on clock difference products and the GNSS original pseudorange and carrier phase observation equations, a non-differential and non-combined precise single-point positioning model is obtained. Then, the GNSS original pseudorange and carrier phase observation equations are solved using the least squares method to determine the oblique ionospheric delay, i.e., the total electron content of the oblique ionospheric.

[0076] Specifically, the GNSS primitive pseudorange and carrier phase observation equations are as follows:

[0077] ;

[0078] ;

[0079] The linearized equation of the above expression is:

[0080] ;

[0081] ;

[0082] In the formula, superscript Indicates satellite; subscript Indicates receiver; subscript Indicates the L1 or L2 frequency band of the carrier signal; These are pseudorange observations; These are carrier phase observations; This is the coefficient matrix of the unknowns after linearization; For receiver position parameters; The geometric distance between the satellite and the receiver; The speed of light; and These are the clock differences between the receiver and the satellite, respectively. This is a tropospheric slant delay; This is the ionospheric slack delay; and These are the pseudorange hardware delays at the receiver and satellite ends, respectively. The carrier phase wavelength; For carrier phase integer ambiguity; and These are the phase hardware delays at the receiver and satellite ends, respectively. and These represent the sum of multipath effects, observation noise, and other unmodeled errors in pseudorange and carrier observations, respectively.

[0083] Specifically, IGS is introduced to release precise ephemeris and clock bias products. IGS uses an ionospheric desaturation combination to estimate satellite clock bias, and the final clock bias product incorporates the satellite pseudorange hardware delay, as shown in the following formula:

[0084] ;

[0085] In the formula, Satellite clock bias estimated for ionospheric combination; This refers to the frequency of the L1 band; This refers to the frequency of the L2 band; For pseudorange observation clock bias of the satellite in the L1 band; This refers to the pseudorange observation clock bias of the satellite in the L2 band.

[0086] Substituting the satellite clock bias estimated by the ionosphere elimination combination into the linearized equations of the original pseudorange and carrier phase observations of GNSS, the final non-differential non-combined precise single-point positioning model is obtained.

[0087] Compared to existing methods for extracting total electron content in the ionosphere using geometrically independent combinations, the proposed method based on a non-differential, non-combined precise single-point positioning model simultaneously corrects for error sources such as precise ephemeris, precise clock error, antenna phase center, and solid tide in the observation equation, thereby improving the extraction accuracy and stability of total electron content in the oblique ionosphere.

[0088] Furthermore, the method for calculating the total electron content of the de-trending oblique ionosphere is as follows:

[0089] The trend term was obtained using the ionospheric STEC sequence from the Savitzky-Golay smoothed ground-based GNSS receiver to each GNSS satellite:

[0090] ;

[0091] In the formula, For in position The total electron content of the oblique ionosphere estimated by the smooth trend at the location (the "trend" of the Savitzky-Golay filter output) is used to characterize the slowly varying background components such as diurnal variation. The center index (integer) for the current trend to be calculated. ; The sequence length; For the relative offset index within the sliding window, from arrive , Half the window length, Indicates the sample on the left side of the center. Indicates the sample on the right side of the center; These are the Savitzky-Golay convolution coefficients, which are fixed weights that vary with... The change is determined by the length of the selected sliding window. ( ) and polynomial order It is uniquely determined by the existing formula;

[0092] The result obtained through a single detrending process is:

[0093] ;

[0094] In the formula, The position obtained using the Savitzky-Golay method The total electron content of the ionosphere at the de-trending angle; For in position The total electron content of the original oblique ionosphere at that location;

[0095] The following frequency band filtering was performed using the Butterworth bandpass filtering method:

[0096] ;

[0097] In the formula, The sampling frequency; The sampling time interval; The lower cutoff frequency; This is the upper cutoff frequency; The upper limit of the target period; The lower limit of the target period; for The normalized cutoff frequency; for The normalized cutoff frequency;

[0098] Employing forward and backward zero-phase filtering Obtain the total electron content of the detrended oblique ionosphere:

[0099] ;

[0100] In the formula, This represents the order of the Butterworth bandpass filter.

[0101] In this invention, Savitzky-Golay smoothing achieves smoothing through local polynomial fitting. While preserving the amplitude and phase characteristics of the ionospheric disturbance waveform, it effectively extracts slowly varying background trends such as diurnal variation, avoiding waveform distortion caused by simple moving averages. The Butterworth bandpass filter has a flat amplitude-frequency response within its passband and a steep transition band, effectively allowing large-scale traveling disturbance signals to pass through within a set period range, while suppressing long-term trends and high-frequency random noise, reducing phase distortion introduced by filtering. The combination of Savitzky-Golay smoothing and the Butterworth bandpass filter can first strip away large-scale background factors such as diurnal ionospheric variation and slowly varying instrument bias, and then perform precise bandpass filtering for the target period range. Compared to single high-pass or bandpass filtering methods, it better balances trend removal and disturbance preservation, improving the signal-to-noise ratio and reliability of the extracted large-scale traveling ionospheric disturbance signals, providing high-quality input data for subsequent 3D inversion. Therefore, this invention can more stably extract typical traveling disturbance signals against a complex ionospheric background, reducing the impact of detrending methods on disturbance morphology and propagation parameter estimation.

[0102] Furthermore, based on the spatiotemporal distribution of multiple stations, satellites, and puncture points at multiple elevation angles, a three-dimensional ionospheric electron density perturbation model is constructed, including:

[0103] Spatial gridding was performed on ionospheric puncture points from different stations, satellites, and elevation angles at the same time. The puncture points were mapped to corresponding voxels according to their geographical latitude, longitude, and reference altitude. By integrating the geographical distribution of multiple stations, the orbital distribution of multiple satellites, and the multi-directional and multi-path constraints formed by multi-elevation angle observations, a set of linear equations between the total electron content of oblique perturbation and the electron density perturbation of each voxel was established. Spatial smoothing constraints and physical rationality constraints were applied to solve the equations, resulting in a three-dimensional ionospheric electron density perturbation model.

[0104] In this invention, ground receivers at different geographical locations provide multi-directional observations of the target area, resulting in a more uniform coverage of ionospheric puncture points on the plane and mitigating the ill-conditioned problems caused by a small number of stations. GNSS satellites on different orbital planes provide multi-directional observation rays for the same station. The staggered crossing of multiple satellite trajectories in the ionosphere effectively increases ray constraints in different azimuths, improving the geometry of the 3D inversion. By utilizing observations at different elevation angles, constraints of different puncture heights and path lengths are formed for the same horizontal position. Low elevation angle observations contribute more to higher-altitude voxels, while high elevation angle observations have stronger constraints on nearby voxels. Thus, when constructing the observation equations, the total electron content of oblique perturbation is allocated to voxels at different altitudes according to the path length, achieving resolution of the vertical structure. For low elevation angle observations, the contribution of the observation path's crossing length in the altitude direction is considered. By calculating the path length of each observation ray in voxels at different altitudes, the total electron content of oblique perturbation is allocated to multiple voxels in the altitude direction to enhance vertical resolution. Building upon the above, the observation area is discretized into three-dimensional voxels. A spatiotemporal distribution consisting of multiple stations, satellites, and elevation-angle penetration points is used to establish a system of linear equations relating the total electron content of the oblique perturbation to the electron density perturbation of each voxel. By applying spatial smoothing and physical rationality constraints, a temporally continuous three-dimensional ionospheric electron density perturbation model can be obtained. This model inherits the advantages of multi-source geometric constraints and significantly improves the horizontal and vertical resolution of three-dimensional ionospheric perturbation monitoring.

[0105] Specifically, the three-dimensional ionospheric electron density perturbation model includes the following steps:

[0106] The total electron content of the de-orbiting ionosphere along the observed ray path can be represented by the integral value of the electron density perturbation along the ray path:

[0107] ;

[0108] In the formula, To observe the total electron content of the de-trending oblique ionosphere along the ray path; This represents the ray path from the GNSS satellite to the ground-based GNSS receiver. For height The electron density perturbation value at that location.

[0109] Based on the distance each ray travels through the grid and The observation equation can be composed as follows:

[0110] ;

[0111] In the formula, To observe the number of rays; This represents the number of grid cells in the tomographic region. For the first The observed ray is at the Intercept values ​​in each grid; For the first Electron density perturbation values ​​in each grid.

[0112] The electron density perturbation value is solved iteratively using the Synchronous Algebraic Reconstruction Algorithm (SART), with the initial background set to 0. The SART algorithm is shown below:

[0113] ;

[0114] In the formula, For the first After the nth iteration Electron density perturbation value of individual elements; For the first The number of rays in a voxel; It is a relaxation factor; For the first The total electron content of the de-orbiting ionosphere on the observed X-ray; For the first The observed ray is at the Intercept values ​​in each grid; The first one provided for the IRI2020 model Background electron density value of individual elements.

[0115] This invention constructs a three-dimensional electron density perturbation field through observations from multiple stations, satellites, and elevation angles, avoiding the loss of altitude information and wave system aliasing problems caused by the simple integral projection of perturbations at different altitude levels in traditional two-dimensional TEC monitoring. Based on the inverted three-dimensional perturbation field, this invention performs spectral and correlation analyses in the time, plane, and altitude dimensions, respectively. It can directly extract the principal period, horizontal propagation direction, horizontal wavelength, and propagation velocity of LSTID at each altitude level, and characterize the coupling relationship between perturbation amplitude and phase with altitude, thus forming a structured three-dimensional description of the LSTID propagation process. Compared with existing two-dimensional methods, this invention integrates the geographical distribution of multiple stations, the orbital distribution of multiple satellites, and the multi-directional and multi-path constraints formed by multi-elevation angle observations to construct a weighted observation matrix. It applies smoothing and physical rationality constraints to the observation matrix, and obtains voxel electron density perturbations through iterative reconstruction, thereby improving the horizontal and vertical resolution and robustness of three-dimensional ionospheric perturbation inversion. This invention not only reveals the vertical coupling characteristics of ionospheric disturbances, but also robustly extracts key propagation elements under multi-wave system and strong noise conditions, providing a more complete and reliable three-dimensional monitoring method for subsequent space weather research and navigation applications.

[0116] Furthermore, methods for obtaining three-dimensional propagation features include:

[0117] On the three-dimensional ionospheric electron density perturbation model, analysis was performed using a sliding time window of fixed duration, and the main period of each sliding time window was determined by simplified spectral analysis.

[0118] Within the local neighborhood of each horizontal position, the propagation direction and horizontal wavelength are extracted using the "spatial spectrum method". A two-dimensional fast Fourier transform is performed on each horizontal slice of the sliding time window to determine the propagation direction and horizontal wavelength of the ionospheric disturbance.

[0119] Given the principal period and propagation direction, the propagation velocity of ionospheric disturbances is estimated using the "time-delay cross-correlation method." Two adjacent profiles are selected along the propagation direction on a horizontal plane, and the cross-correlation function is calculated based on the extracted disturbance time series from these two adjacent profiles. The time lag corresponding to the point of maximum correlation coefficient is obtained, and the horizontal propagation velocity of the disturbance is obtained by dividing the distance between the two adjacent profiles by the time lag. Finally, the three-dimensional propagation parameter characteristics are output for subsequent statistical and comparative analysis.

[0120] Specifically, within each sliding time window, a discrete Fourier transform is performed on the perturbation time series to calculate the power spectral density. The frequency with the highest power within a preset frequency range is selected as the main frequency, and its reciprocal is the main period within that time window.

[0121] Specifically, a two-dimensional fast Fourier transform is performed on the horizontal slice of each sliding time window to determine the propagation direction and horizontal wavelength of the ionospheric disturbance, including:

[0122] Perform a two-dimensional fast Fourier transform on each horizontal slice of the sliding time window:

[0123] ;

[0124] In the formula, The transformed horizontal wavenumber spectrum; and The horizontal wavenumber; For each horizontal slice of the time window; and The horizontal coordinate; The imaginary unit;

[0125] The direction angle of the wavenumber vector where the main spectral peak is located is taken as the propagation direction of the ionospheric disturbance, and the distance from the main spectral peak to the origin is converted into the horizontal wavelength.

[0126] Specifically, horizontal wavelength The calculation formula is:

[0127] ; ;

[0128] In the formula, and These represent the wavenumber components of the main spectral peak in the two-dimensional spatial spectrum in the meridional and latitudinal directions, respectively. The horizontal wavenumber modulus.

[0129] Specifically, the cross-correlation function is:

[0130] ;

[0131] In the formula, and To extract the perturbation time series of two profiles at a given height level or height integral; For the first Samples at each sampling time; For the first There are several time lags. The corresponding time lag is obtained by searching for the maximum value.

[0132] Specifically, the horizontal propagation speed of the disturbance The calculation formula is:

[0133] ;

[0134] In the formula, The distance between the center points of the two cross sections along the propagation direction; This is the propagation time delay of the disturbance signal between the two profiles.

[0135] This invention performs a two-dimensional fast Fourier transform on the horizontal slices of each time window to first determine the main propagation direction and horizontal wavelength of the ionospheric disturbance. After obtaining the propagation direction, a profile can be selected along this direction, and time lag can be calculated using time delay cross-correlation, thus ensuring that the propagation velocity is estimated in the correct direction, resulting in more stable and reliable results. Simultaneously, the horizontal wavelength obtained from the two-dimensional spectrum corresponds to the main period in the time domain, providing a simple physical reference for the propagation velocity, which can be used to check and constrain subsequent velocity inversion.

[0136] This section presents an example of the three-dimensional propagation characteristics of large-scale traveling ionospheric disturbances triggered by GNSS-monitored space weather events (geomagnetic storms).

[0137] Using dual-frequency observation data from ground-based GNSS stations in a certain area, combined with precise ephemeris and clock error products, the observation data of all GNSS stations are preprocessed to remove cycle slips and gross errors. Then, the time series of total ionospheric electron content between each satellite and the receiver is calculated according to the non-differential non-combined precise single-point positioning model described in this invention.

[0138] After obtaining the total electron content of the oblique ionosphere, the Savitzky-Golay smoothing and Butterworth bandpass filtering detrending method described in this invention is used to remove diurnal variation background and high-frequency noise, obtaining large-scale perturbation dSTEC signals for each ray. Further, the study area is discretized into a regular three-dimensional grid along the longitude, latitude, and altitude directions. Observation equations are constructed according to the three-dimensional tomographic inversion steps, and the temporally continuous three-dimensional electron density perturbation field is obtained by solving for it. Based on this, the three-dimensional propagation parameter characteristics are acquired.

[0139] Figure 2 The results of TEC perturbation calculated using the traditional two-dimensional method are presented. These results can only reflect the overall propagation pattern of the perturbation on the plane and cannot distinguish the response differences at different altitudes. Figure 3 The electron density perturbation distributions at multiple height layers obtained by the method of this invention at the same time are presented. The results show that the LSTID amplitude is significantly enhanced near the peak of the F2 layer and exhibits a decay characteristic at higher heights, demonstrating a clear vertical structure. Figure 4 ( Figure 4 The black dashed line in the middle represents the LSTID wavefront. It shows the spatial distribution of propagation direction, horizontal wavelength and propagation speed at different altitudes at the same time, extracted by the method of the present invention. It can be seen that the propagation direction remains consistent at each altitude, while the propagation speed and disturbance amplitude change with altitude, indicating that the event has obvious vertical coupling characteristics.

[0140] See Figure 5 The present invention also provides a large-scale traveling ionospheric disturbance three-dimensional monitoring system based on GNSS, the system being used to implement the above-described large-scale traveling ionospheric disturbance three-dimensional monitoring method based on GNSS, the system comprising:

[0141] The STEC acquisition module for the ionosphere is used to calculate the total electron content of the oblique ionosphere between the satellite and the receiver using non-differential and non-combined precise single-point positioning.

[0142] The detrended ionosphere STEC acquisition module is used to acquire the total electron content of the detrended oblique ionosphere based on Savitzky-Golay smoothing and Butterworth bandpass filtering;

[0143] The three-dimensional propagation feature acquisition module is used to construct a three-dimensional ionospheric electron density perturbation model based on the spatiotemporal distribution of multiple stations, multiple satellites, and multiple elevation angle puncture points, and solve for the three-dimensional propagation features of ionospheric perturbation.

[0144] See Figure 6 The present invention also provides a large-scale traveling three-dimensional monitoring device for ionospheric disturbance based on GNSS, including a memory and a processor;

[0145] The memory is used to store computer program code and to transmit the computer program code to the processor;

[0146] The processor is configured to execute, according to instructions in the computer program code, the large-scale traveling ionospheric disturbance three-dimensional monitoring method based on GNSS as described above.

[0147] The present invention also provides a computer-readable storage medium storing a computer program, which, when executed by a processor, implements the GNSS-based large-scale traveling ionospheric disturbance three-dimensional monitoring method as described above.

[0148] Generally, the computer instructions for implementing the method of the present invention can be carried on any combination of one or more computer-readable storage media. Non-transitory computer-readable storage media can include any computer-readable medium except for the signal itself, which is temporarily propagating.

[0149] Computer-readable storage media can be, for example, but not limited to, electrical, magnetic, optical, electromagnetic, infrared, or semiconductor systems, apparatuses, or any combination thereof. More specific examples (a non-exhaustive list) of computer-readable storage media include: electrical connections having one or more wires, portable computer disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EKROM or flash memory), optical fibers, portable compact disk read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination thereof. In this invention, a computer-readable storage medium can be any tangible medium containing or storing a program that can be used by or in conjunction with an instruction execution system, apparatus, or device.

[0150] Computer program code for performing the operations of this invention can be written in one or more programming languages ​​or a combination thereof. These programming languages ​​include object-oriented programming languages—such as Java, Smalltalk, and C++—as well as conventional procedural programming languages—such as the "C" language or similar programming languages. In particular, Python, suitable for neural network computation, and platform frameworks such as TensorFlow and PyTorch can be used. The program code can be executed entirely on the user's computer, partially on the user's computer, as a standalone software package, partially on the user's computer and partially on a remote computer, or entirely on a remote computer or server. In cases involving remote computers, the remote computer can be connected to the user's computer or to an external computer (e.g., via the Internet using an Internet service provider) through any type of network, including a local area network (LAN) or a wide area network (WAN).

[0151] The aforementioned equipment and non-transitory computer-readable storage media can be found in the detailed description of a GNSS-based large-scale traveling three-dimensional monitoring method for ionospheric disturbances and its beneficial effects, which will not be repeated here.

[0152] Although embodiments of the present invention have been shown and described above, it should be understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those skilled in the art can make changes, modifications, substitutions and variations to the above embodiments within the scope of the present invention.

Claims

1. A method for large-scale, traveling, three-dimensional monitoring of ionospheric disturbances based on GNSS, characterized in that, include: The total electron content of the oblique ionosphere between the satellite and the receiver is calculated using non-differential and non-combined precise single-point positioning. The detrended oblique total electron content of the ionosphere is obtained based on Savitzky-Golay smoothing and Butterworth bandpass filtering; Based on the spatiotemporal distribution of multiple stations, multiple satellites, and multiple elevation angle puncture points, a three-dimensional ionospheric electron density perturbation model is constructed, and the three-dimensional propagation characteristics of ionospheric perturbation are obtained by solving the model.

2. The method for large-scale traveling three-dimensional monitoring of ionospheric disturbances based on GNSS according to claim 1, characterized in that, The three-dimensional ionospheric electron density perturbation model, based on the spatiotemporal distribution of multiple stations, multiple satellites, and multiple elevation angle puncture points, includes: Spatial gridding was performed on ionospheric puncture points from different stations, satellites, and elevation angles at the same time. The puncture points were mapped to corresponding voxels according to their geographical latitude, longitude, and reference altitude. By integrating the geographical distribution of multiple stations, the orbital distribution of multiple satellites, and the multi-directional and multi-path constraints formed by multi-elevation angle observations, a set of linear equations between the total electron content of oblique perturbation and the electron density perturbation of each voxel was established. Spatial smoothing constraints and physical rationality constraints were applied to solve the equations, resulting in a three-dimensional ionospheric electron density perturbation model.

3. The method for large-scale traveling three-dimensional monitoring of ionospheric disturbances based on GNSS according to claim 2, characterized in that, The three-dimensional electron density perturbation model for the ionosphere is as follows: ； ； ； In the formula, For the first After the nth iteration Electron density perturbation value of individual elements; For the first The number of rays in a voxel; It is a relaxation factor; For the first The total electron content of the de-orbiting ionosphere on the observed X-ray; For the first The observed ray is at the Intercept values ​​in each grid; The first one provided for the IRI2020 model Background electron density values ​​of individual elements; To observe the number of rays; This represents the number of grid cells in the tomographic region. For the first The observed ray is at the Intercept values ​​in each grid; For the first Electron density perturbation values ​​in each grid; To observe the total electron content of the de-trending oblique ionosphere along the ray path; This represents the ray path from the GNSS satellite to the ground-based GNSS receiver. For height The electron density perturbation value at that location.

4. The method for large-scale traveling three-dimensional monitoring of ionospheric disturbances based on GNSS according to claim 1, characterized in that, The method for calculating the total electron content of the de-trending oblique ionosphere is as follows: The trend term was obtained using the ionospheric STEC sequence from the Savitzky-Golay smoothed ground-based GNSS receiver to each GNSS satellite: ； In the formula, For in position The total electron content of the ionosphere estimated by the smooth trend at the location; The center index for the current trend to be calculated. ; The sequence length; For the relative offset index within the sliding window, from arrive , Half the window length, Indicates the sample on the left side of the center. Indicates the sample on the right side of the center; These are the Savitzky-Golay convolution coefficients; The result obtained through a single detrending process is: ； In the formula, The position obtained using the Savitzky-Golay method The total electron content of the ionosphere at the de-trending angle; For in position The total electron content of the original oblique ionosphere at that location; The following frequency band filtering was performed using the Butterworth bandpass filtering method: ； In the formula, The sampling frequency; The sampling time interval; The lower cutoff frequency; The upper cutoff frequency; The upper limit of the target period; The lower limit of the target period; for The normalized cutoff frequency; for The normalized cutoff frequency; Employing forward and backward zero-phase filtering Obtain the total electron content of the detrended oblique ionosphere: ； In the formula, This represents the order of the Butterworth bandpass filter.

5. The method for large-scale traveling three-dimensional monitoring of ionospheric disturbances based on GNSS according to claim 1, characterized in that, The method for obtaining the three-dimensional propagation features includes: On a three-dimensional ionospheric electron density perturbation model, analysis was performed using a sliding time window of fixed duration, and the main period of each sliding time window was determined. Perform a two-dimensional fast Fourier transform on the horizontal slice of each sliding time window to determine the propagation direction and horizontal wavelength of the ionospheric disturbance; Two adjacent profiles are selected along the propagation direction on the horizontal plane, and the cross-correlation function is calculated based on the disturbance time series of the two adjacent profiles. The time lag corresponding to the maximum correlation coefficient is obtained, and the horizontal propagation speed of the disturbance is obtained by dividing the distance between the two adjacent profiles by the time lag.

6. The method for large-scale traveling three-dimensional monitoring of ionospheric disturbances based on GNSS according to claim 5, characterized in that, The step of performing a two-dimensional fast Fourier transform on the horizontal slices of each sliding time window to determine the propagation direction and horizontal wavelength of the ionospheric disturbance includes: Perform a two-dimensional fast Fourier transform on each horizontal slice of the sliding time window: ； In the formula, The transformed horizontal wavenumber spectrum; and The horizontal wave number; For each horizontal slice of the time window; and The horizontal coordinate; The imaginary unit; The direction angle of the wavenumber vector where the main spectral peak is located is taken as the propagation direction of the ionospheric disturbance, and the distance from the main spectral peak to the origin is converted into the horizontal wavelength.

7. The method for large-scale traveling three-dimensional monitoring of ionospheric disturbances based on GNSS according to claim 1, characterized in that, The method of calculating the total electron content of the oblique ionosphere between the satellite and the receiver using non-differential, non-combined, precise single-point positioning includes: Collect ground-based multi-system, multi-frequency GNSS observation data and simultaneously acquire precise ephemeris, precise clock error, and differential code deviation products; Construct the GNSS primitive pseudorange and carrier phase observation equations; IGS was introduced to release precise ephemeris and clock bias products. IGS uses an ionospheric desaturation combination to estimate satellite clock bias, and the final clock bias product absorbs the satellite pseudorange hardware delay. Based on clock difference products and the GNSS original pseudorange and carrier phase observation equations, a non-differential and non-combined precise single-point positioning model is obtained. Then, the GNSS original pseudorange and carrier phase observation equations are solved using the least squares method to determine the oblique ionospheric delay, i.e., the total electron content of the oblique ionospheric.

8. The method for large-scale traveling three-dimensional monitoring of ionospheric disturbances based on GNSS according to claim 7, characterized in that, The GNSS primitive pseudorange and carrier phase observation equations are as follows: ； ； In the formula, superscript Indicates satellite; subscript Indicates receiver; subscript Indicates the L1 or L2 frequency band of the carrier signal; These are pseudorange observations; These are carrier phase observations; This is the coefficient matrix of the unknowns after linearization; For receiver position parameters; The geometric distance between the satellite and the receiver; The speed of light; and These are the clock differences between the receiver and the satellite, respectively. This is a tropospheric slant delay; This is the ionospheric slack delay; and These are the pseudorange hardware delays at the receiver and satellite ends, respectively. The carrier phase wavelength; For carrier phase integer ambiguity; and These are the phase hardware delays at the receiver and satellite ends, respectively. and These represent the sum of multipath effects, observation noise, and other unmodeled errors in pseudorange and carrier observations, respectively.

9. A large-scale, traveling, three-dimensional monitoring system for ionospheric disturbances based on GNSS, characterized in that, The system is used to implement the method according to any one of claims 1 to 8, the system comprising: The STEC acquisition module for the ionosphere is used to calculate the total electron content of the oblique ionosphere between the satellite and the receiver using non-differential and non-combined precise single-point positioning. The detrended ionosphere STEC acquisition module is used to acquire the total electron content of the detrended oblique ionosphere based on Savitzky-Golay smoothing and Butterworth bandpass filtering; The three-dimensional propagation feature acquisition module is used to construct a three-dimensional ionospheric electron density perturbation model based on the spatiotemporal distribution of multiple stations, multiple satellites, and multiple elevation angle puncture points, and solve for the three-dimensional propagation features of ionospheric perturbation.

10. A large-scale, traveling, three-dimensional monitoring device for ionospheric disturbances based on GNSS, characterized in that, Including memory and processor; The memory is used to store computer program code and transmit the computer program code to the processor; The processor is configured to execute the method as described in any one of claims 1 to 8 according to instructions in the computer program code.

Images