CN114675307B - Satellite navigation receiver interference detection method based on fractional order pseudo Wigner-Ville distribution - Google Patents
Satellite navigation receiver interference detection method based on fractional order pseudo Wigner-Ville distribution Download PDFInfo
- Publication number
- CN114675307B CN114675307B CN202210248971.2A CN202210248971A CN114675307B CN 114675307 B CN114675307 B CN 114675307B CN 202210248971 A CN202210248971 A CN 202210248971A CN 114675307 B CN114675307 B CN 114675307B
- Authority
- CN
- China
- Prior art keywords
- interference
- signal
- formula
- gnss
- frequency
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000001514 detection method Methods 0.000 title claims abstract description 41
- 238000009826 distribution Methods 0.000 title claims abstract description 33
- 238000000034 method Methods 0.000 claims description 30
- 230000009466 transformation Effects 0.000 claims description 11
- 239000011159 matrix material Substances 0.000 claims description 8
- 238000012545 processing Methods 0.000 claims description 8
- 238000005070 sampling Methods 0.000 claims description 6
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 4
- 238000004364 calculation method Methods 0.000 claims description 4
- 239000000654 additive Substances 0.000 claims description 3
- 230000000996 additive effect Effects 0.000 claims description 3
- 230000004044 response Effects 0.000 claims description 3
- 230000002776 aggregation Effects 0.000 abstract description 5
- 238000004220 aggregation Methods 0.000 abstract description 5
- 238000002474 experimental method Methods 0.000 abstract description 3
- 238000004458 analytical method Methods 0.000 description 14
- 230000007423 decrease Effects 0.000 description 6
- 238000006243 chemical reaction Methods 0.000 description 3
- 230000021615 conjugation Effects 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 230000000116 mitigating effect Effects 0.000 description 3
- 230000035945 sensitivity Effects 0.000 description 3
- 238000001228 spectrum Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 230000002452 interceptive effect Effects 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 238000003491 array Methods 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015556 catabolic process Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000006735 deficit Effects 0.000 description 1
- 238000006731 degradation reaction Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/13—Receivers
- G01S19/23—Testing, monitoring, correcting or calibrating of receiver elements
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/13—Receivers
- G01S19/21—Interference related issues ; Issues related to cross-correlation, spoofing or other methods of denial of service
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Mathematical Physics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Computer Networks & Wireless Communication (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
- Monitoring And Testing Of Transmission In General (AREA)
Abstract
The invention discloses a satellite navigation receiver interference detection method based on fractional order pseudo Wigner-Ville distribution, which is used for interference detection of a global navigation satellite system receiver. The interference detection method is tested through experiments, the tested signal is the BDS-B1I signal of the interfered Beidou satellite navigation system, and the experimental results show that the global navigation satellite system interference detection method provided by the invention effectively eliminates the cross items existing in bilinear time-frequency distribution, simultaneously remarkably improves the time-frequency energy distribution aggregation characteristic and greatly improves the interference detection performance of the global navigation satellite system receiver.
Description
Technical Field
The invention relates to the field of global navigation satellite systems, in particular to a satellite navigation receiver interference detection method based on fractional order pseudo Wigner-Ville distribution (Fractional Pseudo Wigner-Ville Distribution, frPWVD).
Background
Safety issues with global navigation satellite system receivers (e.g., civil aviation and aircraft landing) have attracted considerable attention. Electromagnetic interference in the received signal from the global navigation satellite system (Global Navigation SATELLITE SYSTEM, GNSS) receiver may lead to reduced navigation and positioning performance of the satellite navigation receiver and even to failure to operate properly. Spread spectrum (Direct Sequence Spread Spectrum, DSSS) technology is used in most satellite navigation systems, the principle of which is to spread the received GNSS signal power over a wider bandwidth, ensuring despreading gain in the GNSS receiver, thus reducing the impairments caused by the unwanted interference signals. The DSSS technology makes the satellite navigation system have a certain anti-interference capability, but because the navigation signal power received by the GNSS receiver is very low, even very weak electromagnetic interference signals will cause serious degradation of navigation and positioning performance of the GNSS receiver.
Currently, GNSS interference detection and mitigation techniques (INTERFERENCE DETECTION AND MITIGATION, ID & M) are a very important component in GNSS applications. Different time-frequency representations (Time Frequency Reference, TFR) have been used in GNSS interference detection, such as spectrograms and wegener distribution (Wigner-Ville Distribution, WVD). The spectrogram method has the trade-off problem of Time Frequency (TF) analysis resolution according to the Hessenberg inaccuracy principle, presents poor TF positioning characteristics, and cannot be used for instantaneous Frequency estimation of interference signals. In order to solve the problem of compromise of time-frequency resolution, a WVD method can be used in interference detection of a GNSS receiver. While WVD has many good performances and can provide almost optimal resolution in all time-frequency distributions, it can present severe cross-interference terms in the time-frequency plane due to the interaction of different frequency components, which can lead to severe errors in the GNSS receiver interfering signal instantaneous frequency estimate.
In order to reduce the cross-interference term present in WVD, a reasonable solution is to introduce a window function in the time domain, and therefore, the concept of Pseudo-Wigner-Ville distribution (Pseudo Wigner-Ville Distribution, PWVD) is proposed. The introduction of the window function in PWVD can suppress part of the cross interference term to a certain extent; the disadvantage of this filter window is that the time-frequency analysis energy distribution aggregation characteristic is attenuated, and therefore, the time-frequency positioning accuracy is reduced, and the cross interference term in PWVD can be observed to oscillate parallel to the sweep direction.
In the field of signal processing, fourier transform is widely used as a mature mathematical tool, which is a linear operator for transferring a time signal from a time axis rotation pi/2 to a frequency domain axis, and fractional fourier transform (Fractional Fourier Transform, frFT) is an operator capable of rotating by any angle, which retains the original properties of fourier transform and has new technical advantages, and is considered as a generalized fourier transform. The FrFT can fully embody the transformation characteristic of the signal from the time domain to the frequency domain through continuous transformation of fractional order from 0 to 1. The FrFT can be used as an effective time-frequency analysis tool and can be used for signal estimation instantaneous frequency estimation and phase information recovery.
With the modern progress and rapid development of the global satellite navigation system, the requirements of GNSS interference detection and mitigation technology are higher and higher, and therefore, a more efficient satellite navigation receiver interference detection method based on fractional order pseudo Wigner-Ville distribution needs to be constructed.
Disclosure of Invention
In order to overcome the defects of the prior art, the invention provides a satellite navigation receiver interference detection method based on fractional order pseudo Wigner-Ville distribution, which adopts the combination of FrFT and PWVD to analyze BDS-B1I signals containing sweep frequency interference and extract the characteristics of the sweep frequency interference in GNSS received signals; the FrPWVD method can effectively eliminate the influence of the cross terms existing in the bilinear time-frequency WVD on the GNSS interference detection performance, well reserve the self term components of the GNSS interference signals, and remarkably improve the time-frequency aggregation characteristic of the GNSS interference signals, so that the GNSS receiver interference detection performance is effectively improved.
For the convenience of describing the present invention, firstly, the GNSS received signal will be described, and the sampling frequency of the satellite navigation received signal processed each time is 100MHz, and the digital intermediate frequency is 14.58MHz. The interference signal added in the GNSS effective signal is linear frequency modulation interference. The technical scheme adopted by the invention is as follows:
The satellite navigation receiver interference detection method based on fractional order pseudo Wigner-Ville distribution is characterized by comprising the following steps:
Step 1: and (3) receiving radio frequency signals: receiving a BDS-B1I radio frequency signal r RF (t) through a receiver antenna;
Step 2: an interference signal eta RF (t) is generated by adopting an interference machine and is loaded to a BDS-B1I signal to obtain a signal to be analyzed, and the expression is as follows:
yRF(t)＝rRF(t)+ηRF(t) (1)
Step 3: transmitting the radio frequency signals to a signal collector by using a radio frequency cable for processing, and transmitting the processed signals to a GNSS receiver by using a USB cable for further processing;
Step 4: the GNSS receiver reads a data stream of length n, which can be expressed by the following equation:
y＝[y1,y2,y3,…,yn]T (2)
Step 5: hilbert transform of a data stream y read by a GNSS receiver The received real signal is converted into an analytic signal, and the formula is as follows:
In the formula (3), j is imaginary root units.
Step 6: the FrFT of the analytic signal y a (t) is calculated in the GNSS receiver with the fractional order p=0, and the formula is as follows:
In formula (4), F p (·) represents the FrFT operation, K p (u, t) represents the integral kernel of the FrFT, and the integral kernel function is not only a function of (u, t) but also related to the fractional order p. Wherein K p (u, t) can be expressed as follows:
In the formula (5) of the present invention,
In the formula (5), α=ppi/2, which represents the rotation angle of the time-frequency plane;
In the formula (5), when α=2npi or α= (2n±1) n, K p (u, t) appears as an impulse response;
Obtaining 1 Xn matrix data F p(u)1*n through FrFT;
step 7, adopting a gaussian window function in PWVD, and processing the signal F p(u)1*n after fractional Fourier transformation by PWVD to obtain FrPWVD, wherein the formula is as follows:
In the formula (6) of the present invention, The FrFT conjugation result of the analytic signal y a (t) under the condition that the order is p is shown by the expression conjugate operation;
in the formula (6), h (t) is a time domain window function, and satisfies a generalized time-bandwidth product criterion, and the expression is as follows:
In the formula (7) of the present invention, Representing a gaussian window function of the signal under a minimum Time-bandwidth product (Time-Bandwidth Product, TBP) criterion;
In equation (7), the optimal window function can be described as a window function Is-p-th order FrFT;
In the formula (7) of the present invention, And/>Representing the time width and bandwidth, respectively;
Finally, adopting FrPWVD method to obtain correspondent n×n matrix transformation domain data FrPWVD (u, v) n*n;
Step 8: detecting peak points of matrix transformation domain (u, v) data FrPWVD (u, v) n*n by using a least square method, and fitting the peak points to obtain a connecting line l;
Step 9: calculating an included angle beta between a peak point connecting line l and a coordinate axis of a fractional order transformation domain u, if the deviation E of the included angle beta and a 90-degree angle is larger than sigma, returning to the step 6, increasing the fractional order p of the FrFT by delta p, and repeating the steps 6 to 9; the calculation formula of the deviation E is as follows:
E＝|β-90|/90 (8)
Step 10: if the deviation between the included angle beta and the 90 DEG angle is smaller than epsilon, the result of FrPWVD (u, v) n*n is considered to be an accurate result obtained under the optimal fractional order p, and the intersection point coordinate m of the peak point connecting line l and the fractional transformation domain u axis is determined at the moment;
step 11: and calculating the tuning frequency and the initial frequency of the sweep frequency interference through the optimal fractional order p and the intersection point coordinate m.
GNSS receiver signals are acquired by a receiver single antenna.
The jammer generated GNSS jammer signal η RF (t) may be represented by:
ηRF(t)＝jRF(t)+n(t) (9)
In equation (9), two parts are included, namely a non-stationary radio frequency interference term j RF (t) and a zero-mean additive white gaussian noise term n (t).
Further, the GNSS jammer in step 2 generates an unstable chirp disturbance, which specifically includes the following steps:
The GNSS jammer generates chirp interference with the expression:
In the formula (10), A is the amplitude of the chirp interference, f 0 is the initial frequency of the chirp interference, k is the sweep rate of the chirp interference, Is the initial phase of the disturbance.
In step 3, the length of the data stream processed by the GNSS receiver is calculated, and the formula applied is as follows:
n＝fs×T (11)
In equation (11), f s is the sampling rate of the satellite receiver signal, and T is the scanning period of the interference signal.
The FrFT results of the signals analyzed in step 6 and step 7 are one-dimensional matrix data, and FrPWVD results FrPWVD (u, v) are symmetric arrays.
The beneficial effects of the invention are as follows:
The satellite navigation receiver interference detection method based on fractional order pseudo Wigner-Ville distribution effectively eliminates the influence of cross terms existing in the traditional bilinear time-frequency distribution WVD, reserves the self term energy of GNSS interference signals, and obviously improves the time-frequency energy aggregation characteristic of the GNSS interference signals; in a complex electromagnetic interference environment, the GNSS receiver can accurately identify the GNSS interference type and accurately estimate the corresponding characteristic parameters.
Drawings
Fig. 1 is a schematic block diagram of a satellite navigation interference detection method based on fractional order pseudo-Wigner-Ville distribution according to an embodiment of the present invention.
FIG. 2 is a graph of the results of detection of BDS-B1I signals containing swept interference by an example of the invention at a carrier-to-noise ratio of 46dB-Hz and a dry-to-noise ratio of-2 dB.
FIG. 3 is a graph of the results of detection of BDS-B1I signals containing swept interference by an example of the invention at a carrier-to-noise ratio of 46dB-Hz and a dry-to-noise ratio of-8 dB.
FIG. 4 is a graph showing the analysis of root mean square error obtained by initial frequency estimation of the frequency sweep interference present in the BDS-B1I signal at different interference-to-noise ratios.
FIG. 5 is a graph showing the analysis of root mean square error obtained by performing frequency modulation slope estimation on the frequency sweep interference present in the BDS-B1I signal at different interference-to-noise ratios.
FIG. 6 is a graph showing the analysis of the peak-to-average ratio obtained by interference detection of the swept interference present in the BDS-B1I signal at different interference-to-noise ratios.
Detailed Description
The present invention will be described in further detail with reference to the following examples and the accompanying drawings.
Example 1:
The schematic block diagram of the satellite navigation receiver interference detection method based on fractional order pseudo Wigner-Ville distribution in the embodiment of the invention is shown in figure 1.
Step 1, a single antenna is used for collecting BDS-B1I signals r RF (t), the sampling rate f s is 100MHz, and the digital intermediate frequency f IF of the down-conversion link is 14.58MH.
Step 2: an interference signal eta RF (t) generated by an interference machine is loaded to the BDS-B1I signal to obtain a signal to be analyzed, and the expression is as follows:
yRF(t)＝rRF(t)+ηRF(t) (1)
The interference η RF (t) generated by the software jammer can be expressed by:
ηRF(t)＝jRF(t)+n(t) (2)
In equation (2), two parts are included, namely, a non-stationary interference signal j RF (t) and zero-mean additive white noise n (t).
Step 2.1: the GNSS jammer generates a linear frequency modulation jammer signal with the following expression:
In the formula (3), A is the amplitude of the chirp interference, f 0 is the initial frequency of the chirp interference, k is the sweep rate of the chirp interference, Is the initial phase of the chirp disturbance. The initial frequency f 0 of the constant amplitude chirp interference in this experiment is 13.80MHz, and the sweep rate k is-1.5 x 10 6 MHz/s.
Step 3: and transmitting the GNSS signals to a signal collector by using a radio frequency cable, and transmitting the signals after AD conversion to a GNSS receiver by using a USB cable.
Step 3.1: the length of the data stream processed by the GNSS receiver is calculated, and the following formula is applied:
n＝fs×T (4)
In the formula (4), f s is the sampling rate of the satellite receiver signal, T is the scanning period of the interference signal, and in this experiment, the sampling rate f s of the GNSS receiver signal is 100MHz.
Step 4: the GNSS receiver reads a data stream with length n=2000, and the read data is expressed as:
y(t)＝[y1,y2,y3,…,y2020]T (5)
Step 5: hilbert transform of a data stream y read by a GNSS receiver The received real signal is converted into an analytic signal y a (t), the formula is as follows:
In the formula (6), y a (t) is an analytic signal converted by Hilbert transform, y (t) is a GNSS receiving signal, Representing the Hilbert transform of y (t).
Step 6: the FrFT of y a (t) is calculated with fractional order p=0 and its formula is as follows:
in formula (7), F p (·) represents the FrFT operation, K p (u, t) represents the integral kernel of the FrFT, and the integral kernel function is not only a function of (u, t) but also related to the fraction order p. Wherein K p (u, t) is represented as follows:
In the formula (8), the expression "a",
In the formula (8), α=ppi/2 represents the rotation angle of the time-frequency plane.
In the formula (8), when α=2npi or α= (2n±1) pi, K p (u, t) exhibits an impulse response.
The 1×n matrix-type data F p(u)1*2000 is obtained by the FrFT transform.
Step 7, further processing the signal F p(u)1*2000 of the FrFT with PWVD using a gaussian window function in PWVD, the obtained FrPWVD is represented as follows:
In the formula (9) of the present invention, Is the result of the conjugation of the FrFT of y a (t) when the fractional order is p, and represents the conjugation operation.
H (t) in equation (9) is a window function in the time domain that satisfies the generalized time-bandwidth product criterion, and is expressed as follows:
In the formula (10) of the present invention, Representing the gaussian window function of the signal under the minimum time bandwidth product criterion.
In equation (10), the optimal window function h (t) can be described asIs of the-p-order FrFT.
In the formula (10) of the present invention,And/>Representing time width and bandwidth, respectively.
The corresponding n×n matrix transform domain data FrPWVD (u, v) 2000*2000 is obtained using FrPWVD.
Step 8: and fitting the peak points of matrix-type conversion domain (u, v) data FrPWVD (u, v) 2000*2000 by using a least square method to obtain a fitting straight line l of the peak points.
Step 9: calculating an included angle beta between the peak point fitting straight line l and the coordinate axis of the transformation domain u, if the deviation E between the included angle beta and the 90-degree angle is larger than epsilon, returning to the step 6 to increase the FrFT order p by delta p, and repeating the steps 6 to 9 by delta p=0.01. The calculation formula of the deviation E is as follows:
E＝|β-90|/90 (11)
Step 10: if the deviation between the included angle β and the 90 ° angle is smaller than ε, then FrPWVD (u, v) 2000*2000 is considered to be an accurate result under the condition of the optimal fractional order p, and the intersection point coordinate m of the peak point connecting line l and the transformation domain u axis will be determined.
Step 11: and the frequency modulation slope and the initial frequency of the sweep frequency interference can be calculated through the optimal fractional order p and the intersection point coordinate m.
FIG. 2 (a) shows the result of detection of BDS-B1I signals containing chirp interference by the present embodiment at a dry-to-noise ratio of-2 dB. Since the spectrum of the BDS-B1I useful signal is almost flat, frPWVD analysis makes the variability of the interference signal energy distribution easily visible and distinguishable in the TF plane. This can be further verified in fig. 2 (b), showing the profile of the TF distribution of the present embodiment. For chirp interference, frPWVD presents a tangent plane perpendicular to the u coordinate axis on the transform domain plane with the optimal fractional order, while in other regions the value remains almost zero. The FrPWVD method of the embodiment of the invention effectively suppresses cross-term interference existing in bilinear TF distribution and improves TF energy aggregation characteristics in interference signals. Therefore, the method can be effectively applied to GNSS receiver interference detection in an interference environment, can accurately identify the type of electromagnetic interference and can accurately estimate the characteristic parameters of interference.
FIG. 3 (a) shows the result of detecting BDS-B1I signals containing chirp interference according to this embodiment under the condition that the dry noise ratio is-8 dB. Similar to the influence of the chirp signal interference of fig. 2, under the condition of taking the optimal fractional order, frPWVD of the sweep frequency interference signal presents a tangential plane perpendicular to the u coordinate axis on the transform domain plane, compared with fig. 2, in other areas excluding the sweep frequency interference component, the energy parallel to the interference component is more obvious, and the influence on the GNSS interference detection is very small because the components are far smaller than the contribution of the interference signal, so that the FrPWVD method can effectively detect and identify the GNSS interference signal and accurately estimate the characteristic parameters thereof under the condition of interference.
To fully evaluate the performance of the proposed satellite navigation receiver interference detection method based on fractional order pseudo-Wigner-Ville distribution, a statistical root mean square error (Root Mean Square Error, RMSE) is used for verification analysis. When the BDS-B1I receiving signal has sweep frequency interference, the initial frequency of the sweep frequency interference and the root mean square error of the frequency modulation slope estimation are different under the conditions of different interference to noise ratios (Jamming Noise Ratio, JNR).
Fig. 4 is a quantitative comparison result of initial frequency estimation of frequency sweep interference by different time-frequency analysis methods according to this embodiment, and an RMSE result of initial frequency of chirp interference in BDS-B1I received signals is calculated under different JNR conditions. The RMSE value calculated by adopting the TF distribution method is a function of JNR, and the RMSE of the interference initial frequency decreases with the increase of the JNR value. In fig. 4, the WVD and PWVD methods show very poor frequency estimation accuracy. Wherein, the RMSE of the WVD to the interference initial frequency estimation is almost unchanged with the increase of JNR value, and is maintained at 1×10 -2; RMSE PWVD for the initial frequency estimation of the interfering signal decreases from 1.0×10 -2 to 6.7×10 -3 with increasing JNR value. The RMSE result of the FrFT interferer initial frequency estimate decreases very slowly from 6.9×10 -3 to 2.3×10 -3 with increasing JNR. Along with the change of the JNR value of the frequency sweep interference signal, the RMSE of the frequency sweep interference initial frequency estimation of FrWVD is kept unchanged at the level of 1.9 multiplied by 10 -3, the RMSE of the frequency sweep interference initial frequency estimation of FrPWVD is kept unchanged at the level of 1.9 multiplied by 10 -3, the results of the two interference initial frequency estimation RMSE are the same, and the values of the two interference initial frequency estimation RMSE are smaller than those of the RMSE of other traditional analysis methods, which shows that the method has better stability and effectiveness in the aspect of interference detection parameter estimation compared with the traditional time-frequency analysis method.
Fig. 5 is a comparison analysis result of FrPWVD and the conventional time-frequency analysis method according to the present embodiment for interference frequency modulation slope estimation. Under different JNR conditions, different methods are adopted to estimate the frequency modulation slope of the linear frequency modulation interference in the BDS-B1I receiving signal, and corresponding RMSE is compared and analyzed. The RMSE values of the interference frequency modulation slope estimates calculated by different TF distributions appear as a function of JNR, and it is seen that the RMSE of the frequency sweep interference frequency modulation slope estimates decreases with increasing JNR values as a whole. Wherein, the RMSE of the WVD to the interference frequency modulation slope estimation is almost unchanged with the increase of JNR value, and is maintained at 2.1×10 1; RMSE of PWVD to the jammer chirp rate estimate decreases very slowly with increasing JNR value from 3.1×10 1 to 1.2×10 -1. The RMSE result of FrFT interferer chirp rate estimation slowly decreases from 1x 10 -3 to 3.8 x 10 -4 as JNR increases. The RMSE of FrWVD frequency-sweep interference frequency modulation slope estimation is kept unchanged at the level of 3×10 -4, the RMSE of FrPWVD frequency-sweep interference frequency modulation slope estimation is kept unchanged at the level of 3×10 -4, the results of the two frequency modulation slope estimation are the same, compared with the WVD and PWVD methods, the established FrPWVD method has obvious improvement on the accuracy of the frequency-sweep interference frequency modulation slope estimation, and compared with the FrFT method, the method still has excellent interference parameter estimation performance even under the condition of lower JNR, and the interference detection performance of a GNSS receiver is greatly improved.
To more fully compare the performance of the FrFT, frWVD and FrPWVD methods proposed by the present invention, peak-to-Average Ratio (PAR) will be used to characterize GNSS interference detection sensitivity. PAR is the peak-to-average ratio of the output signal and affects the dynamic range of the power amplifier. The calculation formula of PAR is as follows:
In the formula (12), max { |s (t) | 2 } represents an energy peak value of the signal s (t), and E { |s (t) | 2 } represents an average value of signal energy.
Fig. 6 shows the peak-to-average ratio comparison of FrFT and FrWVD and the proposed FrPWVD method under different JNR conditions. It can be seen that the PAR of all three approaches increases with JNR. The analysis results show that the detection capability of the three methods is improved along with the increase of the JNR value. The PAR results of FrFT range approximately between 1.5 x 10 2 and3 x 10 2; the PAR result of FrWVD appears between 3×10 1 and 9×10 1, which has low interference detection sensitivity. The PAR results of FrPWVD method vary from about 7 x 10 2 to 1.3 x 10 3, which provides a near 470% improvement in interference detection sensitivity compared to the FrFT method. Therefore, compared with the FrFT and FrWVD methods, the FrPWVD method provided by the invention has extremely high interference detection performance, and is expected to solve the interference detection problem of the satellite navigation receiver in the complex electromagnetic environment.
Claims (5)
1. The satellite navigation receiver interference detection method based on fractional order pseudo Wigner-Ville distribution is characterized by comprising the following steps:
step 1: and (3) receiving radio frequency signals: collecting BDS-B1I signals r RF (t) through an antenna;
Step 2: the GNSS interference machine is adopted to generate an interference signal eta RF (t) to be loaded to the BDS-B1I signal, and a signal to be analyzed is obtained, wherein the expression is as follows:
yRF(t)＝rRF(t)+ηRF(t) (1)
step 3: transmitting satellite navigation received signals to a GNSS signal collector for processing by using a radio frequency cable connected with an antenna of the GNSS receiver, and transmitting the collected BDS-B1I signals to the GNSS receiver for processing by using a USB cable;
step 4: the GNSS receiver reads a data stream of length n, the read data being represented by the following formula:
y＝[y1,y2,y3,…,yn]T (2)
Step 5: hilbert transform of a data stream y read by a GNSS receiver Thereby converting the received real signal into an analytic signal, the formula is as follows:
in the formula (3), j is imaginary root units;
Step 6: a fractional fourier transform (Fractional Fourier Transform, frFT) of the analytic signal y a (t) is calculated in the GNSS receiver taking a fractional order p=0, the formula of which is as follows:
In formula (4), F p (·) represents the FrFT operation, K p (u, t) represents the integral kernel of the FrFT, and the integral kernel function is not only a function of (u, t) but also related to the fraction order p; wherein K p (u, t) can be expressed as follows:
In the formula (5) of the present invention,
In the formula (5), α=ppi/2, which represents the rotation angle of the time-frequency plane;
In the formula (5), when α=2npi or α= (2n±1) pi, K p (u, t) appears as an impulse response;
After FrFT, 1 Xn matrix data F p(u)1*n are obtained;
Step 7, using a gaussian window function in a Pseudo-Wigner-Ville distribution (Pseudo-Wigner-Ville Distribution, PWVD), processing the FrFT signal F p(u)1*n with PWVD as follows:
In the formula (6) of the present invention, Is the FrFT of y a (t) under the condition that the fractional order is p, and represents conjugate operation;
in equation (6), h (t) represents a window function in the time domain that satisfies the generalized time-bandwidth product criterion, which can be expressed as follows:
In the formula (7) of the present invention, Representing a gaussian window function of the signal under a minimum Time-bandwidth product (Time-Bandwidth Product, TBP) criterion;
In equation (7), the optimal window function h (t) is expressed as a window function Is-p-th order FrFT;
In the formula (7) of the present invention, And/>Representing the time width and bandwidth, respectively;
Obtaining corresponding n×n-type transform domain matrix data FrPWVD (u, v) n*n by using FrPWVD method;
Step 8: fitting the peak points of the transform domain (u, v) matrix data FrPWVD (u, v) n*n by using a least square method, thereby obtaining a fitting line l connecting the peak points;
Step 9: calculating an included angle beta between a peak point connecting line l and a coordinate axis of a transformation domain u, if the deviation E of the included angle beta and a 90-degree angle is larger than epsilon, returning to the step 6, increasing the fractional order p of the FrFT by delta p, and repeating the steps 6 to 9; the calculation formula of the deviation E is as follows:
E＝|β-90|/90 (8)
step 10: if the deviation of the included angle beta and the 90-degree angle is smaller than epsilon, the result of FrPWVD (u, v) n*n is considered to be an accurate result obtained under the optimal fractional order p, and at the moment, the intersection point coordinate m of the peak point connecting line and the transformation domain u axis is determined;
step 11: and estimating the tuning frequency and the initial frequency of the GNSS sweep frequency interference signal through the optimal fractional order p and the intersection point coordinate m.
2. The satellite navigation receiver interference detection method based on fractional order pseudo-Wigner-Ville distribution according to claim 1, wherein the signals received in the GNSS interference detection process are acquired with a single antenna.
3. The satellite navigation receiver interference detection method based on fractional order pseudo-Wigner-Ville distribution according to claim 1, wherein the GNSS interference signal η RF (t) generated by the jammer can be represented by the following formula:
ηRF(t)＝jRF(t)+n(t) (9)
In equation (9), two parts are included, namely, non-stationary GNSS interference j RF (t) and zero-mean additive white gaussian noise n (t).
4. The method for detecting satellite navigation receiver interference based on fractional order pseudo-Wigner-Ville distribution according to claim 3, wherein the unstable interference j RF (t) generated by the middle scrambler in step 2 is chirped interference, and the specific process is as follows:
The GNSS jammer generates a linear frequency modulation jammer signal with the following expression:
In the formula (10), A is the amplitude of the chirp interference, f 0 is the initial frequency of the chirp interference, k is the sweep rate of the chirp interference, Representing the initial phase of the chirp disturbance.
5. The satellite navigation receiver interference detection method based on fractional order pseudo-Wigner-Ville distribution according to claim 1, wherein the data stream length processed by the GNSS receiver is calculated in step 3, and the following formula is applied:
n＝fs×T (11)
In the formula (11), f s is the sampling rate of the satellite receiver signal, and T is the scanning period of the swept interference signal.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210248971.2A CN114675307B (en) | 2022-03-14 | 2022-03-14 | Satellite navigation receiver interference detection method based on fractional order pseudo Wigner-Ville distribution |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210248971.2A CN114675307B (en) | 2022-03-14 | 2022-03-14 | Satellite navigation receiver interference detection method based on fractional order pseudo Wigner-Ville distribution |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114675307A CN114675307A (en) | 2022-06-28 |
| CN114675307B true CN114675307B (en) | 2024-04-30 |
Family
ID=82074094
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210248971.2A Active CN114675307B (en) | 2022-03-14 | 2022-03-14 | Satellite navigation receiver interference detection method based on fractional order pseudo Wigner-Ville distribution |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114675307B (en) |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9596610B1 (en) * | 2016-07-28 | 2017-03-14 | Universitat Politécnica de Catalunya | System and method for detecting and eliminating radio frequency interferences in real time |
| CN107728166A (en) * | 2017-09-08 | 2018-02-23 | 哈尔滨工程大学 | A kind of more disturbance restraining methods of satellite navigation receiver based on time domain grouping processing |
| CN109375242A (en) * | 2018-10-30 | 2019-02-22 | 重庆邮电大学 | BOC signal capture algorithm under high dynamic based on PMF-FRFT |
| CN110515096A (en) * | 2019-08-20 | 2019-11-29 | 东南大学 | Satellite navigation interference signal identification device and its method based on convolutional neural networks |
| US11194052B1 (en) * | 2018-08-30 | 2021-12-07 | Rockwell Collins, Inc. | Systems and methods for real-time GNSS interference monitoring |
-
2022
- 2022-03-14 CN CN202210248971.2A patent/CN114675307B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9596610B1 (en) * | 2016-07-28 | 2017-03-14 | Universitat Politécnica de Catalunya | System and method for detecting and eliminating radio frequency interferences in real time |
| CN107728166A (en) * | 2017-09-08 | 2018-02-23 | 哈尔滨工程大学 | A kind of more disturbance restraining methods of satellite navigation receiver based on time domain grouping processing |
| US11194052B1 (en) * | 2018-08-30 | 2021-12-07 | Rockwell Collins, Inc. | Systems and methods for real-time GNSS interference monitoring |
| CN109375242A (en) * | 2018-10-30 | 2019-02-22 | 重庆邮电大学 | BOC signal capture algorithm under high dynamic based on PMF-FRFT |
| CN110515096A (en) * | 2019-08-20 | 2019-11-29 | 东南大学 | Satellite navigation interference signal identification device and its method based on convolutional neural networks |
Non-Patent Citations (3)
| Title |
|---|
| A New GNSS Interference Detection Method Based on Rearranged Wavelet–Hough Transform;Kewen Sun 等;Sensors;20210302;1-20 * |
| 基于分数阶傅里叶变换的GNSS接收机扫频干扰检测;赵慧子;孙克文;;信息与电脑(理论版);20170108(01);29-36 * |
| 基于分数阶时频变换的卫星导航干扰检测方法;郑逸飞;中国优秀硕士学位论文全文数据库 信息科技辑 (月刊);20230515;12-68 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114675307A (en) | 2022-06-28 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107678003B (en) | Target detection method under ground wave radar sea clutter background | |
| CN103969635B (en) | The meteorologic signal process IP kernel of low-altitude surveillance radar | |
| CN109387813B (en) | High-precision frequency measurement method based on broadband radar signal reception | |
| WO2015172622A1 (en) | Method for radio-frequency interference suppression of high-frequency ground wave radar | |
| CN111610503B (en) | Linear frequency modulation signal parameter estimation method based on improved LVD | |
| CN101825707A (en) | Monopulse angular measurement method based on Keystone transformation and coherent integration | |
| Ellingson et al. | Mitigation of radar interference in L-band radio astronomy | |
| CN112462356B (en) | Method for improving radar ranging accuracy based on Gaussian interpolation | |
| CN110346772A (en) | A kind of higher-frequency radar significantly ionosphere phase diameter Disturbance Rejection method | |
| CN109975771A (en) | Wideband digital channel method based on three rank phase difference of signal | |
| CN110398722A (en) | Extension target echo detection method based on the limited spectrum of random matrix | |
| CN114089384A (en) | Satellite navigation signal synthesis power spectrum deviation calculation method | |
| CN114675307B (en) | Satellite navigation receiver interference detection method based on fractional order pseudo Wigner-Ville distribution | |
| CN110763930A (en) | Arrester resistive current on-line monitoring system based on Blackman bimodal interpolation harmonic analysis method | |
| Abramovich et al. | Adaptive ionospheric distortion correction techniques for HF skywave radar | |
| CN106646422B (en) | Preprocessing system for enhancing signal-to-noise ratio of Doppler frequency shift signal of coherent wind radar | |
| CN102790650B (en) | Fraction domain odd/even alternation channelized receiving method for broadband LFM signal | |
| CN107769815A (en) | Linear frequency modulation short-range detecting system noise AM interference suppressing method | |
| Ashry et al. | Comparative analysis between SAR pulse compression techniques | |
| CN113406386B (en) | Signal frequency accurate estimation method based on digital down-conversion | |
| CN116087956A (en) | Traffic radar speed measurement method and device, radar and readable storage medium | |
| Sun et al. | A novel interference detection method based on Wigner-Hough transform for GNSS receivers | |
| CN112824927B (en) | Sky wave over-the-horizon radar ionospheric phase pollution correction method based on sparse filtering | |
| Bucci et al. | Validation of pulse compression techniques for meteorological functions | |
| CN114488208B (en) | Beidou signal anti-interference method combining empirical wavelet and SPWVD conversion |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |