A GNSS interference source localization method based on ADS-B data

Abstract

This invention discloses a GNSS interference source localization method based on ADS-B data. The method involves acquiring ADS-B data, saving data points affected by GNSS interference, and, if no NIC is available, constructing ADS-B data with a NIC. The message data is then sorted by ICAO and timestamp. A diagonal matrix W is calculated based on the ICAO, timestamp, NIC, a table showing the correspondence between NIC and received interference power, and the hyperparameter value τ. A residual vector R is calculated based on the message data and the optimization vector, and the optimized vector is then output. This invention analyzes large-scale aviation ADS-B data, sorts and classifies the original message data to determine whether flights are subject to GNSS interference, and abstracts the interference source localization problem into a least-squares problem. An optimization function is constructed, and an iterative optimization algorithm is used to calculate the possible GNSS interference source locations and transmission power, thus more efficiently locating GNSS interference sources.

Description

Technical Field

[0001] The present invention relates to the technical field of radio monitoring. Specifically, it is a method for locating GNSS jamming source based on ADS-B data. Background Art

[0002] With the increasingly wide application fields of GNSS (Global Navigation Satellite System), the GNSS radio electromagnetic environment is becoming more and more complex. There are often interference situations for GNSS signals, which pose a very serious threat to the accuracy, availability, and integrity of satellite navigation. Therefore, a method is needed to efficiently and accurately locate the position of GNSS jamming sources. Since the longitude, latitude, altitude positions of flights, and parameters such as NIC (Navigation Integrity Category) that reflect the GNSS signal quality and GNSS positioning accuracy can be obtained in ADS-B (Automatic Dependent Surveillance-Broadcast) messages, and the ADS-B system has been gradually improved and widely applied to global flights. Therefore, analyzing big data of aviation ADS-B to solve the problem of GNSS jamming source location will be a future trend. Summary of the Invention

[0003] The purpose of the present invention is to provide a method for locating GNSS jamming source based on ADS-B data, which is used to solve the problem that there is no solution to the GNSS jamming source location by analyzing big data of aviation ADS-B in the prior art.

[0004] The present invention solves the above problems through the following technical solutions:

[0005] A method for locating GNSS jamming source based on ADS-B data includes:

[0006] Step S100: Obtain ADS-B data. If there is a navigation integrity category NIC in the ADS-B message data, filter and save the ADS-B message data with NIC < N. When NIC < N is determined that the GNSS is jammed at this position, where N ∈ [0, 11]; then proceed to the next step. Otherwise, first construct the ADS-B data with NIC, and then proceed to the next step;

[0007] Step S200: Sort the message data according to the airport identification code icao and timestamp, and calculate the diagonal matrix W according to the icao number, timestamp, NIC in the message, the correspondence table between NIC and received interference power, and the hyperparameter value τ;

[0008] Step S300: Initialize the optimization vector. The optimization vector includes the longitude, latitude, altitude, and transmit power value of the interference source. Calculate the residual vector R based on the message data and the initially defined optimization vector, using the scalar R... T WR serves as the optimization function. It uses an iterative algorithm to minimize the optimization function and update the optimization vector until the iterative algorithm converges, at which point the optimization vector is output.

[0009] Furthermore, methods for constructing ADS-B data with a NIC include:

[0010] Step A1: Sort the message data according to the airport identification code (icao) and timestamp;

[0011] Step A2: Based on the longitude, latitude, speed, heading, and time difference between the current data and the time in the next data point, calculate the theoretical longitude and latitude of the next data point.

[0012] Step A3: Calculate the distance difference Δd between the theoretical latitude and longitude position and the actual latitude and longitude position of the next point, and construct the NIC of the next data point according to the correspondence table between NIC and containment radius;

[0013] Step A4: Traverse the data and repeat steps SA2 to A3 to construct a NIC from the ADS-B data that does not have a NIC; proceed to step S200.

[0014] Further, step S200 specifically includes:

[0015] Step S210: Sort the message data according to the airport identification code (icao) and timestamp;

[0016] Step S220: Calculate the diagonal matrix W based on the ICA0 number, timestamp, NIC, the NIC-to-received interference power correspondence table, and the hyperparameter value τ in the message. Specifically, this includes:

[0017] Step S221: Calculate the variance of the i-th data of a certain ICA0 number in the message data. Where P rmax P rmin These are the upper and lower bounds of the received interference power corresponding to the NIC, respectively.

[0018] Step S222: Calculate the correlation coefficient of a certain ICA0 number in the message data. Where t i -t j τ is the difference between the timestamps of different message data i and j of the same ICA0 number, and τ is a hyperparameter reflecting the time correlation between different data points of the same ICA0.

[0019] Step S223, according to W ij =ρ ij σ i σ j Calculate the element values ​​of matrix W, where i and j are data with the same icao number;

[0020] Step S224: Traverse the sorted data and repeat steps S221 to S223 to calculate the diagonal matrix W.

[0021] Further, step S300 specifically includes:

[0022] Step S310: Based on the data and the initially defined optimization vector Calculate the residual vector R, where x, y, z are the longitude, latitude, and altitude positions in the ECEF coordinate system, and P t It refers to the transmission power of the transmitter, specifically including:

[0023] Step S311: Calculate the received interference power at the current flight position using the free space propagation model. Where x a y a z a It refers to the location of longitude, latitude, and altitude in the ADS-B message data within the ECEF coordinate system; G t G represents the gain of the transmitting antenna. r λ represents the gain of the receiving antenna; λ is the wavelength.

[0024] Step S312: Calculate the residual element of the i-th data in the message data. in It is the median value of the received interference power range corresponding to the NIC value in message data i;

[0025] Step S313: Traverse the ADS-B message data and repeat steps S141 to S142 to calculate all residual vector elements to form the residual vector R;

[0026] Step S320, using scalar R T WR as an optimization function To optimize the vector, an iterative optimization algorithm is used to minimize the optimization function and update the X vector. Then, the norm difference tol between the new X vector and the previous X vector is calculated.

[0027] Step S330: Repeat steps S311 to S320 until the norm difference tol converges to the threshold x. tol This yields the final optimized vector X, which in turn gives the location and transmission power of the GNSS interference source.

[0028] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0029] (1) By analyzing the aviation ADS-B big data, the present invention sorts and classifies the original message data, determines whether the flight is affected by GNSS interference, abstracts the problem of locating the interference source into a least-squares problem, constructs an optimization function, and uses an iterative optimization algorithm to calculate the possible location and transmission power of the GNSS interference source, more efficiently locating the location of the GNSS interference source, which is of great significance for GNSS interference detection and interference troubleshooting.

[0030] (2) Based on the aviation ADS-B big data, the present invention solves the problem of locating the GNSS interference source more comprehensively, efficiently, time-saving and labor-saving. At the same time, a method for constructing NIC for ADS-B data without NIC is provided. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] Figure 1 is a flow chart of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0032] The present invention will be further described in detail below in conjunction with embodiments, but the embodiments of the present invention are not limited thereto.

[0033] Embodiment:

[0034] Combined with the attached Figure 1 , a method for locating GNSS interference source based on ADS-B data includes:

[0035] Step S100, obtain ADS-B data. If there is a navigation integrity category NIC in the ADS-B message data, filter and save the ADS-B message data with navigation integrity category NIC < N. It is determined that GNSS interference occurs at this location when NIC < N, where N ∈ [0, 11]; in this embodiment, N is selected as 7; then go to the next step, otherwise first construct ADS-B data with NIC and then go to the next step; the method for constructing ADS-B data with NIC includes:

[0036] Step A1, sort the message data according to the airport identification code icao and the timestamp, where the timestamp is in Unix format;

[0037] Step A2, calculate the theoretical longitude and latitude positions of the next data point according to the longitude, latitude, speed, heading in the current message data and the time difference between the next data and the current message data time;

[0038] Step A3, calculate the distance difference Δd between the theoretical longitude and latitude positions and the actual longitude and latitude positions of the next point, and construct the NIC of the next data point according to the NIC and the containment radius correspondence table. The NIC and the containment radius correspondence table is shown in the following table.

[0039] Table 1 Correspondence between NIC and Containment Radius

[0040]

[0041] Step A4: Traverse the data and repeat steps SA2 to A3 to construct a NIC from the ADS-B data that does not have a NIC;

[0042] Step S200: Abstract the problem into a least squares problem and sort the message data according to the airport identification code (icao) and timestamp. Calculate the diagonal matrix W based on the ICAo number, timestamp, NIC, the NIC-to-received interference power correspondence table, and the hyperparameter value τ in the message; specifically including:

[0043] Step S210: Sort the message data according to the airport identification code (icao) and timestamp;

[0044] Step S220: Calculate the diagonal matrix W based on the ICA0 number, timestamp, NIC, the NIC-to-received interference power correspondence table, and the hyperparameter value τ in the message. Specifically, this includes:

[0045] Step S221: Calculate the variance of the i-th data of a certain ICA0 number in the message data. Where P rmax P rmin These are the upper and lower bounds of the received interference power corresponding to the NIC; the correspondence between the NIC and the received interference power is shown in Table 2.

[0046] Table 2. Correspondence between NIC and Received Interference Power

[0047] NIC <![CDATA[P r Range[dBW]]]> 0 [-100,0) 1 [-105,-100) 2 [-110,-105) 3 [-115,-110) 4 [-120,-115) 5 [-125,-120) 6 [-130,-125) 7 (-∞,-130)

[0048] Step S222: Calculate the correlation coefficient of a certain ICA0 number in the message data. Where t i -t j τ is the difference between the timestamps of different message data i and j of the same ICA0 number, and τ is a hyperparameter reflecting the time correlation between different data points of the same ICA0.

[0049] Step S223, according to W ij =ρ ij σ i σ j Calculate the element values ​​of matrix W, where i and j are data with the same icao number;

[0050] Step S224: Traverse the sorted data and repeat steps S221 to S223 to calculate the diagonal matrix W;

[0051] Step S300: Initialize the optimization vector. The optimization vector includes the longitude, latitude, altitude, and transmit power value of the interference source. Calculate the residual vector R based on the message data and the initially defined optimization vector, using the scalar R... T WR serves as the optimization function. An iterative algorithm is used to minimize the optimization function and update the optimization vector until the iterative algorithm converges. The output optimization vector includes:

[0052] Step S310: Based on the data and the initially defined optimization vector Calculate the residual vector R, where x, y, z are the longitude, latitude, and altitude positions (in meters) in the ECEF coordinate system, and P t This refers to the transmitter's power (unit: watts), specifically including:

[0053] Step S311: Calculate the received interference power at the current flight position using the free space propagation model. Where x, y, and z are the elements in the optimization vector X, which represents the location of the interference source. a y a z a It refers to the location of longitude, latitude, and altitude in the ADS-B message data within the ECEF coordinate system; G t G is 1; r =1; λ is 1575.42MHz;

[0054] Step S312: Calculate the residual element of the i-th data in the message data. in It is the median value of the received interference power range corresponding to the NIC value in message data i;

[0055] Step S313: Traverse the ADS-B message data and repeat steps S141 to S142 to calculate all residual vector elements to form the residual vector R;

[0056] Step S320, using scalar R T WR as an optimization function To optimize the vector, use an iterative optimization algorithm (such as the Levenberg-Marquardt algorithm) to minimize the optimization function and update the X vector, and then calculate the norm difference tol between the new X vector and the previous X vector;

[0057] Step S330: Repeat steps S311 to S320 until the norm difference tol converges to the threshold x. tol (e.g., 1×10) -10 This yields the final optimized vector X, which in turn gives the location and transmission power of the GNSS interference source.

[0058] This invention analyzes aviation ADS-B big data, sorts and classifies the raw message data, determines whether a flight is subject to GNSS interference, and abstracts the problem of locating the interference source into a least squares problem. It constructs an optimization function and uses an iterative optimization algorithm to calculate the possible location and transmission power of the GNSS interference source, thus locating the GNSS interference source more efficiently. This is of great significance for GNSS interference detection and troubleshooting.

[0059] Although the present invention has been described herein with reference to illustrative embodiments, the above embodiments are merely preferred embodiments of the present invention, and the implementation of the present invention is not limited to the above embodiments. It should be understood that those skilled in the art can devise many other modifications and implementations, which will fall within the scope and spirit of the principles disclosed in this application.

Claims

1. A GNSS interference source localization method based on ADS-B data, characterized in that, Including: Step S100: Obtain ADS-B data. If there is a navigation integrity category NIC in the ADS-B message data, filter and save the ADS-B message data with NIC < N. When NIC < N is determined that GNSS interference occurs at this location, where N ∈ [0, 11]; then proceed to the next step. Otherwise, first construct ADS-B data with NIC and then proceed to the next step; Step S200: Sort the message data according to the airport identification code icao and the timestamp, and calculate the diagonal matrix W according to the icao number, timestamp, NIC in the message, the NIC and received interference power correspondence table, and the hyperparameter value τ; Step S300: Initialize the optimization vector. The optimization vector includes the longitude, latitude, altitude, and transmit power value of the interference source. Calculate the residual vector based on the message data and the initially defined optimization vector. , in scalar As an optimization function, the optimization function is minimized using an iterative algorithm, and the optimization vector is updated until the iterative algorithm converges, at which point the optimization vector is output. The specific steps of step S200 include: Step S210: Sort the message data according to the airport identification code icao and the timestamp; Step S220: Calculate the diagonal matrix W according to the icao number, timestamp, NIC in the message, the NIC and received interference power correspondence table, and the hyperparameter value τ, specifically including: Step S221: Calculate the variance of the i-th data of a certain ICA0 number in the message data. ,in These are the upper and lower bounds of the received interference power corresponding to the NIC, respectively. Step S222: Calculate the correlation coefficient of a certain ICA0 number in the message data. ,in τ is the difference between the timestamps of different message data i and j of the same ICA0 number, and τ is a hyperparameter reflecting the time correlation between different data points of the same ICA0. Step S223, according to Calculate the element values ​​of matrix W, where i and j are data with the same icao number; Step S224: Traverse the sorted data and repeat steps S221 to S223 to calculate the diagonal matrix W; The specific steps of step S300 include: Step S310: Based on the data and the initially defined optimization vector Calculate the residual vector R, where x, y, and z are the longitude, latitude, and altitude positions in the ECEF coordinate system. It refers to the transmission power of the transmitter, specifically including: Step S311: Calculate the received interference power at the current flight position using the free space propagation model. ,in , , It refers to the position of longitude, latitude, and altitude in the ECEF coordinate system within the ADS-B message data; The gain of the transmitting antenna; The gain of the receiving antenna; Wavelength; Step S312: Calculate the residual element of the i-th data in the message data. ,in It is the median value of the received interference power range corresponding to the NIC value in message data i; Step S313: Traverse the ADS-B message data and repeat steps S141 to S142 to calculate all elements of the residual vector, and form the residual vector R; Step S320, using scalar As an optimization function To optimize the vector, an iterative optimization algorithm is used to minimize the optimization function and update the X vector. Then, the norm difference tol between the new X vector and the previous X vector is calculated. Step S330: Repeat steps S311 to S320 until the norm difference tol converges to the threshold. This yields the final optimized vector X, which in turn gives the location and transmission power of the GNSS interference source.

2. The GNSS interference source localization method based on ADS-B data according to claim 1, characterized in that, The method for constructing ADS-B data with NIC includes: Step A1: Sort the message data according to the airport identification code icao and the timestamp; Step A2: Calculate the theoretical longitude and latitude positions of the next data point according to the longitude, latitude, speed, heading in the current message data, and the time difference between the next data and the current message data time; Step A3: Calculate the distance difference Δd between the theoretical longitude and latitude positions and the actual longitude and latitude positions of the next point, and construct the NIC of the next data point according to the NIC and suppression radius correspondence table; Step A4: Traverse the data and repeat steps SA2 to A3 to construct NIC for the ADS-B data without NIC; proceed to step S200.

Images