A GNSS spoofing interference detection method based on joint SQM square

Abstract

The present invention discloses a GNSS spoofing interference detection method based on combined SQM square, which includes the following steps: obtaining the carrier-to-noise ratio, correlator spacing, and coherent integration time of a GNSS receiver, and calculating the coherent integration values of the early path, prompt path, and late path of the tracking loop; calculating three SQM metrics, namely Delta, Ratio, and ELP, according to the coherent integration values; calculating the statistical characteristics of the SQM metrics based on the carrier-to-noise ratio, correlator spacing, and coherent integration time of the receiver; establishing a combined detection quantity; setting a false alarm rate, and determining a decision threshold according to the statistical distribution of the combined detection quantity; setting a time window length, and making a decision according to the decision threshold within the detection time to achieve spoofing interference detection. The method disclosed by the present invention does not need to change the receiver structure, has a wide range of applicable scenarios, can be applied to static and dynamic receivers, and has good detection effects on small-delay spoofing interference and complex offset spoofing interference.

Description

Technical Field

[0001] The present invention relates to the technical field of satellite signal spoofing interference, and particularly relates to a GNSS spoofing interference detection method based on the square of joint SQM (Signal Quality Monitoring). Background Art

[0002] The Global Navigation Satellite System (GNSS) has been widely used in all walks of life, enabling the quality improvement and upgrading of various industries and deeply integrating into people's daily lives and industrial activities. However, while GNSS devices bring convenience to people's lives, they also pose certain potential threats. Due to the openness of the GNSS signal structure and the weakness of the signal power, it is extremely vulnerable to various intentional and unintentional interferences. Among them, spoofing interference requires little power, has strong concealment, and great harm, and has received extensive attention from the military and experts and scholars in various countries.

[0003] In recent years, many universities and research institutions at home and abroad have conducted a lot of research on spoofing interference and the impact of spoofing signals on receivers. The currently commonly used spoofing interference detection methods have the following defects: (1) Most methods change the structure of the GNSS receiver and cannot be widely promoted under the condition of cost savings. For example, the detection method of navigation information encryption makes it difficult for spoofers to obtain and change the navigation information of satellites by adding signal encryption features, thereby increasing information security. However, it requires the national level to change the satellite signal system, and the receiver side also needs to change the signal decryption algorithm; using the spatial correlation characteristics of spoofing interference signals, the spatial processing technology is used to estimate the spatial characteristics of the received signals to achieve spoofing interference detection. However, most of these methods require an increase in the number of receiving antennas, and even using antennas with special performance requires an increase in the cost of the entire detection system. (2) Some methods do not need to change the structure of the receiver, but the detection applicable scenarios are relatively limited, and the detection effect for some complex and precise spoofing interferences is poor. For example, the detection method based on the time of arrival of the signal is not suitable for dynamic scenarios; the detection method based on Doppler consistency is mostly applied to moving receivers; the detection method based on the number of correlation peaks in the acquisition stage cannot cope with small-delay spoofing interference; the detection method based on the carrier-to-noise ratio (C / N0) has a poor detection effect for gradually offset spoofing interference. Summary of the Invention

[0004] To solve the above technical problems, the present invention provides a GNSS spoofing interference detection method based on the square of joint SQM. This method does not need to change the receiver structure, provides the possibility for commercial receivers to achieve spoofing interference detection under the condition of cost savings, and has a wide range of applicable scenarios. It can be applied to static and dynamic receivers and has a good detection effect for small-delay spoofing interference and complex offset spoofing interference.

[0005] To achieve the above object, the technical solution of the present invention is as follows:

[0006] A GNSS spoofing interference detection method based on combined SQM square, comprising the following steps:

[0007] Step 1, obtain the carrier-to-noise ratio, correlator spacing, and coherent integration time of the GNSS receiver, and calculate the coherent integration values of the early path, prompt path, and late path of the tracking loop;

[0008] Step 2, calculate three SQM metrics, namely Delta, Ratio, and ELP, based on the coherent integration values of the early path, prompt path, and late path;

[0009] Step 3, calculate the statistical characteristics of the SQM metrics according to the carrier-to-noise ratio, correlator spacing, and coherent integration time of the receiver;

[0010] Step 4, set the metric combination method and establish a combined detection quantity;

[0011] Step 5, set the false alarm rate and determine the decision threshold according to the statistical distribution of the combined detection quantity;

[0012] Step 6, set the time window length and make a decision according to the decision threshold within the detection time to achieve spoofing interference detection.

[0013] In the above solution, the specific content of Step 1 is as follows:

[0014] A. Under the condition of no spoofing interference, the intermediate frequency signal model obtained by the GNSS receiver after receiving the signal through the downconverter is:

[0015]

[0016] where i is the satellite number, m is the total number of satellites, c i is the signal power received from the i-th satellite, τ i is the propagation delay of the signal from the i-th satellite, D i (·) is the data code information modulated by the i-th satellite, C i (·) is the C / A code sequence of the i-th satellite, f IF is the intermediate frequency, f d,i is the Doppler frequency shift of the i-th satellite, θ i is the initial carrier phase of the i-th satellite, n fe (t) is the radio frequency front-end noise;

[0017] B. Assume that within the coherent integration time, the data level of the received signal does not change. D(t - τ i ) is denoted as the sampled value D(n), where n is the sampling time, and the integration time T cohFor the high-frequency components in the mixed-frequency signal to be long enough and the integral value to be approximately 0, the in-phase coherent integral value of the r(t) signal in the instant path after passing through the correlator and integral filter can be expressed as:

[0018]

[0019] Among them,

[0020]

[0021] Among them, I P,i (n) represents the in-phase coherent integral value of the instant path of the i-th satellite, t0 is the integral start time, sinc(·) is the sinc function, T coh is the coherent integration time, is the code phase tracking value of the tracking loop for the i-th satellite, is the Doppler frequency tracking value of the tracking loop for the i-th satellite, is the carrier phase tracking value of the i-th satellite, R(·) is the autocorrelation function of the C / A code, is the correlation noise of the in-phase branch of the instant path, represents the code phase tracking error; represents the frequency tracking error; represents the phase tracking error;

[0022] Similarly, the coherent integral outputs of other branches in the tracking loop are expressed as:

[0023]

[0024]

[0025]

[0026]

[0027]

[0028] Among them, d is the correlator spacing, T c is the C / A code chip width, Q P,i (n) represents the quadrature coherent integral value of the instant path of the i-th satellite, I E,i (n) represents the in-phase coherent integral value of the early path of the i-th satellite, Q E,i (n) represents the quadrature coherent integral value of the early path of the i-th satellite, I L,i (n) represents the in-phase coherent integral value of the late path of the i-th satellite, Q L,i (n) represents the quadrature coherent integral value of the late path of the i-th satellite, is the correlation noise of the quadrature branch of the instant path, is the correlated noise of the in-phase branch of the early path, is the correlated noise of the quadrature branch of the early path, is the correlated noise of the in-phase branch of the late path, is the correlated noise of the quadrature branch of the late path.

[0029] In the above solution, step 2 is specifically as follows:

[0030] Calculate the numerical values of each index according to the definition formula of the SQM index. The specific calculation formula is as follows:

[0031]

[0032]

[0033]

[0034] where m Delta,i (n), m Ratio,i (n), m ELP,i (n) respectively represent the Delta, Ratio, and ELP indexes of the i-th satellite signal at time n.

[0035] In the above solution, step 3 is specifically as follows:

[0036] A. Based on the analysis in step 1, for the sampled value of the i-th satellite signal at time n, assume that the sampled value D(n) of the data level is 1, and the tracking loop is operating in a stable state, and the code phase tracking error and the carrier phase error are both 0; Substitute Δτ i = 0, Δf i = 0, Δθ i = 0, D(n) = 1 into equations (4) to (8), and the unified forms of the in-phase coherent integration value and the quadrature coherent integration value of the tracking loop are written as:

[0037]

[0038]

[0039] Since the correlated noise of the in-phase branch and the correlated noise of the quadrature branch are independent of each other and both follow a Gaussian distribution with a mean of 0 and a variance of , the mean and variance of the tracking loop output value are as follows:

[0040]

[0041] E[Q d = 0 (15)

[0042]

[0043] Among them, N0 represents the double-sided power spectral density of Gaussian white noise, and E[I d represents the mean value of the in-phase coherent integration value of the tracking loop, and E[Q d represents the mean value of the quadrature coherent integration value of the tracking loop. D[I d represents the variance of the in-phase coherent integration value of the tracking loop, and D[Q d represents the variance of the quadrature coherent integration value of the tracking loop. When d < 0, I d , Q d respectively represent the in-phase and quadrature coherent integration values of the early path. When d > 0, I d , Q d respectively represent the in-phase and quadrature coherent integration values of the late path. When d = 0, I d , Q d respectively represent the in-phase and quadrature coherent integration values of the prompt path;

[0044] B. According to the definition, the calculation description of the SQM index is a process of first calculating the ratio of the coherent integration values of different branches and then performing a linear combination. Then the general form of the SQM index is expressed as:

[0045]

[0046] Among them,

[0047]

[0048] a j represents the proportionality coefficient, J represents the number of r j proportional expressions, x j , y j respectively represent the coherent integration output values of two different branches of the tracking loop. Assuming that x j and y j are independent of each other, using the Taylor formula expansion and statistical theory analysis, the statistical characteristics of r j are calculated as:

[0049]

[0050]

[0051] Among them, μ x , σ x 2 are respectively the mean value and variance of the variable x j , and μ y , σ y 2 are respectively the mean value and variance of the variable y jThe mean and variance of μ r , σ r 2 are the mean and variance of the variable r j respectively;

[0052] Transform each SQM index into the general form shown in Equation (17), substitute Equations (14) - (16) into Equations (19) - (20), and calculate the theoretical mean and variance of different SQM indices; under the condition of high signal-to-noise ratio, the three SQM indices of Delta, Ratio, and ELP are determined to follow a Gaussian distribution, and the mean and variance are as follows:

[0053] μ Delta = 0 (21)

[0054]

[0055]

[0056]

[0057] μ ELP = 0 (25)

[0058]

[0059] Among them, μ Delta represents the mean of the Delta index, σ Delta 2 represents the variance of the Delta index, μ Ratio represents the mean of the Ratio index, σ Ratio 2 represents the variance of the Ratio index, μ ELP represents the mean of the ELP index, σ ELP 2 represents the variance of the ELP index.

[0060] In the above solution, step 4 is specifically as follows:

[0061] Construct a joint detection quantity M cmb using the statistical property that the SQM index follows a Gaussian distribution, and its specific form is as follows:

[0062]

[0063] Among them, M cmb,i (n) represents the joint detection quantity at the nth moment of the signal of the i-th satellite, m Delta,i (n), m Ratio,i (n), m ELP,i(n) represents the Delta, Ratio, and ELP indicators of the i-th satellite signal at time n, μ Delta represents the mean of the Delta indicator, σ Delta represents the standard deviation of the Delta indicator, μ Ratio represents the mean of the Ratio indicator, σ Ratio represents the standard deviation of the Ratio indicator, μ ELP represents the mean of the ELP indicator, σ ELP represents the standard deviation of the ELP indicator; β1 represents the coefficient of the Delta indicator, β2 represents the coefficient of the Ratio indicator, β3 represents the coefficient of the ELP indicator, β i The value-taking rule is as follows:

[0064]

[0065] Obtain all M cmb,i constituting the sample population M cmb Conforms to the following statistical characteristics:

[0066] M cmb ~χ 2 (k) (29)

[0067] where k = β1 + β2 + β3, representing the number of SQM indicators participating in the joint detection, χ 2 (k) represents the chi-square distribution with k degrees of freedom.

[0068] In the above solution, the specific step 5 is as follows:

[0069] M cmb The probability density function of is:

[0070]

[0071] In the formula, Γ(k / 2) is the Gamma function;

[0072] M cmb The cumulative distribution function of is:

[0073]

[0074] In the formula, is the incomplete Gamma function;

[0075] Thus, the false alarm rate is calculated as:

[0076]

[0077] In the formula, th is the decision threshold, P fa represents the probability that the indicator exceeds the decision threshold in the absence of spoofing interference attacks, that is, the false alarm rate;

[0078] According to the inverse function existence theorem, since the cumulative distribution function F k (·) is strictly monotonically increasing, its inverse function must exist. After rearranging Equation (32), the decision threshold is obtained as follows:

[0079] th = F k -1 (1 - P fa ) (33)

[0080] where F k -1 (·) is the inverse function of F k (·). Given the false alarm rate P fa , there is a unique corresponding decision threshold th; the values of the quantiles corresponding to the χ 2 distribution for different degrees of freedom and right-tail probabilities have been tabulated; the degree of freedom is determined by the index combination method, and the false alarm rate is set according to requirements. The corresponding decision threshold can be obtained by looking up the χ 2 distribution critical value table.

[0081] In the above solution, step 6 is specifically as follows:

[0082] Since the tracking loop outputs a coherent integration value every T coh , the time window length is set to NT coh . Then, there are N combined SQM samples within one detection time interval. P d,i is defined as the ratio of the number of samples where the combined detection quantity M cmb exceeds the threshold to the total number of samples when the i-th satellite is under spoofing interference attack;

[0083] The spoofing interference detection probability P d,i is calculated as follows:

[0084]

[0085] where M cmb,i (n) represents the combined detection quantity of the i-th satellite signal at time n, th is the decision threshold, and I(·) is the indicator function. When M cmb,i (n) > th, it takes the value 1; otherwise, it takes the value 0;

[0086] Based on the value of the spoofing interference detection probability P d,i , it is determined whether the i-th satellite is under spoofing interference.

[0087] Through the above technical solution, a GNSS spoofing interference detection method based on combined SQM square provided by the present invention has the following beneficial effects:

[0088] 1. The present invention calculates using the coherent integration output values of the early path, prompt path, and late path in the receiver tracking loop, without changing the structure of the receiver and without introducing additional devices, providing the possibility for commercial navigation receivers to achieve spoofing interference detection under the condition of saving hardware costs.

[0089] 2. The present invention solves the problem that the detection effect of the joint SQM algorithm based on amplitude combination proposed by scholars on spoofing interference has no improvement compared with the single SQM index. The joint SQM algorithm ignores the situation that different indexes have positive deviation from the upper threshold and negative deviation from the lower threshold, resulting in abnormal cancellation after combination. The method provided by the present invention draws on the idea of least squares, uses the square operation to make all deviations positive and then adds the amplitudes, solving the above abnormal cancellation problem, and greatly improving the detection performance without increasing the algorithm complexity.

[0090] 3. The present invention determines whether the receiver is under spoofing interference attack by identifying the distortion of the correlation peak in the tracking loop, which is applicable to both fixed and moving receivers, and has a good detection effect on small-delay spoofing interference and complex and precise bias spoofing interference. BRIEF DESCRIPTION OF THE DRAWINGS

[0091] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required in the description of the embodiments or the prior art.

[0092] Figure 1 It is a schematic flowchart of a GNSS spoofing interference detection method based on joint SQM square disclosed in the embodiments of the present invention.

[0093] Figure 2 It is a visualization diagram of the Delta index result provided by the embodiments of the present invention;

[0094] Figure 3 It is a visualization diagram of the Ratio index result provided by the embodiments of the present invention;

[0095] Figure 4 It is a visualization diagram of the ELP index result provided by the embodiments of the present invention;

[0096] Figure 5 It is the detection quantity M under the joint mode of selecting three SQM indexes of Delta, Ratio, and ELP provided by the embodiments of the present invention cmb Visualization diagram of the result.

[0097] Figure 6 It is a simulation diagram of the spoofing interference detection result under different joint modes of selecting SQM indexes provided by the embodiments of the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0098] The technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention.

[0099] The present invention provides a GNSS spoofing interference detection method based on combined SQM square, as Figure 1 shown, including the following steps:

[0100] In this embodiment, the second spoofing interference data in the TEXBAT dataset of the University of Texas at Austin in the United States is selected. The data sampling rate is 25 MHz, the C / A code rate is 1.023 MHz, the spoofing signal power is constantly 10 dB higher than the real signal power, the total length of the selected spoofing signal data is 400 s, and the spoofing attack start time is the 100th second. The GNSS software receiver is used to process the data, and the 23rd satellite with the highest signal power in the capture result is selected, that is, i = 23 in the following formula, and the spoofing interference detection is carried out on the tracking result of the 23rd satellite signal.

[0101] Step 1, obtain the carrier-to-noise ratio, correlator spacing, and coherent integration time of the GNSS receiver, and calculate the coherent integration values of the early path, prompt path, and late path of the tracking loop.

[0102] A. Under the condition of no spoofing interference, the intermediate frequency signal model obtained by the GNSS receiver receiving the signal after passing through the downconverter is:

[0103]

[0104] where i is the satellite number, m is the total number of satellites, c i is the signal power received from the i-th satellite, τ i is the propagation delay of the i-th satellite signal, D i (·) is the data code information modulated by the i-th satellite, C i (·) is the C / A code sequence of the i-th satellite, f IF is the intermediate frequency, f d,i is the Doppler frequency shift of the i-th satellite, θ i is the initial carrier phase of the i-th satellite, n fe (t) is the radio frequency front-end noise;

[0105] B. Assume that within the coherent integration time, the data level of the received signal does not change, and D(t - τ i ) is denoted as the sampling value D(n), n is the sampling time, and the integration time T coh is long enough for the high-frequency components in the mixed signal, and the integration value is approximately 0. Then the in-phase coherent integration value of the r(t) signal in the prompt path after passing through the correlator and integration filter can be expressed as:

[0106]

[0107] Among them,

[0108]

[0109] Among them, I P,i (n) represents the in-phase coherent integration value of the prompt path of the i-th satellite, t0 is the integration start time, sinc(·) is the sinc function, and T coh is the coherent integration time, is the code phase tracking value of the tracking loop for the i-th satellite, is the Doppler frequency tracking value of the tracking loop for the i-th satellite, is the carrier phase tracking value of the i-th satellite, R(·) is the autocorrelation function of the C / A code, is the correlation noise of the in-phase branch of the prompt path, represents the code phase tracking error; represents the frequency tracking error; represents the phase tracking error;

[0110] Similarly, the coherent integration outputs of other branches in the tracking loop are expressed as:

[0111]

[0112]

[0113]

[0114]

[0115]

[0116] Among them, d is the correlator spacing, and T c is the C / A code chip width, and Q P,i (n) represents the quadrature coherent integration value of the prompt path of the i-th satellite, and I E,i (n) represents the in-phase coherent integration value of the early path of the i-th satellite, and Q E,i (n) represents the quadrature coherent integration value of the early path of the i-th satellite, and I L,i (n) represents the in-phase coherent integration value of the late path of the i-th satellite, and Q L,i (n) represents the quadrature coherent integration value of the late path of the i-th satellite, is the correlation noise of the quadrature branch of the prompt path, is the correlation noise of the in-phase branch of the early path, is the correlation noise of the quadrature branch of the early path, is the correlation noise of the in-phase branch of the late path, Is the correlated noise of the quadrature branch of the lag path.

[0117] In this embodiment, the carrier-to-noise ratio of the received signal by the receiver is 53 dB·Hz, the correlator spacing is 0.5 chips, the coherent integration duration is 1 ms, and the coherent integration values are stored in a tracking matrix with dimensions of 6×400000.

[0118] Step 2: Calculate three SQM metrics, namely Delta, Ratio, and ELP, based on the coherent integration values of the early path, prompt path, and lag path.

[0119] Calculate the values of each metric according to the definition formulas of the SQM metrics. The specific calculation formulas are as follows:

[0120]

[0121]

[0122]

[0123] where m Delta,i (n), m Ratio,i (n), m ELP,i (n) represent the Delta, Ratio, and ELP metrics of the i-th satellite signal at time n, respectively.

[0124] The visualization results of m Delta,i (n), m Ratio,i (n), m ELP,i (n) calculated in this embodiment are respectively as Figure 2 , Figure 3 , Figure 4 shown.

[0125] Step 3: Calculate the statistical characteristics of the SQM metrics based on the carrier-to-noise ratio, correlator spacing, and coherent integration time of the receiver.

[0126] A. Based on the analysis in Step 1, for the sampled value of the i-th satellite signal at time n, assume that the sampled value D(n) of the data level is 1, and the tracking loop is operating in a stable state, and the code phase tracking error and carrier phase error are both 0; substitute Δτ i = 0, Δf i = 0, Δθ i = 0, D(n) = 1 into equations (4) to (8), and the unified forms of the in-phase coherent integration value and quadrature coherent integration value of the tracking loop are written as:

[0127]

[0128]

[0129] Since the correlated noise of the in-phase branch The correlated noise with the quadrature branch are mutually independent and both follow a Gaussian distribution with a mean of 0 and a variance of . Therefore, the mean and variance of the tracking loop output value are as follows:

[0130]

[0131] E[Q d = 0 (15)

[0132]

[0133] where N0 represents the double-sided power spectral density of Gaussian white noise, E[I d represents the mean of the in-phase coherent integration value of the tracking loop, E[Q d represents the mean of the quadrature coherent integration value of the tracking loop, D[I d represents the variance of the in-phase coherent integration value of the tracking loop, D[Q d represents the variance of the quadrature coherent integration value of the tracking loop. When d < 0, I d , Q d respectively represent the in-phase and quadrature coherent integration values of the early path. When d > 0, I d , Q d respectively represent the in-phase and quadrature coherent integration values of the late path. When d = 0, I d , Q d respectively represent the in-phase and quadrature coherent integration values of the prompt path.

[0134] In the GNSS receiver used in this embodiment, the correlator spacing is 0.5 chip. When d = -0.5, I d , Q d respectively represent the in-phase and quadrature coherent integration values of the early path. When d = 0.5, I d , Q d respectively represent the in-phase and quadrature coherent integration values of the late path. When d = 0, I d , Q d respectively represent the in-phase and quadrature coherent integration values of the prompt path. Substituting different values of d into equations (14) - (16), the means and variances of the coherent integration values of different branches in the tracking loop are calculated when the correlator spacing is 0.5 chip, and the results are as follows:

[0135]

[0136] E[I P = E[Q E = E[Q P = E[Q L = 0

[0137]

[0138] B. According to the definition, the calculation of the SQM index is described as a process of first taking the ratio of the coherent integration values of different branches and then performing a linear combination. For example, the Delta index is written as in the form, then the general form of the SQM index is expressed as:

[0139]

[0140] where

[0141]

[0142] a j represents the proportionality coefficient, J represents the number of j proportional expressions, x j , y j respectively represent the coherent integration output values of two different branches of the tracking loop. Assuming that x j and y j are independent of each other, using Taylor series expansion and statistical theory analysis, the statistical characteristics of r j are calculated as:

[0143]

[0144]

[0145] where μ x , σ x 2 are the mean and variance of the variable x j respectively, μ y , σ y 2 are the mean and variance of the variable y j respectively, μ r , σ r 2 are the mean and variance of the variable r j respectively;

[0146] Transform each SQM index into the general form shown in Equation (17), substitute Equations (14) - (16) into Equations (19) - (20), and calculate the theoretical mean and variance of different SQM indices; under the condition of high signal-to-noise ratio, the three SQM indices of Delta, Ratio, and ELP are determined to follow a Gaussian distribution, and the mean and variance are as follows:

[0147] μ Delta = 0 (21)

[0148]

[0149] μ Ratio = 1 (23)

[0150]

[0151] μ ELP = 0 (25)

[0152]

[0153] Among them, μ Delta represents the mean of the Delta index, and σ Delta 2 represents the variance of the Delta index, μ Ratio represents the mean of the Ratio index, and σ Ratio 2 represents the variance of the Ratio index, μ ELP represents the mean of the ELP index, and σ ELP 2 represents the variance of the ELP index.

[0154] Step 4, set the index combination method and establish the combined measurement.

[0155] Construct the combined measurement M cmb using the statistical property that the SQM index follows a Gaussian distribution. Its specific form is as follows:

[0156]

[0157] Among them, M cmb,i (n) represents the combined measurement at the nth moment of the i-th satellite signal, and m Delta,i (n), m Ratio,i (n), m ELP,i (n) respectively represent the Delta, Ratio, and ELP indexes of the i-th satellite signal at the nth moment. μ Delta represents the mean of the Delta index, and σ Delta represents the standard deviation of the Delta index. μ Ratio represents the mean of the Ratio index, and σ Ratio represents the standard deviation of the Ratio index. μ ELP represents the mean of the ELP index, and σ ELP represents the standard deviation of the ELP index; β1 represents the coefficient of the Delta index, β2 represents the coefficient of the Ratio index, β3 represents the coefficient of the ELP index, and β i has the following value-taking rules:

[0158]

[0159] Obtain all Mcmb,i The sample population M composed of cmb complies with the following statistical characteristics:

[0160] M cmb ~χ 2 (k) (29)

[0161] where k = β1 + β2 + β3, representing the number of SQM indicators participating in the joint detection, and χ 2 (k) represents the chi-square distribution with k degrees of freedom. When β1, β2, and β3 all take the value of 1, the joint detection quantity M cmb The visualization result is as Figure 5 shown.

[0162] Step 5, set the false alarm rate and determine the decision threshold according to the statistical distribution of the joint detection quantity.

[0163] Based on the analysis in Step 4, M cmb follows the χ 2 distribution with k degrees of freedom, and its probability density function is:

[0164]

[0165] In the formula, Γ(k / 2) is the Gamma function;

[0166] M cmb The cumulative distribution function of is:

[0167]

[0168] In the formula, is the incomplete Gamma function;

[0169] Thus, the false alarm rate is calculated as:

[0170]

[0171] In the formula, th is the decision threshold, and P fa represents the probability that the indicator exceeds the decision threshold without spoofing interference attack, that is, the false alarm rate;

[0172] According to the inverse function existence theorem, since the cumulative distribution function F k (·) is strictly monotonically increasing, it must have an inverse function. After arranging Equation (32), the decision threshold is obtained as:

[0173] th = F k -1 (1 - P fa ) (33)

[0174] In the formula, F k -1(·) is F k The inverse function of (·), set the false alarm rate P fa , there is a unique decision threshold th corresponding to it; χ 2 For the χ distribution with different degrees of freedom and right-tail probabilities, the values of the corresponding quantiles have been tabulated; the degree of freedom is determined by the index combination method, the false alarm rate is set according to requirements, and the corresponding decision threshold can be obtained by looking up the χ 2 distribution critical value table.

[0175] In this embodiment, the false alarm rate is set to 0.01. When one SQM index participates in the detection, the threshold value is 6.635; when two SQM indexes participate in the detection, the threshold value is 9.210; when three SQM indexes participate in the detection, the threshold value is 11.345.

[0176] Step 6, set the time window length, and make a decision according to the decision threshold within the detection time to achieve spoofing interference detection.

[0177] Since the tracking loop outputs a coherent integration value every T coh , the time window length is set to NT coh , then there are N combined SQM samples within one detection time interval, and P d,i is defined as the ratio of the number of samples whose combined detection quantity M cmb exceeds the threshold to the total number of samples when there is a spoofing interference attack on the i-th satellite;

[0178] The spoofing interference detection probability P d,i The calculation formula is:

[0179]

[0180] Among them, M cmb,i (n) represents the combined detection quantity of the i-th satellite signal at the n-th moment, th is the decision threshold, and I(·) is the indicator function. When M cmb,i (n) > th, the value is 1, otherwise, the value is 0; according to the spoofing interference detection probability P d,i to judge whether the i-th satellite is under spoofing interference. When P d,i is greater than 0.5, it can be considered that it is under spoofing interference.

[0181] In this embodiment, the time window is selected as 100 ms, there are 100 combined SQM samples within one detection time period, and different combination methods are used to construct the combined SQM detection quantity M cmb for threshold decision, and the spoofing interference detection result is obtained as Figure 6As shown in the figure. The horizontal axis in the figure represents the receiving signal time, and the vertical axis represents the spoofing interference detection rate. Since the spoofing signal is initially phase-aligned with the real signal code, no spoofing is detected from 100 s to 150 s. As time changes, the spoofing signal gradually deviates from the real signal, and the spoofing interference detection probability gradually increases after 150 s, reaching the maximum value at 180 s. After 250 s, the correlation peak of the spoofing signal separates from the real signal, and the spoofing interference detection rate drops to 0. It can be seen from the figure that the detection probability of the proposed joint SQM square algorithm of the present invention is higher than that of the single SQM index, and the detection effect is the best when the three SQM indexes of Delta, Ratio, and ELP are jointly detected.

[0182] The above description of the disclosed embodiments enables those skilled in the art to implement or use the present invention. Various modifications to these embodiments will be obvious to those skilled in the art, and the general principles defined herein can be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to the embodiments shown herein, but will be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A GNSS spoofing interference detection method based on combined SQM square, characterized in that, It includes the following steps: Step 1: Obtain the carrier-to-noise ratio, correlator spacing, and coherent integration time of the GNSS receiver, and calculate the coherent integration values of the early path, prompt path, and late path of the tracking loop; Step 2: Calculate three SQM metrics, namely Delta, Ratio, and ELP, based on the coherent integration values of the early path, prompt path, and late path; Step 3: Calculate the statistical characteristics of the SQM metrics according to the carrier-to-noise ratio, correlator spacing, and coherent integration time of the receiver; Step 4: Set the index combination method and establish a combined detection quantity; Step 5: Set the false alarm rate and determine the decision threshold according to the statistical distribution of the combined detection quantity; Step 6: Set the time window length and make a decision according to the decision threshold within the detection time to achieve spoofing interference detection; The specific content of Step 4 is as follows: Construct the joint detection quantity M by using the statistical property that the SQM index follows a Gaussian distribution cmb , and its specific form is as follows: Among them, M cmb,i (n) represents the joint detection quantity of the i-th satellite signal at the n-th moment, m Delta,i (n), m Ratio,i (n), m ELP,i (n) respectively represent the Delta, Ratio, and ELP indicators of the i-th satellite signal at the n-th moment, μ Delta represents the mean value of the Delta indicator, σ Delta represents the standard deviation of the Delta indicator, μ Ratio represents the mean value of the Ratio indicator, σ Ratio represents the standard deviation of the Ratio indicator, μ ELP represents the mean value of the ELP indicator, σ ELP represents the standard deviation of the ELP indicator; β1 represents the coefficient of the Delta indicator, β2 represents the coefficient of the Ratio indicator, β3 represents the coefficient of the ELP indicator, β i The value-taking rule is as follows: Obtain all of M cmb,i The sample population M composed of (n) cmb Complies with the following statistical characteristics: M cmb ~χ 2 (k) (29) where k = β1 + β2 + β3, representing the number of SQM indicators participating in the joint detection, and χ 2 (k) represents the chi-square distribution with k degrees of freedom.

2. The GNSS spoofing interference detection method based on combined SQM square according to claim 1, characterized in that, The specific content of Step 1 is as follows: A. Under the condition of no spoofing interference, the intermediate-frequency signal model obtained by the GNSS receiver after receiving the signal and passing it through a down-converter is: where i is the satellite number, m is the total number of satellites, and c i is the signal power received from the i-th satellite, τ i is the propagation delay of the signal from the i-th satellite, D i (·) is the data code information modulated by the i-th satellite, C i (·) is the C / A code sequence of the i-th satellite, f IF is the intermediate frequency, f d,i is the Doppler shift of the i-th satellite, θ i is the initial carrier phase of the i-th satellite, n fe (t) is the radio frequency front-end noise; B. Assume that within the coherent integration time, the data level of the received signal does not undergo a jump. Let D(t - τ i ) be denoted as the sampled value D(n), where n is the sampling instant, and the integration time T coh is long enough for the high-frequency components in the mixed-frequency signal, and the integration value is approximately 0. Then the in-phase coherent integration value of the r(t) signal in the instant path after passing through the correlator and integration filter can be expressed as: Where, Among them, I P,i (n) represents the in-phase coherent integration value of the prompt path of the i-th satellite, t0 is the integration start time, sinc(·) is the sinc function, T coh is the coherent integration time, is the code phase tracking value of the tracking loop for the i-th satellite, is the Doppler frequency tracking value of the tracking loop for the i-th satellite, is the carrier phase tracking value of the i-th satellite, R(·) is the autocorrelation function of the C / A code, is the correlation noise of the in-phase branch of the prompt path, represents the code phase tracking error; represents the frequency tracking error; represents the phase tracking error; Similarly, the coherent integration output of other branches in the tracking loop is expressed as: where d is the correlator spacing, T c is the C / A code chip width, Q P,i (n) represents the orthogonal coherent integration value of the prompt path of the i-th satellite, I E,i (n) represents the in-phase coherent integration value of the early path of the i-th satellite, Q E,i (n) represents the orthogonal coherent integration value of the early path of the i-th satellite, I L,i (n) represents the in-phase coherent integration value of the late path of the i-th satellite, Q L,i (n) represents the orthogonal coherent integration value of the late path of the i-th satellite, is the correlation noise of the orthogonal branch of the prompt path, is the correlation noise of the in-phase branch of the early path, is the correlation noise of the orthogonal branch of the early path, is the correlation noise of the in-phase branch of the late path, is the correlation noise of the orthogonal branch of the late path.

3. A GNSS spoofing interference detection method based on combined SQM square according to claim 2, characterized in that, The specific content of Step 2 is as follows: Calculate the values of each metric according to the definition formula of the SQM metric. The specific calculation formulas are as follows: where m Delta,i (n), m Ratio,i (n), m ELP,i (n) respectively represent the Delta, Ratio, and ELP indicators of the i-th satellite signal at time n.

4. A GNSS spoofing interference detection method based on combined SQM square according to claim 3, characterized in that The specific content of Step 3 is as follows: A. Based on the analysis in Step 1, for the sampling value of the i-th satellite signal at time n, assuming that the sampling value D(n) of the data level is 1 and the tracking loop is operating in a steady state, with both the code-phase tracking error and the carrier-phase error being 0; substituting Δτ i = 0, Δf i = 0, Δθ i = 0, and D(n) = 1 into equations (4) to (8), the unified forms of the in-phase coherent integral value and the quadrature coherent integral value of the tracking loop are written as follows: Due to the correlated noise of the in-phase branch and the correlated noise of the quadrature branch are independent of each other and both follow a Gaussian distribution with a mean of 0 and a variance of , the mean and variance of the tracking loop output value are as follows: E[Q d = 0 (15) where N0 represents the two-sided power spectral density of Gaussian white noise, E[I d represents the mean value of the in-phase coherent integration value of the tracking loop, E[Q d represents the mean value of the quadrature coherent integration value of the tracking loop, D[I d represents the variance of the in-phase coherent integration value of the tracking loop, D[Q d represents the variance of the quadrature coherent integration value of the tracking loop. When d < 0, I d , Q d respectively represent the in-phase and quadrature coherent integration values of the early path. When d > 0, I d , Q d respectively represent the in-phase and quadrature coherent integration values of the late path. When d = 0, I d , Q d respectively represent the in-phase and quadrature coherent integration values of the prompt path; B. According to the definition, the calculation of the SQM metric is described as a process of first taking the ratio of the coherent integration values of different branches and then performing a linear combination. Then the general form of the SQM metric is expressed as: Where, a j represents the proportionality coefficient, and J represents r j the number of proportional expressions, x j , y j respectively represent the coherent integration output values of two different branches of the tracking loop. Assume that x j and y j are independent of each other. By using the Taylor formula expansion and statistical theory analysis, the statistical characteristics of r j are obtained as follows: Among them, μ x , σ x 2 are the mean and variance of the variable x j respectively, μ y , σ y 2 are the mean and variance of the variable y j respectively, μ r , σ r 2 are the mean and variance of the variable r j respectively; Transform each SQM metric into the general form shown in Equation (17), substitute Equations (14) - (16) into Equations (19) - (20), and calculate the theoretical mean and variance of different SQM metrics; under the condition of high signal-to-noise ratio, the three SQM metrics, Delta, Ratio, and ELP, are considered to follow a Gaussian distribution, and the mean and variance are as follows: μ Delta = 0 (21) μ ELP = 0 (25) Among them, μ Delta represents the mean of the Delta index, and σ Delta 2 represents the variance of the Delta index, μ Ratio represents the mean of the Ratio index, and σ Ratio 2 represents the variance of the Ratio index, μ ELP represents the mean of the ELP index, and σ ELP 2 represents the variance of the ELP index.

5. A GNSS spoofing interference detection method based on combined SQM square according to claim 1, characterized in that, The specific content of Step 5 is as follows: M cmb The probability density function of: In the formula, Γ(k / 2) is the Gamma function; M cmb The cumulative distribution function of In the formula, is the incomplete Gamma function; Thus, calculate the false alarm rate as: where th is the decision threshold, and P fa represents the probability that the index exceeds the decision threshold without spoofing jamming attack, that is, the false alarm rate; According to the inverse function existence theorem, since the cumulative distribution function F k (·) is strictly monotonically increasing, its inverse function must exist. After rearranging Equation (32), the decision threshold is obtained as follows: th = F k -1 (1 - P fa ) (33) where F k -1 (·) is the inverse function of F k (·), and given the false alarm rate P fa , there is a unique decision threshold th corresponding to it; the values of the quantiles corresponding to the χ 2 distribution for different degrees of freedom and right-tail probabilities have been tabulated; the degree of freedom is determined by the index combination method, the false alarm rate is set according to requirements, and the corresponding decision threshold can be obtained by looking up the χ 2 distribution critical value table.

6. The GNSS spoofing interference detection method based on combined SQM square according to claim 1, wherein The specific content of Step 6 is as follows: Since the tracking loop outputs the coherent integration value every T coh time, the time window length is set to NT coh . Then there are N combined SQM samples within one detection time interval, and P d,i is defined as the ratio of the number of samples whose combined detection quantity M cmb exceeds the threshold to the total number of samples when the i-th satellite is under spoofing interference attack; Probability P of deception jamming detection d,i The calculation formula is as follows: Among them, M cmb,i (n) represents the joint detection quantity at the nth moment of the i-th satellite signal, th is the decision threshold, and I(·) is the indicator function. When M cmb,i (n) > th, the value is 1; otherwise, the value is 0. Determine whether the i-th satellite is under spoofing interference according to the spoofing interference detection probability P d,i value.

Images