A new method for monitoring ionospheric three-dimensional disturbances

Abstract

The present invention discloses a novel ionospheric three-dimensional disturbance monitoring method, which specifically identifies and monitors TID from a three-dimensional spatial level, and is more in line with the spatiality of TID than traditional methods; based on reliable electron density inversion, the present invention identifies and monitors TID by calculating the spatial electron density changes at different times in the three-dimensional space, and can effectively monitor the three-dimensional spatial position, change rate and trend of TID; compared with traditional analytical indexes for identifying and monitoring ionospheric disturbances, the present invention proposes an index for monitoring and identifying ionospheric disturbances in three-dimensional space, which can not only identify ionospheric disturbances, but also analyze the position of disturbances in space; effectively realizes accurate identification and monitoring of space ionospheric disturbances, and helps researchers in the fields of precision positioning, navigation services, and aerospace to better understand the ionospheric disturbance mechanism.

Description

Technical Field

[0001] The invention relates to the technical field of near-Earth space environment detection, and in particular to a novel ionospheric three-dimensional disturbance monitoring method. Background Art

[0002] At present, ionospheric detection technology is mainly divided into ground-based and air-based systems. The ground-based system includes the Global Navigation Satellite Positioning System (GNSS), ionospheric vertical measurement stations and high-power incoherent scattering radars; the air-based system mainly uses borehole balloons and radio occultation as the means. In a comprehensive comparison of operating costs and monitoring efficiency, GNSS technology has become the main means of ionospheric three-dimensional spatial monitoring due to its advantages of all-weather, all-day, wide coverage, good real-time performance and high observation accuracy. The three-dimensional spatial monitoring of the ionosphere mainly relies on computed tomography (CT), that is, the spatial electron density (IED) distribution of the ionosphere is reconstructed through the tilted total electron content (STEC) between the ground station and the satellite, and its temporal and spatial changes are visualized. Commonly used ionospheric tomography methods can be mainly divided into the following according to their solution methods: (1) Iterative algorithms, such as algebraic reconstruction algorithm (ART) and simultaneous iterative algebraic reconstruction algorithm (SIRT). Iterative algorithms require the use of initial values ​​to participate in the inversion, so the inversion process has a certain dependence on the initial values; (2) Non-iterative algorithms, which are basis function methods, are not strongly dependent on the background initial values. Such as principal component analysis (PCA), tomography algorithm based on spherical function station, Kalman filter algorithm, etc. Both types of algorithms have their own advantages and disadvantages. It is particularly important to determine the appropriate inversion algorithm based on the characteristics of the inversion area and the characteristics of the algorithm.

[0003] Ionospheric disturbances can reconstruct the distribution of large-scale electron density in the three-dimensional space of the ionosphere. In mid- and low-latitude regions, common ionospheric disturbances include traveling ionospheric disturbances (TIDs) and equatorial plasma bubbles (EPBs). Among them, TID is the main cause of ionospheric irregularities. TID is generally divided into large-scale TID and medium-scale TID. The horizontal wavelength of large-scale TID is greater than 1000 kilometers, the horizontal speed is greater than 400m / s, and the propagation time interval is between 30 minutes and 3 hours; the horizontal wavelength of medium-scale TID is between hundreds of kilometers, the horizontal speed is between 100-250m / s, and the propagation time interval is tens of minutes. At present, the methods for monitoring TID include high-frequency radar technology, airglow imager technology, and ionospheric sounder technology. Due to the high cost of these equipment, they have not yet been put into use on a large scale, so there are problems of sparse distribution and small coverage. In previous studies, TEC ratio (ROT) and TEC ratio index (ROTI) have been widely used in disturbance research. ROTI is the standard deviation of ROT in a specific period of time. There is a good correspondence between ROTI and EPB, which can be used to effectively identify EPB. However, ROTI is significantly less effective in identifying TID, which means that during TID, the phase change of ionospheric electron density during the periodic deviation of ionospheric background value may be relatively small in a short time; the two-dimensional detrended TEC (DTEC) map can be used to show the time and two-dimensional plane changes of medium-sized traveling ionospheric disturbances; the integral and rate of total electron content index (IROTI) is the integral product of ROTI and DTEC, which monitors ionospheric disturbances based on the amplitude change and phase mutation of IROTI based on the slant electron total content (STEC). Although the existing ionospheric disturbance indices can monitor phenomena such as ionospheric traveling disturbances or plasma bubbles, they often exist on a one-dimensional (altitude) or two-dimensional (longitude × latitude) scale. However, ionospheric disturbances can span hundreds of kilometers on a plane scale and reach tens to hundreds of kilometers in altitude. The use of ionospheric tomography to analyze ionospheric disturbances should provide the location, range, intensity and rate of change of the disturbance. The accurate monitoring of the spatial structure of ionospheric irregularities by tomography technology plays an important role in building a high-precision ionospheric model and improving the accuracy of satellite navigation positioning.

[0004] Therefore, this paper proposes a new method for monitoring ionospheric three-dimensional disturbances. Summary of the invention

[0005] The purpose of the present invention is to design an ionospheric tomography technology based on ground-based GNSS data to detect ionospheric traveling wave disturbances in three-dimensional space, and to provide a new ionospheric three-dimensional disturbance monitoring method that can quickly monitor the location, change rate, range and trend of ionospheric traveling wave disturbances; and to provide a new method for monitoring three-dimensional ionospheric disturbances.

[0006] In order to achieve the above technical effects, the present invention is implemented by the following technical solutions: A novel ionospheric three-dimensional disturbance monitoring method, characterized in that it includes the following steps:

[0007] S1. Discretize the area to be inverted spatially, that is, divide the area to be inverted into three-dimensional grids of a certain resolution in the form of grids according to the latitude, longitude and altitude directions, which are also called voxels, and assign numbers to the discretized voxels;

[0008] S2. Obtain the ground-based GNSS observation value files, GPS navigation message files, precise ephemeris and differential code deviation files in the inversion area;

[0009] S3, solving the tilted total electron content STEC according to the observation file in S2;

[0010] S4, determine the STEC effective rays passing through the voxels in the study area, and calculate the intercepts between the effective rays and the voxels to construct the intercept matrix A;

[0011] S5. Determine the initial value of the electron density of each voxel in the inversion area based on the electron density solved by the IRI2016 model as the initial value;

[0012] S6, using the SART iterative algorithm to solve the three-dimensional ionospheric electron density in the area to be inverted, thereby obtaining the electron density distribution of the ionosphere in three-dimensional space;

[0013] S7, the reconstructed 3D IED was further verified using the electron density profile of the ionospheric radiosonde to ensure the reliability of the inversion results;

[0014] S8. Based on the previous results, further invert multiple time periods in the inversion area and verify the accuracy of the inversion results;

[0015] S9. Arrange the three-dimensional electron density grids of the inversion area in multiple time periods in a certain order, calculate the maximum or minimum value of the ROI in the three directions of XYZ respectively, and finally use the weighted centroid method to determine the center of the TID disturbance.

[0016] Furthermore, in S3, the GNSS signals needed to solve STEC mainly include two common observation quantities: pseudorange and carrier phase, and their observation equations are expressed as:

[0017]

[0018] In the formula, are the carrier phase measurement value and pseudorange measurement value respectively, is the true distance between the satellite and the receiver, and They represent ionospheric delay and tropospheric delay respectively, c represents the speed of light in vacuum, τi and τ j are the satellite clock error and the receiver clock error, is the ambiguity of the carrier phase, d is the instrument bias code delay between the satellite and the receiver, b is the instrument bias phase advance between the satellite and the receiver, ε is the measurement residual; i, k, j represent the satellite serial number, the receiver serial number and frequency respectively.

[0019] Furthermore, in S3, STEC is used to reconstruct the IED distribution during the inversion process of CIT, where STEC is calculated by measuring the modulation of code and carrier phase using a dual-frequency GPS receiver and is expressed as:

[0020]

[0021] In the formula, f represents the frequency of the carrier, c represents the speed of light, and P 4,sm Represents the pseudorange observation value of the phase smoothing code, DCB j and DCB i Represent the differential code deviations of the satellite and the GPS receiver respectively.

[0022] Furthermore, in S5 and S6, IRI2016 and SIRT are selected as the model initial value and inversion model respectively; the SIRT algorithm is expressed as:

[0023]

[0024] In the formula, is the ray-corrected electron density value of the jth voxel after k+1 iterations, P is the number of voxels traversed by the i-th ray, where 0≤P≤n, and λ is the relaxation factor 0≤λ≤1; y i is the STEC of the ith ray, A i,m is the intercept of the i-th ray in the m-th voxel, Δ is the correction of the i-th ray path, and W is the weight assigned to the j-th voxel for the TEC correction of the i-th ray.

[0025] Furthermore, in S9, the ROI index is expressed as:

[0026]

[0027] In the formula, It is the tth i+1 The electron density at time Yes i The electron density at the moment, ΔT is the time interval between the two moments; ROI is the rate of change in the time period, 10≤ΔT≤15 minutes; and ROI≥+μ / ROI≤-μ is defined as the occurrence of TID;

[0028] When inverting the three-dimensional ionospheric electron density, the spatial distribution of the electron density is generally represented by geographic coordinates. In order to better determine the center of the disturbance, the geographic coordinates are converted to geocentric coordinates. After solving the TID through the weighted centroid method, the conversion formula is as follows:

[0029]

[0030] Auxiliary formula 1:

[0031] b=a(1-f),

[0032] Auxiliary formula 2:

[0033]

[0034] Wherein, lon, lat, alt represent longitude, latitude and altitude respectively; X, Y, Z are the corresponding geocentric coordinates; flattening f = 1 / 298.257223563, major semi-axis a = 6378137.0 meters, b is the minor semi-axis; e is the first eccentricity, e′ is the second eccentricity, and N is the radius of the geocentric circle.

[0035] In order to determine the center of the TID disturbance, the weighted centroid method is used to calculate the geometric center of the TID. The calculation method is as follows:

[0036]

[0037] Where XYZ is the location where the disturbance occurs. In this invention, we define ROI as (x) ,ROI (y) ,ROI (z) Equal to max{|ROI (x) |,|ROI (y) |,ROI (z) |}, which means that X, Y, and Z are the positions with the largest absolute values ​​of the corresponding ROI, but this only indicates the size in each direction. The weighted centroid rule is to further determine the geometric center of the disturbance; finally, the geocentric coordinates are converted into geographic coordinates.

[0038] The beneficial effects of the present invention are:

[0039] The monitoring method of the present invention is different from the methods of identifying and monitoring TID proposed by other researchers. The present invention identifies and monitors TID from a three-dimensional spatial level, which is more in line with the spatiality of TID. Based on reliable electron density inversion, the present invention identifies and monitors TID by calculating the changes in spatial electron density at different times, and can effectively monitor the three-dimensional spatial position, change rate and trend of TID.

[0040] Compared with the traditional analytical index for identifying and monitoring ionospheric disturbances, the index proposed in the present invention monitors and identifies ionospheric disturbances in three-dimensional space. It can not only identify ionospheric disturbances, but also analyze the location of the disturbance in space. It can effectively realize the precise identification and monitoring of space ionospheric disturbances, and help researchers in the fields of precision positioning, navigation services, and aerospace to better understand the ionospheric disturbance mechanism. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0042] Figure 1 It is the overall calculation process of the ionospheric disturbance identification and monitoring index ROI of the present invention;

[0043] Figure 2 is the distribution of ground-based GNSS receiving stations in the experimental area of ​​the present invention;

[0044] Figure 3 It is the three-dimensional electron density inversion distribution at 21:30, 22:00, 23:00 and 24:00 on June 22, 2015;

[0045] Figure 4 The experimental area of ​​the present invention has longitude and latitude of 114° and 22°, an altitude range of 90km to 1000km, and a time period of 21:30 on June 22, 2015 to 00:00 on June 23, 2015, and the electron density changes under the ionospheric traveling wave disturbance;

[0046] Figure 5 It is the spatial distribution of the ionospheric traveling wave disturbance index ROI≥+μ / ROI≤-μ at different times of the present invention. Figures a to f are the three-dimensional spatial ROI value distribution of the experimental area from 21:30 on June 22, 2015 to 00:00 on June 23, 2015. DETAILED DESCRIPTION

[0047] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0048] Example 1

[0049] This embodiment describes the principle and technical solution of the present invention in detail with reference to the accompanying drawings and implementation examples; the details are as follows:

[0050] like Figure 1 As shown, the present invention combines the ionospheric three-dimensional electron density inversion method and the spatial characteristics of ionospheric disturbances to provide a three-dimensional ionospheric disturbance identification and monitoring index, and the specific steps include the following:

[0051] Step A, the experimental area is Hong Kong, and the Hong Kong area is discretized, that is, the ionospheric space from 90km to 1000km is divided into 0.5°×0.5°×30km. Note: In order to better invert the three-dimensional electron density in the Hong Kong area, this experiment will appropriately expand the inversion area, and the final inversion longitude and latitude range is 113° to 115° and 20° to 23°.

[0052] Step B, obtain the ground-based GNSS observation value files of the observation area from the Hong Kong Satellite Positioning Reference Station Network and obtain the GPS navigation message files, precise ephemeris and differential code deviation files from the IGS website.

[0053] Step C, calculate the STEC observation value based on the GPS observation data of the ground station. The GNSS signals needed in this step mainly include two common observation quantities: pseudorange and carrier phase. The observation equation can be expressed as:

[0054]

[0055] In the formula, are the carrier phase measurement value and pseudorange measurement value respectively, is the true distance between the satellite and the receiver, and They represent ionospheric delay and tropospheric delay respectively, c represents the speed of light in vacuum, τ i and τ j are the satellite clock error and the receiver clock error, is the carrier phase ambiguity, d is the instrument bias code delay between the satellite and the receiver, b is the instrument bias phase advance between the satellite and the receiver, and ε is the measurement residual. i, k, j represent the satellite serial number, the receiver serial number and the frequency, respectively.

[0056] The STEC of the GPS signal path is used to reconstruct the IED distribution in the inversion process of CIT, where STEC can be calculated by measuring the modulation of code and carrier phase with a dual-frequency GPS receiver and is expressed as

[0057]

[0058] Where f represents the frequency of the carrier, c is the speed of light, P 4,smis the pseudorange observation value of the phase smoothing code, DCB j and DCB i They are the differential code deviations of the satellite and the GPS receiver respectively. The unit of STEC is TECU, 1TECU = 10 16 el / m 2 .

[0059] Step D, the discretized ionospheric three-dimensional grid and STEC rays are combined to solve the intercepts of the rays passing through the voxels in the experimental area to form an intercept matrix A.

[0060] Step E: Extract the electron density of the experimental area from 21:30 on June 22, 2015 to 00:00 on June 23, 2015 from the International Ionospheric Model IRI2016, and use it as the initial value in the iterative algorithm.

[0061] Step F, combining the intercept matrix solved in step D and the initial value of electron density extracted in step E, uses the simultaneous iterative algorithm SIRT to calculate the three-dimensional electron density of the experimental area. The SIRT algorithm is expressed as:

[0062]

[0063] in is the ray-corrected electron density value of the jth voxel after k+1 iterations, P is the number of voxels traversed by the i-th ray, where 0≤P≤n and λ is the relaxation factor 0≤λ≤1. i is the STEC of the ith ray, A i,m is the intercept of the i-th ray in the m-th voxel, Δ is the correction of the i-th ray path, and W is the weight assigned to the j-th voxel for the TEC correction of the i-th ray.

[0064] Step G, analyzing the reliability of the three-dimensional electron density inverted in step F.

[0065] Step H, repeat the above steps to obtain the three-dimensional electron density in multiple time periods.

[0066] Step I: Arrange the three-dimensional electron density grids in the multi-period inversion area in a certain order. First, use the (Rate Of IDE) ROI index calculation in the vertical direction to determine the height position where TID occurs. Similarly, calculate the ROI in the horizontal direction. The calculation of ROI is as follows:

[0067]

[0068] in It is the tth i+1 The electron density at time Yes iThe electron density at the moment, ΔT is the time interval between the two moments; ROI is the rate of change in this time period, generally 10≤ΔT≤15 minutes, and ΔT=15 in this experiment; based on experimental experience, ROI≥3.4562e+07 or ROI≤-3.4562e+07 is considered to be TID.

[0069] When inverting the three-dimensional ionospheric electron density, the spatial distribution of the electron density is generally represented by geographic coordinates. In order to better determine the center of the disturbance, in the present invention, we convert the geographic coordinates into geocentric coordinates; the conversion formula is as follows:

[0070]

[0071] Auxiliary formula 1:

[0072] b=a(1-f),

[0073] Auxiliary formula 2:

[0074]

[0075] Wherein, lon, lat, alt represent longitude, latitude and altitude respectively; X, Y, Z are the corresponding geocentric coordinates; flattening f = 1 / 298.257223563, major semi-axis a = 6378137.0 meters, b is the minor semi-axis; e is the first eccentricity, e′ is the second eccentricity, and N is the radius of the geocentric circle.

[0076] In order to determine the center of the TID disturbance, the present invention uses a weighted centroid method to calculate the geometric center of the TID, and the calculation method is as follows:

[0077]

[0078] Where XYZ is the location where the disturbance occurs. In this invention, we define ROI (x) ,ROI (y) ,ROI (z) Equal to max{|ROI (x) |,|ROI (y) |,|ROI (z) |}, which means that X, Y, and Z are the positions with the largest absolute values ​​of the corresponding ROI, but this only indicates the size in each direction. The weighted centroid rule is to further determine the geometric center of the disturbance (Note: in this experiment, the geocentric coordinates are finally converted into geographic coordinates).

[0079] Figure 2The distribution of ground-based GNSS receiving stations in the experimental area is shown. Since the inversion of three-dimensional electron density will be limited by the number of rays if the area is too small, the inversion area is appropriately expanded in this experiment.

[0080] Figure 3 The three-dimensional electron density inversion distribution at four times on June 22, 2015, 21:30, 22:00, 23:00 and 24:00, is shown. Figure 3 It can be clearly seen from (a) to (d) that the electron density in the three-dimensional ionosphere changes under the influence of TID, and this change is irregular, especially at the altitude of 260km, where the influence of this disturbance is more obvious.

[0081] Figure 4 The figure shows the change of vertical electron density of the ionospheric sounding instrument HKKP from 21:30 on June 22, 2015 to 00:00 on June 23. It can be seen from the figure that the electron density value at the altitude of 230km has been increasing since 21:30, and the change has slowed down until 23:45.

[0082] Figure 5 The ROI value of the experimental area is displayed, that is, the rate of change of the three-dimensional electron density in the area within 15 minutes. Figure 5 From (a) to (c), it can be found that the TID mainly occurs at the altitude of 230 km, and the change is more obvious near longitude 115° and latitude 23°. In addition, from 5 (a) to (c), it can be clearly seen that the change gradually spreads from the altitude of 230 km to the upper and lower layers.

[0083] from Figure 5 (d) The center of the TID gradually moved to the 300km altitude layer, and the range of the disturbance became wider. This disturbance has gradually caused the electron density to change from 120km to 420km from the center of about 230km. Figure 5 From (a) to (d), we find that the center of the TID not only gradually transitions from the 230km altitude layer to the 300km altitude layer, but also gradually transitions to the lower latitudes. Figure 5 From (d) to (f), we find that the influence of TID gradually causes changes in the surrounding electron density, especially from Fig. Figure 5 From (e) to (f), we can find that this disturbance is no longer limited to a certain area, but gradually spreads to the entire ionosphere space.

Claims

1. A novel ionospheric three-dimensional disturbance monitoring method, characterized in that: The following steps are involved: S1. Discretize the area to be inverted spatially, that is, divide the area to be inverted into grids of a certain resolution in the form of grids according to the latitude, longitude and altitude directions, and assign numbers to the discretized grids; S2. Obtain the ground-based GNSS observation value files, GPS navigation message files, precise ephemeris and differential code deviation files in the inversion area; S3, solving the tilted total electron content STEC according to the observation file in S2; S4, determining the STEC effective rays passing through the voxels in the study area, and constructing the intercept matrix A through the intercepts of the effective rays and the voxels; S5. Using the electron density solved by the IRI2016 model as the initial value, determine the initial value of the electron density of each grid in the inversion area; S6. Using the SART iterative algorithm to solve the three-dimensional ionospheric electron density in the area to be inverted, thereby obtaining the electron density distribution in the three-dimensional grid space of the ionosphere; S7, the reconstructed 3D IED was further verified using the electron density profile of the ionospheric radiosonde to ensure the reliability of the inversion results; S8. Based on the previous results, further invert multiple time periods in the inversion area and verify the accuracy of the inversion results; S9, arranging the three-dimensional electron density grids of the multi-period inversion area in a certain order, calculating the maximum or minimum value of the ROI in the three directions of XYZ, and then using the weighted centroid method to determine the center of the TID disturbance; The ROI index is expressed as: In the formula, It is the tth i+1 The electron density at time Yes i The electron density at the moment, ΔT is the time interval between the two moments; ROI is the rate of change in the time period, 10≤ΔT≤15 minutes; and ROI≥+μ / ROI≤-μ is defined to define the occurrence of TID; When inverting the three-dimensional ionospheric electron density, the spatial distribution of the electron density is represented by geographic coordinates; the geographic coordinates are converted to geocentric coordinates, and the conversion formula after solving to TID by the weighted centroid method is as follows: Auxiliary formula 1: Auxiliary formula 2: Where lon, lat, alt represent longitude, latitude and altitude respectively; X, Y, Z are the corresponding geocentric coordinates; flattening f = 1 / 298.257223563, major semi-axis a = 6378137.0 meters, b is the minor semi-axis; e is the first eccentricity, e′ is the second eccentricity, and N is the radius of the azimuth circle; In order to determine the center of the TID disturbance, the weighted centroid method is used to calculate the geometric center of the TID. The calculation method is as follows: Where XYZ is the location where the disturbance occurs, and ROI is defined as (x) ,ROI (y) ,ROI (z) Equal to max{|ROI (x) |,ROI (y) |,ROI (z) |}; At this time, X, Y, and Z correspond to the maximum position of the absolute value of ROI; Finally, the geocentric coordinates are converted into geographic coordinates.

2. A novel ionospheric three-dimensional disturbance monitoring method according to claim 1, characterized in that: In S3, the GNSS signals needed to solve STEC mainly include two common observation quantities: pseudorange and carrier phase. The observation equation is expressed as: In the formula, are the carrier phase measurement value and pseudorange measurement value respectively, is the true distance between the satellite and the receiver, and They represent ionospheric delay and tropospheric delay respectively, c represents the speed of light in vacuum, τ i and τ j are the satellite clock error and the receiver clock error, is the ambiguity of the carrier phase, d is the instrument bias code delay between the satellite and the receiver, b is the instrument bias phase advance between the satellite and the receiver, ε is the measurement residual; i, k, j represent the satellite serial number, the receiver serial number and frequency respectively.

3. The novel ionospheric three-dimensional disturbance monitoring method according to claim 1 is characterized in that: In S3, STEC is used to reconstruct the IED distribution during the inversion process of CIT, where STEC is calculated by measuring the modulation of code and carrier phase using a dual-frequency GPS receiver and is expressed as: In the formula, f represents the frequency of the carrier, c represents the speed of light, and P 4,sm Represents the pseudorange observation value of the phase smoothing code, DCB j and DCB i Represent the differential code deviations of the satellite and the GPS receiver respectively.

4. The novel ionospheric three-dimensional disturbance monitoring method according to claim 1 is characterized in that: In S5 and S6, IRI2016 and SIRT are selected as the model initial value and inversion model respectively; the SIRT algorithm is expressed as: In the formula, is the ray-corrected electron density value of the jth voxel after k+1 iterations, P is the number of voxels traversed by the i-th ray, where 0≤P≤n, λ is the relaxation parameter 0≤λ≤1; y i is the STEC of the ith ray, A i,m is the intercept of the i-th ray in the m-th voxel, Δ is the correction of the i-th ray path, and W is the weight assigned to the j-th voxel for the TEC correction of the i-th ray.

Images