Method for evaluating influence of geomagnetic storm on aviation GNSS service precision based on ADS-B data

Abstract

The invention relates to the field of global satellite navigation system service quality evaluation, and particularly discloses a method for evaluating influence of geomagnetic storm on aviation GNSS service precision based on ADS-B data, which comprises the following steps: S1, acquiring two different content messages in ADS-B, namely, messages of state vectors SV and mode states MS of all ADS-B in a to-be-analyzed range, acquiring flight positions of all aircrafts in the concern area range and navigation precision categories NACp and vertical geometric precision GVA of the aircrafts at the positions through a fusion algorithm; s2, carrying out fusion processing on the two messages in the S1, and carrying out data alignment according to timestamps in the SV message and the MS message; and calculating the position information in the MS message according to the time information and the position information in the SV message through interpolation. According to the method for evaluating the influence of the geomagnetic storm on the precision of the aviation GNSS service based on the ADS-B data, the flight position and navigation precision parameters of the aircraft are adopted, and then the precision change of the GNSS service under different geomagnetic conditions is evaluated.

Description

Technical Field

[0001] The present invention relates to the field of global navigation satellite system service quality assessment, and particularly to a method for evaluating the influence of geomagnetic storms on the service accuracy of aviation GNSS based on ADS-B data. Background Art

[0002] With the development of space technology, the significant impact of space weather events on the air transportation system has attracted increasing attention, especially their interference with the aviation safety critical systems - Communication, Navigation and Surveillance (CNS). The International Civil Aviation Organization (ICAO) has clearly pointed out that space weather events such as geomagnetic storms can cause an increase in GNSS positioning errors, signal loss of lock or even interruption by disturbing the Earth's electromagnetic environment and ionospheric structure, thereby threatening aviation navigation and surveillance services that rely on GNSS. Research shows that there is a correlation between space weather events and accidents such as flight delays, runway excursions and air crashes, highlighting the urgency of evaluating the performance of airborne GNSS during geomagnetic storms.

[0003] Traditional research is mostly based on ground monitoring data. Although ground-based monitoring can reflect the influence of geomagnetic changes on GNSS, it cannot directly characterize the performance of airborne GNSS during actual flight. Although the airborne quick access recorder can record GNSS parameters, its installation rate is currently low globally and the data is not publicly available, and the broadband ground-air real-time data link with the ability to transmit navigation performance indicators has not been popularized, resulting in a lack of airborne GNSS quality data.

[0004] In the prior art, the Automatic Dependent Surveillance-Broadcast (ADS-B) system, as an open aviation surveillance system, broadcasts flight parameters through the 1090 MHz frequency band. Its positioning information depends on GNSS, providing the possibility for large-scale monitoring of airborne GNSS performance. Existing research has used ADS-B data to achieve GNSS radio frequency interference detection and positioning, but the evaluation research on the influence of geomagnetic storms is still in its infancy. Existing work such as Hospodka et al. analyzing solar flare events based on the MS reports of local area short-period ADS-B, and Xue et al. simulating the influence of geomagnetic storms through historical data, all have the defects of limited spatio-temporal range or methods being divorced from the actual scenario. Summary of the Invention

[0005] The object of the present invention is to provide a method for evaluating the influence of geomagnetic storms on the accuracy of aviation GNSS services based on ADS-B data, extracting the ADS-B state vectors and mode status reports within a specified time and geographical range and performing effective fusion, as well as dividing the area to be analyzed into grids and segmenting the data within the grids according to their timestamps; establishing a mathematical model for the ADS-B fusion data in each grid and each time period and calculating its accuracy index; considering that some aircraft use satellite-based augmentation navigation services and conducting grouped evaluations, and still being able to perform a relatively accurate evaluation in the case where there are differences in the positioning services adopted by different aircraft.

[0006] To achieve the above object, the present invention provides a method for evaluating the influence of geomagnetic storms on the accuracy of aviation GNSS services based on ADS-B data, comprising the following steps:

[0007] S1. Obtain two different types of message reports in ADS-B, namely, the state vector (SV) and mode status (MS) message reports of all ADS-B within the range to be analyzed, and obtain the flight positions of all aircraft within the area of interest and the navigation accuracy category (NACp) and vertical geometric accuracy (GVA) of the aircraft at that position through a fusion algorithm;

[0008] S2. Perform fusion processing on the two types of message reports in S1, align the data according to the timestamps in the SV message report and the MS message report; interpolate and calculate the position information in the MS message report according to the time information and position information in the SV message report;

[0009] S3. Divide the area to be analyzed into grids, fill the message reports obtained after fusion processing into the specified grids, and group the message reports in each grid according to time in groups of 5 minutes;

[0010] S4. Classify the aircraft using the satellite-based augmentation system (SBAS) and the standard position service (SPS) based on an unsupervised classification method;

[0011] S5. Based on the assumption that the horizontal dilution of precision (HFOM) follows a Rayleigh distribution, construct a cost function, and for the aircraft appearing within 5 minutes in a single grid, count the number N of ADS-B reports with NACp = k when k = {8, 9, 10, 11} k , according to N k estimate the short-term average value of HFOM

[0012] S6. Based on the assumption that the vertical dilution of precision (VFOM) follows a semi-normal distribution, construct a cost function, and for the aircraft appearing within 5 minutes in a single grid, count the number M of ADS-B reports with GVA = k when k = {1, 2} k , according to M k estimate the short-term average value of VFOM

[0013] Preferably, in S2, the interpolation method selects first-order linear interpolation or second-order interpolation for calculation.

[0014] Preferably, in S3, during the gridding process, due to the non-uniformity of ADS-B data in geographical distribution, the grid division should adopt a dynamic division method to ensure that each grid has sufficient data at each moment. The basic principle is to ensure that, at intervals of five minutes, the number of ADS-B data contained in each network in each time period is not less than 1000 during peak hours to ensure the accuracy of the parameter estimation process.

[0015] Preferably, in S4, the unsupervised classification method selects the K-means clustering algorithm to classify the data. The input of the clustering process is the historical NACp mean value of each aircraft. For the massive data across days, the average value is calculated day by day, and the mean value with a larger data volume is selected and input into the clustering process; aircraft with a mean value less than 8 are excluded in advance. For areas where ADS-B data is sparse or at times when there are fewer flights in the night airspace, unified parameter estimation can be performed for multiple five-minute intervals.

[0016] Preferably, the specific steps for classifying aircraft are as follows:

[0017] S41. Expand the clustering input dimension to multi-parameter joint features, including NACp, GVA, NIC, and SIL. Each parameter eliminates the dimension difference through normalization in advance;

[0018] S42. Replace the K-means clustering method with a Gaussian mixture model, set the covariance matrix type to a diagonal matrix, and constrain the number of clusters to 2 to distinguish SBAS / SPS service types;

[0019] S43. Perform feature weighting on the multi-dimensional joint features. Among them, the weight of the NACp parameter is 0.4 - 0.5, the weight of GVA is 0.25 - 0.35, and NIC and SIL are each 0.1 - 0.15; the sum of all weights is equal to 1;

[0020] S44. Dynamically generate SBAS / SPS discrimination rules based on the clustering results: when the cluster center satisfies NACp≥9.5, GVA≥1.8, SIL≥2.9, and NIC≥9, it is determined as the SBAS service-dominated cluster; the rest are determined as the SPS service-dominated cluster;

[0021] S45. Filter the boundary samples using the probability membership threshold, and only retain the samples with a membership greater than 80% to participate in the geomagnetic storm impact analysis.

[0022] Preferably, in S5, the cost function is established as the sum of four parts with the same form. When k = {1, 2, 3, 4}, the cost function l H (σ) is calculated as follows:

[0023]

[0024] Frequency N 12-k and the integral value of the theoretical Rayleigh distribution probability density function in the to interval form the log-likelihood function as follows:

[0025]

[0026] wherein, and are two sets of constants, representing the upper and lower bounds when the ADS-B airborne equipment maps the continuous parameter HFOM output by the GNSS receiver to different values of the discrete parameter NACp in the MS report of ADS-B.

[0027] Preferably, in S6, the cost function is established as the sum of two parts with the same form. The cost function l V (σ) is calculated as follows:

[0028]

[0029] Frequency M k and the integral value of the theoretical semi-normal distribution probability density function in the to interval form the log-likelihood function as follows:

[0030]

[0031] wherein, and respectively represent the upper and lower bounds when the ADS-B airborne equipment maps the continuous parameter VFOM output by the GNSS receiver to different values of the discrete parameter GVA in the MS report of ADS-B.

[0032] Therefore, the present invention adopts the above method for evaluating the influence of geomagnetic storms on the accuracy of aviation GNSS services based on ADS-B data, and the beneficial effects are as follows:

[0033] (1) By utilizing the navigation service accuracy parameters NACp and GVA broadcast by the Automatic Dependent Surveillance - Broadcast (ADS - B) system, the present invention can evaluate the impact of geomagnetic storms on the service accuracy of the airborne Global Navigation Satellite System (GNSS). By fusing the State Vector (SV) and Mode Status (MS) reports of ADS - B, the flight position and navigation accuracy parameters of the aircraft are extracted, and then the accuracy changes of GNSS services under different geomagnetic conditions are evaluated.

[0034] (2) Through ADS - B data, the present invention can directly obtain the navigation accuracy parameters of the aircraft during actual flight, thereby more accurately evaluating the impact of geomagnetic storms on aviation GNSS services.

[0035] (3) By using unsupervised classification methods such as K - means clustering or Gaussian mixture models, the aircraft using the Satellite - Based Augmentation System (SBAS) and Standard Positioning Service (SPS) are classified, and their navigation accuracies are evaluated respectively. This method can more accurately reflect the performance differences of different navigation services under geomagnetic storms. Description of the Drawings

[0036] Figure 1 is the overall flowchart of an embodiment of the method for evaluating the impact of geomagnetic storms on the accuracy of aviation GNSS services based on ADS - B data of the present invention;

[0037] Figure 2 is the change trend of the HFOM and VFOM estimated values near Calgary from April 22 to 25, 2023 in an embodiment of the method for evaluating the impact of geomagnetic storms on the accuracy of aviation GNSS services based on ADS - B data of the present invention;

[0038] Figure 3 is the global geomagnetic index K p 、ASY, SYM change trend between April 21 and 29, 2023. Detailed Embodiment

[0039] The technical solutions of the present invention will be further described below with reference to the drawings and embodiments.

[0040] Unless otherwise defined, the technical terms or scientific terms used in the present invention should have the ordinary meanings understood by those of ordinary skill in the field to which the present invention belongs.

[0041] Embodiment

[0042] As Figure 1 shown, a method for evaluating the impact of geomagnetic storms on the accuracy of aviation GNSS services based on ADS - B data includes the following steps:

[0043] S1. Obtain the geomagnetic index and ADS - B data. In this embodiment, Kaiyuan ADS - B data and the global geomagnetic index K are used as an example. During the implementation of the fusion algorithm, first align the timestamps of the SV report and the MS report in ADS - B, and perform one - dimensional linear interpolation based on the position information in the SV report and the timestamps of the MS and SV reports respectively. After interpolation, each fused MS report will correspond to a position information. p For example. During the implementation of the fusion algorithm, first align the timestamps of the SV report and the MS report in ADS - B, and perform one - dimensional linear interpolation based on the position information in the SV report and the timestamps of the MS and SV reports respectively. After interpolation, each fused MS report will correspond to a position information.

[0044] S2. Data pre - processing includes cleaning the ADS - B data broadcast by low - MACp aircraft, and marking the ADS - B data during abnormal geomagnetic index (K > 5) periods to avoid their participation in navigation service classification. p >5) during the ADS - B data for marking to avoid its participation in the navigation service classification.

[0045] S3. Perform grid processing on the area to be analyzed. Since the ADS - B data has non - uniformity in geographical distribution, the grid division should be carried out in a dynamic division method so that each individual network has sufficient data at each moment to implement parameter estimation. The grid division uses quadtree segmentation, and the specific steps are as follows:

[0046] S31. At the beginning of the algorithm, regard the entire area to be analyzed as 1 grid.

[0047] S32. If the number N of ADS - B messages contained in each grid in the current round is greater than the threshold N max , then the current rectangle is divided into 4 sub - regions along the mid - points of longitude and latitude.

[0048] S33. When the actual grid coverage area corresponding to the longitude - latitude network is less than 100 km 2 , stop the quadtree segmentation.

[0049] S34. If the data volume of adjacent sub - networks is less than the threshold N min , merge them into the parent network.

[0050] S4. Use 5 minutes as the basic interval for time segmentation. Depending on the data density of the analysis area, the segmentation interval can be appropriately adjusted.

[0051] S5. For the ADS - B data covered by each network in each time period, use K - means clustering. Use the mean value of the historical NACp of each aircraft as the input. For data across days, it can be calculated day by day, and select the larger mean value as the input for the clustering process. Aircraft with a mean value less than 8 are regarded as aircraft not using GNSS and are pre - excluded in advance. At the same time, based on the pre - input geomagnetic index and the change of the total electron content in the ionosphere, try to avoid inputting the NACp data corresponding to the time period with large geomagnetic or ionospheric changes during the classification process to avoid contaminating the clustering process.

[0052] S6. When the navigation accuracy parameters NACp and GVA, and the navigation integrity parameters NIC and SIL can all be obtained, the clustering input dimension can be extended to a multi-parameter joint feature, and GMM can be used to replace K-means clustering. The specific steps are as follows:

[0053] S61. Set the covariance matrix type to a diagonal matrix and constrain the number of clusters to 2 to distinguish SBAS / SPS service categories;

[0054] S62. Perform feature weighting on the multi-dimensional data, where the weight of the NACp parameter is 0.45, the weight of GVA is 0.25, and the weights of NIC and SIL are each 0.15; ensuring that the sum of all weights is equal to 1;

[0055] S63. Dynamically generate SBAS / SPS discrimination rules based on the clustering results: when the cluster center satisfies NACp≥9.5, GVA≥1.8, SIL≥2.9, and NIC≥9, it is determined as the SBAS service dominant cluster; the rest are determined as the SPS service dominant cluster;

[0056] S64. Filter the boundary samples using a probability membership threshold, and only retain the samples with a membership greater than 90% to participate in the geomagnetic storm impact analysis.

[0057] S7. Before estimating the HFOM and VFOM parameters, it is necessary to first check whether the distributions of NACp and GVA are pathological. If NACp or GVA is completely concentrated in one interval, parameter estimation cannot be performed. When estimation is not possible, the cumulative frequency is carried over to the next moment for estimation.

[0058] S8. Estimate HFOM and VFOM using the method of optimizing the maximum likelihood function. Assume that the positioning error of each aircraft follows an independent and identically distributed zero-mean Gaussian distribution in all directions. Then, the HFOM of all aircraft in a certain area over a period of time follows a Rayleigh distribution, and the VFOM follows a semi-normal distribution. The parameters to be estimated are the average HFOM of all aircraft in a certain area over a period of time and the average VFOM is The input parameter is N k represents the number of messages with NACp = k in the messages broadcast by all aircraft in a single grid within a short period of time through ADS-B, and M k represents the number of messages with GVA = k in the messages broadcast by all aircraft in a single grid within a short period of time through ADS-B.

[0059] The cost function of HFOM is calculated as follows:

[0060]

[0061] Frequency N 12-kThe log-likelihood function formed by the integral value of the theoretical Rayleigh distribution probability density function in the range from to is as follows:

[0062]

[0063] where and are two sets of constants, representing the upper and lower bounds when the ADS-B airborne equipment maps the continuous parameter HFOM output by the GNSS receiver to different values of the discrete parameter NACp in the MS report of ADS-B, as shown in Table 1:

[0064] Table 1 NACp mapping parameters

[0065]

[0066] The cost function of VFOM is calculated as follows:

[0067]

[0068] Frequency M k The log-likelihood function formed by the integral value of the theoretical semi-normal distribution probability density function in the range from to is as follows:

[0069]

[0070] where and respectively represent the upper and lower bounds when the ADS-B airborne equipment maps the continuous parameter VFOM output by the GNSS receiver to different values of the discrete parameter GVA in the MS report of ADS-B, as shown in Table 2:

[0071] Table 2 VFOM mapping parameters

[0072]

[0073] The effects of the present invention can be further illustrated by the following actual data analysis results. The data uses the open-source ADS-B data during the period from April 22 to 25, 2023 as input. There was a geomagnetic storm of G4 level during this period. In this example, the effectiveness of the present invention is verified by analyzing the relationship between the change of navigation accuracy broadcast by ADS-B of flights in the southern region of Canada during this geomagnetic storm and the change of the geomagnetic index K p of.

[0074] Data description: The geomagnetic data includes the K p , AYM, and SYM indices from April 21 to April 29, 2023, as Figure 2As shown. It can be seen that there were obvious changes in the geomagnetic index during the period from April 23 to 24. The ADS-B data included approximately 30 million ADS-B MS reports and SV reports in southern Canada (within the range of latitude 49°N to 60°N and longitude 141°W to 90°W) from April 22 to 25, 2023. After being processed by the fusion algorithm, approximately 10 million MS messages with position information were obtained. After being classified by the navigation service, approximately 55% of the messages were classified as SBAS services, and the remaining 45% were SPS. After the parameter estimation algorithm, the horizontal and vertical accuracy indicators HFOM and VFOM of the Calgary area grid SBAS and SPS are as Figure 3 shown.

[0075] From Figure 3 the comparison with Figure 2 it can be found that at about 21:00 in the evening of April 23, the accuracy of the aviation navigation service increased with the increase of the K p index, indicating that the aviation navigation accuracy is indeed affected by geomagnetic changes, and the method proposed in the present invention is effective.

[0076] Therefore, the present invention adopts the above-mentioned method for evaluating the influence of geomagnetic storms on the accuracy of aviation GNSS services based on ADS-B data. By fusing the status vector (SV) and mode status (MS) reports of ADS-B, the flight position and navigation accuracy parameters of the aircraft are extracted, and then the accuracy changes of GNSS services under different geomagnetic conditions are evaluated; through ADS-B data, the navigation accuracy parameters of the aircraft during actual flight can be directly obtained, so as to more accurately evaluate the influence of geomagnetic storms on aviation GNSS services; it can more accurately reflect the performance differences of different navigation services under geomagnetic storms.

[0077] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit them. Although the present invention has been described in detail with reference to the preferred embodiments, those of ordinary skill in the art should understand that they can still modify or equivalently replace the technical solutions of the present invention, and these modifications or equivalent replacements cannot make the modified technical solutions deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A method for evaluating the impact of geomagnetic storms on aviation GNSS service accuracy based on ADS-B data, characterized in that: The following steps are involved: S1. Obtain two different content messages in ADS-B, namely, the state vector SV and mode state MS messages of all ADS-Bs within the range to be analyzed, and obtain the flight positions of all aircraft within the area of ​​interest and the navigation accuracy category NACp and vertical geometric accuracy GVA of the aircraft at that position through a fusion algorithm; S2, merge the two messages in S1, align the data according to the timestamps in the SV message and the MS message; interpolate and calculate the location information in the MS message according to the time information and location information in the SV message; S3, gridding the area to be analyzed, filling the messages obtained after fusion processing into the specified grid, and grouping the messages in each grid according to time order and processing them in groups of 5 minutes; S4. Classify aircraft using satellite-based augmentation system (SBAS) and standard position service (SPS) based on unsupervised classification method; S5. Based on the assumption that the precision factor HFOM follows the Rayleigh distribution, a cost function is constructed. For aircraft that appear in a single grid within 5 minutes, the number of ADS-B reports N with NACp = k when k = {8, 9, 10, 11} is counted. k , according to N k Estimate the short-term mean of HFOM S6. Based on the assumption that the precision factor VFOM follows a half-normal distribution, a cost function is constructed. For aircraft that appear in a single network within 5 minutes, the number of ADS-B reports M with GVA = k when k = {1, 2} is counted. k , according to M k Estimate the short-term average value of VFOM 2. The method for evaluating the impact of geomagnetic storms on aviation GNSS service accuracy based on ADS-B data according to claim 1, characterized in that: In S2, the interpolation method uses first-order linear interpolation or second-order interpolation for calculation.

3. The method for evaluating the impact of geomagnetic storms on aviation GNSS service accuracy based on ADS-B data according to claim 1, characterized in that: In S3, due to the non-uniform geographical distribution of ADS-B data, dynamic grid division should be adopted in the gridding process so that a single grid has sufficient data at every moment. The basic principle is to ensure that at a five-minute interval, the ADS-B data contained in each network in each time period is not less than 1,000 at peak times to ensure the accuracy of the parameter estimation process.

4. The method for evaluating the impact of geomagnetic storms on aviation GNSS service accuracy based on ADS-B data according to claim 3 is characterized in that: In S4, the unsupervised classification method uses the K-means clustering algorithm to classify the data. The input of the clustering process is the historical NACp mean of each aircraft. The average value is calculated daily for the massive data across days, and the data with larger means are selected to input the clustering process; aircraft with means less than 8 are pre-excluded. For areas with sparse ADS-B data or times when there are fewer flights in the night airspace, unified parameter estimation can be performed for multiple time periods with an interval of five minutes.

5. The method for evaluating the impact of geomagnetic storms on aviation GNSS service accuracy based on ADS-B data according to claim 4, characterized in that: The specific steps to classify an aircraft are as follows: S41, the clustering input dimension is expanded to multi-parameter joint features, including NACp, GVA, NIC, and SIL, and each parameter is normalized in advance to eliminate dimensional differences; S42, using Gaussian mixture model instead of K-means clustering method, setting the covariance matrix type to diagonal matrix, and constraining the number of clusters to 2 to distinguish SBAS / SPS service types; S43, weighting the multi-dimensional joint features, wherein the NACp parameter weight is 0.4-0.5, the GVA weight is 0.25-0.35, and the NIC and SIL weights are 0.1-0.15 respectively; and the sum of the weights is equal to 1; S44. Dynamically generate SBAS / SPS discrimination rules based on clustering results: when the cluster center satisfies NACp≥9.5, GVA≥1.8, SIL≥2.9 and NIC≥9, it is determined to be the SBAS service dominant cluster; the rest are determined to be the SPS service dominant cluster; S45. The boundary samples are filtered using a probability membership threshold, and only samples with a membership greater than 80% are retained to participate in the geomagnetic storm impact analysis.

6. The method for evaluating the impact of geomagnetic storms on aviation GNSS service accuracy based on ADS-B data according to claim 5, characterized in that: In S5, the cost function is established as the sum of four parts of the same form. When k = {1, 2, 3, 4}, the cost function l H (σ) is calculated as follows: Frequency N 12-k The probability density function of the theoretical Rayleigh distribution is arrive The log-likelihood function formed by the integral value on the interval is as follows: in, and These are two sets of constants, representing the upper and lower bounds when the ADS-B airborne equipment maps the continuous parameter HFOM output by the GNSS receiver to different values ​​of the discrete parameter NACp in the ADS-B MS report.

7. The method for evaluating the impact of geomagnetic storms on aviation GNSS service accuracy based on ADS-B data according to claim 6, characterized in that: In S6, the cost function is established as the sum of two parts of the same form, the cost function l V (σ) is calculated as follows: Frequency M k The theoretical half-normal distribution probability density function is arrive The log-likelihood function formed by the integral value on the interval is as follows: in, and They respectively represent the upper and lower bounds when the ADS-B airborne equipment maps the continuous parameter VFOM output by the GNSS receiver to different values ​​of the discrete parameter GVA in the ADS-B MS report.

Images