A satellite navigation deception jamming recognition method and system based on multi-feature fusion

Abstract

The application discloses a satellite navigation spoofing interference recognition method and system based on multi-feature fusion and relates to the technical field of satellite navigation security. The application comprises the following steps: estimating mutual information of channel carrier-to-noise ratio (CNR) gradient change and carrier motion parameters to determine physical layer confidence; modeling the distribution form of pseudo-range residual on line to extract bimodal and thick tail significance and determine information layer confidence, and generating environment state quantity according to the over-limit state of bimodal and thick tail significance of each channel; reconstructing a multi-dimensional error observation vector phase space of a tracking loop to calculate a dynamics complexity index, and determining dynamics layer confidence and triggering a hierarchical early warning according to the relationship between the dynamics complexity index and an early warning threshold. The application extracts three heterogeneous features, i.e., mutual information, bimodal / thick tail significance and maximum Lyapunov exponent, from the signal power domain, pseudo-range residual domain and loop state domain respectively, realizes essential distinction between spoofing and multipath, and reduces false alarm and missed alarm rates.

Description

Technical Field

[0001] This invention relates to the field of satellite navigation security technology, specifically to a method and system for identifying satellite navigation deception and interference based on multi-feature fusion. Background Technology

[0002] Satellite navigation systems are widely used in critical applications such as time synchronization in power systems, unmanned aerial vehicles (UAVs) for transmission line inspections, automated inspection robots for substation equipment, and power emergency communications. However, satellite navigation signals are weak when they reach the ground, making them highly vulnerable to deception and interference attacks. Deception interference generates fake signals with structures similar to real satellite signals, gradually pulling the receiver's tracking loop away from the true signal. This causes the receiver to output incorrect positioning or timing results, potentially leading to serious consequences such as power system time synchronization failure, PMU (Phasor Measurement Unit) data errors, UAV inspection deviations, and abnormal positioning of substation automation equipment. Therefore, accurate and timely identification of deception interference in satellite navigation applications within the power industry is an urgent requirement for ensuring the safe and stable operation of the power system.

[0003] In typical power system application environments, such as densely packed metal structures in substations, mountainous environments along transmission line corridors, urban underground utility tunnels, and power tunnels, metal structures, hillsides, and tunnel walls generate dense multipath reflections. The received signal is simultaneously affected by multipath fading and potential spoofing attacks. Current technologies typically employ the following methods to detect spoofing interference: monitoring the absolute level of the carrier-to-noise ratio (CNR) of each channel and comparing it with a fixed threshold; statistically monitoring pseudorange residuals to detect anomalies; or using signal quality monitoring indicators to analyze signal snapshots at individual moments. These methods have a certain detection capability against spoofing signals with sudden power changes in open environments.

[0004] However, the aforementioned existing technologies have the following shortcomings in the complex and dynamic environment of the power industry: relying on single or a small number of signal features makes it difficult to distinguish between strong multipath and generative deception features, which can easily lead to false alarms and missed alarms; multi-feature fusion uses fixed weights or static thresholds, which makes it difficult to adapt to the dynamic changes in the signal environment and results in poor recognition performance; and snapshot-style discrimination is performed on the signal state at a single moment, which does not make sufficient use of the temporal evolution features of multiple parameters in the deception process and lacks early warning capabilities. Summary of the Invention

[0005] The purpose of this invention is to provide a satellite navigation spoofing and interference identification method and system based on multi-feature fusion, so as to solve the problems mentioned in the background art.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a satellite navigation deception interference identification method based on multi-feature fusion, comprising: estimating the mutual information between the carrier-to-noise ratio change gradient of each channel and the carrier motion parameters, and determining the physical layer confidence level; modeling the distribution pattern of the pseudorange residual online, extracting the bimodal saliency and the heavy-tailed saliency, determining the information layer confidence level, and generating an environmental state quantity based on the out-of-limit state of the bimodal saliency and the heavy-tailed saliency; reconstructing the phase space of the multidimensional error observation vector of the tracking loop, calculating the dynamic complexity index, and determining the dynamic layer confidence level and triggering a graded warning based on the relationship between the dynamic complexity index and the warning threshold; adjusting the warning threshold and the fusion weights of each layer according to the environmental state quantity, weighting and fusing the confidence levels of each layer according to the adjusted fusion weights, and performing graded decision-making based on the fusion result to identify satellite navigation deception interference.

[0007] Furthermore, the carrier motion parameter is selected by default as velocity amplitude. When the variance of the velocity amplitude within the sliding window is lower than a preset variance threshold, the carrier motion parameter is switched from velocity amplitude to displacement increment. The estimation of mutual information includes: establishing a mapping model between carrier noise ratio and satellite elevation angle, subtracting the slowly varying component caused by elevation angle change from the carrier noise ratio time series to obtain the rapidly varying component, obtaining the carrier noise ratio change gradient by sliding window difference of the rapidly varying component, and estimating the mutual information between the carrier noise ratio change gradient and the carrier motion parameter based on the K-nearest neighbor method.

[0008] Furthermore, when the speed amplitude within the sliding window is lower than the preset speed threshold and remains so for a preset duration, active micro-motion excitation is initiated.

[0009] Furthermore, the online modeling involves fitting a bi-Gaussian mixture model to the pseudorange residuals using an online expectation-maximization algorithm with a forgetting factor; the bimodal significance is the ratio of the difference between the mean values ​​of the Gaussian components to the mean standard deviation, and the heavy-tailed significance is the ratio of the proportion of samples whose residual tails exceed a preset threshold to the probability of the tails in the standard normal distribution.

[0010] Furthermore, the generation of the environmental state quantity includes: each channel performs single-channel determination based on the excessive states of bimodal significance and heavy-tailed significance; combining the determination results of all channels, a strong multipath environment is generated when the number of channels with strong multipath determination exceeds a preset proportion threshold; a deception competition lock environment is generated when the number of channels with deception competition lock determination exceeds a preset proportion threshold; otherwise, a normal environment is generated; the switching of the environmental state quantity is processed by hysteresis anti-jitter; after combining the determination results of each channel, the candidate state of the environmental state quantity must be continuously maintained at a preset number of sampling points before the switching is confirmed; and a hysteresis interval is set when exiting.

[0011] Furthermore, the dynamic complexity index is the maximum Lyapunov exponent: after removing the mean and normalizing the multidimensional error observation vector of the tracking loop, the phase space is reconstructed using the delayed coordinate embedding method. In the reconstructed phase space, the divergence distance of each neighboring point pair within a preset number of evolution steps is tracked. The logarithm of the divergence distance is taken and the average value is calculated. Then, a linear fit is performed with the number of evolution steps, and the resulting slope is the maximum Lyapunov exponent. The warning threshold includes a zero-reaching threshold and a negative-turning threshold. When the maximum Lyapunov exponent is lower than the zero-reaching threshold, a first-level warning is triggered, and when it is lower than the negative-turning threshold, a second-level warning is triggered.

[0012] Furthermore, the adjustment of the environmental state quantity to the warning threshold is as follows: under strong multipath environment, it is adjusted to the reference threshold minus the preset correction amount; under deceptive contention locking environment, it is adjusted to the reference threshold plus the preset correction amount; under normal environment, the reference threshold is taken.

[0013] Furthermore, the fusion weights of each layer are highest in a strong multipath environment (physical layer), highest in a deception-competitive locking environment (information layer), and equal in a normal environment. The hierarchical decision-making includes: when the fusion result is in a first preset interval, it is determined to be suspicious, the monitoring frequency is increased to a preset multiple, and the decision period is shortened to a preset duration; when the fusion result is in a second preset interval, it is determined to be highly suspicious, a first-level alarm is issued, and preparations are made to switch to a backup navigation source; when the fusion result is in a third preset interval, it is determined to confirm deception, a second-level alarm is issued, a backup navigation source is switched, and a security handling procedure is triggered.

[0014] Furthermore, when the receiver is configured with dual antennas, the carrier motion parameter is the relative displacement increment of the two antennas due to the change in carrier attitude, the carrier-to-noise ratio change gradient is the difference between the carrier-to-noise ratio change gradients of the two antennas, and the reference threshold of mutual information is the position-related threshold indexed by the position estimated by the carrier track. The position-related threshold is determined based on the mutual information statistics of each grid in the fixed route when the historical decision is that there is no deception.

[0015] A satellite navigation spoofing and interference identification system based on multi-feature fusion includes: The physical layer feature extraction module is used to estimate the mutual information between the carrier-to-noise ratio change gradient of each channel and the carrier motion parameters, and to determine the physical layer confidence level. The information layer feature extraction module is used to model the distribution pattern of pseudorange residuals online, extract bimodal significance and heavy-tailed significance, determine the information layer confidence, and generate environmental state quantities based on the out-of-limit state of the bimodal significance and the heavy-tailed significance. The dynamic layer feature extraction module is used to reconstruct the phase space of the multidimensional error observation vector of the tracking loop, calculate the dynamic complexity index, and determine the dynamic layer confidence and trigger graded warnings based on the relationship between the dynamic complexity index and the warning threshold. The fusion decision module is used to adjust the early warning threshold and the fusion weight of each layer according to the environmental state quantity, weight the confidence of each layer according to the adjusted fusion weight, and perform hierarchical decision-making based on the fusion result to identify satellite navigation deception interference.

[0016] Compared with the prior art, the beneficial effects of the present invention are: A satellite navigation deception and interference identification method and system based on multi-feature fusion solves the problem of physical layer features being unable to be extracted under extremely low dynamic conditions by initiating active micro-motion excitation when the speed is below a preset threshold and introducing controllable motion when the self-motion excitation is insufficient to ensure that the mutual information is estimable. By setting a hysteresis anti-shake mechanism, requiring the feature quantity to continuously maintain a preset number of sampling points when the limit is exceeded, and setting a hysteresis interval when exiting, the problem of threshold and weight oscillation caused by frequent switching of environmental states is avoided, thereby improving the stability of the decision.

[0017] Meanwhile, by extracting three heterogeneous features—mutual information, bimodal / fat-tailed saliency, and maximum Lyapunov exponent—from the signal power domain, pseudorange residual domain, and loop state domain, respectively, the essential distinction between deception and multipath propagation is achieved, reducing false alarm and missed alarm rates. An environmental state variable is generated by integrating the out-of-limit states of bimodal and fat-tailed saliency across all channels, and this environmental state variable simultaneously drives the adaptive adjustment of the fusion weights and decision thresholds of each layer, maintaining stable recognition performance under different environmental conditions. By reconstructing the phase space of the tracking loop error state and calculating the maximum Lyapunov exponent, the dynamic phase transition law in the deception traction process is captured, and two-level warning thresholds are set. A graded warning is issued before the positioning error exceeds the limit, reserving response time for the navigation system. Attached Figure Description

[0018] Figure 1 This is a schematic diagram of the overall method flow of the present invention; Figure 2 This is a schematic diagram illustrating the evolution of the pseudorange residual distribution morphology in this invention; Figure 3 This is a schematic diagram of the dual adaptive adjustment logic of the environmental state quantity E driving threshold and weight in this invention; Figure 4 This is a schematic diagram of a deception interference identification scenario according to Embodiment 1 of the present invention; Figure 5(a) is a schematic diagram comparing the flight trajectories of the power transmission line inspection drone in Embodiment 1 of the present invention; Figure 5(b) is a schematic diagram of the evolution curve of the comprehensive confidence level C over time in Embodiment 1 of the present invention; Figure 6(a) is a schematic diagram of the substation inspection area layout in Example 2; Figure 6(b) is a schematic diagram of the substation inspection robot's navigation and multipath environment in Example 2; Figure 7(a) is a schematic diagram comparing the physical layer differential mutual information and the location fingerprint threshold in Embodiment 2 of the present invention; Figure 7(b) is a schematic diagram of the evolution of fusion decision confidence and alarm triggering in Embodiment 2 of the present invention; Figure 8 This is a schematic diagram comparing the running trajectories of the substation inspection robot in Embodiment 2 of the present invention. Detailed Implementation

[0019] 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, and 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.

[0020] Example 1: This example uses a UAV for power transmission line inspection in a densely multipath urban canyon environment as its application. The UAV flies between urban buildings 50m to 150m high, where glass curtain walls and metal facades generate dense multipath reflections. Its flight speed is 5 to 15m / s, allowing it to perform hovering tasks. The GNSS receiver is configured with multiple constellations and frequencies (GPS L1 / L5 and BDS B1 / B2), and auxiliary sensors include an IMU (100Hz update frequency), a barometric altimeter, and a magnetometer. For example... Figure 1 As shown, the present invention provides a technical solution: a satellite navigation spoofing and interference identification method based on multi-feature fusion, comprising the following steps: Step 1: Define a multidimensional heterogeneous feature set.

[0021] Each channel constitutes a multidimensional feature vector F, denoted as F=[I,B,T,λ]. max ].

[0022] I represents the signal power domain characteristics. The mutual information between the carrier-to-noise ratio (CNR) gradient and motion parameters is used to determine whether there is a physical causal coupling relationship between the received signal and the carrier motion: the fading spectrum structure of real multipath fading has a deterministic coupling with the carrier motion, while the spoofing signal is actively injected from a single source, and its CNR shows a smooth change along a single path, without the above coupling with the motion. This distinguishes between real signals and spoofing signals.

[0023] B and T are pseudorange residual domain features. B is the bimodal significance, defined as the ratio of the difference between the means of the two Gaussian components in the pseudorange residual bi-Gaussian mixture fitting to the average standard deviation, indicating the degree of deceptive competitive locking. T is the heavy-tailed significance, defined as the ratio of the proportion of samples whose tails in the pseudorange residual distribution exceed a preset threshold to the theoretical value of the standard normal distribution, indicating the degree of strong multipath environment. B and T are used to distinguish between deceptive competitive locking and strong multipath: in deceptive competitive locking, the residuals exhibit a bimodal distribution, and B increases significantly; in strong multipath, the residuals exhibit a heavy-tailed distribution, and T increases significantly.

[0024] The B and T values ​​of each channel together generate the environmental state variable E, which further drives the adaptive adjustment of the threshold and weights.

[0025] λ max The loop state domain characteristics are defined by the maximum Lyapunov exponent, calculated after reconstructing the phase space of the multidimensional error observation vectors of the receiver tracking loop. This exponent measures the phase transition of the system from high-complexity behavior to deterministic low-complexity behavior: λ during normal tracking. max When λ is positive, the deception signal pulls the signal. max The signal collapses towards zero, and the deception signal lock is completed when λ is λ. max Turning to a negative value can provide an early warning before the positioning error exceeds the safety envelope.

[0026] The heterogeneous features of the above three observation domains (signal power domain, pseudorange residual domain, and loop state domain) are converged and then adaptively adjusted by E-driven threshold and weights to finally output a hierarchical comprehensive decision.

[0027] Step 2: Physical layer spoofing detection, detrending preprocessing of carrier-to-noise ratio and calculation of mutual information.

[0028] Carrier-to-noise ratio (CNR) detrending preprocessing. Obtain the CNR time series C / N0(t) for each channel from the receiver correlator output (typical update frequency 1Hz), and establish a third-order polynomial mapping model between CNR and elevation angle:

[0029] Where θ is the satellite elevation angle in degrees; C / N0(θ) is the fitted carrier-to-noise ratio in dB-Hz; and a0, a1, a2, and a3 are fitting coefficients. This model characterizes the slowly varying deterministic components determined by the atmospheric attenuation path and antenna gain pattern. The model is fitted in segments based on the elevation angle within a preset interval (e.g., every 5 degrees), with the fitting data taken from long-term historical observation data over a preset time period (e.g., the past hour). The rapidly varying component δC / N0(t) is obtained by subtracting the slowly varying component caused by changes in elevation angle from the original C / N0.

[0030] It is important to note that the C / N0 trend caused by elevation angle changes is a slowly varying deterministic component; C / N0 increases as the elevation angle increases and decreases as the elevation angle decreases. If this component is not subtracted, the slowly varying trend of the carrier-to-noise ratio with elevation angle will be included in the rapidly varying component, introducing a trend component unrelated to multipath fading into the mutual information estimation. This leads to an overestimation of the mutual information, weakening the ability to detect spoofing signals. The residual δC / N0 after detrending mainly contains the rapidly varying component caused by motion-coupled multipath fading and motion-independent thermal noise.

[0031] Calculate the carrier-to-noise ratio gradient. Perform sliding window differencing on the δC / N0 sequence. The sliding window length W1 is determined based on the signal fading characteristics to cover multiple complete fading periods, and the difference interval is... The coherence time is greater than that of multipath fading. For example, in a typical urban canyon power transmission line inspection drone flight scenario, the drone speed is 5-15 m / s, the carrier-to-noise ratio is high, and the fast fading period is about 2-5 seconds. A W1=20-second window can cover at least 4 complete fading periods, ensuring statistical representativeness; the coherence time of multipath fading is usually 1-3 seconds, and the differential interval is... The time interval can effectively extract the trend-based changes in carrier-to-noise ratio caused by motion, rather than instantaneous noise fluctuations. Carrier-to-noise ratio change gradient. (Unit: dB / second) Calculate using the following formula:

[0032] Among them, δC / N0(t), t and The carrier-to-noise ratio (CNR) fast-changing component after time-detrended processing, in dB; The difference interval is expressed in seconds. This formula converts the absolute level of the carrier-to-noise ratio (CNR) into a trend, eliminating the impact of CNR differences between different satellite channels.

[0033] Extract the time series of motion parameters. The velocity vector magnitude V(t) (unit: m / s) is selected by default as the motion parameter. After being acquired by the IMU, it is low-pass filtered (cutoff frequency 0.5Hz) and downsampled to align with the carrier-to-noise ratio timestamp. When the variance of V(t) within the sliding window is lower than a preset variance threshold (for constant speed cruise), the motion parameter is switched to a displacement increment sequence. (Unit: m):

[0034] Where P(t) and t and The position vector of the time carrier in the navigation coordinate system is obtained by IMU track estimation or wheel speed measurement.

[0035] It should be noted that under uniform speed conditions, the velocity amplitude is constant and the mutual information tends to zero. Switching to displacement increment can ensure that the mutual information can be continuously estimated under any non-static conditions.

[0036] Mutual information is estimated quickly using the K-nearest neighbor method. Mutual information I( The KSG estimator can capture the nonlinear coupling between multipath fading and motion, outperforming the Pearson correlation coefficient, which can only capture linear relationships. The nearest neighbor parameter k is selected based on the bias-variance trade-off principle under sample size constraints; in this embodiment, k=4. The window slides at a preset step size (e.g., every 5 seconds), and a moving average is performed on the estimates of several consecutive (e.g., 3) adjacent windows to suppress small sample estimation variance.

[0037] The specific calculation process for mutual information I is as follows: The N=20 sample points within the sliding window form a two-dimensional dataset {( V i )}, V represents the gradient of the carrier-to-noise ratio change at the i-th sample point. i Let be the velocity amplitude of the i-th sample point.

[0038] For each sample point i, find its k-th nearest neighbor in two-dimensional space and record the Chebyshev distance ε. i .exist The maximum difference between the two dimensions, V and ε, is taken as ε. i This ensures that both dimensions have equal weight in distance calculation.

[0039] Statistics based on ε i Within a neighborhood with a radius of 2, respectively in Number of samples n in dimension x,i and the number of samples n in the V dimension y,i . Specifically, n x,i To satisfy | |<ε i / 2 is the number of sample points j; n y,i To satisfy |V i -V j |<ε i The number of sample points j is 2 / 2.

[0040] The mutual information I is estimated using the following formula:

[0041] in, is the digamma function; k is the K-nearest neighbor parameter; N is the total number of sample points within the sliding window. The first term is the baseline term determined by the K-nearest neighbor parameter k; the summation average term is the average of the digamma function values ​​of the number of neighborhood samples in both dimensions for each sample point, reflecting the degree of dependence between the two variables in the actual data. The term is a finite sample correction term determined by the total number of samples N. This formula estimates mutual information by comparing the statistical properties of the number of neighboring samples in two dimensions: if... If n is independent of V, then x,i and n y,i The statistical properties of n are consistent with the expected value determined solely by the marginal distribution, and the mutual information is close to zero; if there is a statistical dependence between the two, then n x,i and n y,i The joint distribution deviates from the independence assumption, the mutual information is positive, and the larger the value, the stronger the coupling.

[0042] The above mutual information estimation is based on the natural logarithm. The original unit of the estimation result is kilotons. In this embodiment, it is divided by ln2 to convert it to bits before being used in the threshold comparison.

[0043] Adaptive threshold comparison. Adaptive threshold T I (Unit: bits) Calculated using the following formula:

[0044] Among them, T I,base The threshold is determined based on the quantile of the mutual information distribution in typical scenarios. For example, in a typical urban canyon environment, the 3% quantile of the mutual information distribution in 100 real signal scenarios is 0.15 bits. Using this value as the threshold can reserve a detection margin while ensuring a real signal pass rate of over 95%. The correction amount is determined by the environmental state quantity E, and the specific value is issued in step 3.

[0045] If I < T I If the channel outputs a suspicious flag, then the channel will be flagged. I =1 indicates that the coupling relationship between the channel signal and the carrier motion is weaker than the expected level of a normal multipath signal, suggesting possible deception; if I≥T I Then flag I =0 indicates that the signal characteristics of this channel are normal.

[0046] Calculate the physical layer confidence score. Physical layer confidence score c phy The distance between mutual information and the threshold is calculated using a Sigmoid-type mapping function:

[0047] Where, α phyThe physical layer kurtosis parameter is set to a preset value (e.g., 5), with a transition bandwidth approximately ±30% of the threshold value; I represents the currently calculated mutual information; T I For adaptive mutual information threshold; T I,base This is the baseline threshold. phy The value ranges from [0,1], with larger values ​​indicating a stronger likelihood of being flagged as deception. The formula uses T... I,base Normalization makes the mapping scale uniform across different threshold settings. When I=T I c phy =0.5, corresponding to the decision boundary.

[0048] Active excitation compensation under low dynamic operating conditions. When V(t) is detected to be less than a preset speed threshold V... th When the speed reaches 2 m / s and continues for a preset duration (e.g., 30 seconds), the carrier is determined to be in a hovering or extremely low-speed condition, and active micro-motion excitation is initiated: a small-amplitude sinusoidal position disturbance signal is injected into the flight controller through the flight control system, causing the antenna to oscillate with amplitude A in the horizontal plane. exc =8cm, frequency f exc The parameters were selected based on the following: an amplitude of 8 cm, approximately 0.4 times the GPS L1 carrier wavelength, which can cause a significant change in carrier phase exceeding 1 / 4 wavelength, resulting in detectable fluctuations in C / N0; and a frequency of 0.3 Hz, corresponding to a period of approximately 3.3 seconds, with approximately 3 sampling points per period at a sampling rate of 1 Hz, satisfying the Nyquist sampling theorem and falling within the typical multipath fading frequency range of 0.2 to 0.5 Hz, thus exhibiting an excitation effect similar to natural multipath fading.

[0049] Optionally, if the receiver supports a carrier-to-noise ratio output rate of 5 Hz or higher, the excitation frequency can be increased to 1.5 Hz; if the carrier is equipped with dual antennas, the relative motion of the two antennas due to attitude changes can be used as the excitation source to calculate the differential mutual information of the carrier-to-noise ratio fluctuation of the two antennas.

[0050] Step 3: Information layer deception detection, online fitting of pseudorange residuals using a dual Gaussian mixture model and extraction of environmental state variables E.

[0051] This step involves online modeling of the statistical distribution of pseudorange residuals, extracting two features: bimodal significance (B) and heavy-tailed significance (T). This distinguishes between deceptive competitive locking and strong multipath propagation. The environmental state variable E is generated from B and T and fed back to the threshold adjustments in steps 2 and 4, as well as the fusion weight adjustment in step 5. The specific operation steps are as follows.

[0052] Online fitting of a dual Gaussian mixture model. A dual Gaussian mixture model is applied to the pseudorange residual r(t) for each channel. The pseudorange residual r(t) is defined as the difference between the pseudorange observation and the pseudorange prediction calculated based on the current positioning result. The probability density function p(r) of the dual Gaussian mixture model is:

[0053] Where π1 and π2 are the mixing weights of the first and second Gaussian components, respectively, with values ​​ranging from [0,1], and π2 = 1 - π1; N(r|μ,σ 2 ) is a expression with mean μ and variance σ. 2 The Gaussian probability density function; μ1 and μ2 are the means of the two Gaussian components, in units of m; σ1 2 σ2 2 The variance of the two Gaussian components, in m. 2 This dual-Gaussian mixture model uses a mixture of two Gaussian distributions to simulate the empirical distribution of pseudorange residuals, simultaneously characterizing three distribution patterns: normal unimodal (π1≈1 or π2≈1), deceptive competitive locking bimodal (π1 and π2 are both significantly greater than 0 and μ1 and μ2 are separated), and strong multipath with thick tails (increased variance of the two components and increased probability density of the tail). The evolution of the residual distribution pattern is tracked in real time through online estimation of model parameters.

[0054] The parameters, including the forgetting factor λ, are estimated using the online expectation-maximization (EM) algorithm. EM The effective memory length is set according to the deception traction characteristic time. The effective memory length covers the key stages of the deception traction process (power increases at a rate of 0.5 dB / s, and the process from initial intervention to complete lock-up typically lasts 10 to 20 seconds). In this embodiment, λ EM =0.92.

[0055] The specific calculation process of the online expectation-maximization algorithm is as follows: Initialization: flag I When π = 0, it runs as a single Gaussian model, π1 = 1, π2 = 0, μ1 is the mean of the current residual sliding window, and σ1 2 Take the variance of the sliding window. (flag) I When μ = 1, the double Gaussian model is activated, and the current residual samples are divided into two groups according to their positive and negative signs. The mean and variance of the positive group are used as μ1 and σ1, respectively. 2 The initial values, the mean and variance of the negative group are taken as μ2 and σ2. 2 The initial values ​​are π1=π2=0.5.

[0056] Expected step: For each new residual sample r t Calculate the posterior probability γ of it belonging to the j-th Gaussian component (j=1,2). j (t):

[0057] Where, π j μ is the mixing weight of the j-th Gaussian component; jσ is the mean of the j-th Gaussian component; j 2 Let be the variance of the j-th Gaussian component.

[0058] Maximize the step: with the forgetting factor λ EM Update the parameters of each component recursively.

[0059] Valid sample count update:

[0060] Where, n j (t) represents the number of effective cumulative samples of the j-th component at time t; n j (t-1) represents the number of effective cumulative samples of the j-th component at time t-1; λ EM =0.92. This recursive formula means that the contribution of each new sample to the number of valid samples is its posterior probability γ. j (t), the contribution of historical samples is λ EM The decay factor decays exponentially.

[0061] Hybrid weight update:

[0062] Mean update:

[0063] Where, μ j (t), μ j (t-1) are the mean values ​​of the j-th Gaussian component at times t and t-1, respectively; Indicates the new sample r t The weight of the contribution to the mean update varies with the number of valid samples n. j (t) decreases as γ increases, and decreases with increasing γ. j (t) increases as it increases.

[0064] Variance update:

[0065] Where, σ j 2 (t) and σ j 2 (t-1) represents the variance of the j-th Gaussian component at times t and t-1, respectively.

[0066] Extract the distribution morphological features B and T.

[0067] The significance of a bimodal distribution (B) is defined as the ratio of the difference between the means of the two Gaussian components to their average standard deviation.

[0068] This formula normalizes the bimodal distance to a dimensionless measure B in units of the mean standard deviation. B > 1 indicates that the bimodal distance exceeds one standard deviation, and the two components are statistically distinguishable; the larger B is, the higher the significance of the bimodal distance, and the higher the possibility of deceptive competitive locking.

[0069] The significance of the fat tail, T, is defined as the actual tail probability, P. e The tail probability P of standard normal theory th The ratio:

[0070] The tail boundary threshold is ±k t ·σ s ;k t The tail-bound coefficient is k. t =2.5; σ s This represents the standard deviation obtained by fitting the current residual sequence using a single Gaussian model. The actual tail probability P... e This represents the proportion of samples in the current residual sequence that exceed this threshold; since the threshold varies with σ s Scaling, the theoretical probability of actually falling into this interval is always equal to the probability of exceeding ±k under the standard normal distribution. t The probability, therefore P th P is a constant. th =2×[1-Φ(k t )], where Φ(·) is the cumulative distribution function of the standard normal distribution. k t When P = 2.5, th ≈0.0124. T=1 indicates that the tail thickness is consistent with the expected Gaussian distribution; T>1 indicates the presence of a thick tail feature, and the larger the value, the more significant the thick tail and the higher the probability of strong multipath.

[0071] like Figure 2 The diagram shown illustrates the evolution of the pseudorange residual distribution. Figure 2 In the middle (a), the residuals exhibit a unimodal distribution centered at zero during normal tracking; μ is the mean of the pseudorange residuals under normal tracking conditions, which are distributed in a single Gaussian pattern. Figure 2 In (b), after the deception signal is introduced, the tracking loop alternately locks between the true and false correlation peaks, and the residual splits into a bimodal distribution with two Gaussian components superimposed. The distance between the two peaks increases with the increase of traction, and the significance B of the bimodal distribution increases. μ1 is the mean of the first Gaussian component in the bi-Gaussian mixture model, corresponding to the residual peak formed by the real signal. μ2 is the mean of the second Gaussian component in the bi-Gaussian mixture model, corresponding to the residual peak formed by the deception signal. Figure 2 After the deception in (c) completes the locking, the loop is completely pulled to the deception signal, and the residual converges back to a single peak, but its peak center has shifted from the zero point (i.e., the theoretical mean position of the true residual) to... This offset represents the equivalent pseudorange error caused by the deceptive signal. Figure 2 As can be seen from the three-stage evolution in (a), (b), and (c), the residuals after normal tracking and deception locking are both unimodal. The unimodal shape or amplitude of the residuals alone cannot distinguish between the two. However, the bimodal characteristics of the deception competition locking stage can serve as a priori signs of deception identification and provide a basis for graded early warning.

[0072] Asymmetric threshold linkage decision. Setting a deception competition lock detection threshold B. th =1.2, strong multipath detection threshold T th =2.0. The threshold selection is based on the following: when the spoofing signal power is approximately 3dB higher than the real signal, the tracking loop begins to alternately lock onto the two correlation peaks of the true and spoof signals. The mean difference between the residual double peaks is approximately 0.3 to 0.5 times the GPS C / A chip width. After thermal noise and multipath broadening, the standard deviation of each component is approximately 0.1 to 0.2 times the chip width, corresponding to a B value between 1.5 and 5.0. In this embodiment, B is chosen. th =1.2 is below the lower limit of this range, and detection can be triggered when the deception signal power advantage is about 2dB. This allows for early warning in the early stage of traction when the deception signal power is still below 3dB, reserving more response time for subsequent fusion decision and navigation source switching. th It is determined based on measured data from extreme multipath environments in densely packed metal structures of substations, ensuring reliable identification in dense multipath conditions while avoiding false triggering in moderate multipath scenarios.

[0073] After extracting B and T independently for each channel, a single-channel determination is made based on the out-of-limit status: T > T th And B≤B th When determining strong multipath, B > B th When determining deceptive competition lock-in, B≤B th And T≤T th This is a normal judgment; when B > B th And T > T th At that time, based on the principle of safety and conservatism, it is preferentially judged as a deceptive competition lock. A global environment state variable E is generated by combining the judgment results of all channels: when the number of channels with strong multipath determination exceeds a preset proportion threshold (50% in this embodiment), E is set to E0. mp When the number of channels involved in the deception competition lockout exceeds a preset proportional threshold, E is set to E. sp In all other cases, E is set to E. norm .

[0074] It is important to note that B th and T th As the criterion for determining the environmental state variable E, it is kept at a fixed value in this embodiment and is not included in the E-driven adaptive adjustment loop. Optionally, when long-term changes in the application environment cause an overall shift in the normal residual statistical characteristics, it can be adjusted by setting E=E. normThe state continues for longer than the preset long window T long (Recommended value 300 seconds) Update threshold: Set T th Updated to max(T) th,initial ,μ T +2σ T ), where T th,initial As the initial value, μ T and σ T These represent the mean and standard deviation of T within the long window, respectively; the updated threshold is applied in subsequent T... long The window remains unchanged, and an evaluation is only conducted at the end of the next long window to determine if an update is needed. th It does not participate in adaptive adjustment to avoid security vulnerabilities caused by the B value being continuously too high during deception attacks, which could lead to the threshold being incorrectly raised and resulting in missed alarms. It also prevents the system from being unable to exit on its own in a contention-locked state.

[0075] The value of the environmental state variable E directly determines the subsequent threshold correction and fusion weight configuration. If E frequently jumps between different states due to instantaneous fluctuations in feature variables B or T, the threshold and weights will oscillate accordingly, affecting the stability of the decision. Therefore, a hysteresis anti-jitter mechanism is introduced into the state switching logic of E. Specifically, after synthesizing the decision results from all channels, the candidate state of E needs to be continuously maintained at N... h Switching is confirmed only after 10 sampling points (corresponding to 10 seconds, approximately 2 to 4 times the fast fading period of the carrier-to-noise ratio), filtering out transient noise or short-term multipath anomalies that cause occasional single-channel overshooting and resulting E candidate state fluctuations; a hysteresis interval is set upon exiting, from E... mp Exit requirement T < T th -h T (h) T Take T th 15%, or 0.3), from E sp Exit requirement B < B th -h B (h) B Take B th The hysteresis value is set at 15% of the threshold (i.e., 0.18) to prevent the characteristic from repeatedly crossing the boundary when it fluctuates slightly near the threshold. Setting the hysteresis value to 15% of the threshold value strikes a balance between preventing oscillations and maintaining sensitivity during state switching.

[0076] Based on the value of E, a threshold correction amount is generated and distributed to the mutual information threshold calculation formula in step 2 and the dynamic complexity threshold calculation formula in step 4, respectively. The generation of the threshold correction amount specifically includes: when E = E... mp At that time, the threshold is lowered to reduce false alarms in strong multipath environments, such as... = -0.05 bits, =-0.003. When E=E spAt the same time, the threshold is raised to improve the detection sensitivity of deceptive competition locking phase and trigger early warnings to prevent missed alarms, such as... =+0.05 bits, =+0.003. When E=E norm At that time, each threshold value remains unchanged from the baseline value.

[0077] Calculate the information layer confidence level c based on the judgment result. info .

[0078] When classified as a strong multipath environment:

[0079] When classified as a competitive locked environment:

[0080] Where, α info This is the information layer kurtosis parameter, with a value of 3. When the classification is normal, c... info The value is fixed at 0.1. This value is an empirical calibration value, representing the baseline confidence level when the information layer is in a normal state. This ensures that the information layer contributes a negligible basic component in the fusion weighting, avoiding the influence of c. info =0 results in the complete disconnection of that branch during the fusion operation. The information layer confidence level characterizes the degree of anomaly in the residual domain: c in a strong multipath environment info Although it increases with T, E mp The state reduces its fusion weight to 0.15, and at this point, the physical layer confidence level c... phy confidence level c of the dynamic layer dyn All are at low levels, and the overall confidence level C will not exceed the deception judgment threshold, thus avoiding multipath misjudgment as deception from the fusion mechanism.

[0081] The environmental state variable E drives the dual adaptive adjustment logic of threshold and weight as follows: Figure 3 As shown.

[0082] Step 4: Dynamics layer deception detection, reconstructing the phase space of the tracking loop error observation vector and providing graded early warning based on the dynamics complexity index.

[0083] This step performs nonlinear dynamic analysis on the multidimensional error observation vector of the receiver tracking loop. Through phase space reconstruction and complexity index calculation, it monitors the phase transition process of the system from a high-complexity state to a deterministic low-complexity state, and achieves graded early warning before the positioning error exceeds the safety envelope. The specific operation steps are as follows.

[0084] A multidimensional error observation vector is constructed and preprocessed. A three-dimensional observation vector S(t) is extracted from the receiver tracking loop, including: code phase error (output from the delay-locked loop discriminator, unit: chip); carrier frequency error (output from the frequency-locked loop or phase-locked loop discriminator, unit: Hz); and carrier-to-noise ratio (CNR), unit: dB-Hz. The rationale for selecting these three components is that code phase error and carrier frequency error reflect the tracking loop's locking state to the signal in the time and frequency domains, respectively; CNR is an observation quantity of the loop rather than an internal state quantity. Including these components in the observation vector allows the phase space to simultaneously contain information from the time, frequency, and power domains, improving the multidimensional observability of spoofing signals. S(t) is then subjected to mean subtraction (subtracting the preset duration W2 = 60 seconds moving average) and normalization (dividing by the moving standard deviation of each dimension) to eliminate the influence of dimensional differences on the phase space geometry. W2 is only used for preprocessing statistics and is related to λ. max The lengths L of the observation sequences on which the calculations are based are independent of each other.

[0085] The phase space of the preprocessed S(t) is reconstructed using the delayed coordinate embedding method (based on Takens embedding theorem).

[0086] The embedding dimension *m* is adaptively determined by the spurious nearest neighbor method within a preset range (e.g., 4 to 8). The original observation vector dimension is 3. According to Takens' embedding theorem, *m* ≥ 2d + 1 (where *d* is the system's degrees of freedom). The tracking loop has at least 2 degrees of freedom, therefore *m* is at least 5. The search range of 4 to 8 can cover different operating conditions from low degrees of freedom (pure signal tracking without multipath) to high degrees of freedom (multipath reflection introducing additional degrees of freedom in a strong multipath environment). The spurious nearest neighbor method's determination process is as follows: For candidate embedding dimension *m*, find the nearest neighbor of each point in the reconstructed phase space and calculate the distance; after increasing the embedding dimension to *m + 1*, recalculate the distance of the same pair of nearest neighbors in the new dimension. If the ratio of the new distance to the original distance exceeds the threshold *R*, the spurious nearest neighbor method will be used to determine the nearest neighbor. fnn If the percentage is 15%, then the nearest neighbor is considered a false neighbor. The percentage of false neighbors for all points is calculated, and when the percentage drops below 5%, the current m is considered to be the minimum sufficient embedding dimension.

[0087] The delay time τ is determined by the mutual information method. The mutual information I(τ) between the original multidimensional sequence S(t) and its delayed sequence S(t+τ) is calculated to measure the self-dependency of the same sequence under different time delays. The optimal delay time is taken as the value of τ when I(τ) first drops to 1 / e of I(τ=0). The physical meaning of this selection criterion is that when the delay time is τ, a suitable amount of statistical correlation is still retained between the sequences before and after the delay. It is neither highly redundant (too small τ leads to almost overlapping coordinate axes) nor completely independent (too large τ leads to loss of deterministic structure). Thus, in the reconstructed phase space, each coordinate axis can provide independent information while jointly characterizing the dynamic structure of the system. For a typical receiver tracking loop with a 1Hz error output, τ is usually taken as 3 to 8 seconds.

[0088] The reconstructed phase space state vector X(t) = [S(t), S(t-τ), S(t-2τ), ..., S(t-(m-1)τ)] has a dimension of m×3 and contains complete information on the current and historical states of the tracking loop.

[0089] Calculate the maximum Lyapunov exponent λ max . λ max >0 indicates exponential separation between adjacent orbits, λ max <0 indicates that the orbit is converged.

[0090] It should be noted that the normal tracking loop is driven by random factors such as thermal noise, and its λ max A positive value is not strictly equivalent to chaos, but rather a representation of highly complex stochastic behavior; this invention uses λ max The relative evolutionary trend of collapsing from positive to zero to negative values, rather than the absolute sign, is used as the criterion for deception traction. When the deception signal is completely locked, the tracking loop degenerates into a deterministic system, λ. max It turns into a negative value.

[0091] λ max The observation sequence length L on which the calculation is based is 120 seconds, and the sliding step size is 5 seconds. This length satisfies L > (m-1) × τ max +p, i.e., 120 > 7 × 8 + 10 = 66, and retains approximately 1 statistical margin for nearest neighbor search and linear fitting. λ is calculated using the Rosenstein small data method. max In the reconstructed phase space, for each reference point X(t) i Find its nearest neighbor X(t) j Set a minimum time interval constraint |t i -t j |>T mean / 4(T) mean(The average period is estimated using the autocorrelation function), ensuring that the selected nearest neighbor pairs come from different orbital segments in phase space. Track the divergence distance d of each pair of nearest neighbors after p steps of evolution. k (p) (e.g., Euclidean distance, k=1,2,...,K,K=5); the average of all K pairs of nearest neighbors is obtained.<lnd(p)> ,right<lnd(p)> Perform a least-squares linear fit with the number of evolution steps p:

[0092] Where d0 is the initial average distance; The time interval between adjacent steps is 1 second; the slope of the fitted line is λ. max Setting the number of nearest neighbor pairs K to 5 can achieve a balance between averaging the effects of local noise and preserving divergent statistical characteristics. Setting the number of evolution steps p to 10 corresponds to an evolution window of about 10 seconds, which is smaller than the characteristic time of deception traction (10 to 20 seconds) to ensure that changes in orbital divergence characteristics can be captured before traction is completed.

[0093] Two-level threshold for tiered early warning. The adaptive threshold is calculated using the following formula:

[0094]

[0095] Among them, T λ1 This is the first-level warning threshold (zero-approach judgment threshold); T λ2 This is the second-level warning threshold (the threshold for determining whether the situation turns negative). The correction amount issued in step three; T λ1,base and T λ2,base As the baseline threshold, based on λ in normal scenarios max Distribution statistics determine: the first-level benchmark threshold corresponds to the low-end quantile of the normal distribution; values ​​below this threshold indicate a significant decrease in system complexity from the normal high-value range. The second-level benchmark threshold is slightly below zero, indicating that the system has moved from the high-complexity region to the deterministic region, and the deception signal is about to be locked or has already been locked. In this embodiment, T is taken as... λ1,base =0.005, T λ2,base =-0.003.

[0096] The rules for tiered early warning are: λ max <T λ1 A level one warning is triggered at this time, i.e., λ max It is in the stage of collapsing from a positive value to a zero value, and the sampling interval of the observation vector is shortened from 1 second to 0.5 seconds; λ max <T λ2 A level 2 warning is triggered at this time, i.e., λ maxThe value has turned negative, indicating that the tracking loop has converged to a deterministic state, the deception lock-on is basically complete, and the positioning error is about to exceed the safety envelope. Warning lead time Δt w The second-level warning trigger time t T2 The positioning error first exceeds the safe envelope time t e The time interval, as determined by calibration experiments, is typically 5 to 15 seconds.

[0097] Calculate the confidence level c of the dynamic layer dyn c dyn From λ max The distance to the two-level threshold is mapped via a piecewise function: When λ max >T λ1 hour:

[0098] Where, λ max,norm λ represents the normal tracking state in the current scenario. max The nominal upper bound is determined based on the statistical characteristics of calibration data in normal scenarios; in this embodiment, λ is used. max,norm It is 0.035. When λ max More than λ max,norm The value calculated using the above formula will be lower than 0.1, so the max operation truncates it to 0.1 to avoid negative confidence scores. This piecewise function ensures that when λ... max As c approaches the first threshold from the normal high-value range, dyn The value smoothly increases from 0.1 to 0.5, providing progressive early warning information even within the safe range.

[0099] When T λ2 ≤λ max ≤T λ1 hour:

[0100] When λ max <T λ2 hour:

[0101] The above three piecewise functions together constitute c dyn Complete mapping: λ max When it is in the high range of normal values, c dyn It increases linearly as it decreases; λ max After falling below the first-level threshold and entering the abnormal range, c dyn It grows linearly from 0.5 to 1.0; λ max When the value is below the second-level threshold, c dyn A value of 1.0 indicates that the deception has been highly confirmed.

[0102] It should be noted that in step 4, the dynamic complexity index can also be replaced by permutation entropy PE instead of λ. max With an embedding dimension of 5 and a delay of 4 seconds, the two-level thresholds are set at 60% and 40% of the median value of the sliding window PE under normal conditions, respectively, resulting in lower computational complexity and suitability for resource-constrained scenarios. Alternatively, a recursive quantitative analysis of the deterministic parameter DET can be used, with a recursive threshold set at 15% of the maximum diameter of the phase space and a minimum diagonal length of 5 points. The two-level thresholds are DET > 75% and DET > 90%, offering fast response and serving as a first-level coarse detection indicator.

[0103] Step 5: Multi-feature fusion and hierarchical decision-making.

[0104] This step fuses the heterogeneous features and confidence scores of the three observation domains output from steps 2 to 4. The environmental state variable E drives the adaptive adjustment of the fusion weights, outputting a unified hierarchical comprehensive decision. The specific operation steps are as follows.

[0105] Feature aggregation. Physical layer upload of I and C. phy Information layer uploads B, T, c info and E; λ uploaded from the dynamic layer max and c dyn The confidence levels for each layer have been uniformly mapped to the [0,1] interval.

[0106] Environmental adaptive weight adjustment. The fusion weights are dynamically adjusted by E, and the three domains satisfy w. phy +w info +w dyn =1, where w phy w info w dyn These are the weights for the physical layer, information layer, and dynamics layer, respectively.

[0107] When E=E mp At that time, the physical layer has the highest weight, and the information layer has the lowest weight, such as w phy =0.55、w info =0.15、w dyn =0.30. At this point, the physical layer weight is the highest because the heavier the multipath, the stronger the witnessing effect of mutual information, the higher the mutual information between carrier-to-noise ratio fluctuation and motion, and the greater the discrimination margin for spoofing signals; the information layer weight is the lowest because pseudorange residuals are severely contaminated in multipath-dense environments, and the reliability of B and T features decreases.

[0108] When E=E sp At that time, the information layer has the highest weight, and the physical layer has the lowest weight, such as w phy =0.25、w info =0.40、w dyn=0.35. At this point, the information layer weight is the highest because the direct evidence of deceptive contention locking comes from the bimodal distribution characteristics of the pseudorange residual; the dynamics layer weight is the second highest because contention locking is a precursor to the completion of deceptive traction, and the dynamic complexity is collapsing; the physical layer weight is the lowest because the deceptive signal has partially intervened, and the mutual information may have begun to decrease.

[0109] When E=E norm At that time, the weights of the three domains are basically balanced, such as w phy =0.33, w info =0.34, w dyn =0.33. At this point, the weights of the three domains are basically balanced, serving as backups for each other. Any anomaly in any domain will be reflected in the fusion result.

[0110] The weights switch with E through a linear smooth transition, with 5 transition frames. During the transition, the weights are linearly interpolated from the current value to the target value to avoid instantaneous changes in the overall confidence C caused by a step change in the weights.

[0111] The overall confidence level C is calculated using the following formula:

[0112] The value of C ranges from [0,1]. A tiered decision is made based on the C value: if C < the first preset threshold (e.g., 0.3), it is determined that there is no deception, and the system operates normally; if the first preset threshold ≤ C < the second preset threshold (e.g., 0.5), it is determined to be suspicious, the monitoring frequency is increased (e.g., increased to twice the original value), the decision period is shortened (e.g., shortened to 0.5 seconds), and the update rate of each layer of features is simultaneously increased; if the second preset threshold ≤ C < the third preset threshold (e.g., 0.7), it is determined to be highly suspicious, a first-level alarm is issued, and the navigation system is notified to prepare to switch to the backup navigation source; if C ≥ the third preset threshold, it is determined to be confirmed deception, a second-level alarm is issued, the system immediately switches to the backup navigation source, and a safe docking or return-to-home procedure is triggered. The shortening of the sampling interval and the increase in the update frequency must be implemented under the condition that the receiver supports the corresponding output frequency.

[0113] Alternatively, in addition to the decision-level weighted fusion described above, a feature-level fusion method can also be used: the feature vector F=[I,B,T,λ] is... max The data can be directly input into a lightweight classifier for fusion. A linear SVM classifier or a shallow fully connected network can be used; training data can be collected and labeled in a calibrated scenario, or online weakly supervised training can be performed using the independent decision results of each layer.

[0114] When a coordinated suppression and spoofing attack is detected (manifested as a sudden drop in the real signal C / N0 exceeding a preset threshold, such as 5 dB, accompanied by abnormal spoofing signals), an enhancement strategy is executed. Specifically, the enhancement strategy is as follows: (The code is issued in step 3...) On top of that, the weights are further reduced by 0.03 bits to avoid triggering false alarms due to reduced mutual information after the real signal is suppressed; the weights of the dynamics layer are temporarily increased by 0.1, while the weights of the physical and information layers are scaled proportionally to maintain normalization. This is because the dynamics layer is based on the internal state of the tracking loop and is least likely to be directly interfered with by deceptive signals when dealing with coordinated attacks; based on the remaining available channels, continuous monitoring is performed, and a partial capability degradation status flag is output to the system.

[0115] Figure 4 This is a schematic diagram of a deception interference identification scenario in Example 1. In an urban canyon environment, a power transmission line inspection drone simultaneously receives real satellite signals, multipath reflection signals reflected by building glass curtain walls, and deception signals emitted by ground deception sources. Multipath fading is coupled with the carrier's motion, while the deception signal is not. The invention achieves accurate detection of deception interference through the joint identification of the three observation domains.

[0116] As shown in Figure 5(a), in an urban canyon scenario, a power line inspection drone flies along a preset route between buildings. Without deception, the GNSS positioning trajectory fluctuates slightly around the actual flight trajectory, always remaining within a 40m safety envelope (point-line), and the system does not generate an alarm. After the deception signal is injected from the deception initiation point at t=60s, the deceived GNSS positioning trajectory is gradually pulled away from the actual route along the traction deviation direction. During this process, the present invention sequentially issues a first-level warning (t=80s), a second-level warning (t=87s), and triggers an alarm and switches to backup navigation at t=91s. All of these responses are earlier than the positioning error exceeding the limit point (t=102s), indicating that the present invention can complete deception identification and handling before the power line inspection drone deviates from the safe range.

[0117] As shown in Figure 5(b), before the deception signal injection (t < 60s), the overall confidence level C remains low and is in the normal range (C < 0.3), and multipath disturbances do not cause abnormal accumulation of confidence. After the deception injection at t = 60s, the physical layer, information layer, and dynamic layer features respond successively, and C continues to rise monotonically: at t = 80s, entering the suspicious zone triggers a first-level warning; at t = 87s, entering the highly suspicious zone triggers a second-level warning; and at t = 91s, crossing the alarm threshold of 0.7, entering the confirmed deception zone and triggering an alarm. In contrast, if the positioning error exceeds the limit (t = 102s) as the detection criterion, the response is significantly delayed. This invention gains approximately 15s of warning advance time from the second-level warning to the error exceeding the limit, providing sufficient margin for switching to backup navigation or performing safe maneuvers, demonstrating the timeliness advantage of multi-feature fusion decision over single signal quality detection and fixed-weight fusion methods.

[0118] Example 2: This example uses an automated inspection robot for substation equipment in a densely packed metal structure area. The substation inspection robot operates in this area, surrounded by densely packed metal structures and equipment (metal walls with a reflectivity greater than 0.9). Its operating speed is 1 to 5 km / h, and its path is fixed and repetitive. The GNSS receiver uses a dual-antenna configuration, with both antennas mounted laterally along the roof of the vehicle, with a baseline length of 0.5 m. Auxiliary sensors include a wheel speedometer and an IMU. Figure 6(a) shows a schematic diagram of the substation inspection area layout; Figure 6(b) shows a schematic diagram of the substation inspection robot's navigation and multipath environment elevation. The main difference between this example and Example 1 lies in step 2, which will be explained in detail below.

[0119] In this embodiment, the sliding window length W1 is adaptively scaled according to the substation inspection robot's running speed V:

[0120]

[0121] Among them, W1 is adaptively scaled inversely proportional to the carrier speed, and a lower limit value is set to ensure consistency with Example 1; The differential interval is dynamically adjusted according to W1(V) / 4 to maintain the constraint that the differential interval is greater than the fading coherence time.

[0122] This embodiment employs a dual-antenna differential excitation method for active excitation. Because the substation inspection robot operates at speeds below 2 m / s for most of the time, active micro-motion excitation might interfere with its operational stability; therefore, a dual-antenna differential excitation scheme is adopted. The specific implementation of dual-antenna differential excitation is as follows: The slight changes in the robot's attitude during movement and turning cause relative motion between the two antennas. When the robot turns, the change in the yaw angle causes the two antennas to displace in opposite directions on the horizontal plane. The difference in the multipath environment experienced by the two antennas creates a detectable differential carrier-to-noise ratio fluctuation. The relative displacement increment of the two antennas is calculated from the IMU attitude angle. :

[0123] Where, d base Baseline length This represents the change in the yaw angle. V(t) is used as a motion parameter. After performing carrier-to-noise ratio (CNR) detrending preprocessing and calculating the CNR change gradient for both antenna channels, the differential gradient is taken. To suppress interference components shared by both antennas:

[0124] in, , These are the carrier-to-noise ratio gradients calculated for antenna channel 1 and antenna channel 2, respectively.

[0125] Then, the mutual information I is estimated based on the K-nearest neighbor method. ; The subsequent threshold comparison, confidence mapping, and steps 3 to 5 are the same as in Example 1.

[0126] The reference threshold T in this embodiment I,base A fixed-route multipath fingerprint threshold is adopted. Using the position calculated from wheel speedometer and IMU track as an index (divided into 2m grids), the differential mutual information I at each time a no-spoofing decision (C < first preset threshold) is made is stored in a baseline database according to the position grid. The mean value μ of I within each grid is calculated. I and standard deviation σ I Reference threshold T I,base Change to position-dependent threshold T I,base (P):

[0127] Instead of the fixed reference threshold of 0.15 bits in Example 1, the adaptive threshold becomes:

[0128] To prevent deceptive data from polluting the benchmark database, a pre-verification condition is added before data is entered into the database: the difference between the IMU track-estimated position and the GNSS positioning result is less than a preset threshold (e.g., 0.5m), and the information layer environmental state variable E≠E. sp If GNSS has been spoofed, its positioning results will deviate significantly from the IMU's estimated location. This pre-processing can effectively intercept contaminated data from being stored in the database; simultaneously, E=E sp The status clearly indicates the presence of deceptive competition locking; during this period, characteristic data will not be entered into the database, preventing deceptive signal characteristics from polluting the benchmark database at the source. For grids not covered by the benchmark database (such as first-time access), a fixed benchmark threshold will be used for rollback. This adjustment method changes the decision benchmark from a uniform calibration value across the entire field to a location-specific measured benchmark, eliminating threshold mismatch caused by multipath intensity differences in different areas of the substation.

[0129] The operating configuration of this embodiment is as follows: in the dense metal structure area of ​​the substation, T is consistently in a high range (2.5 to 4.0), B is usually in a low range (less than 0.5), and according to the rule of step 3 in embodiment 1, E is determined to be stable. mp Threshold correction amount = -0.05 bits, =-0.003, fusion weight w phy =0.55、w info =0.15、w dyn=0.30; After the deceptive signal was injected, B exceeded B. th At that time, switch to Arbitration Rule E according to Arbitration Rule E. sp , =+0.05 bits, =+0.003, fusion weight w phy =0.25、w info =0.40、w dyn =0.35, switching is constrained by the hysteresis anti-shake mechanism.

[0130] As shown in Figure 7(a), in the high multipath section of the substation's dense metal structure area (approximately 40s to 70s, during which the substation inspection robot traverses the dense metal structure area), the differential mutual information I increases significantly due to multipath reflection. If a fixed reference threshold (0.15 bits) is used, multiple false alarms will occur in this section (marked with × in the figure); however, the fingerprint threshold established by the present invention based on 2m grid statistics increases adaptively with the grid, and the mutual information in the entire high multipath section never exceeds the fingerprint threshold, thus preventing false alarms. After approximately 90s of spoofing signal injection, the mutual information rapidly increases and exceeds the fingerprint threshold, and the physical layer outputs a valid abnormal response.

[0131] As shown in Figure 7(b), the environmental state E = E mp During this period, the overall confidence level C remained low, and the C level remained low in the strong multipath segment, without triggering the alarm threshold of 0.7, indicating that the present invention has low false alarm characteristics in a strong multipath environment. After the spoofing signal was injected, the bimodal salience B of the information layer exceeded the threshold at approximately 94 seconds, triggering the environmental state to change from E. mp Switch to E sp The fusion weights are adjusted accordingly, the overall confidence level C rises rapidly, and crosses the alarm threshold at approximately 102.7 seconds, triggering a spoofing alarm.

[0132] like Figure 8 As shown, under non-deception conditions, the GNSS positioning trajectory fluctuates slightly within the safety envelope (5m) around the true trajectory, and the system does not issue false alarms. After the deception signal is injected, the deceived positioning trajectory gradually deviates from the true trajectory and is pulled towards the dense metal structure area. This invention issues an alarm at the alarm trigger point, at which time the trajectory deviation is less than 1m, earlier than the position where the deceived trajectory exceeds the 5m safety envelope. This provides time for switching to wheel speedometer / IMU integrated navigation or safe stopping, indicating that this invention can provide timely warning before the deception causes substantial harm.

[0133] A satellite navigation spoofing and interference identification method and system based on multi-feature fusion is provided, which uses any of the above-mentioned satellite navigation spoofing and interference identification methods based on multi-feature fusion.

[0134] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended embodiments and their equivalents.

Claims

1. A satellite navigation spoofing and interference identification method based on multi-feature fusion, characterized in that, include: Estimate the mutual information between the carrier-to-noise ratio change gradient of each channel and the carrier motion parameters to determine the physical layer confidence level; The distribution pattern of pseudorange residuals is modeled online, bimodal significance and heavy-tailed significance are extracted, the information layer confidence is determined, and environmental state variables are generated based on the over-limit states of bimodal significance and heavy-tailed significance of each channel; the phase space of the multidimensional error observation vector of the tracking loop is reconstructed, the dynamic complexity index is calculated, and the dynamic layer confidence is determined and a graded warning is triggered based on the relationship between the dynamic complexity index and the warning threshold; The warning threshold and the fusion weights of each layer are adjusted according to the environmental state quantity. The confidence of each layer is weighted and fused according to the adjusted fusion weights. A hierarchical decision is made based on the fusion result to identify satellite navigation deception interference.

2. The satellite navigation spoofing and interference identification method based on multi-feature fusion according to claim 1, characterized in that: The default selection of the carrier motion parameter is the velocity amplitude. When the variance of the velocity amplitude within the sliding window is lower than the preset variance threshold, the carrier motion parameter is switched from the velocity amplitude to the displacement increment. The estimation of mutual information includes: establishing a mapping model between carrier noise ratio and satellite elevation angle; subtracting the slowly varying component caused by elevation angle change from the carrier noise ratio time series to obtain the rapidly varying component; obtaining the carrier noise ratio change gradient by sliding window difference on the rapidly varying component; and estimating the mutual information between the carrier noise ratio change gradient and the carrier motion parameters based on the K-nearest neighbor method.

3. The satellite navigation spoofing and interference identification method based on multi-feature fusion according to claim 2, characterized in that: When the speed amplitude within the sliding window is lower than the preset speed threshold and remains so for a preset duration, active micro-motion excitation is initiated.

4. The satellite navigation spoofing and interference identification method based on multi-feature fusion according to claim 1, characterized in that: The online modeling involves fitting a bi-Gaussian mixture model to the pseudorange residuals using an online expectation-maximization algorithm with a forgetting factor; the bimodal significance is the ratio of the difference between the mean of the Gaussian components to the mean standard deviation, and the heavy-tailed significance is the ratio of the proportion of samples whose residual tails exceed a preset threshold to the probability of the tails in the standard normal distribution.

5. The satellite navigation spoofing and interference identification method based on multi-feature fusion according to claim 1, characterized in that: The generation of the environmental state quantity includes: each channel performs single-channel judgment based on the excessive states of bimodal significance and heavy-tailed significance; combining the judgment results of all channels, a strong multipath environment is generated when the number of channels with strong multipath judgment exceeds a preset proportion threshold; a deception competition lock environment is generated when the number of channels with deception competition lock judgment exceeds a preset proportion threshold; otherwise, a normal environment is generated; the switching of the environmental state quantity is processed by hysteresis anti-jitter; after combining the judgment results of all channels, the candidate state of the environmental state quantity must be continuously maintained at a preset number of sampling points before the switching is confirmed; and a hysteresis interval is set when exiting.

6. The satellite navigation spoofing and interference identification method based on multi-feature fusion according to claim 1, characterized in that: The dynamic complexity index is the maximum Lyapunov exponent: after removing the mean and normalizing the multidimensional error observation vector of the tracking loop, the phase space is reconstructed using the delayed coordinate embedding method. In the reconstructed phase space, the divergence distance of each neighboring point pair within a preset number of evolution steps is tracked. The logarithm of the divergence distance is taken and the average value is calculated. Then, a linear fit is performed with the number of evolution steps, and the resulting slope is the maximum Lyapunov exponent. The warning threshold includes a zero-reaching threshold and a negative-turning threshold. When the maximum Lyapunov exponent is lower than the zero-reaching threshold, a first-level warning is triggered, and when it is lower than the negative-turning threshold, a second-level warning is triggered.

7. The satellite navigation spoofing and interference identification method based on multi-feature fusion according to claim 1, characterized in that: The adjustment of the environmental state quantity to the warning threshold is as follows: under strong multipath environment, it is adjusted to the baseline threshold minus the preset correction amount; under deceptive contention locking environment, it is adjusted to the baseline threshold plus the preset correction amount; under normal environment, the baseline threshold is taken.

8. The satellite navigation spoofing and interference identification method based on multi-feature fusion according to claim 1, characterized in that: The fusion weights of each layer are highest in the physical layer under strong multipath environment, highest in the information layer under deceptive competition locking environment, and equal in the fusion weights of each layer under normal environment. The hierarchical decision-making includes: when the fusion result is in the first preset range, it is determined to be suspicious, the monitoring frequency is increased to a preset multiple and the decision period is shortened to a preset duration; when the fusion result is in the second preset range, it is determined to be highly suspicious, a first-level alarm is issued and preparations are made to switch to the backup navigation source; when the fusion result is in the third preset range, it is determined to confirm deception, a second-level alarm is issued, the backup navigation source is switched and the security handling procedure is triggered.

9. The satellite navigation spoofing and interference identification method based on multi-feature fusion according to claim 1, characterized in that: When the receiver is configured with dual antennas, the carrier motion parameter is the relative displacement increment of the two antennas due to the change in carrier attitude, and the carrier-to-noise ratio change gradient is the difference between the carrier-to-noise ratio change gradients of the two antennas; the reference threshold of mutual information is the position-related threshold indexed by the position estimated from the carrier track, and the position-related threshold is determined based on the mutual information statistics of each grid in the fixed route when the historical decision is that there is no deception.

10. A satellite navigation spoofing and interference identification system based on multi-feature fusion, characterized in that, include: The physical layer feature extraction module is used to estimate the mutual information between the carrier-to-noise ratio change gradient of each channel and the carrier motion parameters, and to determine the physical layer confidence level. The information layer feature extraction module is used to model the distribution pattern of pseudorange residuals online, extract bimodal significance and heavy-tailed significance, determine the information layer confidence, and generate environmental state quantities based on the out-of-limit state of the bimodal significance and the heavy-tailed significance. The dynamic layer feature extraction module is used to reconstruct the phase space of the multidimensional error observation vector of the tracking loop, calculate the dynamic complexity index, and determine the dynamic layer confidence and trigger graded warnings based on the relationship between the dynamic complexity index and the warning threshold. The fusion decision module is used to adjust the early warning threshold and the fusion weight of each layer according to the environmental state quantity, weight the confidence of each layer according to the adjusted fusion weight, and perform hierarchical decision-making based on the fusion result to identify satellite navigation deception interference.

Images