A method for identifying global small and medium scale ionospheric disturbances using GPS observation data

Abstract

The application discloses a method for identifying global small and medium scale ionospheric disturbances by using GPS observation data, which comprises the following steps: step 1, reading GPS satellite ephemeris file data, calculating satellite position information, satellite elevation angle and azimuth angle, and calculating the piercing point longitude and latitude according to the satellite elevation angle, azimuth angle, ground receiving point longitude and latitude and ionospheric piercing point height; step 2, calculating the slant path relative phase TEC according to the GPS double-frequency receiver carrier phase observation data; step 3, calculating the difference of the slant path relative phase TEC adjacent observation time, eliminating abnormal values, and converting the slant path difference value into the vertical direction difference value; and step 4, calculating the average value of the absolute value of the vertical direction difference value of each position point in the world, and taking the average value as the representation index of the global small and medium scale ionospheric disturbance degree. The method disclosed by the application is used as a guarantee index for satellite navigation positioning and other space applications, and provides a reference for performance evaluation of related electronic information system equipment.

Description

Technical Field

[0001] This invention belongs to the field of ionospheric physics research, and specifically relates to a method for identifying small-scale disturbances in the global ionosphere using GPS observation data. Background Technology

[0002] For users of GPS, BeiDou, and other radio satellite navigation and positioning systems, ionospheric delay error is one of the main sources of error affecting navigation and positioning accuracy. Providing users with real-time early warning information about ionospheric anomalies during periods of abnormal ionospheric changes is crucial. Total Electron Content (TEC) is an important parameter describing the morphology and structure of the ionosphere. Theoretically, TEC and its spatiotemporal variations can reflect the main characteristics of the ionosphere. Therefore, by detecting and analyzing TEC parameters, we can study ionospheric physical processes at different spatiotemporal scales, including various ionospheric variation characteristics at synoptic and climatological scales. In engineering applications, ionospheric TEC parameters are closely related to the time delay and phase delay of radio waves propagating through the ionosphere and can be used for radio wave propagation correction in space applications such as satellite positioning and navigation. Similarly, ionospheric TEC is also an important environmental factor of concern for systems such as ground-based telemetry and control radar and space-based SAR radar.

[0003] During periods of intense ionospheric disturbance, electron density undergoes significant fluctuations, substantially impacting radio information systems such as communication, navigation, radar, and telemetry. Increases or decreases in the ionospheric TEC (Electrode Temperature Coefficient) can increase the error in ionospheric delay correction models for navigation and positioning systems. Increased electron density gradients disrupt the spatiotemporal correlation of the ionosphere, thus affecting systems like GPS to varying degrees. Besides large-scale disturbances during ionospheric storms, the ionosphere also contains medium- and small-scale inhomogeneous structures. The generation and development of these inhomogeneities cause rapid fluctuations in radio signals passing through them, resulting in short-term irregular changes in the amplitude and phase of radio signals—a phenomenon known as ionospheric scintillation. Ionospheric scintillation often leads to deep fading and distortion of radio signals received by ground receivers, causing satellite communication failures, increased bit error rates, multipath effects, and reduced satellite navigation accuracy. Summary of the Invention

[0004] The technical problem to be solved by this invention is to provide a method for identifying global small- and medium-scale ionospheric disturbances using GPS observation data. In view of the shortcomings of existing technologies in identifying global small- and medium-scale ionospheric disturbances, this invention introduces the Global Ionospheric TEC Disturbance Index, which represents small- and medium-scale disturbances, based on globally distributed IGS GPS observation data. This can provide a simple and intuitive measure of global ionospheric disturbance changes.

[0005] The present invention adopts the following technical solution:

[0006] An improved method for identifying global small- and medium-scale ionospheric disturbances using GPS observation data includes the following steps:

[0007] Step 1: Read GPS satellite ephemeris file data, calculate satellite position information, satellite elevation angle, and azimuth angle. Based on the satellite elevation angle, azimuth angle, latitude and longitude of the ground receiving point, and the ionospheric puncture point height, calculate the latitude and longitude of the puncture point:

[0008] Step 11: Calculate the satellite position based on the ephemeris data. In the geocentric coordinate system, the satellite position coordinates are:

[0009]

[0010] In the above formula, r is the distance from the Earth's center to the satellite, Ω er ω is the angle between the Greenwich Meridian and the ascending node direction, v is the horizontal anomaly angle, ω is the angle between the perigee direction and the ascending node direction of the satellite orbital plane, and i is the angle between the satellite orbital plane and the equatorial plane.

[0011] Step 12: For a point P on the Earth's surface, let P be the origin of the station-centered coordinate system. Draw a plane tangent to the Earth's ellipsoid through P. Take due north as the Y-axis, due east as the X-axis, and the Z-axis pointing towards the normal direction. In the station-centered coordinate system, the formulas for calculating the satellite elevation angle α and azimuth angle β are:

[0012]

[0013] Step 13, the direction cosine DC from the user to the satellite in the geocentric coordinate system is expressed as:

[0014]

[0015] Where (x) s ,y s ,z s ), (x u ,y u ,z u The positions of the satellite and the ground receiving point are respectively in the geocentric coordinate system. Transform them to the station-centered coordinate system:

[0016] DC ENU =R e2t DC (4)

[0017] R e2t The rotation matrix from the geocentric coordinate system to the station-centered coordinate system is expressed as:

[0018]

[0019] In the above formula, λ is the longitude of the ground receiving point, and φ is the latitude of the ground receiving point, in the DC coordinate system. ENU Represented as (κ) e ,κ n ,κ u ), calculate (κ) according to formula (4) e ,κ n ,κ u Then, according to formula (2), the satellite elevation angle α and azimuth angle β are calculated;

[0020] Step 14: Calculate the latitude and longitude (λ) of the puncture point using the satellite elevation angle α, azimuth angle β, latitude and longitude of the ground receiving point, and the altitude of the ionospheric puncture point. PP ,φ PP The formula for ) is:

[0021] λ PP =λ+arcsin(sin(angle)sinβ / cosφ)

[0022] φ PP =arcsin(sinφcos(angle)+cosφsin(angle)cosβ)

[0023] In the above formula, R = 6370 km is the Earth's radius, and h is the height of the ionospheric puncture point;

[0024] Step 2: Calculate the slant path relative phase (TEC) based on the GPS dual-frequency receiver carrier phase observation data.

[0025]

[0026]

[0027]

[0028] ρ j It is the sum of the geometric distance from the GPS dual-frequency receiver to the satellite, the clock bias of the GPS dual-frequency receiver, the clock bias of the satellite, and the tropospheric delay error;

[0029] k = 1, 2, which are the sum of the phase ambiguity of the two frequencies, the hardware phase delay of the GPS dual-frequency receiver, and the hardware phase delay of the satellite, respectively;

[0030] The ionospheric delay along the path from the GPS dual-frequency receiver to the satellite;

[0031] The expression for the ionospheric TEC obtained from the above formula is as follows:

[0032] Wherein, λ1 and λ2 are the wavelengths of the two carrier waves of the GPS signal;

[0033] Step 3: Calculate the difference between adjacent observation times of the relative phase TEC along the oblique path, remove outliers, and convert the oblique path difference values ​​into vertical difference values:

[0034] Step 31: The time resolution of the GPS observation data is 15 seconds. Calculate the difference in phase TEC between two consecutive sets of data. The calculation formula is as follows:

[0035]

[0036] Step 32: Remove outliers caused by phase ambiguity. If |dtec|>=0.3 is considered an outlier, remove this outlier.

[0037] Step 33: Convert the difference values ​​of the oblique path into difference values ​​in the vertical direction. The calculation formula is as follows:

[0038] Vdtecj=dtecj*cos(α)

[0039] Step 4: Calculate the average of the absolute values ​​of the vertical differences at all locations globally. Use this average as a characterization index of the degree of ionospheric disturbance at small and medium scales globally. The calculation formula is as follows:

[0040]

[0041] In the above formula, N represents the number of puncture points.

[0042] The beneficial effects of this invention are:

[0043] This invention discloses a method for identifying small-scale disturbances in the global ionosphere using GPS observation data. Based on globally distributed IGS GPS observation data, it introduces the Global Ionospheric TEC Disturbance Index (GTEC), which represents both large-scale and small-scale disturbances. This provides a simple and intuitive measure of global ionospheric disturbance changes. Users of different electronic information systems can compare the impact on their system performance with the GTEC. Based on the comparison results, the degree of impact of ionospheric disturbances on information systems can be obtained through the GTEC. Therefore, the GTEC can serve as an important indicator for studying ionospheric changes and predicting ionospheric disturbance states. It can also be used as a guarantee indicator for space applications such as satellite navigation and positioning, providing a reference for the performance evaluation of related electronic information system equipment. Attached Figure Description

[0044] Figure 1 This is a schematic diagram showing the positional relationship between the satellite and the ground receiving point;

[0045] Figure 2 This is a flowchart of the method of the present invention;

[0046] Figure 3 This is a graph showing the ionospheric storm index and the geomagnetic activity index Kp during the severe geomagnetic storm event on March 19, 2015.

[0047] Figure 4 This is a graph showing the ionospheric storm index and the geomagnetic activity index Kp during a moderate geomagnetic storm on July 16, 2017. Detailed Implementation

[0048] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0049] Example 1 discloses a method for identifying global small- and medium-scale ionospheric disturbances using GPS observation data, such as... Figure 2 As shown, it includes the following steps:

[0050] Step 1: Read GPS satellite ephemeris file data, calculate satellite position information, satellite elevation angle, and azimuth angle. Based on the satellite elevation angle, azimuth angle, latitude and longitude of the ground receiving point, and the ionospheric puncture point height, calculate the latitude and longitude of the puncture point:

[0051] Step 11: Calculate the satellite position based on the ephemeris data. In the geocentric coordinate system, the satellite position coordinates are:

[0052]

[0053] In the above formula, r is the distance from the Earth's center to the satellite, Ω er ω is the angle between the Greenwich Meridian and the ascending node (the satellite's orbital plane intersects the Earth's equatorial plane on a straight line; the ascending node is defined as the point along this line where the satellite crosses the equatorial plane to the north), v is the mean anomaly angle, ω is the angle between the perigee direction and the ascending node direction, and i is the angle between the satellite's orbital plane and the equatorial plane; the satellite's position relative to the Earth is as follows: Figure 1 As shown.

[0054] Step 12: For a point P on the Earth's surface, let P be the origin of the station-centered coordinate system ENU. Draw a plane tangent to the Earth's ellipsoid through point P. Take due north as the Y-axis, due east as the X-axis, and the Z-axis pointing in the direction of the normal. In the station-centered coordinate system ENU, the formulas for calculating the satellite elevation angle α and azimuth angle β are:

[0055]

[0056] Step 13, the direction cosine DC from the user to the satellite in the geocentric coordinate system is expressed as:

[0057]

[0058] Where (x) s ,y s ,z s ), (x u ,y u ,z u The positions of the satellite and the ground receiving point are respectively in the geocentric coordinate system (GCS) and converted to the station-centered coordinate system (ENU).

[0059] DC ENU =R e2t DC (4)

[0060] R e2t The rotation matrix from the geocentric coordinate system ECEF to the station-centered coordinate system ENU is expressed as:

[0061]

[0062] In the above formula, λ is the longitude of the ground receiving point, and φ is the latitude of the ground receiving point, in the ENU coordinate system DC. ENU Represented as (κ) e ,κ n ,κ u ), calculate (κ) according to formula (4) e ,κ n ,κ u Then, according to formula (2), the satellite elevation angle α and azimuth angle β are calculated;

[0063] Step 14, the positional relationship between the satellite and the ground receiving point is as follows: Figure 1 As shown, the latitude and longitude of the puncture point (λ) are calculated from the satellite elevation angle α, azimuth angle β, latitude and longitude of the ground receiving point, and the height of the ionospheric puncture point. PP ,φ PP The formula for ) is:

[0064] λ PP =λ+arcsin(sin(angle)sinβ / cosφ)

[0065] φ PP =arcsin(sinφcos(angle)+cosφsin(angle)cosβ)

[0066] In the above formula, R = 6370 km is the Earth's radius, and h is the height of the ionospheric puncture point;

[0067] Step 2: Calculate the slant path relative phase (TEC) based on the GPS dual-frequency receiver carrier phase observation data.

[0068] The GPS dual-frequency receiver provides carrier phase information for the j-th satellite at two frequencies, f1 and f2. The observation equation is:

[0069]

[0070]

[0071] ρ j It is the sum of the geometric distance from the GPS dual-frequency receiver to the satellite, the clock bias of the GPS dual-frequency receiver, the clock bias of the satellite, and the tropospheric delay error;

[0072] k = 1, 2, which are the sum of the phase ambiguity of the two frequencies, the hardware phase delay of the GPS dual-frequency receiver, and the hardware phase delay of the satellite, respectively;

[0073] Ionospheric delay (in meters) along the path from the GPS dual-frequency receiver to the satellite.

[0074] The expression for the ionospheric TEC obtained from the above formula is as follows:

[0075] Wherein, λ1 and λ2 are the wavelengths of the two carrier waves of the GPS signal;

[0076] Step 3: Calculate the difference between adjacent observation times of the relative phase TEC along the oblique path, remove outliers, and convert the oblique path difference values ​​into vertical difference values:

[0077] Step 31: The time resolution of the GPS observation data is 15 seconds. Calculate the difference in phase TEC between two consecutive sets of data. The calculation formula is as follows:

[0078]

[0079] Step 32: Remove outliers caused by phase ambiguity. Values ​​with |dtec| >= 0.3 are considered outliers and are removed. Values ​​with |dtec| >= 0.3 are not considered when calculating the characterization index of global small- and medium-scale ionospheric disturbance.

[0080] Step 33: Convert the difference values ​​of the oblique path into difference values ​​in the vertical direction. The calculation formula is as follows:

[0081] Vdtecj=dtecj*cos(α)

[0082] Step 4: Calculate the average of the absolute values ​​of the vertical differences at all locations globally. Use this average as a characterization index of the degree of ionospheric disturbance at small and medium scales globally. The calculation formula is as follows:

[0083]

[0084] In the above formula, N represents the number of puncture points.

[0085] To verify the effectiveness of this patent in detecting ionospheric storm events, this embodiment provides a global distribution map of small-scale ionospheric disturbances during the strong geomagnetic storm event on March 19, 2015. During this geomagnetic storm event, the minimum geomagnetic activity index Dst reached -228 nT, and the kp index reached 9. Figure 3 The variation of the small-scale perturbation index of the global ionosphere with the geomagnetic activity index kp is given.

[0086] This embodiment also provides a global distribution map of small-scale ionospheric disturbances during the moderate geomagnetic storm event on July 16, 2017. During this geomagnetic storm event, the minimum geomagnetic activity index Dst reached -69nT, and the kp index reached 6. Figure 4 The variation of the small-scale perturbation index of the global ionosphere with the geomagnetic activity index kp is given.

Claims

1. A method for identifying global small- and medium-scale ionospheric disturbances using GPS observation data, characterized in that, Includes the following steps: Step 1: Read GPS satellite ephemeris file data, calculate satellite position information, satellite elevation angle, and azimuth angle. Based on the satellite elevation angle, azimuth angle, latitude and longitude of the ground receiving point, and the ionospheric puncture point height, calculate the latitude and longitude of the puncture point: Step 11: Calculate the satellite position based on the ephemeris data. In the geocentric coordinate system, the satellite position coordinates are: In the above formula, r is the distance from the Earth's center to the satellite, Ω er ω is the angle between the Greenwich Meridian and the ascending node direction, v is the horizontal anomaly angle, ω is the angle between the perigee direction and the ascending node direction of the satellite orbital plane, and i is the angle between the satellite orbital plane and the equatorial plane. Step 12: For a point P on the Earth's surface, let P be the origin of the station-centered coordinate system. Draw a plane tangent to the Earth's ellipsoid through P. Take due north as the Y-axis, due east as the X-axis, and the Z-axis pointing towards the normal direction. In the station-centered coordinate system, the formulas for calculating the satellite elevation angle α and azimuth angle β are: Step 13, the direction cosine DC from the user to the satellite in the geocentric coordinate system is expressed as: Where (x) s ,y s ,z s ), (x u ,y u ,z u These represent the positions of the satellite and the ground receiving point in the Earth-centered, Earth-fixed coordinate system, respectively. Transform it to the station-centered coordinate system: DC ENU =R e2t ·DC (4) R e2t The rotation matrix from the geocentric coordinate system to the station-centered coordinate system is expressed as: In the above formula, λ is the longitude of the ground receiving point, and φ is the latitude of the ground receiving point, in the DC coordinate system. ENU Represented as (κ) e ,κ n ,κ u ), calculate (κ) according to formula (4) e ,κ n ,κ u Then, according to formula (2), the satellite elevation angle α and azimuth angle β are calculated; Step 14: Calculate the latitude and longitude (λ) of the puncture point using the satellite elevation angle α, azimuth angle β, latitude and longitude of the ground receiving point, and the altitude of the ionospheric puncture point. PP ,φ PP The formula for ) is: l PP =λ+arcsin(sin(angle)sinβ / cosφ) f PP =arcsin(sinφcos(angle)+cosφsin(angle)cosβ) In the above formula, R = 6370 km is the Earth's radius, and h is the height of the ionospheric puncture point; Step 2: Calculate the slant path relative phase (TEC) based on the GPS dual-frequency receiver carrier phase observation data. The GPS dual-frequency receiver provides carrier phase information for the j-th satellite at two frequencies, f1 and f2. The observation equation is: ρ j It is the sum of the geometric distance from the GPS dual-frequency receiver to the satellite, the clock bias of the GPS dual-frequency receiver, the clock bias of the satellite, and the tropospheric delay error; These are the sum of the phase ambiguity of the two frequencies, the hardware phase delay of the GPS dual-frequency receiver, and the hardware phase delay of the satellite, respectively. The ionospheric delay along the path from the GPS dual-frequency receiver to the satellite; The expression for the ionospheric TEC obtained from the above formula is as follows: Wherein, λ1 and λ2 are the wavelengths of the two carrier waves of the GPS signal; Step 3: Calculate the difference between adjacent observation times of the relative phase TEC along the oblique path, remove outliers, and convert the oblique path difference values ​​into vertical difference values: Step 31: The time resolution of the GPS observation data is 15 seconds. Calculate the difference in phase TEC between two consecutive sets of data. The calculation formula is as follows: Step 32: Remove outliers caused by phase ambiguity. If |dtec|>=0.3 is considered an outlier, remove this outlier. Step 33: Convert the difference values ​​of the oblique path into difference values ​​in the vertical direction. The calculation formula is as follows: Vdtec j =dtec j *cos(α) Step 4: Calculate the average of the absolute values ​​of the vertical differences at all locations globally. Use this average as a characterization index of the degree of ionospheric disturbance at small and medium scales globally. The calculation formula is as follows: In the above formula, N represents the number of puncture points.

Images