# Interference Mitigation Based on Low Complexity Subspace Tracking Method for Vehicle Positioning

Authors: Ping Zhang; Qiang Wu; Binbin Jing
Journal: Discrete Dynamics in Nature and Society
 (2022)

---

## Abstract

Advanced driver assistance systems and automatic driving have drawn lots of attentions due to user’s multiapplication needs. Among different technologies, the integration of Inertial Navigation System (INS) and Global Navigation Satellite System (GNSS) is a key technique with the decreasing costs of Inertial Measurement Unit (IMU) and the increasing deployment of GNSS for navigation applications. Although the integration may have immune ability to interference in a certain extent, interference suppression strategy applied on the integration system can greatly improve the performance. A complete framework including interference monitor, detection, and mitigation is provided jointly using subspace tracking method and Kalman filter measurements innovation of INS and GNSS. The proposed monitor is simple, efficient, and can decide when to start the interference mitigation procedure, which can reduce the computation complexity significantly; Subspace-based method is introduced to track the interference basis matrix that can be feed to RootMUSIC to detect the frequency of interference, followed by the interference mitigation with a notch filter. The proposed framework has superior performance than the well-known time-frequency algorithm, such as short-time Fourier transform (STFT) and smoothed pseudo Wigner-Ville distribution (SPWVD), and can improve the performance of INS/GNSS system in the presence of interference with simulation verification.


---

## Body

## 1. Introduction

Advanced driver assistance systems (ADAS) and automatic driving play an important role and evolve gradually with different technologies [1]. Precise and stable positioning and time synchronization are key techniques for connectivity and collaboration of vehicle-to-vehicle, vehicle-to-devices, and vehicle-to-infrastructure. Among different positioning techniques, such as Global Navigation Satellite System (GNSS), Inertial Measurement Unit (IMU), Lidar [2], and cameras [3], GNSS and IMU are the most commonly used since Lidar and cameras perform poorly in low-visibility conditions such as in the dense fog and heavy rain situation. With the decreasing costs of IMU and more GNSS deployment, the integration of Inertial Navigation System (INS) and GNSS has drawn a lot of attentions for navigation applications, such as airlines, unmanned plane, driverless cars, and pedestrian navigation [4]. INS can provide continuous output (usually at least 50 Hz) of velocity, position, and attitude for vehicles with less hardware fault after a proper initial value, as well as the immunity to interference. However, due to its accumulated error with time, the position accuracy of INS starts to deteriorate after navigation and eventually becomes unacceptable, which is common for the low-cost IMU [5]. Therefore, INS needs to be corrected regularly. Compared with INS, with more cheaper device cost and the increasing satellite system deployment, real time kinematic (RTK) based on GNSS can provide position accuracy about centimeter level. However, GNSS may become unstable due to signals outage especially in urban environment. So the combination of INS/GNSS has become a common way to utilize its advantages and provide long- or short-term high accuracy navigation solution. In general, according to the type of measurements used in GNSS, three kinds named loosely coupled, tightly coupled, and ultra-tightly coupled INS/GNSS are used. Here the tightly coupled INS/GNSS approach is adopted as the main architecture since its flexible operation and medium complexity. For this scheme, pseudorange and pseudorange-rate measurements of GNSS are applied to aid INS to provide a corrected navigation solution.The performance of INS/GNSS system can be improved by different methods such as building a proper model for the INS errors and bias terms, adopting a robust Kalman filter [6] or denoising the measurements of IMU sensors [7]. However, these methods are mainly considered from the aspect of INS. As the ever-increasing risk for GNSS infrastructure, methods are less considered from the aspect of GNSS in the presence of interference. As is well known, interference has become one of the major threats to GNSS. Although the integration of INS/GNSS has the ability to suppress the interference or jamming in some extent, as the description above-mentioned, the errors from IMU sensors, such as bias, scale factor, and white noise, can cause an unacceptable output of INS without the GNSS aiding for a long time. So the performance of INS/GNSS system can be improved significantly with the interference detection and mitigation (D&M) technique. In this paper, we focus on the interference D&M technique for GNSS using low complexity method, while how to improve the performance of GNSS signal outages in urban environment is out of scope of this paper.The problem of interference D&M for the single GNSS receiver is widely studied in the literature. Generally, interference can be dealt with from two aspects: hardware configuration and signal processing. From the point of hardware strategies, automatic gain control (AGC) [8] and ADC [9] are commonly used methods to mitigate the interference. The main weakness of hardware settings is the lack of flexibility. From the aspect of signal processing, interference detection and mitigation technology can be applied in time domain [10, 11], time-frequency domain [12, 13], and spatial domain [14, 15]. In this paper, methods for interference detection in time-frequency domain are our main concern. In practice, transforming the received data to the frequency domain is a feasible way to detect the interference. Here we classify these approaches into three types. The first and well-known type of method is named short-time Fourier transform (STFT), based on the Fourier transform (FT) [16], and applied for the analysis of interference spectrum in time-frequency domain. The second type of methods, based on Wigner-Ville distribution (WVD) or Smooth Pseudo WVD (SPWVD), can also offer time and frequency information of interference for GNSS receiver [12]. However, the main drawback of methods based on WVD is its cross-interfering term, which can degrade the performance of frequency analysis seriously. In addition, Wavelet-Packet-Transform method is also used for the chirp-type interference D&M of GNSS receiver [17]. The similar method based on simplified Welch algorithm and notch filtering is also used to reduce the noise [18].In comparison with the above-mentioned methods, a competitive solution based on subspace decomposition (SD) of covariance matrix of sample data, named Karhunen-Loève transform (KLT), is proposed to detect very weak RF signals in noisy environment with the signal-to-noise ratio (SNR) as low as -23 decibels for GNSS receiver [19], which is a level that the art-of-the-state approaches cannot reach. Our proposed algorithm, which is also based on SD, is similar to KLT. However, the framework is mainly used to track the interference space basis matrix, which can be applied to estimate the interference frequencies. In addition, based on subspace tracking method, the computation complexity can be reduced.In this paper, we focus on the low-complexity subspace-based technique to mitigate the interference. Span basis matrix of sample covariance matrix, whether signal subspace or noise subspace, should be estimated by all the subspace-based methods. Conventional batch singular value decomposition (SVD), such as used in KLT, is not fit for online interference detection for INS/GNSS receiver because of its high computation complexity requiringON3 operations (we assume that the observation vector r∈ℝN and the rank of interference subspace is L. In general case, N is much larger than L. In the literature, some low complexity algorithms have been provided with computation as low as ONL operations, Depending on tracking types of the subspace, these algorithms can be categorized into two types [20]. By applying recursive least squares (RLS) method, projection approximation subspace tracking (PAST) [21] and fast approximated power iteration (FAPI) [22] are well-known subspace tracking algorithms. The main drawback is that this kind of algorithms can only track a single subspace. Another kind of methods can track signal or noise subspace simultaneously with a simple sign change using the same code structure. Two popular methods of this kind algorithm are FDPM [23] and the algorithm proposed in [24]. Comparing with the first type, these algorithms have more simple code design. In this paper, we adopt the second kind of algorithms as our detection method due to their simple code structure.A framework for interference monitor, detection, and mitigation is provided based on low complexity subspace tracking method for vehicle positioning. The whole work consists of two phases: (1) interference monitor phase and (2) interference detection and mitigation phase. In the first phase, consistency checks of innovation measurement between GNSS and INS are used to decide when the interference is present; in the second phase, once the interference occurs, an algorithm based on the low-complexity subspace tracking algorithm in time-frequency domain is presented to detect and mitigate the interference. Proposed framework belongs to the pre-correlation class of receiver and need not modify the existing structure of system. The overview of the proposed framework is described, as shown in Figure1.Figure 1
Overview of the proposed framework.The paper is organized as follows. In Section2, signal and system model in the presence of interference is provided. We propose a complete framework including interference monitor, detection, and mitigation in INS/GNSS system in Section 3. In Section 4, simulations and verification for the proposed framework are provided. Finally, conclusion and future work are presented in Section 5.


## 2. Signal and System Model

In this section, the signal and system model is described first in the presence of interference environment. Here we notice that although our proposed scheme is applied on the INS/GNSS system, INS is an autonomous system and do not receive any RF signals. Therefore, the interference detection and mitigation method are mainly used for GNSS signals. The signals received from GNSS front-end can be modeled as(1)rt=∑k=1Mxkt+∑k=1Likt+υt,where ∑k=1Mxkt represent useful signals received from M satellites in view, ∑k=1Likt are interference signals including L components, and υt is a zero-mean Gaussian process caused by the thermal noise. For the kth single signal of satellite, xkt can be represented as(2)xkt=Akdkt−τkckt−τcos2πfkt+ϕk,where Ak, fk, τk are the estimated components of amplitude, time delay, and carrier frequency including the intermediate frequency plus Doppler shift frequency, respectively. The above-mentioned parameters are all time varying and need to be estimated. dkt is the transmitted data, ckt is the spreading code, and ϕk is the carrier phase. For the lth single signal of interference, ilt is usually in the form of a periodic signal and has a general model as(3)ilt=Ailtexpjωltt+φl,where Ailt, ωlt, and φl denote the time-varying amplitude, instantaneous frequency, and phase of interference, respectively. In general, according to the type of jammers, ωlt may have different forms such as continuous wave (CW) interference when ωlt is a constant, or sweep-like interference when ωlt regularly changes over time slowly. In addition, we assume that Ailt is a constant and Ail for simplicity. In this paper, instantaneous estimation of the interference frequency, ωt, is our main concern. Equation (3) can be further written as(4)ilt=Blelt,where Bl=Ailexpφl, elt=expjωltt.Due to the utilization of direct-sequence spread spectrum technique in most GNSS systems, the signal characteristic of the GNSS is different from that of the interference; meanwhile, the power density of received GNSS signals at the input of front-end is far lower than that of thermal noise about 20 dB. So satellites signals can be viewed as wideband Gaussian white noise. Equation (1) can be rewritten as(5)rt=∑k=1Likt+Υt,where Υt is the zero-mean Gaussian white noise including term ∑k=1lxkt plus term υt as the above description.After down-conversion, sampling, and quantization in the front-end of GNSS device, the discrete form ofrt can be written as(6)rn=∑k=1Likn+Υn.Thus, different digital signal processing method can be used to detect and mitigate the interference term∑k=1Likn. After obtaining N data samples, rn, eln, and Υn can be rewritten in the vector form as(7)r=r0r1⋯rN−1T,(8)el=1expjωlexpj2ωl⋯expjN−1ωlT,(9)Υ=Υ0Υ1⋯ΥN−1T.Using (4), (7)–(9),we derive the covariance matrix of r(10)Cr=ErrH=∑l=1LPelelH+σ2IN,where E⋅ represents the expectation operator, Pl=Ail2 is the power of interference signal, and IN is the identity matrix with rank N. Applying eigenvalue decomposition (EVD) to Cr, we have(11)Cr=UiΣiUiT+UΥΣΥUΥT=UiUΥΣi00ΣΥUiUΥT,where ∑i=diagλ1,λ2,…,λL mainly contains the interference eigenvalues, ∑Y=diagλL+1,λL+2,…,λN is the noise eigenvalues, and Ui and UΥ are corresponding interference and noise subspace basis matrix, which are both orthonormal and estimated by the subspace tracking algorithms. Once Ui is estimated, RootMusic method can be used to detect the frequency of interference [25].


## 3. Framework for the Interference Detection and Mitigation Using Data Fusion

In this section, a complete framework for interference detection and mitigation by subspace tracking method is presented in a tightly coupled INS/GNSS system. The innovation of Kalman filter is used to provide the interference monitor. Then subspace tracking method is used to estimate the interference basis matrix, which can be used to feed to the frequency estimation algorithm, such as RootMusic, to estimate the frequencies of interference. Finally, with a simple filter method, notch filter, to mitigation the interference.


### 3.1. Tightly-Coupled INS/GNSS with Kalman Filter

Data fusion based on Kalman filter is a common method for vehicle positioning, and consistent check using the innovation of INS/GNSS system can provide a monitor of interference. By tightly coupled mechanization, the pseudo-range and pseudo-range rate measurements from GNSS equipment are used to aid INS. Meanwhile, INS can also aid GNSS receiver to track the satellites signals more fastly. Extended Kalman filter (EKF) is used for data fusion due to its simpleness and stability, which provides an optimal estimation of the states. A general design frame for the tightly coupled INS/GNSS system is studied, as shown in Figure2.Figure 2
A common design for tightly coupled INS/GNSS.Corrected measurements from the GNSS code tracking loop in earth centered earth fixed (ECEF) frame can be written as(12)ρ˜=ρ˜R−δρ^I−δρ^T+δρ^C,ρ˙˜=ρ˙˜R+δρ˙^C,where ρ˜R, ρ˙˜R denote the pseudo-range and pseudo-range rate from GNSS receiver, respectively. δρ^I, δρ^T are the ionospheric and tropospheric refraction errors, respectively. δρ^C, δρ˙^C, respectively, represent the clock error and clock error shift. An approximate estimate of ρ˜R can be written as(13)ρ˜R=xa−xs2+ya−ys2+za−zs2=ts−ta×c,where xa,ya,za, xs,ys,zs are the receiver and satellite position in ECEF frame, ta−ts represents the signal transmission time from satellite to GNSS receiver, and c is the speed of light. It should be noticed that some superscript or subscript are omitted for simplicity. For example, the signal from lth satellite in ECEF frame may be written as ρ˜Re,l.Since the states of INS and GNSS are independent, the state vector and state transformation matrix are given by(14)x=xINSxGNSS,ϕ=ϕINS00ϕGNSS.The estimated 17 components of state vector are written as(15)x=δΨδυδrbgbaδρδρ˙,where δp, δv, and δA are the position, velocity, and attitude state errors, respectively, bg and ba are the gyroscope and accelerometer biases, respectively, which should be estimated for any applications; δρ is the receiver clock error; and δρ˙ is the receiver clock drift, which are processed and corrected by GNSS. The component of x may have different forms in different coordinate frame. For instant, when in local navigation frame δr can be represented as latitude, longitude, and height error. After corrected, INS output can be written as(16)C^αβγ=δC^αβγC˜αβγ,r^=δr^−r˜,υ^=δυ˜−υ˜,where α, β, γ represent the body, reference, and projection coordinate system, respectively.The predicted statex using inertial navigation equation is expressed as(17)x^k−=ϕkx^k−1+εk,where x^k− denotes the predicted state vector at epoch k, ϕk is the state transformation matrix, x^k−1 is the updated state vector at epoch k − 1, andεk is the process noise. ϕk with first-order approximate in ECEF frame is presented as(18)ϕke=I3−Ωτs030303Cbeτs00F21eτsI3−2ΩτsF23eτsCbeτs030003I3τsI3030300030303I3030001,301,301,301,3I30001,301,301,301,3031τs01,301,301,301,30301,where Ω is the skew symmetric matrix of the earth rotation angle velocity ωie in ECEF frame(19)Ω=0−ωie0ωie00000,τs is the time epoch of IMU processing unit. F21e and F23e are, respectively, defined as(20)F21e=−C^bef^ibb∧,F23e=−2γ^ibereSeL^br^eber^ebe.The measurement model can be written as(21)zk=Hkxk+νk.where zk is the measurements vector of pseudo-range and pseudo-range rate from GNSS, νk is the measurement noise at epoch k, and Hk is the geometry matrix including line-of-sight vectors from vehicle to satellites in the ECEF frame and can be written as(22)δzk=δzρ,kδzρ˙,k,Hk=01,301,3ua1γT01,301,31001,301,3ua2γT01,301,310⋮⋮⋮⋮⋮⋮⋮01,301,3uamγT01,301,310⋯⋯⋯⋯⋯⋯⋯01,3ua1γT01,301,301,30101,3ua2γT01,301,301,301⋮⋮⋮⋮⋮⋮⋮01,3uamγT01,301,301,301,where u is the line of sight, defined as(23)u=∂ρ∂x,∂ρ∂y,∂ρ∂x.


### 3.2. Monitor Based on the Consistent Check of Measurements

Some alarm mechanisms should be provided first when a fault occurred. In this work, we assume that the fault is caused by the interference only for simplicity, and an interference monitor algorithm is proposed to check whether the interference is present or not. In general, consistent check of the measurements is a common method used in receivers. Similar to receiver autonomous integrity monitoring (RAIM) derived from the least squares estimation for GNSS receiver and can be viewed as a hypothesis testing problem [26], we use the innovation of Kalman filter as a monitor in the INS/GNSS system. The main concept can be explained as follows: when there is no interference, the measurements of GNSS and INS should be with higher correlation. While in the presence of interference, due to the GNSS signals degrade seriously, the deviation of INS and GNSS measurement become unacceptable and inconsistent. So by checking sudden changes of the measurement residual compared with a predefined threshold, we can decide whether the interference is present or not. The innovation formula is(24)ηk=zk−Hkx^k−,where xk− is a prior of state vector at time k from the inertial navigation equation defined in (17).From the hypothesis testing theory, a test statisticrk as the quadratic form of the normalized innovation at time k is defined as(25)rk=ηkTCk−1ηk,where Ck is the covariance matrix of the innovation ηk at time k and can be calculated as(26)Ck=EηkηkT=HkPk−HkT+Rk,where Rk is the covariance of measurement noise and Pk− is the predicted covariance of state error vector at time k(27)Pk−=Ex^k−−xkx^k−−xkT.For a properly selected thresholdΓk, interference occurrence or not can be inferred by comparing the value of rk with Γk. Let m be the number of components of ηk, then rk in equation (25) is a chi-square distributed with m−1 degrees of freedom. With a fixed fault alarm value Γk is obtained from the inverse chi-square cumulative distribution function, which can be done by the offline Monte Carlo simulation. Additionally, Ck is the denominator of Kalman filter gain so as to compute rk using equation (25) only, the proposed monitor algorithm is simple and needs less computation burden, and is flexible without modifying the existing integration structure. In order to make a more clear investigation, the output of the proposed interference monitor in/without the presence of interference is studied, as shown in Figure 3. Figure 3(a) presents a relatively smooth line without any interference, which exhibits better correlation of the monitor output; the performance of the monitor in the presence of interference is provided; and jamming-to-noise ratio (JNR) is used as a metric, as shown in Figure 3(b). It is obviously seen that an abrupt change happens when the interference occurs, and this sudden change can be acquired by the monitor if the abrupt is higher than the predefined Γk.Figure 3
Output plot of the interference monitor. (a) Without the interference. (b) In the presence of interference, JNR = 18 dB.
(a)(b)


### 3.3. Subspace Tracking Algorithm

Once the interference is present and has been monitored, subspace tracking algorithm is used to track the interference signal basis so as to estimate the instantaneous frequency of interference. In this paper, fast data projection method (FDPM) is chose as the main detection method, which can track either signal or noise subspace weighting basis matrix simultaneously [23].


#### 3.3.1. Estimate of the Covariance Matrix

Sample covariance matrix should be first estimated since its eigenvalues are used by the subspace-based detection algorithm and the number estimate of interference source in this paper. In general, two types of methods that are called exponential window and truncated window estimate are used, respectively. For the exponential window, we can write:(28)Ct=∑i=−∞tγt−ixixiH,where γ is the scale factor 0<γ≤1. For the truncated window, we can write:(29)Ct=∑i=t−l+1tγt−ixixiH,where l is the truncated length and sliding window is the case of γ=1.


#### 3.3.2. Determining the Number of Interference Sources

For any subspace tracking algorithm, the number of interference sources, which means the rank of signal subspace basis matrix, should be first detected. In this work, the interference number is assumed to be known for simplicity. In fact, the number can be accomplished by information theoretic criteria [27], such as Akaike information criterion (AIC) and minimum description length (MDL), or sources dependence previously (SDC) [28].


#### 3.3.3. FDPM Subspace Tracking Method

Orthogonal iteration is a common design principle for any subspace tracking method, which can be applied to compute the dominant subspace basis matrix of sample covariance iteratively, and can be written as(30)Uk=orthonormCUk−1,where “orthonorm” represents orthonormalization that can be implemented by means of QR-decomposition or Gram-Schmidt. For the known number of signals L, orthogonal iteration will converge to the signal basis matrix UN×L with an exponential rate λL+1/λLk, which means(31)limk⟶∞Uk=u1u2⋯uL,where uii=1L is the column basis vector.FDPM starts from the well-known data projection method (DPM) and is a fast implementation version [29]. With the simple estimation of covariance Ck=rkrkH, DPM can be implemented by the following formula(32)yk=Uk−1Trk,Sk=Uk−1±μrkykT,Uk=orthonormSk,where “orthonorm” is performed by Gram-Schmidt process, μ>0 is a step size and usually small for numerical stability, “+” stands for tracking the signal subspace, and while “−” for the noise subspace. Since the utilization of Gram-Schmidt is time consuming requiring ONL2 operations that is not suitable for online application such as GNSS equipment, we can use other algorithms to replace Gram-Schmidt in “orthonorm” step in order to reduce the computation load. FPDM is one of these algorithms.For orthnormalizing (32), Householder transformation matrix, defined as Hk, is used(33)ak=yk−ykξ1,Hk=IL−2ak2akakH,where ξ1=1,0,…,0T and Hk is used to perform orthnormalization Uk=TkHk satisfying(34)UkHUk=IL.The steps of FDPM are shown in Table1, and the computation complexity is ONL. More details in terms of complexity, convergence rate, and pseudocode can be found in [23]. The tracking performance of FDPM algorithm is studied, as shown in Figure 4. To evaluate sudden change of frequency, a sine signal with frequency change at sample 100 is used. As Figure 4 described, FDPM algorithm has the similar performance with the algorithm using SVD while lower complexity.Table 1
Pseudocode and cost of FDPM subspace tracking method.
Initial valueU0 should be orthonormal.Input sequential data vector:rkCost1.yk=Uk−1HrkNL2.Tk=Uk−1+μrkykHL3.ak=yk−ykξ1L4.Zk=Tk−2/ak2TkakakHNL5.Uk=normalizeZkLFigure 4
Frequency estimate using FDPM subspace tracking method with exponential widow, a sine signal with 0.3 Hz is used to test,μ=0.008, SNR = 15 dB.In the end of this section, we simply consider the interference mitigation method. For narrow-band interference suppression, notch filter is an effective and commonly used technique [30]. With the filter, the part of spectrum held by the interference will be removed. The weakness of notch filtering is that only narrow frequency band can be mitigated. In this paper, notch filter is adopted as the mitigation method for its simplicity and easy-to-use.


## 3.1. Tightly-Coupled INS/GNSS with Kalman Filter

Data fusion based on Kalman filter is a common method for vehicle positioning, and consistent check using the innovation of INS/GNSS system can provide a monitor of interference. By tightly coupled mechanization, the pseudo-range and pseudo-range rate measurements from GNSS equipment are used to aid INS. Meanwhile, INS can also aid GNSS receiver to track the satellites signals more fastly. Extended Kalman filter (EKF) is used for data fusion due to its simpleness and stability, which provides an optimal estimation of the states. A general design frame for the tightly coupled INS/GNSS system is studied, as shown in Figure2.Figure 2
A common design for tightly coupled INS/GNSS.Corrected measurements from the GNSS code tracking loop in earth centered earth fixed (ECEF) frame can be written as(12)ρ˜=ρ˜R−δρ^I−δρ^T+δρ^C,ρ˙˜=ρ˙˜R+δρ˙^C,where ρ˜R, ρ˙˜R denote the pseudo-range and pseudo-range rate from GNSS receiver, respectively. δρ^I, δρ^T are the ionospheric and tropospheric refraction errors, respectively. δρ^C, δρ˙^C, respectively, represent the clock error and clock error shift. An approximate estimate of ρ˜R can be written as(13)ρ˜R=xa−xs2+ya−ys2+za−zs2=ts−ta×c,where xa,ya,za, xs,ys,zs are the receiver and satellite position in ECEF frame, ta−ts represents the signal transmission time from satellite to GNSS receiver, and c is the speed of light. It should be noticed that some superscript or subscript are omitted for simplicity. For example, the signal from lth satellite in ECEF frame may be written as ρ˜Re,l.Since the states of INS and GNSS are independent, the state vector and state transformation matrix are given by(14)x=xINSxGNSS,ϕ=ϕINS00ϕGNSS.The estimated 17 components of state vector are written as(15)x=δΨδυδrbgbaδρδρ˙,where δp, δv, and δA are the position, velocity, and attitude state errors, respectively, bg and ba are the gyroscope and accelerometer biases, respectively, which should be estimated for any applications; δρ is the receiver clock error; and δρ˙ is the receiver clock drift, which are processed and corrected by GNSS. The component of x may have different forms in different coordinate frame. For instant, when in local navigation frame δr can be represented as latitude, longitude, and height error. After corrected, INS output can be written as(16)C^αβγ=δC^αβγC˜αβγ,r^=δr^−r˜,υ^=δυ˜−υ˜,where α, β, γ represent the body, reference, and projection coordinate system, respectively.The predicted statex using inertial navigation equation is expressed as(17)x^k−=ϕkx^k−1+εk,where x^k− denotes the predicted state vector at epoch k, ϕk is the state transformation matrix, x^k−1 is the updated state vector at epoch k − 1, andεk is the process noise. ϕk with first-order approximate in ECEF frame is presented as(18)ϕke=I3−Ωτs030303Cbeτs00F21eτsI3−2ΩτsF23eτsCbeτs030003I3τsI3030300030303I3030001,301,301,301,3I30001,301,301,301,3031τs01,301,301,301,30301,where Ω is the skew symmetric matrix of the earth rotation angle velocity ωie in ECEF frame(19)Ω=0−ωie0ωie00000,τs is the time epoch of IMU processing unit. F21e and F23e are, respectively, defined as(20)F21e=−C^bef^ibb∧,F23e=−2γ^ibereSeL^br^eber^ebe.The measurement model can be written as(21)zk=Hkxk+νk.where zk is the measurements vector of pseudo-range and pseudo-range rate from GNSS, νk is the measurement noise at epoch k, and Hk is the geometry matrix including line-of-sight vectors from vehicle to satellites in the ECEF frame and can be written as(22)δzk=δzρ,kδzρ˙,k,Hk=01,301,3ua1γT01,301,31001,301,3ua2γT01,301,310⋮⋮⋮⋮⋮⋮⋮01,301,3uamγT01,301,310⋯⋯⋯⋯⋯⋯⋯01,3ua1γT01,301,301,30101,3ua2γT01,301,301,301⋮⋮⋮⋮⋮⋮⋮01,3uamγT01,301,301,301,where u is the line of sight, defined as(23)u=∂ρ∂x,∂ρ∂y,∂ρ∂x.


## 3.2. Monitor Based on the Consistent Check of Measurements

Some alarm mechanisms should be provided first when a fault occurred. In this work, we assume that the fault is caused by the interference only for simplicity, and an interference monitor algorithm is proposed to check whether the interference is present or not. In general, consistent check of the measurements is a common method used in receivers. Similar to receiver autonomous integrity monitoring (RAIM) derived from the least squares estimation for GNSS receiver and can be viewed as a hypothesis testing problem [26], we use the innovation of Kalman filter as a monitor in the INS/GNSS system. The main concept can be explained as follows: when there is no interference, the measurements of GNSS and INS should be with higher correlation. While in the presence of interference, due to the GNSS signals degrade seriously, the deviation of INS and GNSS measurement become unacceptable and inconsistent. So by checking sudden changes of the measurement residual compared with a predefined threshold, we can decide whether the interference is present or not. The innovation formula is(24)ηk=zk−Hkx^k−,where xk− is a prior of state vector at time k from the inertial navigation equation defined in (17).From the hypothesis testing theory, a test statisticrk as the quadratic form of the normalized innovation at time k is defined as(25)rk=ηkTCk−1ηk,where Ck is the covariance matrix of the innovation ηk at time k and can be calculated as(26)Ck=EηkηkT=HkPk−HkT+Rk,where Rk is the covariance of measurement noise and Pk− is the predicted covariance of state error vector at time k(27)Pk−=Ex^k−−xkx^k−−xkT.For a properly selected thresholdΓk, interference occurrence or not can be inferred by comparing the value of rk with Γk. Let m be the number of components of ηk, then rk in equation (25) is a chi-square distributed with m−1 degrees of freedom. With a fixed fault alarm value Γk is obtained from the inverse chi-square cumulative distribution function, which can be done by the offline Monte Carlo simulation. Additionally, Ck is the denominator of Kalman filter gain so as to compute rk using equation (25) only, the proposed monitor algorithm is simple and needs less computation burden, and is flexible without modifying the existing integration structure. In order to make a more clear investigation, the output of the proposed interference monitor in/without the presence of interference is studied, as shown in Figure 3. Figure 3(a) presents a relatively smooth line without any interference, which exhibits better correlation of the monitor output; the performance of the monitor in the presence of interference is provided; and jamming-to-noise ratio (JNR) is used as a metric, as shown in Figure 3(b). It is obviously seen that an abrupt change happens when the interference occurs, and this sudden change can be acquired by the monitor if the abrupt is higher than the predefined Γk.Figure 3
Output plot of the interference monitor. (a) Without the interference. (b) In the presence of interference, JNR = 18 dB.
(a)(b)


## 3.3. Subspace Tracking Algorithm

Once the interference is present and has been monitored, subspace tracking algorithm is used to track the interference signal basis so as to estimate the instantaneous frequency of interference. In this paper, fast data projection method (FDPM) is chose as the main detection method, which can track either signal or noise subspace weighting basis matrix simultaneously [23].


### 3.3.1. Estimate of the Covariance Matrix

Sample covariance matrix should be first estimated since its eigenvalues are used by the subspace-based detection algorithm and the number estimate of interference source in this paper. In general, two types of methods that are called exponential window and truncated window estimate are used, respectively. For the exponential window, we can write:(28)Ct=∑i=−∞tγt−ixixiH,where γ is the scale factor 0<γ≤1. For the truncated window, we can write:(29)Ct=∑i=t−l+1tγt−ixixiH,where l is the truncated length and sliding window is the case of γ=1.


### 3.3.2. Determining the Number of Interference Sources

For any subspace tracking algorithm, the number of interference sources, which means the rank of signal subspace basis matrix, should be first detected. In this work, the interference number is assumed to be known for simplicity. In fact, the number can be accomplished by information theoretic criteria [27], such as Akaike information criterion (AIC) and minimum description length (MDL), or sources dependence previously (SDC) [28].


### 3.3.3. FDPM Subspace Tracking Method

Orthogonal iteration is a common design principle for any subspace tracking method, which can be applied to compute the dominant subspace basis matrix of sample covariance iteratively, and can be written as(30)Uk=orthonormCUk−1,where “orthonorm” represents orthonormalization that can be implemented by means of QR-decomposition or Gram-Schmidt. For the known number of signals L, orthogonal iteration will converge to the signal basis matrix UN×L with an exponential rate λL+1/λLk, which means(31)limk⟶∞Uk=u1u2⋯uL,where uii=1L is the column basis vector.FDPM starts from the well-known data projection method (DPM) and is a fast implementation version [29]. With the simple estimation of covariance Ck=rkrkH, DPM can be implemented by the following formula(32)yk=Uk−1Trk,Sk=Uk−1±μrkykT,Uk=orthonormSk,where “orthonorm” is performed by Gram-Schmidt process, μ>0 is a step size and usually small for numerical stability, “+” stands for tracking the signal subspace, and while “−” for the noise subspace. Since the utilization of Gram-Schmidt is time consuming requiring ONL2 operations that is not suitable for online application such as GNSS equipment, we can use other algorithms to replace Gram-Schmidt in “orthonorm” step in order to reduce the computation load. FPDM is one of these algorithms.For orthnormalizing (32), Householder transformation matrix, defined as Hk, is used(33)ak=yk−ykξ1,Hk=IL−2ak2akakH,where ξ1=1,0,…,0T and Hk is used to perform orthnormalization Uk=TkHk satisfying(34)UkHUk=IL.The steps of FDPM are shown in Table1, and the computation complexity is ONL. More details in terms of complexity, convergence rate, and pseudocode can be found in [23]. The tracking performance of FDPM algorithm is studied, as shown in Figure 4. To evaluate sudden change of frequency, a sine signal with frequency change at sample 100 is used. As Figure 4 described, FDPM algorithm has the similar performance with the algorithm using SVD while lower complexity.Table 1
Pseudocode and cost of FDPM subspace tracking method.
Initial valueU0 should be orthonormal.Input sequential data vector:rkCost1.yk=Uk−1HrkNL2.Tk=Uk−1+μrkykHL3.ak=yk−ykξ1L4.Zk=Tk−2/ak2TkakakHNL5.Uk=normalizeZkLFigure 4
Frequency estimate using FDPM subspace tracking method with exponential widow, a sine signal with 0.3 Hz is used to test,μ=0.008, SNR = 15 dB.In the end of this section, we simply consider the interference mitigation method. For narrow-band interference suppression, notch filter is an effective and commonly used technique [30]. With the filter, the part of spectrum held by the interference will be removed. The weakness of notch filtering is that only narrow frequency band can be mitigated. In this paper, notch filter is adopted as the mitigation method for its simplicity and easy-to-use.


## 3.3.1. Estimate of the Covariance Matrix

Sample covariance matrix should be first estimated since its eigenvalues are used by the subspace-based detection algorithm and the number estimate of interference source in this paper. In general, two types of methods that are called exponential window and truncated window estimate are used, respectively. For the exponential window, we can write:(28)Ct=∑i=−∞tγt−ixixiH,where γ is the scale factor 0<γ≤1. For the truncated window, we can write:(29)Ct=∑i=t−l+1tγt−ixixiH,where l is the truncated length and sliding window is the case of γ=1.


## 3.3.2. Determining the Number of Interference Sources

For any subspace tracking algorithm, the number of interference sources, which means the rank of signal subspace basis matrix, should be first detected. In this work, the interference number is assumed to be known for simplicity. In fact, the number can be accomplished by information theoretic criteria [27], such as Akaike information criterion (AIC) and minimum description length (MDL), or sources dependence previously (SDC) [28].


## 3.3.3. FDPM Subspace Tracking Method

Orthogonal iteration is a common design principle for any subspace tracking method, which can be applied to compute the dominant subspace basis matrix of sample covariance iteratively, and can be written as(30)Uk=orthonormCUk−1,where “orthonorm” represents orthonormalization that can be implemented by means of QR-decomposition or Gram-Schmidt. For the known number of signals L, orthogonal iteration will converge to the signal basis matrix UN×L with an exponential rate λL+1/λLk, which means(31)limk⟶∞Uk=u1u2⋯uL,where uii=1L is the column basis vector.FDPM starts from the well-known data projection method (DPM) and is a fast implementation version [29]. With the simple estimation of covariance Ck=rkrkH, DPM can be implemented by the following formula(32)yk=Uk−1Trk,Sk=Uk−1±μrkykT,Uk=orthonormSk,where “orthonorm” is performed by Gram-Schmidt process, μ>0 is a step size and usually small for numerical stability, “+” stands for tracking the signal subspace, and while “−” for the noise subspace. Since the utilization of Gram-Schmidt is time consuming requiring ONL2 operations that is not suitable for online application such as GNSS equipment, we can use other algorithms to replace Gram-Schmidt in “orthonorm” step in order to reduce the computation load. FPDM is one of these algorithms.For orthnormalizing (32), Householder transformation matrix, defined as Hk, is used(33)ak=yk−ykξ1,Hk=IL−2ak2akakH,where ξ1=1,0,…,0T and Hk is used to perform orthnormalization Uk=TkHk satisfying(34)UkHUk=IL.The steps of FDPM are shown in Table1, and the computation complexity is ONL. More details in terms of complexity, convergence rate, and pseudocode can be found in [23]. The tracking performance of FDPM algorithm is studied, as shown in Figure 4. To evaluate sudden change of frequency, a sine signal with frequency change at sample 100 is used. As Figure 4 described, FDPM algorithm has the similar performance with the algorithm using SVD while lower complexity.Table 1
Pseudocode and cost of FDPM subspace tracking method.
Initial valueU0 should be orthonormal.Input sequential data vector:rkCost1.yk=Uk−1HrkNL2.Tk=Uk−1+μrkykHL3.ak=yk−ykξ1L4.Zk=Tk−2/ak2TkakakHNL5.Uk=normalizeZkLFigure 4
Frequency estimate using FDPM subspace tracking method with exponential widow, a sine signal with 0.3 Hz is used to test,μ=0.008, SNR = 15 dB.In the end of this section, we simply consider the interference mitigation method. For narrow-band interference suppression, notch filter is an effective and commonly used technique [30]. With the filter, the part of spectrum held by the interference will be removed. The weakness of notch filtering is that only narrow frequency band can be mitigated. In this paper, notch filter is adopted as the mitigation method for its simplicity and easy-to-use.


## 4. Experimental Results

In this section, the performance of the proposed framework based on INS/GNSS monitor and subspace tracking method in the presence of CW interference is evaluated by experimental data [5]. A tightly coupled integration architecture is applied and extended Kalman filter is used as data fusion. The frequencies of two CW interference are 9.25 MHz and 11.25 MHz, respectively. In order to illustrate that the subspace tracking method has excellent performance for frequency estimation, other two well-known time-frequency methods, STFT and SPWVD, are used to make a comparison. GPS L1 signal is used with the carrier-to-noise ratio CN0=45 dB-Hz. The trajectory data of a vehicle are obtained from [5] with 60 seconds time duration, as shown in Figure 5. More detailed configurations are provided in Table 2.Figure 5
Trajectory of the vehicle for 60 seconds in NED frame.Table 2
Configurations for the tightly coupled INS/GNSS integration.
ParameterValueGyro bias0.1°/hGyro random walk0.01°/h1/2Accelerometer bias100μgAccelerometer random walk10μgGPS intermediate frequency12.45 MHzGPS sampling frequency25 MHzGPS position measurement noise10 mGPS velocity measurement noise0.1 m/sGPS initial clock drift uncertainty0.1 m/sGPS initial clock uncertainty10 m


### 4.1. Detection Probability and Receiver Operating Characteristic

Receiver operating characteristic (ROC) [31] is commonly utilized to assess the detection performance of a receiver. For a given JNR, ROC is a function of the detection probability and the probability of false alarm. Since jamming signals mainly affect the GNSS performance, we focus on ROC analysis from the GNSS receiver side. In this paper, CW signals are used as the jamming signals, and Monte Carlo method is considered as offline experiment. GNSS signals are treated as white Gaussian noise (WGN) due to its low power. We, respectively, use JNR = 3 dB, 5 dB, and 7 dB to analyze the ability of jamming detection of the proposed algorithm, as shown in Figure 6. We can see that for a given false arm ratio, Pfa, detection probability, Pd, increases with the lager JNR.Figure 6
ROC with the JNR = 3 dB, 5 dB, and 7 dB, respectively.


### 4.2. D&M Analysis for the Integrated GNSS/INS

The contour plots of GPS L1 signal in the presence of two CW interference are presented, respectively, using the proposed algorithm, STFT, and SPWVD in time-frequency domain, as shown in Figure7. In Figure 7(a), two clear frequency lines are obtained by proposed algorithm without any disturbance frequencies, and it is easy to estimate the frequency of interference. In Figure 7(b), the window size is 80, two relatively clear lines present while some distortions also appear due to cross-interfering terms, so in this range of CN0, SPWVD has a poor performance, and the frequency estimation is not so accurate as that by our provided method. The worst case is that in Figure 7(c) by STFT, where we hardly find any frequency information in the time-frequency plane. Although the larger number of data sample can make a better resolution of frequency, the degradation of time resolution and computation complexity also increase. In general, FDPM algorithm outperforms other two time-frequency algorithms.Figure 7
Contour plot for GPS L1 signal in the presence of CW interference.C/N0=45 dB-Hz, JNR = 18 dB. (a) Contour plot using FDPM with samples 80. (b) Contour plot using SPWVD with a Kaiser window of samples 80. (c) Contour plot using STFT with a Kaiser window of samples 80.
(a)(b)(c)Position errors from the north, east, and down directions in local navigation coordinate are analyzed, respectively, as shown in Figure8. For simplicity, only CW interference (JNR = 45 dB) is added and the interference signal lasts for 50 seconds, which leading to about 1 kilometer pseudo-range deviation for GNSS receiver. Here we should notice that the RF interference only affects the GNSS receiver, and IMU errors are ignored due to antijamming interference. Mean square error (MSE) is used as a criterion. From Figure 8, we can observe that the interference can affect the position accuracy seriously (red line) for the INS/GNSS system. When the proposed framework is applied, the interference can be mitigated and the performance of the integrated system can be improved (blue line). In addition, Figure 8(c) shows a relative larger error than those in Figures 8(a) and 8(b), respectively. The reason is that before corrected by GNSS, the low-cost INS commonly has larger position error than that of GNSS, so position errors are mainly caused by GNSS for coupled GNSS/INS system, and we know that down (up) direction error of GNSS is larger than that of horizontal direction for GNSS receiver.Figure 8
Position errors comparison between using interference mitigation and without mitigation. (a) North position error. (b) East position error. (c) Down position error.
(a)(b)(c)


## 4.1. Detection Probability and Receiver Operating Characteristic

Receiver operating characteristic (ROC) [31] is commonly utilized to assess the detection performance of a receiver. For a given JNR, ROC is a function of the detection probability and the probability of false alarm. Since jamming signals mainly affect the GNSS performance, we focus on ROC analysis from the GNSS receiver side. In this paper, CW signals are used as the jamming signals, and Monte Carlo method is considered as offline experiment. GNSS signals are treated as white Gaussian noise (WGN) due to its low power. We, respectively, use JNR = 3 dB, 5 dB, and 7 dB to analyze the ability of jamming detection of the proposed algorithm, as shown in Figure 6. We can see that for a given false arm ratio, Pfa, detection probability, Pd, increases with the lager JNR.Figure 6
ROC with the JNR = 3 dB, 5 dB, and 7 dB, respectively.


## 4.2. D&M Analysis for the Integrated GNSS/INS

The contour plots of GPS L1 signal in the presence of two CW interference are presented, respectively, using the proposed algorithm, STFT, and SPWVD in time-frequency domain, as shown in Figure7. In Figure 7(a), two clear frequency lines are obtained by proposed algorithm without any disturbance frequencies, and it is easy to estimate the frequency of interference. In Figure 7(b), the window size is 80, two relatively clear lines present while some distortions also appear due to cross-interfering terms, so in this range of CN0, SPWVD has a poor performance, and the frequency estimation is not so accurate as that by our provided method. The worst case is that in Figure 7(c) by STFT, where we hardly find any frequency information in the time-frequency plane. Although the larger number of data sample can make a better resolution of frequency, the degradation of time resolution and computation complexity also increase. In general, FDPM algorithm outperforms other two time-frequency algorithms.Figure 7
Contour plot for GPS L1 signal in the presence of CW interference.C/N0=45 dB-Hz, JNR = 18 dB. (a) Contour plot using FDPM with samples 80. (b) Contour plot using SPWVD with a Kaiser window of samples 80. (c) Contour plot using STFT with a Kaiser window of samples 80.
(a)(b)(c)Position errors from the north, east, and down directions in local navigation coordinate are analyzed, respectively, as shown in Figure8. For simplicity, only CW interference (JNR = 45 dB) is added and the interference signal lasts for 50 seconds, which leading to about 1 kilometer pseudo-range deviation for GNSS receiver. Here we should notice that the RF interference only affects the GNSS receiver, and IMU errors are ignored due to antijamming interference. Mean square error (MSE) is used as a criterion. From Figure 8, we can observe that the interference can affect the position accuracy seriously (red line) for the INS/GNSS system. When the proposed framework is applied, the interference can be mitigated and the performance of the integrated system can be improved (blue line). In addition, Figure 8(c) shows a relative larger error than those in Figures 8(a) and 8(b), respectively. The reason is that before corrected by GNSS, the low-cost INS commonly has larger position error than that of GNSS, so position errors are mainly caused by GNSS for coupled GNSS/INS system, and we know that down (up) direction error of GNSS is larger than that of horizontal direction for GNSS receiver.Figure 8
Position errors comparison between using interference mitigation and without mitigation. (a) North position error. (b) East position error. (c) Down position error.
(a)(b)(c)


## 5. Conclusions

This paper provides a low complexity interference D&M method for vehicle positioning. A framework of interference monitor, detection, and mitigation for the INS/GNSS system is proposed in a tightly coupled architecture, with a Kalman filter for data fusion. The framework can significantly improve the performance of the integration system in the presence of interference. A simple and effective monitor based on the INS/GNSS innovation measurements is proposed to make an alarm when the interference is present, which can decide when to start the interference detection procedure so as to reduce the computation burden. Interference detection using subspace tracking method outperforms well-known time-frequency algorithms such as STFT and SPWVD in time-frequency domain. The performance of proposed framework is also verified by the computer simulation.Future research of this work will include the implementation of interference monitor in the different coupled architecture; the implementation and performance analysis in the presence of various interference types, such as chirp-type interference, should also be considered.



---

*Source: 7502191-2022-05-20.xml*