A method for fitting GNSS interference sources using civil aviation ADSB data

Abstract

This invention provides a method for fitting GNSS interference sources using civil aviation ADSB data, comprising: receiving raw messages and preprocessing them to obtain standardized data; extracting interference segments based on the interference state, the interference segments including interfered points and boundary points; integrating these to obtain a set of interference segments; cleaning the set of interference segments to obtain a cleaned result, which serves as the interference and boundary point set for cone fitting; using the interference and boundary point set for cone fitting to optimize a three-dimensional shell cone model to solve for the optimal cone parameters, the optimal cone parameters including the coordinates and half-vertex angle of the interference source; and outputting the optimal cone parameters. This method for fitting GNSS interference sources using civil aviation ADSB data results in higher physical rationality and more accurate fitting results, and it relies entirely on existing airborne data, making deployment more convenient.

Description

Technical Field

[0001] This invention belongs to the field of air traffic management, and specifically relates to a method for fitting GNSS interference sources using civil aviation ADSB data. Background Technology

[0002] GNSS interference has become a major threat to civil aviation safety. GPS interference has already led to aircraft go-arounds, air traffic control, and airport closures. With the rapid development of the low-altitude economy, GNSS interference has become a major hidden danger affecting air traffic safety, and the International Civil Aviation Organization (ICAO) has even required its member states to establish interference monitoring and protection capabilities.

[0003] Currently, traditional GNSS interference detection methods include:

[0004] 1) Dedicated radio frequency monitoring stations monitor GNSS interference frequency band signals. The cost of a single station is >500,000 yuan and a spectrum license is required. Monitoring stations detect GNSS interference through spectrum analysis. The spectrum signal is combined with the cross-location of multiple spectrum stations to determine the target. This solution is more accurate, but the cost is higher. The amount of spectrum data is large and more difficult to analyze, making it difficult to achieve comprehensive real-time monitoring.

[0005] 2) Manually piloting the calibration aircraft with a receiver to fly around it costs more than 200,000 yuan per flight and has poor real-time performance.

[0006] In addition, the impact of GNSS interference can be reduced through proactive identification and intervention, including:

[0007] 1) Deploying a carrier phase differential (TDOA / FDOA) target identification system to resolve interference requires synchronization of ≥3 stations, with an inter-station baseline of <30 km. Deployment in mountainous areas is difficult and too costly.

[0008] 2) Deploy airborne equipment with anti-interference or communication encryption capabilities. However, due to the influence of international certification, security and quality on the upgrading of airborne equipment, the corresponding application process is lengthy and costly.

[0009] ADS-B (Automatic Dependent Surveillance-Broadcast System) reverse positioning is a mandatory air-to-ground data link for civil aircraft, broadcasting position, altitude, and NUC (Navigation Uncertainty Category) fields every second. According to European control rules, NUC=0~15 represents the GPS error level, with each level representing an error distance range. Data is transmitted in the downlink ADSB during flight. NUC=0 indicates unmeasurable error, which is the most common case. 1-5 indicates excessive error leading to poor signal, while 6 and above generally indicate a normal error level. When an aircraft enters a GNSS interference zone, NUC<5, and a systematic position shift occurs, forming an "interference track cluster."

[0010] Currently, existing technologies disclose the use of this NUC field to achieve ADS-B reverse positioning, thereby enabling GNSS interference detection. For example, the interference source determination method disclosed in the literature ["Research on GNSS Interference Source Localization Based on ADS-B", Acta Aeronautica Sinica, 2021] uses the planar intersection method. Its core idea is to project the interference segment track onto the ground and obtain the planar coordinates by intersecting two-dimensional straight lines. However, it ignores the height dimension, resulting in large errors; it lacks half-vertex angle estimation, making it impossible to assess the interference beamwidth; manual mapping and positioning are time-consuming and have poor real-time performance; the results lack 3D visualization verification; and ground detection is easily misled by buildings or terrain obstruction. Summary of the Invention

[0011] The purpose of this invention is to provide a method for fitting GNSS interference sources using civil aviation ADSB data, which makes the fitting results more physically reasonable and more accurate, and is entirely dependent on existing airborne data, making deployment more convenient.

[0012] To achieve the above objectives, the present invention provides a method for fitting GNSS interference sources using civil aviation ADSB data, comprising:

[0013] S1: Receive the original message and preprocess it to obtain standardized data. The message has a field related to navigation accuracy as the NUC field.

[0014] S2: Divide the NUC field into different interference states according to the NUC threshold, extract interference segments according to the interference states, the interference segments include the interfered points and boundary points; integrate to obtain a set of interference segments;

[0015] S3: Robust removal of outliers: The set of interference segments is cleaned to obtain the cleaned result, which serves as the set of interference and boundary points for cone fitting;

[0016] S4: The three-dimensional shell cone model is fitted and optimized using the interference and boundary point set of the cone fitting to solve for the optimal cone parameters, which include the coordinates of the interference source and the half-apex angle.

[0017] Step S5: Output the optimal cone parameters.

[0018] In step S1, the original message is a civil aviation ADS-B message in ASTERIX format. The preprocessing includes message parsing, field mapping and anomaly filtering. The standardized data is a CSV format data table, which includes the aircraft's timestamp, latitude, longitude, altitude, NUC field and S-mode address code.

[0019] In step S2, the NUCs in the standardized data are classified into different interference states according to the NUC threshold, specifically including:

[0020] S21: Perform the time series sorting for a single aircraft to obtain the time series of the NUC of the single aircraft;

[0021] S22: Perform NUC binarization, divide it into different interference states, and generate an interference state sequence;

[0022] When the NUC < NUC threshold, the interference state is marked as 1, and when the NUC >= NUC threshold, the interference state is marked as 0.

[0023] In the step S2, extract the interference segments according to the interference states, specifically including:

[0024] S23: Perform a difference operation on the interference state sequence to obtain the start position and end position of the interference segment;

[0025] S24: Perform screening processing according to the start position and end position of the interference segment to obtain the final interference segment and boundary points; process the jump boundaries of the interference states, including length filtering and boundary point capture.

[0026] The step S3 specifically includes;

[0027] S31: Perform local plane projection on the coordinates of all points in the set of interference segments to achieve the coordinate conversion from the longitude and latitude coordinate system to the plane coordinate system;

[0028] S32: Execute an outlier algorithm with fixed parameters to identify outliers and valid clusters;

[0029] S33: Select the valid cluster with the largest number of points as the main interference cluster, and further use it as the interference and boundary point set for conical fitting.

[0030] In the step S31, perform local plane projection with the data center as the origin on the coordinates of all points in the set of interference segments to achieve the coordinate conversion from the longitude and latitude coordinate system to the plane coordinate system, and the data center is the average latitude and average longitude of all points in the set of interference segments;

[0031] In the step S32, the fixed parameters include the neighborhood radius and the minimum number of clustering samples.

[0032] In the step S3, set at least one outlier determination condition. If the outlier determination condition is met, perform a downgrade process, and directly use the set of interference segments output by the step S2 as the interference and boundary point set for conical fitting;

[0033] The outlier determination condition includes at least one of the following: the number of points in the set of interference segments output by the step S2 < the minimum number of clustering samples, no valid clusters, and the number of points in the main interference cluster is less than the cleaning posterior sample quantity threshold.

[0034] In step S4, the three-dimensional shell cone model is fitted and optimized using the interference and boundary point set of the cone fitting, specifically including:

[0035] S41: Initial parameter estimation;

[0036] S42: Establish an objective function to evaluate the interference between the current 3D shell-cone model and the cone fitting, and the degree of fitting of the boundary point set;

[0037] S43: Optimize the coordinates and half-vertex angle of the interference source using an optimization algorithm to obtain the optimal cone parameters when the objective function is minimized.

[0038] In step S4, the three-dimensional shell cone model is a rotationally symmetric structure, with its vertex being the ground coordinates of the interference source and the cone axis pointing vertically upwards;

[0039] The objective function includes interference point conical shell loss, boundary point loss, excessive angle penalty term, and excessive angle penalty term. The conical shell constraint term is based on the deviation of the interference point from the conical shell layer, penalizing points that deviate from the conical shell layer and rewarding points that are within the conical shell layer, thereby guiding the interference points to be distributed within the conical shell layer centered on the target conical surface and with a specified conical shell thickness, the conical shell thickness being proportional to the radius of the conical cross-section at the same height. The boundary point constraint term sets strong constraints on the boundary points, making them close to the target conical surface. The excessive angle penalty term and the excessive angle penalty term penalize the half-apex angle that exceeds the range.

[0040] In step S43, the optimization algorithm is a combination of a global optimization algorithm and a local optimization algorithm. During optimization, starting from the initial coordinates of the interference source and the half-vertex angle, the global optimization algorithm is first used to perform a global search to determine the parameter combination that is close to the global optimum. Then, using this combination as the initial value, the local optimization algorithm is used to refine the local parameters to further reduce the objective function value and finally solve for the optimal cone parameters.

[0041] This invention employs a three-dimensional shell conical model for fitting and optimization, resulting in higher physical plausibility and greater accuracy in the fitting results. In particular, the optimization of interference point conical shell loss, boundary point loss, excessive angle penalty terms, and excessively small angle penalty terms as objective functions further reduces the error in the fitting results. Furthermore, this invention utilizes only the NUC field and location field from existing airborne ADS-B data, relying entirely on existing airborne data such as ADS-B, eliminating the need for new monitoring stations or flight calibration, thus simplifying deployment. Additionally, this invention establishes a three-dimensional shell conical model with latitude, longitude, and altitude to simultaneously output the three-dimensional coordinates of the interference source and the beam half-apex angle for power / antenna type inference, resulting in low error. It also includes a built-in interactive 3D conical view with keyboard interaction for zooming, panning, and resetting. During fitting and optimization, this invention first performs geographic denoising, followed by global-local two-level optimization to ensure meter-level convergence, improve angle estimation accuracy, and better reflect actual physical conditions. An automated pipeline completes minute-level three-dimensional localization and visualization of the interference source in less than 3 minutes; airborne detection coverage is relatively comprehensive. Therefore, the method of the present invention reduces the time required to analyze interference, establishes a visual and reasonable model, and assists ground personnel in detecting and eliminating interference sources.

[0042] Therefore, this invention forms a differentiated competitive advantage in three dimensions: data requirements, physical rationality, and ease of deployment, making it particularly suitable for application scenarios that require rapid on-site deployment, beam characteristic analysis, and strong interpretability of results. Attached Figure Description

[0043] Figure 1 This is a flowchart of a method for fitting GNSS interference sources using civil aviation ADSB data according to an embodiment of the present invention.

[0044] Figure 2A and Figure 2B This is a schematic diagram of the three-dimensional shell cone model used in this invention; Figure 2A The case of a thin shell layer is shown. Figure 2B The case of a thick shell is shown.

[0045] Figure 3 This is a plan view output by the method of fitting GNSS interference sources using civil aviation ADSB data according to the present invention.

[0046] Figure 4 This is a stereoscopic view output by the method of fitting GNSS interference sources using civil aviation ADSB data according to the present invention. Detailed Implementation

[0047] The embodiments of the present invention are described in detail below. These embodiments are implemented based on the technical solution of the present invention, and provide detailed implementation methods and specific operation processes. However, the scope of protection of the present invention is not limited to the following embodiments.

[0048] The core principle of the method for fitting GNSS interference sources using civil aviation ADSB data in this invention lies in fully utilizing the physical meaning of the NUC index in ADS-B data. If the NUC index is less than a certain threshold, it directly indicates that the aircraft navigation system has detected a decline in GNSS signal quality. This information is more reliable than traditional signal strength measurements because it originates from the aircraft receiver's self-diagnosis rather than indirect observations from ground stations. Subsequently, the NUC change is converted into the entry / exit time of the interference event using a time-series differential algorithm. Combined with the spatial distribution of multiple aircraft, the three-dimensional location and radiation characteristics of the interference source can be inferred. Therefore, the method for fitting GNSS interference sources using civil aviation ADSB data in this invention processes the original ADS-B message into interference source localization results.

[0049] When GNSS interference occurs, it is generally classified into deceptive interference and suppression interference. Deceptive interference simulates a real GNSS transmitter to mislead the target's GNSS reception, causing the target's position to be redirected to another location, typically the location of the interference source. Suppression interference uses a high-power full pulse to prevent the target from acquiring real GNSS data, resulting in the target's GNSS position being an all-1 signal. The method of fitting GNSS interference sources using civil aviation ADSB data in this invention is mainly applicable to suppression interference, i.e., interference where the GNSS position is an all-1 signal, but it also has some ability to detect deceptive interference.

[0050] The method of fitting GNSS interference sources using civil aviation ADSB data according to the present invention specifically includes:

[0051] Step S1: Receive the raw ADS-B message and preprocess it to obtain standardized data;

[0052] The preprocessing includes CSV parsing of the original ADS-B messages, field mapping, and anomaly filtering. The original ADS-B messages are contained in ASTERIX format civil aviation surveillance data. CSV parsing involves parsing the ASTERIX format civil aviation surveillance data and saving the original ADS-B messages as CSV files. Field mapping involves renaming each field in the CSV file. The anomaly filtering rule is to automatically remove abnormal records with altitude <500m and airspeed <20 or >1200. Therefore, when judging GNSS interference, the method of this invention judges based on whether there are abnormal records with altitude <500m and airspeed <20 or >1200 at the existing target location, thus achieving the identification of deceptive interference. Data with altitude <500m and airspeed <20 or >1200 often indicates that the target has been deceived to other locations; the location of the interference source can be directly determined and investigated based on the location of the abnormal records.

[0053] In other embodiments, the ADS-B message can be replaced with other types of messages, as long as they have fields related to navigation accuracy to serve as the NUC field.

[0054] Among them, the data items of the ADS-B message in the civil aviation surveillance data in ASTERIX format are parsed and converted into a CSV format data table through CSV parsing as the standardized data. The CSV format data table includes the timestamp, latitude, longitude, altitude, NUC field, and S-mode address code of the aircraft. Each core field and its semantics are shown in Table 1.

[0055] Table 1: Each core field and its semantics of the converted ADS-B message

[0056]

[0057] Step S2: Intelligent extraction of interference segments: The NUC field in the standardized data is divided into different interference states according to the NUC threshold, and the interference segments are extracted according to the interference states. The interference segments include the interfered points and the boundary points; the set of interference segments is integrated.

[0058] Thus, the present invention obtains the change position of the interference state through the NUC field, and further determines the interference range through the clustering of the change positions.

[0059] In the step S2, dividing the NUC in the standardized data into different interference states specifically includes:

[0060] Step S21: Perform the time series sorting of a single aircraft to obtain the time series of the NUC of a single aircraft;

[0061] Among them, group by the S-mode address code to distinguish different aircraft, and sort in ascending order of the timestamp within the group to reconstruct a strict time series relationship.

[0062] The algorithm complexity of the present invention mainly comes from the grouping and sorting operations, and can be completed in seconds for typical data sets (hundreds of thousands of points).

[0063] Step S22: Perform NUC binarization, divide into different interference states, and generate an interference state sequence;

[0064] Among them, the interference state sequence is a binary sequence, the interference state is represented as nuc, and the interference state sequence is represented as affected_mask = (nuc_vals> 0). In this embodiment, the NUC threshold is 6. Therefore, when NUC < NUC threshold (i.e., NUC < 6), the interference state nuc is marked as 1, that is, marked as True (interference, affected), and when NUC >= NUC threshold, the interference state is marked as 0, that is, marked as False (normal).

[0065] In the existing technology, the NUC=0 state usually indicates interference, but it is a false indicator. NUC=1~5 corresponds to interference with a higher degree of directionality. NUC>=6 is a normal state without interference. NUC=-1 is an interpolation state.

[0066] Therefore, the present invention can transform continuous NUC values ​​into discrete interference states, and then locate the transition boundaries of the interference states through differential operations.

[0067] In step S2, the interference segment is extracted based on the interference state, specifically including:

[0068] Step S23: Perform a difference operation on the interference state sequence to obtain the start and end positions of the interference segment;

[0069] Step S23 specifically includes: padding the beginning and end of the interference state sequence with zeros to obtain an extended sequence; performing a first-order difference on the extended sequence to obtain the interference entry marker and interference exit marker of the difference result; the positions corresponding to the interference entry marker and interference exit marker are respectively used as the start position and end position of the interference segment, that is, as the transition boundary of the interference state.

[0070] The expression for the first-order difference is diff_mask = np.diff([0] + affected_mask +[0]), where diff_mask is the difference result of the first-order difference, affected_mask is the sequence of disturbance states, np.diff represents the first-order difference, and nuc represents the disturbance state.

[0071] In the difference result of the first-order difference, the +1 transition of the interference state sequence (i.e., diff_mask=+1) is the interference entry indicator, and the -1 transition of the interference state sequence (i.e., diff_mask=-1) is the interference exit indicator.

[0072] Step S24: Filter the interference segments based on their start and end positions to obtain the final interference segments and boundary points.

[0073] The transition boundaries of the disturbance state are processed, including length filtering and boundary point capture.

[0074] Step S24 specifically includes:

[0075] Step S241: Based on the start and end positions of the interference segments, filter out interference segments with a continuous length ≥ the length threshold, thereby achieving length filtering and eliminating instantaneous noise;

[0076] In this embodiment, the length threshold is 2.

[0077] Step S243: The start and end positions of the filtered interference segments are taken as boundary points, and the points outside the boundary points of the interference segments are the interference points, so as to achieve boundary point capture.

[0078] The boundary points are a key innovative design element of this technical solution. Their physical significance lies in the fact that the transition point from nuc=0 to >0 in the interference state corresponds to the moment the aircraft enters the boundary of the interference beam, and the transition point from nuc>0 to =0 in the interference state corresponds to the moment the aircraft leaves the boundary of the interference beam. The starting and ending positions of the interference segment serve as the "near-end" and "far-end" constraints of the conical surface during subsequent fitting, respectively.

[0079] The spatial location of these boundary points should theoretically be near the surface of the cone, providing strong constraints for subsequent half-apex angle estimation. Typical interference signals have two layers in a practical sense. The first layer is distance-related; the radiation itself has power, and attenuation results in interference with a distance range. The second layer involves a blind zone in the overhead area; the radiated waveform generally points towards the distance, creating a blind zone directly above the interference device. The principle of the interference device is the same as that of low-altitude detection equipment. Low-altitude radar can rotate (especially with this type of installation, the overhead blind zone is quite large), forming a cone-shaped interference ring. Therefore, relative to a plane, this invention can be fitted using a three-dimensional shell cone.

[0080] Therefore, this invention first determines the size of the interference range by the jump of the interference state, then captures the boundary points of the interference segment to help determine the boundary points and thus the size of the cone, and finally fits a shell cone based on the boundary points to ensure the interference range of the outer circle and the blind zone range of the inner circle.

[0081] After filtering, all the boundary points of the interference segments and the interference points together form the set of interference segments combined_df. The boundary point identifier is_boundary to distinguish whether it is a boundary point (the boundary point identifier is_boundary is a boolean field), so that different penalty terms of the objective function can be input in subsequent steps.

[0082] Step S3: Robust removal of outliers: Clean the set of interference segments to obtain the cleaned result, which serves as the set of interference and boundary points for cone fitting;

[0083] The cleaning of the set of interference segments specifically includes coordinate transformation from latitude and longitude coordinate system to plane coordinate system, DBSCAN clustering to remove outliers, and retaining the largest cluster.

[0084] Step S3 specifically includes:

[0085] Step S31: Perform local planar projection on the coordinates of all points in the set of interference segments to achieve coordinate transformation from latitude and longitude coordinate system to plane coordinate system;

[0086] Since ADS-B data uses a spherical latitude and longitude coordinate system, while the DBSCAN algorithm requires Euclidean distance measurement, in step S31, the coordinates of all points in the interference segment set are subjected to a local planar projection with the data center as the origin, to achieve coordinate transformation from the latitude and longitude coordinate system to the planar coordinate system. The data center is the average latitude and average longitude of all points in the interference segment set.

[0087] After performing local planar projection, the north and east coordinates of the points in the set of interference segments are:

[0088] North coordinates:

[0089] Eastward coordinates:

[0090] Where y is the north coordinate in meters, x is the east coordinate in meters; lat is the latitude of the current point, and lon is the longitude of the current point, which are obtained from the latitude and longitude values ​​in the ADS-B message, respectively. The latitude of the data center is the average latitude of all points in the set of interference segments. This refers to the radian value corresponding to the latitude of the data center. Used to compensate for the radius contraction effect of latitudinal circles; The longitude of the data center is the average longitude of all points in the set of interference segments; 111000 is the standard meridian length in the latitude direction, in meters per degree (the circumference of an Earth's meridian is approximately 40,008,000 meters, divided by 360 ≈ 111,133 meters / degree ≈ 111,000 meters / degree). This approximation has an error of <0.5% at a regional scale (<100km), meeting the clustering accuracy requirements.

[0091] Step S32: Execute the outlier algorithm using fixed parameters to identify outliers and valid clusters;

[0092] In this embodiment, the fixed parameters include the neighborhood radius and the minimum number of clustered samples min_samples. Points in the interference point set that satisfy all the fixed parameters are considered to be valid clusters; otherwise, they are judged as outliers.

[0093] In this embodiment, the outlier algorithm is DBSCAN spatial clustering. In other embodiments, DBSCAN spatial clustering can be replaced by other outlier algorithms, such as Isolation Forest or LOF, as long as the purpose of geographic denoising is achieved.

[0094] Step S32 finally outputs the clustering label of each point in the interference point set; if the clustering label is -1, it indicates an outlier point; otherwise, the clustering label indicates a valid cluster.

[0095] The settings and engineering basis for the fixed parameters are shown in Table 2.

[0096] Table 2: Fixed Parameter Settings and Engineering Basis

[0097]

[0098] Step S33: Select the effective cluster with the most points as the main interference cluster, and then use it as the interference and boundary point set for cone fitting.

[0099] In this process, the number of sample points in each effective cluster is counted, and the effective cluster with the most points is selected as the main interference cluster, while all others are removed.

[0100] This achieves the maximum cluster retention strategy, where the main interference cluster represents the real interference area. The rationale for this strategy is that real interference events typically affect multiple aircraft, forming a spatially dense main interference cluster; while outliers are sparsely and randomly distributed, making them easy to identify.

[0101] In step S3, multiple outlier detection conditions are set. If an outlier detection condition is met, a downgrade process is performed, and the set of interference segments output in step S2 is directly used as the interference and boundary point set for cone fitting. This ensures sufficient subsequent fitting data and avoids algorithm failure.

[0102] In this embodiment, the abnormal point determination conditions include: the number of points in the set of interference segments output in step S2 is less than the minimum clustering sample number min_samples, there are no valid clusters, and the number of points in the main interference cluster is less than the sample number threshold after cleaning (e.g., less than 10).

[0103] Step S3 of this invention uses spatial clustering to remove discrete points and clutter generated by other interference sources, retaining only the main interference cluster. This can eliminate 5%-15% of outliers while retaining more than 95% of the valid data. The center of the main interference cluster obtained after clustering is basically consistent with the results of traditional interference center estimation algorithms (such as the bisection method), and can be used as a reliable initial position for subsequent three-dimensional cone fitting.

[0104] The code example for step S3 is as follows:

[0105] # Use DBSCAN to remove outliers on a 2D plane. Inputs include latitude array lats, longitude array lons, cluster neighborhood radius eps_km (kilometers), and minimum number of cluster samples min_samples.

[0106] def remove_outliers_2d(lats, lons, eps_km=15, min_samples=3):

[0107] #If the number of data points is less than the minimum cluster size, return all data directly (without performing clustering).

[0108] if len(lats) <min_samples:

[0109] return lats, lons, np.ones(len(lats), dtype=bool)

[0110] # Convert latitude and longitude to an approximate planar coordinate system (for clustering only)

[0111] # Calculate the latitude (center_lat) and longitude (center_lon) of the data center

[0112] center_lat = np.mean(lats)

[0113] center_lon = np.mean(lons)

[0114] #Using coordinate transformation formulas, we obtain the north coordinates (y) and east coordinates (x) of points in the set of interference segments (for clustering only).

[0115] x = (lons - center_lon) * 111000 * np.cos(np.radians(center_lat))

[0116] y = (lats - center_lat) * 111000

[0117] # Combine into a two-dimensional coordinate point set [x, y]

[0118] points = np.column_stack([x, y])

[0119] # DBSCAN spatial clustering, obtaining cluster labels

[0120] clustering = DBSCAN(eps=eps_km*1000, min_samples=min_samples).fit(points)

[0121] # Count the cluster labels and number of all clusters to obtain the effective clusters.

[0122] unique_labels, counts = np.unique(clustering.labels_, return_counts=True)

[0123] if -1 in unique_labels: # Remove outliers based on cluster labels

[0124] noise_idx = np.where(unique_labels == -1)[0][0]

[0125] counts = np.delete(counts, noise_idx)

[0126] unique_labels = np.delete(unique_labels, noise_idx)

[0127] # If no valid clusters are found, return all data.

[0128] if len(unique_labels) == 0:

[0129] return lats, lons, np.ones(len(lats), dtype=bool)

[0130] # Select the effective cluster with the most points as the main interference cluster, representing the actual interference area, and mark the points belonging to the main interference cluster.

[0131] main_cluster_label = unique_labels[np.argmax(counts)]

[0132] inliers = clustering.labels_ == main_cluster_label

[0133] # Returns: Latitude, longitude, and mask of all points in the main interference cluster

[0134] return lats[inliers], lons[inliers], inliers.

[0135] Step S4: The three-dimensional shell cone model is fitted and optimized using the interference and boundary point set obtained from the cone fitting to solve for the optimal cone parameters (lat, lon, θ). The optimal cone parameters include the coordinates of the interference source and the half-apex angle. Thus, the three-dimensional cone optimization fitting of the GNSS interference source is realized.

[0136] The three-dimensional shell-cone model is a rotationally symmetric structure. The vertices of the three-dimensional shell-cone model represent the ground coordinates of the interference source, and the cone's axis points vertically upwards. In this embodiment, the generatrix of the three-dimensional shell-cone model during rotation is a straight line, and the equation of the conical surface is r = h•tanθ, where h is the aircraft altitude and r is the cross-sectional radius at altitude h. In other embodiments, the three-dimensional shell-cone model with a straight generatrix can be replaced with other types of cone models, such as three-dimensional shell-cone models with parabolic, elliptical, or piecewise linear generatrixes, also fitted using nonlinear optimization.

[0137] The three-dimensional radiation characteristics of the interference source are described using a three-dimensional shell-cone model, with the interference source as the vertex of the cone. The geometric structure of the three-dimensional shell-cone model is shown in Figure 2, and the geometric parameters are shown in Table 3.

[0138] Table 3: Geometric parameters of the three-dimensional shell cone model

[0139]

[0140] When optimizing a 3D shell cone model using the interference and boundary point set obtained from cone fitting, for any observation point P(lat, lon, h), its radial deviation ε from the candidate cone is defined as: ε = d - h·tanθ, where d is the horizontal distance between the horizontal coordinates of the observation point and the ground coordinates of the interference source, θ is the half-apex angle, and h is the height. That is, d = geodesic((lat, lon), (lat0, lon0)), (lat, lon) are the horizontal coordinates of the observation point, (lat0, lon0) are the ground coordinates of the interference source, and geodesic is the horizontal distance from the ground.

[0141] The specific steps for fitting and optimizing a 3D shell conical model using interference and boundary point sets include:

[0142] Step S41: Initial parameter estimation: Use the horizontal coordinates of the data center of the interference and boundary point set as the initial coordinates of the interference source, and determine the initial half-apex angle based on the 90th percentile distance / height ratio of the interference and boundary point set.

[0143] In addition, step S41 may also include estimating the shell thickness: setting an initial shell thickness t based on the average deviation from the interference point to the conical surface, providing a benchmark for subsequent reward items.

[0144] In this embodiment, the search boundary and initial values ​​of the parameters to be estimated are shown in Table 4.

[0145] Table 4: Search Boundaries and Initial Values ​​of Parameters

[0146]

[0147] Specifically, extending the boundary by 50% of the data span ensures that the true solution is included; the initial angle value uses the 90th quantile distance / height ratio × 1.5, which both suppresses the influence of extreme values ​​and reserves adjustment space for optimization. Here, lat is latitude, and lon or lng is longitude.

[0148] Step S42: Establish an objective function to evaluate the interference between the current 3D shell cone model and the cone fitting and the degree of fitting of the boundary point set;

[0149] The objective function is used to evaluate the degree of fit between the current three-dimensional shell conical model and the cone fitting interference and boundary point set. All points of the cone fitting interference and boundary point set should be inside the three-dimensional shell conical model and have a headspace blind zone.

[0150] The objective function includes the interference point conical shell loss, boundary point loss, excessive angle penalty term, and excessive angle penalty term. The conical shell constraint term is based on the deviation of the interference point from the conical shell layer, penalizing points that deviate from the shell layer and rewarding points within it. This guides the interference points to be distributed within a conical shell layer centered on the target cone surface and with a specified thickness, where the shell thickness is proportional to the radius of the cone section at the same height. The boundary point constraint term imposes strong constraints on boundary points, ensuring they are close to the target cone surface. The excessive angle penalty term and the excessive angle penalty term penalize half-apex angles that exceed reasonable ranges to ensure the physical feasibility of the model.

[0151] The expression for the objective function is:

[0152] + + max + max ,

[0153] in, For the point of interference, conical shell layer loss, Let P be the set of disturbance points in the disturbance and boundary point set, B be the set of boundary points in the disturbance and boundary point set, i be the ordinal number of the disturbance point, and j be the ordinal number of the boundary point. Let be the horizontal distance between the i-th interference point and the ground line of the interference source. Let be the horizontal distance between the j-th boundary point and the ground line of the interference source. Let be the radius of the cone at the height of the i-th interference point. Let be the radius of the cone at the height of the j-th boundary point. , These are the weighting coefficients for the penalty terms for excessively large and small angles. The weight coefficients for the boundary point loss are max. As a penalty for excessive angle, max Penalty for angles that are too small. It is a semi-vertical angle. The upper limit of the semi-vertical angle, As the lower limit of the half-vertex angle, in this embodiment, , .

[0154] In this embodiment, the interference point conical shell layer loss for:

[0155] ,

[0156] in, Let be the horizontal distance between the i-th interference point and the ground line of the interference source. Let t be the radius of the cone at the height of the i-th interference point, and t be the thickness of the cone shell. The inside-the-cone point is rewarded, and the rest are outside-the-cone point penalties.

[0157] Boundary point loss for:

[0158] ,

[0159] in, Let be the horizontal distance between the i-th interference point and the ground line of the interference source. Let t be the radius of the cone at the height of the i-th interference point, and t be the thickness of the cone shell.

[0160] It should be noted that 0.8, 1.2, and 0.2 in the above formula are actually adjustable parameters based on the shell thickness t. In this embodiment, using the boundary distance * 1 as the standard, the shell thickness t allows for an allowable error range of ±20%. Therefore, 0.8 to 1.2 is considered to be the target on the boundary surface of the three-dimensional shell conical model. If it is less than 0.8, it is within the blind zone, and the penalty weight for loss needs to be increased; if it is greater than 1.2, the penalty weight for loss is generally lower because the receiver may not be deactivated. If it is inside the three-dimensional shell conical model, a reward weight for loss is set and a reward is given.

[0161] In other embodiments, many parameters of the objective function have a certain degree of variability. The present invention aims to protect a scheme that utilizes the interference point conical shell loss and boundary point loss objective functions to jointly obtain the objective function, and then uses it to solve the problem. For example, the cone thickness can be described as thickness or blind zone angle; the specific values ​​of the parameters are arbitrary. The scheme mainly involves simulating the function through a three-dimensional inverse conical graph using fitted graphics and known data sources. In other embodiments, when the interference source device is a gun-type directional reflector, the side-conical interference caused by a directional antenna can also be fitted by changing the objective function (increasing the cone's start and end angles).

[0162] Weighting coefficients of boundary point loss Employing high weights, preferably 100, for the weighting coefficients of penalty terms for excessively large and small angles. , It adopts an ultra-high weight, preferably 1000.

[0163] The implementation code for the objective function is as follows:

[0164] #Optimized cone fitting function

[0165] def fit_single_cone_optimized(lats, lons, alts, boundary_lats=None, boundary_lons=None, boundary_alts=None):

[0166] # Improved initial estimate

[0167] center_lat_init = np.mean(lats)

[0168] center_lon_init = np.mean(lons)

[0169] # Calculate the initial half-vertex angle: based on all points, not just the maximum point.

[0170] dists = []

[0171] for lat, lon, alt in zip(lats, lons, alts):

[0172] if alt > 0:

[0173] dist = calculate_distance(lat, lon, center_lat_init,center_lon_init)

[0174] dists.append(dist / alt)

[0175] if dists:

[0176] # Use the 90th percentile instead of the maximum value

[0177] theta_init = np.arctan(np.percentile(dists, 90)) * 1.5 # Magnify by 50%

[0178] theta_init = min(theta_init, np.radians(85)) # Limit the maximum angle

[0179] theta_init = max(theta_init, np.radians(5)) # Limit the minimum angle

[0180] else:

[0181] theta_init = np.radians(45) # Default 45 degrees

[0182] print(f"Initial estimates: Center({center_lat_init:.6f}, {center_lon_init:.6f}), Angles{np.degrees(theta_init):.2f}°")

[0183] # Improved objective function

[0184] def objective(params):

[0185] center_lat, center_lon, theta, shell_thickness = params

[0186] # Main penalty: The incentive point is inside the cone

[0187] cost = 0

[0188] n_points = len(lats)

[0189] for i in range(n_points):

[0190] dist = calculate_distance(lats[i], lons[i], center_lat,center_lon)

[0191] cone_radius = alts[i] * np.tan(theta)

[0192] # Calculate the distance to the conical shell

[0193] dist_to_surface = abs(dist - cone_radius)

[0194] # Annular region logic: The point should be close to the conical shell surface and far from the central axis

[0195] if dist < cone_radius * 0.8: # Enter the blind zone (too deep inside the cone)

[0196] cost += (cone_radius * 0.8 - dist) ** 2 * 10 # Penalty

[0197] elif dist > cone_radius * 1.2: # Beyond the interference range (too far outside the cone)

[0198] cost += (dist - cone_radius * 1.2) ** 2 # Penalty

[0199] else:

[0200] # Near the shell (±20%), reward

[0201] cost += -0.5 * (shell_thickness - dist_to_surface) ** 2

[0202] # Boundary point handling: Encourage near the conical surface

[0203] if boundary_lats is not None and len(boundary_lats) > 0:

[0204] for i in range(len(boundary_lats)):

[0205] dist = calculate_distance(boundary_lats[i], boundary_lons[i], center_lat, center_lon)

[0206] cone_radius = boundary_alts[i] * np.tan(theta)

[0207] # The boundary point should be within ±10% of the cone surface.

[0208] surface_dist = abs(dist - cone_radius)

[0209] Ideally, the boundary point should be on the surface of the cone.

[0210] ideal_surface_ratio = 1.0

[0211] surface_ratio = dist / cone_radius if cone_radius> 0else 1.0

[0212] # Increase penalty if the boundary point is too far from the surface.

[0213] If abs(surface_ratio - ideal_surface_ratio) > 0.2: # Allow a 20% deviation

[0214] cost += (surface_ratio - ideal_surface_ratio) **2 * 100

[0215] # Angle regularization: Prevents angles from being too large or too small

[0216] cost += max(0, theta - np.radians(85)) ** 2 * 1000 # Penalty for angles greater than 85 degrees

[0217] cost += max(0, np.radians(5) - theta) ** 2 * 1000 # Penalty for angles less than 5 degrees

[0218] Return cost.

[0219] Step S43: Use an optimization algorithm to optimize the coordinates and half-vertex angle of the interference source, and solve for the optimal cone parameters when the objective function is minimized.

[0220] In this embodiment, the optimization algorithm is a combination of global and local optimization algorithms. During optimization, starting with the initial coordinates and half-apex angle of the interference source, a global search is first performed using a global optimization algorithm (such as differential evolution) to determine a parameter combination close to the global optimum. Then, using this combination as the initial value, a local optimization algorithm (such as the quasi-Newton method L-BFGS-B) is used for local refinement to further reduce the objective function value, ultimately solving for the optimal cone parameters, including the coordinates of the interference source and the half-apex angle of the beam. In other embodiments, the optimization algorithm can be replaced with one or more combinations of Particle Swarm Optimization (PSO), Simulated Annealing (SA), Bayesian optimization, differential evolution, and quasi-Newton methods.

[0221] The execution of the Differential Evolution (DE) algorithm specifically includes:

[0222] Step A1: Initialize the population: Using the initial parameters as the center, generate multiple candidate parameter combinations (individuals) within a reasonable range.

[0223] Step A2: Iterative Evolution: Through mutation, crossover, and selection operations, individuals with lower objective function values ​​are continuously selected, gradually approaching the global optimal region.

[0224] Step A3: Termination condition: When the preset number of iterations is reached or the population converges, output the current optimal parameter combination as the initial value for subsequent local optimization algorithms.

[0225] This enables global search, avoiding sensitivity to initial values ​​and local optima.

[0226] When executing the differential evolution algorithm, the core parameters are popsize=20, maxiter=200, and seed=42. The population size of 20 corresponds to 3 optimization variables × approximately 7 individuals / variable, allowing for parallel exploration within a broad parameter space.

[0227] The quasi-Newton method L-BFGS-B utilizes gradient information (numerical difference approximation) to achieve superlinear convergence. Two-stage collaboration results in total function evaluations requiring less than 7000 iterations and taking less than 5 seconds in typical scenarios. The core parameters are maxiter=500 and ftol=1e. -10 .

[0228] Step S5: Result Output and Visualization: Output the optimal cone parameters.

[0229] The optimal cone parameters include the coordinates of the interference source and the half-apex angle. In some embodiments, the optimal cone parameters may also include the shell thickness. The output optimal cone parameters are in CSV parameter form.

[0230] Step S5 may also include:

[0231] Fit quality assessment: Quality indices are calculated to verify the fitting effect. These indices include, for example, the proportion of points within the cone, the proportion near the surface, the average distance error, the maximum distance error, and the true position deviation, quantifying the reliability of the fitting results. The five-dimensional quality indices are shown in Table 5.

[0232] Abnormal result filtering: If the optimal cone parameters exceed the physically reasonable range (e.g., the half-apex angle is outside the range, or the positional deviation is too large), it is judged as a fitting failure and a prompt is returned.

[0233] Rendering of Planar / Stereoscopic Views: Based on the optimal conic parameters and aircraft trajectory data, corresponding planar / stereoscopic views are drawn and displayed. The planar / stereoscopic views are 300dpi PNG images, meeting the printing quality requirements of technical reports and academic publications. The planar view is implemented using the `visualize_fitting_results_3d` function, and its output is shown in Table 6; the stereoscopic view is implemented using the `create_detailed_cone_visualization` function, and its output is shown in Table 7.

[0234] Table 5: Five-Dimensional Quality Indicators

[0235]

[0236] These metrics are automatically calculated after optimization and encapsulated in a results dictionary, supporting automated filtering and manual review.

[0237] Table 6: Output Contents of the Plan View

[0238]

[0239] Table 7: Output Contents of Stereoscopic Views

[0240]

[0241] The four-dimensional encoding system for aircraft track points is the core innovation of this invention for rendering planar / stereoscopic views: spatial position (x, y, z) determines geometric coordinates; NUC value (>0 / =0) determines color (red / green family); cone-inside / outside determination determines marker shape (× or o / ^ or .); aircraft identifier determines specific color (tab20 color wheel). This encoding makes the interference state, geometric relationship, and attribution of each point immediately clear.

[0242] The final output results are shown in Table 8.

[0243] Table 8: Output Results

[0244]

[0245] Based on the above-mentioned method for fitting GNSS interference sources using civil aviation ADSB data, the system for fitting GNSS interference sources using civil aviation ADSB data of the present invention includes a data preprocessing module, an interference segment extraction module, a hybrid optimization module, and a result output module.

[0246] The data preprocessing module is configured to receive raw ADS-B messages and preprocess them to obtain standardized data for the ADS-B messages. For specific steps, please refer to step S1 above.

[0247] The interference segment extraction module is configured to: divide the NUC field in the standardized data into different interference states according to the NUC threshold; extract interference segments according to the interference states, wherein the interference segments include disturbed points and boundary points; and integrate them to obtain a set of interference segments. Optionally, the interference segment extraction module can also be configured to: clean the set of interference segments to obtain the cleaned result, which serves as the interference and boundary point set for cone fitting; or, directly use the set of interference segments as the interference and boundary point set for cone fitting. Specific steps can be found in steps S2 and S3 above.

[0248] The hybrid optimization module is configured to: use the interference and boundary point set obtained from cone fitting to fit and optimize the three-dimensional shell cone model in order to solve for the optimal cone parameters (lat, lon, θ), which include the coordinates of the interference source and the half-apex angle; thus, the three-dimensional cone optimization fitting of the GNSS interference source is realized. For specific steps, please refer to step S4 above.

[0249] The output module is set to output the optimal cone parameters.

[0250] Optionally, the results output module can also be configured as follows: Fit quality assessment: Calculate quality indicators to verify the fitting effect. Quality indicators include, for example, the proportion of points inside the cone, the proportion near the surface, the average distance error, the maximum distance error, and the true position deviation, to quantitatively assess the reliability of the fitting results. Abnormal result filtering: If the optimal cone parameters exceed the physically reasonable range (e.g., the half-apex angle exceeds the interval, or the position deviation is too large), it is determined that the fitting has failed and a prompt is returned.

[0251] For specific steps, please refer to step S5 above.

[0252] These modules implement data flow through standard Python data structures (NumPy arrays, pandas DataFrames), with low coupling between modules, facilitating independent testing and iterative optimization. In particular, the hybrid optimization module and the result output module can form a closed-loop feedback—optimization results directly drive the calculation of quality indicators, while quality indicators can serve as the basis for convergence judgment and result selection.

[0253] The innovations of this invention in fitting GNSS interference sources using civil aviation ADSB data include improvements to the objective function and the optimization strategy. Unlike traditional geometric fitting that aims to minimize the point-to-surface distance, the objective function of this technical solution deeply integrates the physical generation mechanism of GNSS interference.

[0254] Internal-external asymmetry: Points inside the cone receive a negative penalty (reward), encouraging the "inclusive" expansion of the solution; points outside the cone receive a positive penalty, with the penalty increasing as the distance increases. This design avoids the suboptimal behavior of shrinking the cone to accommodate individual edge points.

[0255] Boundary point constraints: Transforming the temporal information of NUC jumps into spatial surface constraints significantly improves the physical plausibility of half-vertex angle estimation. Experimental results show that introducing boundary constraints can improve surface_ratio by 10-20 percentage points.

[0256] Multi-scale penalty: Distance-related weight adjustment (1→5 times) enables adaptive suppression of severe outliers without the need for a preset hard threshold.

[0257] Furthermore, the optimization algorithm of this invention combines the DE algorithm with L-BFGS-B, and the improvement in performance is shown in Table 9.

[0258] Table 9: Improvement effect of DE algorithm combined with L-BFGS-B

[0259]

[0260] The hybrid strategy employs a two-stage division of labor: DE is responsible for "coarse localization" to a high-quality solution neighborhood, while L-BFGS-B is responsible for "fine refinement" to machine accuracy. This achieves an optimal balance between efficiency and accuracy in typical scenarios.

[0261] Furthermore, the innovations in the specific implementation process of this invention are shown in Table 10.

[0262] Table 10: Innovations in Engineering Implementation

[0263]

[0264] Actual running results:

[0265] The validation results based on a typical dataset (Nanjing Lukou Airport 30nm range scene, multiple interference sources) are as follows: Figures 3-4 As shown in Tables 11-13. Among them, Figure 4 The diagram shows NUC anomalies, fitting interference sources, the lower conical blind zone and the upper spherical interference boundary, as well as the near-shell boundary and the far-shell boundary. Figure 4The red area represents the lower cone-shaped blind zone, while the blue area represents the upper spherical interference attenuation boundary. In other words, the interference region is the portion of the three-dimensional shell cone model whose distance from the interference source is less than the interference attenuation threshold.

[0266] Table 11: Result Error

[0267]

[0268] Table 12: Fitting Quality Indicators

[0269]

[0270] Table 13: Processing Efficiency

[0271]

[0272] As shown in Tables 10-12, the method of fitting GNSS interference sources using civil aviation ADSB data of the present invention has achieved excellent results in terms of result error, fitting quality index and processing efficiency.

[0273] This invention employs a three-dimensional shell conical model for fitting and optimization, resulting in higher physical plausibility and greater accuracy in the fitting results. In particular, the optimization of interference point conical shell loss, boundary point loss, excessive angle penalty terms, and excessively small angle penalty terms as objective functions further reduces the error in the fitting results. Furthermore, this invention utilizes only the NUC field and location field from existing airborne ADS-B data, relying entirely on existing airborne data such as ADS-B, eliminating the need for new monitoring stations or flight calibration, thus simplifying deployment. Additionally, this invention establishes a three-dimensional shell conical model with latitude, longitude, and altitude to simultaneously output the three-dimensional coordinates of the interference source and the beam half-apex angle for power / antenna type inference, resulting in low error. It also includes a built-in interactive 3D conical view with keyboard interaction for zooming, panning, and resetting. During fitting and optimization, this invention first performs geographic denoising, followed by global-local two-level optimization to ensure meter-level convergence, improve angle estimation accuracy, and better reflect actual physical conditions. An automated pipeline completes minute-level three-dimensional localization and visualization of the interference source in less than 3 minutes; airborne detection coverage is relatively comprehensive. Therefore, the method of the present invention reduces the time required to analyze interference, establishes a visual and reasonable model, and assists ground personnel in detecting and eliminating interference sources.

[0274] Therefore, this invention forms a differentiated competitive advantage in three dimensions: data requirements, physical rationality, and ease of deployment, making it particularly suitable for application scenarios that require rapid on-site deployment, beam characteristic analysis, and strong interpretability of results.

[0275] The above description is merely a preferred embodiment of the present invention and is not intended to limit the scope of the invention. Various variations can be made to the above embodiments of the present invention. All simple and equivalent changes and modifications made in accordance with the claims and description of this application fall within the protection scope of the claims of this patent. All aspects not described in detail in this invention are conventional technical content.

Claims

1. A method for fitting GNSS interference sources using civil aviation ADSB data, characterized in that, Including: Step S1: Receive the original message and preprocess it to obtain standardized data. The message has a field related to navigation accuracy as the NUC field. Step S2: Divide the NUC field into different interference states according to the NUC threshold, extract the interference segments according to the interference states. The interference segments include interference points and boundary points; integrate to obtain a set of interference segments. Step S3: Clean the set of interference segments to obtain a cleaning result as the interference and boundary point set for conical fitting. Step S4: Use the interference and boundary point set for conical fitting to fit and optimize the three-dimensional shell conical model to solve the optimal conical parameters. The optimal conical parameters include the coordinates and semi-apex angle of the interference source. Step S5: Output the optimal conical parameters.

2. The method for fitting GNSS interference sources using civil aviation ADSB data according to claim 1, characterized in that, In step S1, the original message is an ASTERIX format civil aviation ADS-B message. The preprocessing includes message parsing, field mapping, and exception filtering. The standardized data is a CSV format data table. The CSV format data table includes the timestamp, latitude, longitude, altitude, NUC field, and S-mode address code of the aircraft.

3. The method for fitting GNSS interference sources using civil aviation ADSB data according to claim 1, characterized in that, In the step S2, divide the NUC in the standardized data into different interference states according to the NUC threshold, specifically including: Step S21: Perform single-aircraft time series sorting to obtain the time series of the NUC of a single aircraft; Step S22: Perform NUC binarization, divide it into different interference states, and generate an interference state sequence; When NUC < NUC threshold, the interference state is marked as 1; when NUC >= NUC threshold, the interference state is marked as 0.

4. The method for fitting GNSS interference sources using civil aviation ADSB data according to claim 3, characterized in that, In the step S2, extract the interference segments according to the interference states, specifically including: Step S23: Perform a difference operation on the interference state sequence to obtain the start and end positions of the interference segments; Step S24: Perform screening processing according to the start and end positions of the interference segments to obtain the final interference segments and boundary points; process the jump boundaries of the interference states, including length filtering and boundary point capture.

5. The method for fitting GNSS interference sources using civil aviation ADSB data according to claim 1, characterized in that, Step S3 specifically includes; Step S31: Perform local plane projection on the coordinates of all points in the set of interference segments to achieve coordinate conversion from the longitude-latitude coordinate system to the plane coordinate system; Step S32: Execute an outlier algorithm with fixed parameters to identify outliers and valid clusters; Step S33: Select the valid cluster with the largest number of points as the main interference cluster, and further use it as the interference and boundary point set for conical fitting.

6. The method for fitting GNSS interference sources using civil aviation ADSB data according to claim 5, characterized in that, In step S31, perform local plane projection with the data center as the origin on the coordinates of all points in the set of interference segments to achieve coordinate conversion from the longitude-latitude coordinate system to the plane coordinate system. The data center is the average latitude and average longitude of all points in the set of interference segments. In the step S32, the fixed parameters include the neighborhood radius and the minimum number of clustering samples.

7. The method for fitting GNSS interference sources using civil aviation ADSB data according to claim 5, characterized in that, In the step S, set at least one outlier determination condition. If the outlier determination condition is met, perform downgrading processing, and directly use the set of interference segments output by step S2 as the interference and boundary point set for conical fitting; The anomaly determination criteria include at least one of the following: the number of points in the set of interference segments output in step S2 is less than the minimum clustering sample number, there are no valid clusters, and the number of points in the main interference cluster is less than the threshold of the number of samples after cleaning.

8. The method for fitting GNSS interference sources using civil aviation ADSB data according to claim 1, characterized in that, In step S4, the three-dimensional shell cone model is fitted and optimized using the interference and boundary point set of the cone fitting, specifically including: Step S41: Initial parameter estimation; Step S42: Establish an objective function to evaluate the interference between the current 3D shell cone model and the cone fitting and the degree of fitting of the boundary point set; Step S43: Use an optimization algorithm to optimize the coordinates and half-vertex angle of the interference source, and solve for the optimal cone parameters when the objective function is minimized.

9. The method for fitting GNSS interference sources using civil aviation ADSB data according to claim 8, characterized in that, In step S4, the three-dimensional shell cone model is a rotationally symmetric structure, with its vertex being the ground coordinates of the interference source and the cone axis pointing vertically upwards; The objective function includes interference point conical shell loss, boundary point loss, excessive angle penalty term, and excessive angle penalty term. The conical shell constraint term is based on the deviation of the interference point from the conical shell layer, penalizing points that deviate from the conical shell layer and rewarding points that are within the conical shell layer, thereby guiding the interference points to be distributed within the conical shell layer centered on the target conical surface and with a specified conical shell thickness, the conical shell thickness being proportional to the radius of the conical cross-section at the same height. The boundary point constraint term sets strong constraints on the boundary points, making them close to the target conical surface. The excessive angle penalty term and the excessive angle penalty term penalize the half-apex angle that exceeds the range.

10. The method for fitting GNSS interference sources using civil aviation ADSB data according to claim 8, characterized in that, In step S43, the optimization algorithm is a combination of a global optimization algorithm and a local optimization algorithm. During optimization, starting from the initial coordinates of the interference source and the half-vertex angle, the global optimization algorithm is first used to perform a global search to determine the parameter combination that is close to the global optimum. Then, using this combination as the initial value, the local optimization algorithm is used to refine the local parameters to further reduce the objective function value and finally solve for the optimal cone parameters.

Images