Probability model influence assessment method of unmanned aerial vehicle countering device on flight GPS

Abstract

The invention relates to a probability model influence assessment method of unmanned aerial vehicle countering equipment on a flight GPS. The method comprises the following steps: determining an interference probability model expression of the unmanned aerial vehicle countering equipment on the flight GPS; acquiring flight broadcast type automatic related monitoring system data; calculating a probability function of the GPS interference probability changing along with the distance; solving constants in the interference probability model; probability function conversion; converting a probability distribution function through variable substitution; solving a probability density function corresponding to the changed probability distribution function; fitting the parameterized function to a probability density function; solving a function in the interference probability model and obtaining the interference probability model; and obtaining an interfered probability distribution diagram according to the obtained interference probability model. According to historical data of an unmanned aerial vehicle countering equipment peripheral flight broadcast type automatic correlation monitoring system, under the condition that the layout mode is similar to the surrounding environment, the interference combination probability of equipment of the same model to flight GPSs at different surrounding positions is quantitatively evaluated.

Description

Technical Field

[0001] This invention belongs to the field of interference source impact assessment, specifically relating to a probabilistic model impact assessment method for unmanned aerial vehicle (UAV) countermeasures equipment on flight GPS. Background Technology

[0002] To address the security risks posed by drones, drone countermeasures are being widely used. Most existing drone countermeasures rely on radio interference with the drone's GPS (Global Positioning System) or remote control signals. However, while these countermeasures work, they may also interfere with normally operating equipment.

[0003] GPS is also a crucial data source for positioning information on civil aviation flights. Radio interference signals emitted by drone countermeasures equipment can also interfere with the acquisition of GPS positioning information for civil aviation flights, causing errors or loss of GPS positioning data and affecting flight safety. Therefore, it is necessary to assess the potential interference of drone countermeasures equipment on the GPS of civil aviation flights.

[0004] However, existing assessments of the potential interference of drone countermeasures equipment with GPS on civil aviation flights are mostly qualitative. There is currently a lack of effective quantitative assessment methods that combine probability to evaluate the potential interference of drone countermeasures equipment with GPS on civil aviation flights. Summary of the Invention

[0005] The purpose of this invention is to overcome the shortcomings of the prior art and provide a quantitative and reliable method for assessing the probabilistic model impact of unmanned aerial vehicle (UAV) countermeasures equipment on flight GPS.

[0006] To solve the above problems, the technical solution of the present invention is as follows:

[0007] This invention provides a method for assessing the probabilistic impact of unmanned aerial vehicle (UAV) countermeasures devices on flight GPS, characterized by comprising:

[0008] Step 1: Establish a probability model of interference between drone countermeasures and flight GPS. The expression for the interference probability model is defined as follows:

[0009] P(d) = K·L(d);

[0010] In the formula, d is the distance between the location to be analyzed and the UAV countermeasure device, P(d) is the probability that the GPS of the flight at a distance d from the UAV countermeasure device is interfered with, L(d) is the probability that the typical received power of the flight at a distance d from the UAV countermeasure device exceeds the interference threshold power, and K is the probability that the GPS of the flight is interfered with under the condition that the typical received power of the flight receiving the UAV countermeasure device exceeds the interference threshold power.

[0011] Step 2: Obtain historical data from the Automatic Dependent Surveillance System (ADS) of flights during the period of disruption caused by drone countermeasures equipment in cases of drone disruption.

[0012] Step 3: Based on the historical data of the broadcast automatic dependent surveillance system, calculate the proportion of the number of flights whose GPS is interfered with in the total number of flights within different distance intervals from the UAV countermeasures device, and construct a probability function S(d) of the probability of flight GPS interference as a function of distance based on the proportion.

[0013] Step four: Obtain the maximum value of the probability function S(d), and use the maximum value as the estimate of K in the expression of the interference probability model.

[0014] Step 5: Convert the probability function S(d) into a first probability distribution function T(d). The first probability distribution function T(d) corresponds to the probability that the typical received power of the signal from the UAV countermeasure device at a distance d from the UAV countermeasure device is less than the interference threshold power. The formula used for converting to the first probability distribution function T(d) is as follows:

[0015]

[0016] In the formula, The expression for the interference probability model is the estimated result of K, where d is the distance between the location to be analyzed and the UAV countermeasure device, T(d) is the first probability distribution function, and S(d) is the probability function.

[0017] Step 6: Perform variable substitution using the formula x = v(d), where d is the distance between the location to be analyzed and the UAV countermeasure device, v(d) is the variable substitution function, and x is the independent variable reflecting the distance after the variable substitution. The first probability distribution function T(d) about the independent variable d is transformed into the second probability distribution function Q(x) about the independent variable x.

[0018] Step 7: Solve for the corresponding probability density function q(x) based on the second probability distribution function Q(x);

[0019] Step 8: Fit the probability density function q(x) using the parameterized function r(x) and estimate the parameters of the parameterized function r(x);

[0020] Step nine, using the parameterized function r(x), and employing the formula... The estimation result of L(d) in the expression of the interference probability model is obtained, which represents the probability that the typical received power of the UAV countermeasures device signal received by a flight at a distance d from the UAV countermeasures device exceeds the interference threshold power. In the formula, v(d) is the variable substitution formula used in step six by performing variable substitution through the formula x = v(d), combined with the estimation result of K in the expression of the interference probability model obtained in step four. Substituting the values ​​into the expression of the interference probability model, we obtain the estimation result of the interference probability model. for

[0021] Step 10: For the countermeasures device against the same type of UAV to be evaluated, using the device as the center, within a horizontal plane at the altitude to be analyzed, calculate the results based on the estimation results of the interference probability model. The probability distribution map of GPS interference affecting flights around the drone countermeasure equipment is obtained as the output result of the impact assessment.

[0022] In an optional implementation, in step six, variable substitution is performed using the formula x = v(d), where d is the distance from the location to be analyzed to the UAV countermeasure device, v(d) is the variable substitution function, and x is the independent variable reflecting the distance after the variable substitution. The first probability distribution function T(d) with respect to the independent variable d is transformed into a second probability distribution function Q(x) with respect to the independent variable x. Specifically, the variable substitution formula x = v(d) is x = 20log 10 d.

[0023] In one optional implementation, in step one, an interference probability model of the UAV countermeasure device on the flight GPS is established, and the expression of the interference probability model is set as: P(d)=K·L(d), specifically including: where d is the distance between the location to be analyzed and the UAV countermeasure device, and the distance is the straight-line distance.

[0024] In an optional implementation, in step three, based on historical data from the automatic dependent surveillance broadcast system, the proportion of flights whose GPS is interfered with within different distance intervals from the UAV countermeasures device is calculated out of the total number of flights. A probability function S(d) is then constructed from this proportion, which specifically includes:

[0025] The historical data of the broadcast automatic dependent surveillance system is preprocessed. The preprocessing includes dividing the track segments with normal time intervals between adjacent track points, removing track points that deviate significantly from the normal track position, and supplementing track points between track points with abnormal time intervals between adjacent track points to obtain the preprocessed broadcast automatic dependent surveillance system data.

[0026] The navigation quality index values ​​in the preprocessed broadcast automatic dependent surveillance system data are used to determine whether each waypoint is an interfered waypoint.

[0027] Within the range where the distance to the UAV countermeasure device does not exceed the maximum analysis distance, starting from the minimum analysis distance, the analysis distance from the minimum analysis distance to the maximum analysis distance is divided into different distance intervals;

[0028] For a distance interval (a, b) with a minimum distance of a and a maximum distance of b from the UAV countermeasures device, all trackpoints in the preprocessed Automatic Dependent Surveillance-Broadcast (ADBBS) data within the distance interval are analyzed. The total number of flights corresponding to all trackpoints is the total number of flights in the distance interval. If the number of interfered trackpoints of a flight within the distance interval is not less than a predetermined threshold, then the flight is a flight with GPS interference. Within the distance interval, the ratio of the number of flights with GPS interference to the total number of flights is the proportion of the number of flights with GPS interference in the distance interval to the total number of flights. The proportion of the number of flights with GPS interference in the distance interval to the total number of flights is calculated sequentially for distance intervals at different distances from the UAV countermeasures device.

[0029] The proportion of the number of flights with GPS interference within the distance interval (a,b] out of the total number of flights is taken as the function value of the probability function S(d) of the probability of GPS interference as a function of distance at the midpoint of the distance interval. By calculating the proportion of the number of flights with GPS interference within the distance interval at different distances out of the total number of flights, the function value of the probability function S(d) at different distances is obtained, thus obtaining the probability function S(d) of the probability of GPS interference as a function of distance.

[0030] In an optional implementation, in step seven, the probability density function q(x) is solved according to the second probability distribution function Q(x), specifically including: taking the derivative of the second probability distribution function Q(x) using the formula q(x)=Q'(x) to obtain the probability density function q(x).

[0031] In an optional implementation, in step eight, the probability density function q(x) is fitted using a parameterized function r(x), and the parameters of the parameterized function r(x) are estimated. Specifically, the parameterized function r(x) used is: In the formula, the parameterized function r(x) is the Gaussian probability density function, μ is the mean of the Gaussian probability density function, σ is the standard deviation of the Gaussian probability density function, x is the independent variable of the function reflecting distance after variable substitution, and the parameters of the estimated parameterized function r(x) are the mean μ and standard deviation σ of the estimated Gaussian probability density function.

[0032] In one optional implementation, in step ten, the countermeasures device for the same type of UAV to be evaluated is analyzed in a horizontal plane at an altitude of the altitude to be analyzed, centered on the device, based on the estimation results of the interference probability model. The probability distribution map of GPS interference affecting flights around the UAV countermeasure equipment is obtained as the output of the impact assessment. Specifically, this includes: calculating the distance d from different locations to the same type of UAV countermeasure equipment to be assessed within a horizontal plane at the altitude to be analyzed, and based on... Calculate the corresponding interference probability, draw a two-dimensional data distribution map with different colors to represent different interference probability value ranges, and overlay this two-dimensional data distribution map on the map of the corresponding area in a semi-transparent manner to obtain the GPS interference probability distribution map of flights around the UAV countermeasure equipment, and use the GPS interference probability distribution map of flights as the output result of the impact assessment.

[0033] In an optional implementation, in step three, within the range where the distance to the UAV countermeasure device does not exceed the maximum analysis distance, starting from the minimum analysis distance, the analysis distance from the minimum analysis distance to the maximum analysis distance is divided into different distance intervals, specifically including: dividing the analysis distance from the minimum analysis distance to the maximum analysis distance into different distance intervals with equal lengths by equal division.

[0034] In one optional implementation, in step three, determining whether each waypoint is an interfered waypoint based on the navigation quality index values ​​in the preprocessed Automatic Dependent Surveillance-Broadcast (ADS-B) data specifically includes: judging waypoints in tracks whose navigation integrity category index values ​​are not all invalid based on the navigation integrity category index values ​​in the ADS-B data; if the navigation integrity category index value of a waypoint is less than the normal threshold, then the waypoint is judged to be an interfered waypoint.

[0035] In one optional implementation, in step three, the division of the track segment with normal time intervals between adjacent track points specifically includes: for the acquired historical data of the Automatic Dependent Surveillance-Broadcast System, for each flight's track point data, based on whether the "Ground Station Received Position Time" data item contains valid data, dividing the track point data into track point data with timestamps and track point data without timestamps; for the track point data with timestamps of the same flight, sorting them in ascending order of time; for the sorted track point data, calculating the time interval between adjacent track point data; if the time interval between adjacent track points exceeds the normal time interval threshold, it is considered that the time interval between adjacent track points is abnormal; setting the two adjacent track points with abnormal time intervals as the start and end points of the track break, respectively; and dividing the track point data with timestamps of the same flight into multiple different track segments using the first and last track points of the same flight and the start and end points of the track break.

[0036] In an optional implementation, in step three, removing track points that significantly deviate from the normal track position specifically includes: after dividing the track segments with normal time intervals between adjacent track points, performing track filtering on the track point data within each track segment of each flight. The track filtering, based on the three-dimensional position vector filter output value and the three-dimensional velocity vector filter output value from the previous moment, and the current moment's observed three-dimensional position vector value, yields the current moment's predicted three-dimensional position vector value, the three-dimensional position vector filter output value, and the three-dimensional velocity vector filter output value; the deviation between the observed three-dimensional position vector value and the predicted three-dimensional position vector value is calculated: e i =z i -p i / i-1 In the formula, i is an integer not less than 2, and z i Let p be the three-dimensional position vector observation value at time i. i / i-1 Let e ​​be the predicted value of the three-dimensional position vector at time i. i Let be the deviation between the observed and predicted 3D position vector values ​​at time i. Calculate the standard deviation of the magnitude of the deviation for all times other than time i within the same track segment. If the magnitude of the deviation at any time is greater than three times the standard deviation, then the track point corresponding to that time is a track point that significantly deviates from the normal track position, and the track point corresponding to that time is removed.

[0037] In one optional implementation, in step three, the step of supplementing track points between track points with abnormal time intervals to obtain preprocessed broadcast automatic dependent surveillance system data specifically includes:

[0038] After dividing the track segments with normal time intervals between adjacent track points, each track segment contains timestamped track point data for the same flight. Between adjacent track segments of the same flight, the last track point of the previous track segment and the first track point of the next track segment are obtained. Starting from the last track point of the previous track segment, M track points are sequentially taken backward, where M is a non-negative integer. The average distance interval Δd between the M+1 adjacent track points is calculated. M and the average value of the time interval Δt M Starting from the first trackpoint of the next track segment, take N trackpoints sequentially, where N is a non-negative integer. Calculate the average distance Δd between the N+1 adjacent trackpoints. N and the average value of the time interval Δt N The average distance interval Δd between adjacent waypoints is obtained by weighted averaging:

[0039]

[0040] The average time interval Δt between adjacent waypoints is obtained by weighted averaging:

[0041]

[0042] Where M and N are both 10. If the number of track points before the last track point of the previous track segment is less than 10, then M is the total number of track points before the last track point of the previous track segment. If the number of track points after the first track point of the next track segment is less than 10, then N is the total number of track points after the first track point of the next track segment.

[0043] Between the last track point of the previous track segment and the first track point of the next track segment, supplementary track points are added. The number of supplementary track points is:

[0044]

[0045] In the formula, D is the distance between the last track point of the previous track segment and the first track point of the next track segment, Δd is the average distance interval between adjacent track points obtained by the weighted average, and W is the number of supplementary track points. Indicates rounding down;

[0046] The three-dimensional position vector of the w-th supplementary track point is:

[0047]

[0048] In the formula, z s z is the three-dimensional position vector of the last track point in the previous track segment. eLet z be the three-dimensional position vector of the first track point in the next track segment, w be the index of the supplementary track point, where 1≤w≤W, and W be the number of the supplementary track points. w The three-dimensional position vector of the w-th supplementary track point;

[0049] The time for the wth supplementary waypoint is:

[0050]

[0051] In the formula, t s t represents the time of the last waypoint of the previous waypoint segment. e The time of the first track point in the next track segment is given, w is the sequence number of the supplementary track point, where 1 ≤ w ≤ W, W is the number of the supplementary track points, and t is the time of the first track point in the next track segment. w The time for the wth additional waypoint;

[0052] All additional waypoints have their navigation integrity category index value set to -2 in the navigation quality metrics.

[0053] Add waypoints sequentially between all adjacent track segments of all flights;

[0054] After dividing the track segments with normal time intervals between adjacent track points and removing track points that deviate significantly from the normal track position, all track segment data for all flights are track point data with timestamps. This data, together with all supplementary track points corresponding to the flights, constitutes the preprocessed broadcast automatic dependent surveillance system data.

[0055] In one optional implementation, in step three, the track filtering specifically includes:

[0056] Track filtering using the α-β filtering method:

[0057]

[0058] In the formula, i is an integer not less than 2, j is an integer determined according to i, and p i-1 / i-1 Let v be the three-dimensional position vector filter output value at time i-1. i-1 / i-1 Let t be the three-dimensional velocity vector filter output value at time i-1. i / i-1 p represents the time interval from the (i-1)th time point to the ith time point. i / i-1 Let z be the predicted value of the three-dimensional position vector at time i. i Let p be the three-dimensional position vector observation value at time i. i / i Let v be the three-dimensional position vector filter output value at time i. i / i Let α be the three-dimensional velocity vector filter output value at time i, where α is the position filter parameter and β is the velocity filter parameter.

[0059] The three-dimensional position vector filter output value p at time i-1=1 is... 1 / 1 and the three-dimensional velocity vector filter output value v 1 / 1 The value can be:

[0060]

[0061] In the formula, z1 is the three-dimensional position vector observation value at the first time step, p 1 / 1 v is the three-dimensional position vector filter output value at time i-1=1. 1 / 1 This is the three-dimensional velocity vector filter output value at time i-1=1.

[0062] Compared with the prior art, the present invention has the following beneficial effects:

[0063] The present invention provides a method for assessing the impact of unmanned aerial vehicle (UAV) countermeasures equipment on flight interference based on the probability of GPS interference with flights. This method can quantitatively assess the interference of the same type of UAV countermeasures equipment on the GPS of flights at different locations in the vicinity of the surrounding area, based on historical data from the automatic dependent surveillance system broadcast to flights around the UAV countermeasures equipment, under conditions where the deployment method and the surrounding environment are similar. Attached Figure Description

[0064] Figure 1 This is a flowchart of the method for assessing the probabilistic model impact of the UAV countermeasure equipment on flight GPS according to the present invention;

[0065] Figure 2 This is a diagram illustrating the probability function S(d) of flight GPS interference probability as a function of distance, obtained from training data, according to an embodiment of the present invention.

[0066] Figure 3 This is a diagram illustrating a transformed first probability distribution function T(d) obtained from training data according to an embodiment of the present invention;

[0067] Figure 4 This is a diagram illustrating a second probability distribution function Q(x) of the independent variable x obtained from training data and transformed through variable substitution, according to an embodiment of the present invention.

[0068] Figure 5 This is a diagram illustrating the interpolated probability density function q(x) and the fitted parameterized function r(x) corresponding to the second probability distribution function Q(x) obtained from training data in an embodiment of the present invention.

[0069] Figure 6 This is an estimation result of an interference probability model obtained from training data according to an embodiment of the present invention. Comparison with the probability function S(d) of the flight GPS interference probability as a function of distance in the training data;

[0070] Figure 7 This is a probability distribution map of GPS interference affecting flights around a drone countermeasure device, obtained from training data, according to an embodiment of the invention.

[0071] Figure 8 This is an estimation result of an interference probability model obtained from training data according to an embodiment of the present invention. Comparison with the probability function S(d) of the GPS interference probability of flight 1 as a function of distance;

[0072] Figure 9 This is an estimation result of an interference probability model obtained from training data according to an embodiment of the present invention. Comparison with the probability function S(d) of the GPS interference probability of flight 2 as a function of distance;

[0073] Figure 10 This is an estimation result of an interference probability model obtained from training data according to an embodiment of the present invention. Comparison with the probability function S(d) of the GPS interference probability of flight 3 as a function of distance. Detailed Implementation

[0074] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0075] The following is combined Figures 1 to 10 The following describes embodiments of the present invention.

[0076] like Figure 1 As shown in the embodiment of the present invention, the method for assessing the probabilistic model impact of unmanned aerial vehicle (UAV) countermeasures equipment on flight GPS includes the following steps:

[0077] Step S101: Determine the probability model expression for the interference of the UAV countermeasures equipment on the flight's GPS;

[0078] Step S103: Obtain data from the Automatic Dependent Surveillance-Broadcast (ADS-B) system for flights;

[0079] Step S105: Calculate the probability function of GPS interference probability as a function of distance;

[0080] Step S107: Solve for the constant K in the disturbance probability model;

[0081] Step S109: Probability function transformation;

[0082] Step S111: Transform the probability distribution function through variable substitution;

[0083] Step S113: Solve for the probability density function corresponding to the changed probability distribution function;

[0084] Step S115: Parameterize the function to fit the probability density function;

[0085] Step S117: Solve for the function L(d) in the interference probability model and obtain the interference probability model;

[0086] Step S119: Obtain the interference probability distribution map based on the obtained interference probability model.

[0087] Step S101 involves determining the probability model expression for the interference of the UAV countermeasures device on the flight's GPS, specifically including:

[0088] An interference probability model of unmanned aerial vehicle (UAV) countermeasures equipment on flight GPS is established, and the expression of the interference probability model is defined as follows:

[0089] P(d) = K·L(d);

[0090] In the formula, d is the distance between the location to be analyzed and the UAV countermeasure device, P(d) is the probability that the GPS of the flight at a distance d from the UAV countermeasure device is interfered with, L(d) is the probability that the typical received power of the flight at a distance d from the UAV countermeasure device exceeds the interference threshold power, and K is a constant representing the probability that the GPS of the flight is interfered with under the condition that the typical received power of the flight receiving the signal from the UAV countermeasure device exceeds the interference threshold power. The flight GPS refers to the onboard GPS receiver of the flight.

[0091] In the above formula, d is the distance between the location to be analyzed and the UAV countermeasure device, and the distance is a straight-line distance.

[0092] Regarding step S103: This involves acquiring data from the Automatic Dependent Surveillance-Broadcast (ADS) system, specifically including:

[0093] In cases of unmanned aerial vehicle (UAV) countermeasures equipment disrupting air traffic, historical data from the Automatic Dependent Surveillance System (ADS) of flights during the period of disruption caused by the UAV countermeasures equipment were obtained.

[0094] Among them, the area surrounding the drone countermeasure equipment can be measured at a straight-line distance not exceeding the maximum analysis distance d. max The range; therefore, in this step, in cases of drone countermeasures disrupting flight, the straight-line distance to the drone countermeasures device does not exceed the maximum analysis distance d. maxHistorical data from the Automatic Dependent Surveillance-Broadcast (ADS) system during the period of disruption to air traffic within the specified range will be used. Subsequently, the percentage of GPS-interfered flights within the total number of flights will be calculated at different distances from the UAV countermeasures equipment, assuming the distance does not exceed the maximum analysis distance d. max The length of the distance interval from the drone countermeasure device within the specified range is Δl; since the interference range of the drone countermeasure device usually does not exceed 100km, for ease of calculation, 100km is set at no more than the maximum analysis distance d. max The midpoint of the range of distances from the drone countermeasure device within the specified range is the maximum analysis distance d. max You can take d max =100km+Δl / 2, where Δl can be taken as Δl = 5km, so d max Take d max =102.5km.

[0095] Step S105 involves calculating the probability function of GPS interference probability as a function of distance, specifically including:

[0096] Based on the historical data of the broadcast automatic dependent surveillance system, the proportion of the number of flights whose GPS is interfered with in the total number of flights is calculated within different distance intervals from the UAV countermeasures device. The probability function S(d) of the probability of flight GPS interference changing with distance is constructed from the proportion, where d is the distance of the location to be analyzed from the UAV countermeasures device.

[0097] The process in step S105 can be implemented through the following specific steps:

[0098] Step S301: Preprocess the historical data of the broadcast-type automatic dependent monitoring system.

[0099] Preprocessing includes dividing the track segments into those with normal time intervals between adjacent track points, removing track points that deviate significantly from the normal track positions, and supplementing track points between those with abnormal time intervals between adjacent track points to obtain preprocessed broadcast automatic dependent surveillance system data.

[0100] The method for dividing track segments with normal time intervals between adjacent waypoints is as follows:

[0101] For the acquired historical data from the Automatic Dependent Surveillance-Broadcast System, for each flight's track data, the track point data is divided into track point data with timestamps and track point data without timestamps based on whether the "Ground Station Received Position Time" data item contains valid data.

[0102] For the time-stamped track point data of the same flight, sort them in ascending order of time. For the sorted track point data, calculate the time interval between adjacent track point data. If the time interval between adjacent track points exceeds the normal time interval threshold, it is considered that the time interval between adjacent track points is abnormal. Set the two adjacent track points with abnormal time intervals as the start and end points of the track break. Through the first and last track points of the same flight and the start and end points of the track break, the time-stamped track point data of the same flight is divided into multiple different track segments.

[0103] The normal time interval threshold used is 5 seconds;

[0104] The method for removing waypoints that significantly deviate from the normal track position is as follows:

[0105] After dividing the track segments with normal time intervals between adjacent track points, track point data within each track segment of each flight is subjected to track filtering. The track filtering is based on the three-dimensional position vector filter output value, the three-dimensional velocity vector filter output value, and the three-dimensional position vector observation value at the current moment to obtain the three-dimensional position vector prediction value, the three-dimensional position vector filter output value, and the three-dimensional velocity vector filter output value at the current moment.

[0106] The deviation between the observed and predicted 3D position vector values ​​is calculated using the following formula:

[0107] e i =z i -p i / i-1 ,

[0108] In the formula, i is an integer not less than 2, z i Let p be the three-dimensional position vector observation value at time i. i / i-1 Let e ​​be the predicted value of the three-dimensional position vector at time i. i The deviation between the observed and predicted three-dimensional position vector values ​​at time i is given.

[0109] The deviation e is calculated for all times other than the first time point within the same track segment. i The modulus ||e i The standard deviation σ of || e If the modulus value at one of the time intervals is ||e i || Greater than three times the stated standard deviation 3σ e If the track point corresponding to the stated time is significantly deviated from the normal track position, then the track point corresponding to the stated time will be removed.

[0110] The method for track filtering is as follows:

[0111] Track filtering is performed using the α-β filtering method with the following formula:

[0112]

[0113] in

[0114]

[0115] In the formula, i is an integer not less than 2, j is an integer determined according to i, and p i-1 / i-1 Let v be the three-dimensional position vector filter output value at time i-1. i-1 / i-1 Let t be the three-dimensional velocity vector filter output value at time i-1. i / i-1 p represents the time interval from the (i-1)th time point to the ith time point. i / i-1 Let z be the predicted value of the three-dimensional position vector at time i. i Let p be the three-dimensional position vector observation value at time i. i / i Let v be the three-dimensional position vector filter output value at time i. i / i Let α be the three-dimensional velocity vector filter output value at time i, where α is the position filter parameter and β is the velocity filter parameter.

[0116] The three-dimensional position vector filter output value p at time i-1=1 is... 1 / 1 and the three-dimensional velocity vector filter output value v 1 / 1 Values

[0117]

[0118] In the formula, z1 is the three-dimensional position vector observation value at the first time step, p 1 / 1 v is the three-dimensional position vector filter output value at time i-1=1. 1 / 1 The output value of the three-dimensional velocity vector filter at time i-1 = 1;

[0119] The method for supplementing track points between track points with abnormal time intervals between adjacent track points to obtain preprocessed automatic dependent surveillance system (ADS) data is as follows:

[0120] After dividing the track segments with normal time intervals between adjacent track points, each track segment contains timestamped track point data for the same flight. Between adjacent track segments of the same flight, the last track point of the previous track segment and the first track point of the next track segment are obtained. Starting from the last track point of the previous track segment, M track points are sequentially taken backward, where M is a non-negative integer. The average distance interval Δd between the M+1 adjacent track points is calculated. M and the average value of the time interval Δt MSimilarly, starting from the first waypoint of the next waypoint segment, take N waypoints sequentially forward, where N is a non-negative integer, and calculate the average distance Δd between the N adjacent waypoints among these N+1 waypoints. N and the average value of the time interval Δt N Through the formula

[0121]

[0122] The weighted average is used to obtain the average distance interval Δd between adjacent waypoints, which is obtained through the formula.

[0123]

[0124] The average time interval Δt between adjacent waypoints is obtained by weighted averaging; where M and N are both 10. If the number of waypoints before the last waypoint of the previous waypoint is less than 10, then M is the total number of waypoints before the last waypoint of the previous waypoint. If the number of waypoints after the first waypoint of the next waypoint is less than 10, then N is the total number of waypoints after the first waypoint of the next waypoint.

[0125] Between the last track point of the previous track segment and the first track point of the next track segment, supplementary track points are added. The number of supplementary track points is:

[0126]

[0127] In the formula, D is the distance between the last track point of the previous track segment and the first track point of the next track segment, Δd is the average distance interval between adjacent track points obtained by the weighted average, and W is the number of supplementary track points. Indicates rounding down;

[0128] The three-dimensional position vector of the w-th supplementary track point is:

[0129]

[0130] In the formula, z s z is the three-dimensional position vector of the last track point in the previous track segment. e Let z be the three-dimensional position vector of the first track point in the next track segment, w be the index of the supplementary track point, where 1≤w≤W, and W be the number of the supplementary track points. w The three-dimensional position vector of the w-th supplementary track point;

[0131] The time for the wth supplementary waypoint is:

[0132]

[0133] In the formula, t s t represents the time of the last waypoint of the previous waypoint segment. e The time of the first track point in the next track segment is given, w is the sequence number of the supplementary track point, where 1 ≤ w ≤ W, W is the number of the supplementary track points, and t is the time of the first track point in the next track segment. w The time for the wth additional waypoint;

[0134] All additional waypoints have their navigation integrity category index value set to -2 in the navigation quality metrics.

[0135] Add waypoints sequentially between all adjacent track segments of all flights;

[0136] After dividing the track segments with normal time intervals between adjacent track points and removing track points that deviate significantly from the normal track position, all track segment data for all flights are track point data with timestamps. This data, together with all supplementary track points corresponding to the flights, constitutes the preprocessed broadcast automatic dependent surveillance system data.

[0137] Step S303: Determine whether each waypoint is an interfered waypoint by using the navigation quality index values ​​in the preprocessed automatic dependent surveillance system data.

[0138] Based on the navigation integrity category index value in the navigation quality index of the Automatic Dependent Surveillance-Broadcast system data, the track points in the track whose navigation integrity category index values ​​are not all invalid are judged. If the navigation integrity category index value of the track point is less than the normal threshold, the track point is judged to be an interfered track point.

[0139] The normal threshold for the navigation integrity category can be 7. If the navigation integrity category index value of a waypoint is less than 7, then the waypoint is judged to be an interfered waypoint.

[0140] Step S305: Within the distance of the UAV countermeasures device, not exceeding the maximum analysis distance d max Within the range, from the minimum analytical distance d min To begin, we will start with the minimum analysis distance d. min to the maximum analysis distance d max The analysis distance is divided into different distance intervals.

[0141] For ease of processing, the minimum analysis distance d will be used. min to the maximum analysis distance d max The analysis distance is obtained by dividing the distance interval into different distance intervals of equal length;

[0142] From the minimum analysis distance d min to the maximum analysis distance d maxWhen the analysis distance is divided into different distance intervals, the midpoint of each distance interval can be taken as an integer distance for ease of calculation.

[0143] In step S103, the maximum analysis distance d has been taken. max The length of the distance interval from the UAV countermeasures device within the range is Δl = 5km, and the maximum analysis distance is d. max =102.5km; therefore, the length of each of the different equally divided distance intervals can be Δl = 5km, and the minimum analytical distance d can be taken. min = 2.5km, so the distance intervals for different distances are (2.5km, 7.5km], (7.5km, 12.5km], ..., (97.5km, 102.5km], and the midpoints of the distance intervals for different distances are 5km, 10km, ..., 100km.

[0144] Step S307: Calculate the percentage of flights whose GPS is interfered with in the total number of flights within different distance ranges from the drone countermeasure device.

[0145] For a distance interval (a, b) between a minimum distance and a maximum distance of b from the UAV countermeasures device, all track points in the preprocessed automatic dependent surveillance system data within the distance interval are analyzed. The total number of flights corresponding to all track points is the total number of flights in the distance interval. If the number of interfered track points of one flight within the distance interval is not less than a predetermined threshold, then the flight is a flight whose GPS is interfered with.

[0146] The predetermined threshold used is 3;

[0147] Within the specified distance range, the ratio of the number of flights whose GPS is interfered with to the total number of flights is the percentage of the number of flights whose GPS is interfered with within the specified distance range to the total number of flights. Using this method, the percentage of the number of flights whose GPS is interfered with within the total number of flights is calculated sequentially within different distance ranges from the UAV countermeasure device.

[0148] An approximate calculation method for the number of flights affected by GPS interference and the total number of flights within the analyzed distance range is as follows:

[0149] Since the flight altitude is relatively small compared to the horizontal distance between the flight and the drone countermeasures equipment, the horizontal distance between the flight and the drone countermeasures equipment is used directly to replace the straight-line distance between the flight and the drone countermeasures equipment.

[0150] Within a coordinate system using latitude and longitude, the area surrounding the interference source is divided into grids based on latitude and longitude, with the interference source as the center. A uniform grid division can be used. The size of each grid can be calculated as follows: with the interference source as the center point, calculate the longitude deviation |Δl0| between the location d0 east of the interference source and the location d0 north of the interference source and the location d0 latitude. Then, |Δl0| and |Δb0| can be used as the magnitudes of each grid in the longitude and latitude directions to uniformly divide the area surrounding the interference source; where d0 can be taken as 5km.

[0151] Calculate the horizontal distance d from the center of each uniformly divided latitude and longitude grid to the interference source. c If d c Within the distance interval (a, b), all track points within the latitude and longitude grid are analyzed as track points within the distance interval (a, b) to calculate the proportion of the number of flights with GPS interference to the total number of flights. In the previous steps, the distance intervals for different distances were set as (2.5km, 7.5km], (7.5km, 12.5km], ..., (97.5km, 102.5km]. For example, if the horizontal distance from the center of a certain latitude and longitude grid to the interference source is 11.5km, then all track points within that latitude and longitude grid are analyzed as track points within the distance interval (7.5km, 12.5km], and the total number of flights and the number of flights with GPS interference are calculated according to the aforementioned method, and then the proportion of the number of flights with GPS interference to the total number of flights is calculated.

[0152] Step S309: Construct a probability function S(d) for the probability of GPS interference of a flight changing with distance, based on the proportion of the number of flights whose GPS is interfered with in the total number of flights, where d is the distance between the location to be analyzed and the UAV countermeasure device.

[0153] The proportion of flights with GPS interference within the distance interval (a, b) out of the total number of flights is taken as the probability function S(d) representing the change in the probability of GPS interference with distance, at the midpoint d of the distance interval (a, b). m The function value at a given distance is obtained by calculating the proportion of the number of flights whose GPS is interfered with in the total number of flights within different distance intervals. The function value of the probability function S(d) at different distances is thus obtained, and the probability function S(d) of the probability of GPS interference with flight changes with distance is obtained.

[0154] In one example, the probability function S(d) of the flight GPS interference probability varying with distance is shown in the training data. Figure 2 As shown in the figure, the dots represent the function values ​​of the probability function S(d) for the midpoint position in different distance intervals.

[0155] Step S107 involves solving for the constant K in the interference probability model, specifically including:

[0156] Obtain the maximum value of the probability function S(d) representing the probability of GPS interference affecting the flight as a function of distance, and use this maximum value as an estimate of the constant K in the expression P(d)=K·L(d) of the interference probability model.

[0157] Step S109 involves probability function transformation, specifically including:

[0158] The probability function S(d) representing the probability of GPS interference affecting a flight as a function of distance is converted into a first probability distribution function T(d). This first probability distribution function T(d) corresponds to the probability that the typical received power of the signal from the UAV countermeasures device at a distance d from the UAV countermeasures device is less than the interference threshold power. The formula used for this conversion is as follows:

[0159]

[0160] In the formula, The expression P(d) = K·L(d) in the interference probability model is the estimated result of K, where d is the distance between the location to be analyzed and the UAV countermeasure device, T(d) is the first probability distribution function, and S(d) is the probability function.

[0161] The transformed first probability distribution function T(d) obtained from one training data example is as follows: Figure 3 As shown in the figure, the dots represent the function values ​​of the midpoint positions of different distance intervals.

[0162] For step S111, the probability distribution function is transformed through variable substitution, specifically including:

[0163] By substituting variables using the formula x = v(d), where d is the distance between the location to be analyzed and the UAV countermeasure device, v(d) is the variable substitution function, and x is the independent variable reflecting the distance after the variable substitution, the first probability distribution function T(d) about the independent variable d is transformed into the second probability distribution function Q(x) about the independent variable x.

[0164] The variable substitution formula x = v(d) is specifically x = 20log 10 d.

[0165] An example of a second probability distribution function Q(x) obtained from training data is as follows: Figure 4 As shown in the figure, the horizontal axis represents the distance value corresponding to the variable x.

[0166] Step S113 involves solving for the probability density function corresponding to the changed probability distribution function, specifically including:

[0167] Based on the second probability distribution function Q(x), solve for the corresponding probability density function q(x);

[0168] The specific method is as follows: differentiate the second probability distribution function Q(x) using the formula q(x)=Q'(x);

[0169] The probability density function q(x) is obtained by interpolating the differentiated function at equal intervals according to a certain interval, where cubic spline interpolation is used.

[0170] For step S115, which involves fitting the probability density function to the parameterized function, specifically includes:

[0171] The probability density function q(x) is fitted using the parameterized function r(x), and the parameters of the parameterized function r(x) are estimated.

[0172] The parameterized function r(x) used is

[0173]

[0174] The parameterized function r(x) is a Gaussian probability density function, where μ is the mean of the Gaussian probability density function, σ is the standard deviation of the Gaussian probability density function, and x is the independent variable of the distance function after variable substitution. The parameters of the estimated parameterized function r(x) are the mean μ and standard deviation σ of the estimated Gaussian probability density function.

[0175] The fitting method is as follows:

[0176] Fitting is performed according to the criterion of minimizing mean square error, and the result is selected that minimizes the mean square error.

[0177]

[0178] The parameterized function is constructed using the minimum mean parameter μ and the standard deviation parameter σ. In the above formula for calculating the mean square error, x i Let I be the i-th sampling point of the function's independent variable, and let I be the total number of sampling points.

[0179] The optimal mean parameter μ and standard deviation parameter σ can be found using an exhaustive search method;

[0180] For the range of values ​​for the mean parameter μ, the maximum analytical distance d maxWhen the distance is 102.5km, the previous steps have set the distance intervals for different distances as (2.5km, 7.5km], (7.5km, 12.5km], ..., (97.5km, 102.5km], and the maximum value of the midpoint of the distance intervals for different distances is d'. max =100km, we can take x max =20log10(d' max In a distance interval where the number of flights is not zero, the distance with the smallest distance to the midpoint of the interval is denoted as d'. min , can take x min =20log10(d' min ), the value of the mean parameter μ during exhaustive search can be found in [x min ,x max Select at equal intervals within the range, with an interval of Δμ;

[0181] Regarding the range of values ​​for the standard deviation parameter σ, it is determined that when the mean μ is x max When, its corresponding μ-3σ is located at x min At that point, select the maximum value of σ. max =(x max -x min The minimum value of σ can be taken as σ / 3. min =Δμ / 3, the value of σ during exhaustive search can be found in [σ min ,σ max Select at equal intervals within the range, with an interval of Δσ;

[0182] When calculating the mean squared error, for a specific μ and σ, the mean squared error of the data in the range [μ-3σ, μ+3σ] should be calculated at least. During the calculation, the negative portions of the probability density function q(x) and the portions exceeding x should be excluded. max The part is assigned a value of 0 for calculation.

[0183] The interpolated probability density function q(x) and the fitted parameterized function r(x) corresponding to the second probability distribution function Q(x) obtained from the training data in the example are as follows: Figure 5 As shown in the figure, the dashed line represents the interpolated probability density function q(x) corresponding to the second probability distribution function Q(x), and the solid line represents the fitted parameterized function r(x).

[0184] Step S117 involves solving for the function L(d) in the interference probability model and obtaining the interference probability model, specifically including:

[0185] Using the parameterized function r(x), and employing the formula The expression P(d) = K·L(d) of the interference probability model is obtained, which estimates the probability L(d) that the typical received power of a flight at a distance d from the UAV countermeasures device exceeds the interference threshold power when receiving the signal from the UAV countermeasures device. In the formula, v(d) is the variable substitution formula used in step S111 to perform variable substitution through the formula x = v(d), combined with the estimation result of the constant K in the expression P(d) = K·L(d) of the interference probability model obtained in step four. Substituting the values ​​into the expression of the interference probability model, we obtain the estimation result of the interference probability model. for

[0186] When the variable substitution formula x = v(d) in step S111 is specifically x = 20log 10 When d, the formula The specific form is

[0187] An estimation result of the interference probability model obtained from training data in one embodiment. A comparison with the probability function S(d) of the flight GPS interference probability as a function of distance in the training data is shown below. Figure 6 As shown. The estimation results of the interference probability model obtained from the training data. like Figure 6 As shown by the solid line, for ease of comparison, Figure 6 The diagram also uses dashed lines to show the probability function S(d) of the flight GPS interference probability as a function of distance, obtained from the training data. Figure 6 The dotted line in the middle is Figure 2 The curve in the figure, where the dots are still the midpoints of different distance intervals, represents the function value of the probability function S(d).

[0188] For step S119, the interference probability distribution map is obtained based on the obtained interference probability model, specifically including:

[0189] For the countermeasures against the same type of UAV to be evaluated, with the device as the center, within a horizontal plane at the altitude to be analyzed, the estimation results of the interference probability model are used. The probability distribution map of GPS interference affecting flights around the drone countermeasure equipment is obtained as the output result of the impact assessment.

[0190] Within a horizontal plane at the height to be analyzed, calculate the distance d from different locations to the countermeasures device of the same type of UAV to be evaluated, and based on... Calculate the corresponding interference probability, draw a two-dimensional data distribution map with different colors to represent different interference probability value ranges, and overlay this two-dimensional data distribution map on the map of the corresponding area in a semi-transparent manner to obtain the GPS interference probability distribution map of flights around the UAV countermeasure equipment, and use the GPS interference probability distribution map of flights as the output result of the impact assessment.

[0191] An estimation result of the interference probability model obtained from training data in one embodiment. The generated probability distribution map of GPS interference affecting flights around the drone countermeasure equipment is as follows: Figure 7 As shown in the figure. Only the maximum analysis distance d is given in the figure. max =The results within a range of 102.5km. Areas not overlaid with semi-transparent methods in the figure are areas beyond the maximum analysis distance range, as well as areas within the maximum analysis distance range that have no evaluation results, such as areas very close to drone countermeasures devices. Because of the lack of supporting data, these areas have no evaluation results, thus avoiding giving incorrect results due to insufficient data support.

[0192] As described in step S101, an interference probability model of the UAV countermeasures device on the flight's GPS is established, and the expression of the interference probability model is set as follows:

[0193] P(d)=K·L(d),

[0194] In the formula, d is the distance between the location to be analyzed and the UAV countermeasures device, P(d) is the probability that the flight's GPS at a distance d from the UAV countermeasures device is interfered with, L(d) is the probability that the typical received power of the flight at a distance d from the UAV countermeasures device exceeds the interference threshold power, and K is a constant representing the probability that the flight's GPS is interfered with under the condition that the typical received power of the flight's received signal exceeds the interference threshold power. The flight's GPS refers to the onboard GPS receiver of the flight. The estimation result of the constant K is then estimated in subsequent steps. And using the parameterized function r(x), using the formula The estimation result of L(d) in the expression P(d)=K·L(d) of the interference probability model is obtained. The estimation results of the interference probability model are then obtained. The principle is:

[0195] Whether a flight's GPS is interfered with is closely related to the power of the received signal from the UAV countermeasures device. However, when the typical received power of the UAV countermeasures signal exceeds the interference threshold, interference does not necessarily occur due to factors such as the flight attitude and the receiver's anti-interference capabilities. Therefore, let parameter K be the probability of flight GPS interference under the condition that the typical received power of the UAV countermeasures signal exceeds the interference threshold. Parameter K should generally not fluctuate significantly, so it is simplified to a constant in the model. Thus, the probability of flight GPS interference at a distance d from the UAV countermeasures device is obtained as follows:

[0196] P(d)=K·L(d)

[0197] In the formula, L(d) is the probability that the typical received power of the signal from the UAV countermeasure device at a distance d from the UAV countermeasure device exceeds the interference threshold power.

[0198] Since the probability represented by L(d) should be a value between 0 and 1, the observed value of the interference probability P(d) and the maximum value of the probability function S(d) representing the change in the flight GPS interference probability with distance are used as the estimate of K in the model P(d) = K·L(d). then The value is between 0 and 1.

[0199] The main problem below is the estimation of the function L(d).

[0200] Typical received power P of a flight's GPS receiving signal from a drone countermeasure device at a distance d from the device. r can be derived from formula

[0201]

[0202] We obtain P in the formula t For the signal transmission power of drone countermeasure equipment, G t G represents the gain of the transmitting antenna of the UAV countermeasures device in the flight direction, λ represents the wavelength of the transmitted signal of the UAV countermeasures device, and G represents the gain of the transmitting antenna of the UAV countermeasures device in the flight direction. r0 Here, l represents the typical receiving antenna gain for a flight GPS antenna to receive signals from a drone countermeasure device, and l represents the signal loss during transmission.

[0203] Due to the gain G of the transmitting antenna of the drone countermeasure equipment in the flight direction t It will be affected by the pitch angle of the flight relative to the drone countermeasure equipment and the surrounding environment, thus making

[0204] G t =G t0 k t

[0205] In the formula G t0 k represents a typical value for the transmitting antenna gain of a drone countermeasure device. t To describe G t With G t0 The coefficient of difference.

[0206] G t =G t0 k t Substitution

[0207]

[0208] And by organizing, we can obtain

[0209]

[0210] Taking the logarithm of both sides of the above equation, we get

[0211]

[0212] In the above formula, the left side of the equal sign represents the received power value in dB at a distance d from the UAV countermeasure device, denoted as P. r(dB) (d); The first term on the right side of the equation is a constant value, which can be denoted as C; The second term on the right side of the equation is related to the distance d between the flight and the drone countermeasure equipment; The third term on the right side of the equation is because k t Both l and are affected by random factors such as the environment, and are therefore random variables. This random variable is denoted as R = R0 + ΔR, where R0 is the mean of the random variable R, which is a constant, and ΔR is a random variable with zero mean. Thus, the above equation simplifies to:

[0213] P r(dB) (d)=C+R0-20log 10 d+Δ R

[0214] Let P be the interference threshold power expressed in dB. th(dB) Meanwhile, let the random variable ΔR take the value of 0 and the received power P r(dB) (d)=P th(dB) At that time, the distance between the corresponding flight and the drone countermeasure device is d0, therefore we have

[0215] P th(dB) =C + R0 - 20log 10 d0

[0216] Let Y = ΔR + 20log 10 If d0, then Y is a vector with a mean of 20log. 10Let d0 be a random variable with the same variance as ΔR. Let the probability density function of the random variable Y be f(y) and its probability distribution function be F(y). Then, at a distance d from the UAV countermeasures device, the probability, expressed in dB, that the typical received power (in dB) of the signal received by the flight's GPS from the UAV countermeasures device exceeds the interference threshold power is:

[0217]

[0218] Let x = 20log 10 d, the above equation becomes

[0219]

[0220] And thus obtain

[0221]

[0222] Differentiating both sides with respect to x, we can obtain

[0223]

[0224] because

[0225]

[0226] And P(P) r(dB) (d)>p th(dB) This is also equal to the probability L(d) that the typical received power of the signal from the UAV countermeasure device received by a flight at a distance d from the UAV countermeasure device exceeds the interference threshold power. Therefore, the probability density function f(x) can be represented by a parameterized function r(x). Determining the parameters of this parameterized function determines r(x), and r(x) can be used as the estimate of the probability density function f(x).

[0227] Then utilize

[0228]

[0229] P(P) can then be obtained r(dB) (d)>p th(dB) The estimated value of ) That is, the estimated value of L(d)

[0230]

[0231] In specific implementation, the probability function S(d) representing the change in the probability of GPS interference with distance is the observed value of the interference probability P(d). The observed value corresponding to the probability that the typical received power of a flight at a distance d from the UAV countermeasures device exceeds the interference threshold power when receiving the signal from the UAV countermeasures device is also equal to P(Pr(dB) (d)>p th(dB) The value of ).

[0232] Using Relations

[0233]

[0234] Knowable formula

[0235] Let f(x) be the observed value of the probability density function f(x). Therefore, in step S109, the equation is used... To obtain T(d), in step S111, by variable substitution x = 20log 10 d yields Q(x), and then in step S113, the probability density function q(x) is obtained by differentiation. According to the following formula, the calculation shows that...

[0236]

[0237] q(x) is the observed value of the probability density function f(x).

[0238] In step S115, the probability density function q(x) is fitted using the parameterized function r(x), and the fitted result r(x) is used as an estimate of f(x) according to the following formula.

[0239]

[0240] Thus, the estimated result of the required function L(d) is obtained. for

[0241]

[0242] The probability density function f(x) corresponds to the random variable Y = ΔR + 20log 10 Since the random variable Y is affected by a large number of factors, according to the central limit theorem, it can be assumed that it follows a Gaussian probability distribution, and the Gaussian probability density function is used as the parameterization function to model f(x). Thus, in step S115, the probability density function q(x) is fitted with the Gaussian probability density function r(x).

[0243] The effectiveness of this application can be further illustrated by the following experiments.

[0244] Experiments were conducted using historical data from the Automatic Dependent Surveillance System (ADS-S) during the period of disruption caused by drone countermeasures equipment in real-world flight disruption cases.

[0245] Because existing measured data makes it difficult to find historical ADS-B data corresponding to the interference of flight GPS by the same type of UAV countermeasures equipment in different regions, and because different devices of the same model have similar impacts on surrounding flights under similar deployment methods and surrounding environments, we use the interference probability of the same device on the GPS of flights at different locations in the surrounding area at different time periods to replace the interference probability of the same type of device at different locations on the GPS of flights at different locations in the surrounding area. This is to verify the applicability of the interference probability obtained by the method of this invention to the same type of device.

[0246] In a case study of unmanned aerial vehicle (UAV) countermeasures equipment causing air traffic disruption, historical data from the Automatic Dependent Surveillance-Broadcast (ADS-B) system during one time period of disruption was used as training data to calculate the interference probability model for flight GPS. Historical ADS-B data from other time periods of the same disruption were then used as test data to verify the performance of the proposed interference probability model. Three hours of training data were used, and three sets of three different test data were used, resulting in three sets of different test data.

[0247] According to the method provided in this invention, the probability function S(d) obtained from the training data, which represents the change in the probability of GPS interference affecting a flight as a function of distance, is as follows: Figure 2 As shown in the figure, the dots represent the function values ​​of the probability function S(d) for the midpoint position in different distance intervals; this training data is used to obtain... Figure 2 The first probability distribution function T(d) after the transformation of the probability function S(d) is as follows: Figure 3 As shown in the figure, the dots represent the function values ​​at the midpoints of different distance intervals; the second probability distribution function Q(x) with respect to the independent variable x after variable substitution is as follows: Figure 4 As shown in the figure, the horizontal axis represents the distance value corresponding to the variable x; Figure 5 The dashed line represents the interpolated probability density function q(x) corresponding to the second probability distribution function Q(x), and the solid line represents the fitted parameterized function r(x). Figure 6 The solid line represents the estimation result of the disturbance probability model obtained from the training data. For ease of comparison, Figure 6 The diagram also uses dashed lines to show the probability function S(d) of the flight GPS interference probability as a function of distance, obtained from the training data. Figure 6 The dotted line in the middle is Figure 2 The curve in the figure shows that the dots represent the function values ​​of the probability function S(d) for the midpoint position of different distance intervals. From Figure 6 The results of the estimation of the perturbation probability model obtained from the training data can be seen in the image. It has good consistency with the probability function S(d).

[0248] Using a method similar to that used for the training data, the probability function of flight GPS interference probability as a function of distance was calculated for different test data. The estimation results of the interference probability model obtained from the training data were then used. Compared with the probability functions of flight GPS interference probability as a function of distance corresponding to different test data, the comparison results of the three different test data are as follows: Figure 8 , Figure 9 , Figure 10 As shown in the figure, the solid line represents the estimation result of the interference probability model obtained from the training data. The dashed line represents the probability function of flight GPS interference probability as a function of distance, obtained from the test data. As can be seen from the figure, regardless of the test data, the estimation results of the interference probability model obtained from the training data... The model shows good agreement with the probability function of the probability of interference in the test data as a function of distance, indicating that the estimation results of the interference probability model obtained through training are satisfactory. It can effectively describe the probability of interference in the test data. This also demonstrates that, under conditions of similar deployment and surrounding environment, the estimation results of the interference probability model obtained from the training data are satisfactory. The system can accurately describe the probability of interference to GPS of flights at different locations in the vicinity of the same type of UAV countermeasure equipment, and can effectively quantify the interference situation of GPS of flights at different locations by combining the probability.

Claims

1. A method for evaluating the influence of a UAV countermeasure device on a flight GPS probability model, characterized in that, The method comprises the following steps: Step one, establishing a drone countermeasure device interference probability model for flight GPS, setting the expression of the interference probability model as: P(d) = K·L(d); In the formula, d is the distance from the position to be analyzed to the drone countermeasure device, P(d) is the probability of flight GPS interference at the position with a distance d from the drone countermeasure device, L(d) is the probability that the typical receiving power of the flight receiving the signal of the drone countermeasure device exceeds the interference threshold power at the position with a distance d from the drone countermeasure device, and K is the flight GPS interference probability under the condition that the typical receiving power of the flight receiving the signal of the drone countermeasure device exceeds the interference threshold power. Step two, obtaining the broadcast automatic dependent surveillance system historical data of flights in the drone countermeasure device flight interference case in the flight interference period around the drone countermeasure device; Step three, calculating the proportion of the number of flight GPS interference flights in the distance interval with different distances from the drone countermeasure device in the total number of flights according to the broadcast automatic dependent surveillance system historical data, and constructing a probability function S(d) of flight GPS interference probability changing with distance from the proportion; Step four, obtaining the maximum value of the probability function S(d), and taking the maximum value as the estimation result of K in the expression of the interference probability model Step five, converting the probability function S(d) into a first probability distribution function T(d), the first probability distribution function T(d) corresponding to the probability that the typical receiving power of the flight receiving the signal of the drone countermeasure device is less than the interference threshold power at the position with a distance d from the drone countermeasure device, and the formula used for converting into the first probability distribution function T(d) is In the formula, K is an estimation result in an expression of the interference probability model, d is a distance from the position to be analyzed to the unmanned aerial vehicle countermeasure device, T(d) is the first probability distribution function, and S(d) is the probability function. Step six, performing variable substitution through the formula x = v(d), in which d is the distance from the position to be analyzed to the drone countermeasure device, v(d) is a variable substitution function, and x is a function independent variable reflecting distance after variable substitution, changing the first probability distribution function T(d) related to the independent variable d into a second probability distribution function Q(x) related to the independent variable x; Step seven, solving the corresponding probability density function q(x) according to the second probability distribution function Q(x); Step eight, fitting the probability density function q(x) using a parameterized function r(x) and estimating the parameters of the parameterized function r(x); Step nine, using the parameterized function r(x), and employing the formula... The estimation result of L(d) in the expression of the interference probability model is obtained, which represents the probability that the typical received power of the UAV countermeasures device signal received by a flight at a distance d from the UAV countermeasures device exceeds the interference threshold power. In the formula, v(d) is the variable substitution formula used in step six by performing variable substitution through the formula x = v(d), combined with the estimation result of K in the expression of the interference probability model obtained in step four. Substituting the values ​​into the expression of the interference probability model, we obtain the estimation result of the interference probability model. for Step ten, for the same type of UAV countermeasure equipment to be evaluated, taking the equipment as the center, in the horizontal plane with the height to be analyzed, according to the estimation result of the interference probability model The GPS interference probability distribution diagram of the flight around the UAV countermeasure equipment is obtained as the output result of the influence evaluation.

2. The method of claim 1, wherein the probability model of the flight GPS affected by the UAV countermeasure device is evaluated, and In step six, the first probability distribution function T(d) related to the independent variable d is changed into the second probability distribution function Q(x) related to the independent variable x through variable substitution through the formula x = v(d), in which d is the distance from the position to be analyzed to the drone countermeasure device, v(d) is a variable substitution function, and x is a function independent variable reflecting distance after variable substitution, and the specific process comprises the following steps: The variable substitution formula x = v(d) is specifically x = 20 log 10 d.

3. The method of claim 1, wherein the probability model of the flight GPS affected by the UAV countermeasure device is evaluated, and In step one, the expression of the interference probability model of the drone countermeasure device for flight GPS is set as P(d) = K·L(d), and the specific process comprises the following steps: In the formula, d is the distance from the position to be analyzed to the drone countermeasure device, and the distance is a straight-line distance.

4. The method of claim 1 to 3, wherein, In the third step, according to the broadcast automatic dependent surveillance system historical data, the proportion of the number of flight GPS jammed flights in the total number of flights in the distance interval at different distances from the unmanned aerial vehicle countermeasure equipment is calculated, and a probability function S(d) of flight GPS jam probability changing with distance is constructed from the proportion. The broadcast automatic dependent surveillance system historical data is preprocessed, and the preprocessing includes adjacent track point time interval normal track segment division, removal of track points obviously deviating from the normal track position, and supplement of preprocessed broadcast automatic dependent surveillance system data between adjacent track point time interval abnormal track points. Whether each track point is a jammed track point is judged by a navigation quality index value in the preprocessed broadcast automatic dependent surveillance system data. In the range of a distance not more than a maximum analysis distance from the unmanned aerial vehicle countermeasure equipment, the analysis distance from the minimum analysis distance to the maximum analysis distance is divided into different distance intervals starting from the minimum analysis distance. For a distance interval (a, b] with a minimum distance a and a maximum distance b from the unmanned aerial vehicle countermeasure equipment, all track points in the distance interval in the preprocessed broadcast automatic dependent surveillance system data are analyzed, the total number of flights corresponding to the track points is the total number of flights in the distance interval, and if the number of jammed track points of one flight in the distance interval is not less than a predetermined threshold, the flight is a flight GPS jammed flight. The ratio of the number of flight GPS jammed flights to the total number of flights in the distance interval is the proportion of the number of flight GPS jammed flights to the total number of flights in the distance interval. The proportion of the number of flight GPS jammed flights to the total number of flights in the distance interval is calculated for different distances from the unmanned aerial vehicle countermeasure equipment. The proportion of the number of flight GPS jammed flights to the total number of flights in the distance interval (a, b] is taken as the function value of the probability function S(d) of flight GPS jam probability changing with distance at the midpoint of the distance interval. The function value of the probability function S(d) at different distances is obtained by calculating the proportion of the number of flight GPS jammed flights to the total number of flights in the distance interval, so as to obtain the probability function S(d) of flight GPS jam probability changing with distance.

5. The method of claim 1 to 3, wherein, In the seventh step, according to the second probability distribution function Q(x), the corresponding probability density function q(x) is solved, specifically including: The second probability distribution function Q(x) is differentiated by Q'(x) to obtain the probability density function q(x).

6. The method of claim 1 to 3, wherein, In the eighth step, the probability density function q(x) is fitted using a parameterized function r(x), and the parameters of the parameterized function r(x) are estimated, specifically including: The parameterized function r(x) used is: In the formula, the parameterized function r(x) is a Gaussian distribution probability density function, μ is the mean of the Gaussian distribution probability density function, σ is the standard deviation of the Gaussian distribution probability density function, x is a variable-substituted function argument reflecting distance, and the parameters of the estimated parameterized function r(x) are the estimated mean μ and standard deviation σ of the Gaussian distribution probability density function.

7. The method of claim 1 to 3, wherein, In the step ten, for the same type of UAV countermeasure equipment to be evaluated, taking the equipment as the center, in the horizontal plane with the height to be analyzed, according to the estimation result of the interference probability model The GPS interference probability distribution diagram of the flight around the UAV countermeasure equipment is obtained as the output result of the influence evaluation, specifically including: In the horizontal plane with the height being the height to be analyzed, the distances d of different positions to the same type of UAV countermeasure equipment to be evaluated are calculated, and the corresponding interference probability is calculated according to The corresponding interference probability is calculated, a two-dimensional data distribution diagram represented by different colors in different interference probability value intervals is drawn, and the two-dimensional data distribution diagram is superimposed on the map of the corresponding area in a semi-transparent manner to obtain a UAV countermeasure equipment surrounding flight GPS interference probability distribution diagram, and the flight GPS interference probability distribution diagram is taken as the output result of the influence evaluation.

8. The method of claim 4, wherein the method further comprises: The range from the minimum analysis distance to the maximum analysis distance is divided into different distance intervals, starting from the minimum analysis distance, and the maximum analysis distance is not more than the maximum analysis distance from the unmanned aerial vehicle countermeasure equipment, and the method comprises the following steps: The analysis distance from the minimum analysis distance to the maximum analysis distance is divided into different distance intervals with equal length by equal division.

9. The method of claim 4, wherein the method further comprises: The navigation quality index value in the pre-processed broadcast automatic dependent surveillance system data is used to determine whether each track point is a disturbed track point, and the method comprises the following steps: According to the navigation integrity category index value in the navigation quality index in the broadcast automatic dependent surveillance system data, the track points in the track with navigation integrity category index values other than invalid values are determined, and if the navigation integrity category index value of the track point is less than the normal value threshold, the track point is determined to be a disturbed track point.

10. The method of claim 4, wherein the method further comprises: The adjacent track point time interval normal track segment division comprises the following steps: The obtained broadcast automatic dependent surveillance system historical data is divided into time stamp track point data and non-time stamp track point data according to whether the "ground station receiving position time" data item contains valid data. The time interval between adjacent track point data is calculated for the sorted track point data, and if the time interval between adjacent track points exceeds the normal time interval threshold, it is considered that the time interval between adjacent track points is abnormal, and the two adjacent track points with abnormal time interval are set as the start point and end point of track breakage.

11. The method of claim 10, wherein the method further comprises: The track filtering is performed on the track point data in each track segment of each flight obtained after the adjacent track point time interval normal track segment division, and the track filtering obtains the three-dimensional position vector prediction value, three-dimensional position vector filtering output value and three-dimensional velocity vector filtering output value according to the three-dimensional position vector filtering output value, three-dimensional velocity vector filtering output value and three-dimensional position vector observation value at the current time. The deviation between the three-dimensional position vector observation value and the three-dimensional position vector prediction value is calculated. The standard deviation of the modulus value of the deviation at all times except the first time in the same track segment is calculated, and if the modulus value at a time is greater than three times the standard deviation, the track point corresponding to the time is a track point obviously deviating from the normal track position, and the track point corresponding to the time is removed. e i = z i - p i / i-1 ; where i is an integer not less than 2, z i is a three-dimensional position vector observation value at the i-th moment, p i / i-1 is a three-dimensional position vector prediction value at the i-th moment, e i is a deviation of the three-dimensional position vector observation value from the three-dimensional position vector prediction value at the i-th moment; ​ 12. The method of claim 10, wherein the method further comprises: The preprocessed ADS-B data obtained by supplementing supplementary track points between adjacent track points with abnormal time interval, specifically comprises: After dividing the normal flight path section with the adjacent flight path point time interval, the timestamped flight path point data in each flight path section is of the same flight, between the adjacent flight path sections of the same flight, the last flight path point of the previous flight path section and the first flight path point of the next flight path section are obtained, M flight path points are sequentially taken from the last flight path point of the previous flight path section, the M is a non-negative integer, the average value Δd of M adjacent flight path point distance intervals among the M+1 flight path points is calculated M and the average value Δt of the time interval M N flight path points are sequentially taken from the first flight path point of the next flight path section, the N is a non-negative integer, the average value Δd of N adjacent flight path point distance intervals among the N+1 flight path points is calculated N and the average value Δt of the time interval N , and the average distance interval Δd between the adjacent flight path points is obtained by weighted average The average time interval Δt between adjacent track points is obtained by weighted average: Wherein M and N are 10, if the number of track points before the last track point of the previous track segment is less than 10, M takes the number of all track points before the last track point of the previous track segment, if the number of track points after the first track point of the next track segment is less than 10, N takes the number of all track points after the first track point of the next track segment; Supplement supplementary track points between the last track point of the previous track segment and the first track point of the next track segment, the number of supplementary track points is: wherein D is the distance between the last track point of a previous track segment and the first track point of a subsequent track segment, Ad is the average distance interval between adjacent track points resulting from the weighted average, and W is the number of additional track points added as a supplement, denotes the floor function. The three-dimensional position vector of the wth supplementary track point is: wherein z s is a three-dimensional position vector of the last track point of the previous track segment, z e is a three-dimensional position vector of the first track point of the subsequent track segment, w is the serial number of the additional track point, wherein 1≤w≤W, W is the number of the additional additional track points, z w is a three-dimensional position vector of the wth additional track point; The time of the wth supplementary track point is: wherein t s is the time of the last track point of the previous track segment, t e is the time of the first track point of the following track segment, w is the number of the additional track point, wherein 1≤w≤W, W is the number of the additional additional track points, t w is the time of the wth additional track point; The navigation integrity category index value of the navigation quality index of all supplementary track points is set to -2; Supplement track points between all adjacent track segments of all flights in turn; After division and removal of track points obviously deviating from normal track positions, all track segment data of all flights are time-stamped track point data, and the data and all supplementary track points corresponding to the flights jointly constitute the preprocessed ADS-B data.

13. The method of claim 11, wherein the method further comprises: The track filtering, specifically comprises: Using α-β filtering method for track filtering: where i is an integer not less than 2, j is an integer determined according to i, p i-1 / i-1 is a three-dimensional position vector filtered output value at the i-1th moment, v i-1 / i-1 is a three-dimensional velocity vector filtered output value at the i-1th moment, t i / i-1 is a time interval from the i-1th moment to the ith moment, p i / i-1 is a three-dimensional position vector predicted value at the ith moment, z i is a three-dimensional position vector observed value at the ith moment, p i / i is a three-dimensional position vector filtered output value at the ith moment, v i / i is a three-dimensional velocity vector filtered output value at the ith moment, α is a position filtering parameter, and β is a velocity filtering parameter. i-1 = 1, the first time of three-dimensional position vector filter output value p 1 / 1 and three-dimensional velocity vector filter output value v 1 / 1 The value is: In the formula, z1 is a three-dimensional position vector observation value at the first time, p 1 / 1 is a three-dimensional position vector filter output value at the first time when i-1 = 1, v 1 / 1 is a three-dimensional velocity vector filter output value at the first time when i-1 = 1.

Images