Detection and Mitigation of Spoofing Attacks Against Time Synchronization and Positioning

Global Navigation Satellite System (GNSS) receivers are vulnerable to intentional spoofing attacks which can manipulate position, velocity, and time (PVT) measurements. Previous work has demonstrated that Time Synchronization Attacks (TSAs) can be detected and mitigated using sparse optimization techniques which reveal spoofers’ presence in inflicted signals’ derivative domains. Aiming to provide an efficient protection algorithm against spoofing, this paper expands the scope of the sparse signal processing framework to address more complicated attacks for stationary and low-dynamic receivers. In particular, TSAs against stationary receivers should be addressed differently if the position and velocity are manipulated by the spoofer at the same time. A new sparse processing method is presented employing a novel linearization of the measurement equation that includes attacks against time, position, and velocity. The method is assessed against both authentically and synthetically spoofed signals to verify its robustness in two test beds: 1) a lab-based software defined GPS receiver; and 2) a commercial hand-held device.


I. INTRODUCTION
The emergence of Global Navigation Satellite Systems (GNSS) in a wide range of infrastructures, such as the smart power grid, financial transactions, and aviation mapping, upholds its growing popularity in contemporary applications [1].Despite its wide use of applications, GPS is vulnerable to various types of external disturbances.
While some perturbations are attributed to unavoidable environmental noise from the ionosphere and troposphere, other disturbances originate from a harmful third-party hacker.Particular mechanisms used by hackers against GNSS receivers include jamming and spoofing [2].Both types of malicious GPS attacks intend to induce receivers to compute incorrect position, velocity, and timing (PVT) The associate editor coordinating the review of this manuscript and approving it for publication was Tianhua Xu .solutions.The jamming attacks completely block one or more communication channels from satellites.Smart spoofing occurs when the actor injects counterfeit signals to all open GPS channels uniformly with deliberate deviations.With more than five deliberately deceived channels, a sophisticated spoofer has the capability of bypassing the traditional receiver autonomous integrity monitoring (RAIM) method [3].The common spoofing schemes are categorized depending on the impacted PVT domain into time synchronization attacks (TSAs) and position spoofing.
The purpose of TSA is to disrupt timing estimations in GNSS receivers, while position and velocity estimates may remain intact.Thus, such attacks are frequently targeting stationary receivers.One well-known example is the power grid which employs Phasor Measurement Units (PMU) to read synchrophasor measurements from network nodes.The measurements are globally synchronized with timing provided by GPS receiver clocks [4].Consequently, the impact of TSAs is a significant concern [5].In the work of [6], power system state estimation with TSA-impacted PMU measurements is developed.TSAs can also affect the localization of transmission line faults and other grid monitoring routines [7], [8].
Position attacks aim to corrupt spatial-related estimates in GNSS receivers [9].Intermediate position attacks are capable of deceiving a victim receiver by forcing them to compute deviated user location and velocity.To implement position attacks on a commercial off-the-shelf (COTS) receiver, counterfeit signals shift the code phase in baseband correlators [10], and the aftermath of attack is manifested in navigation domain data, such as pseudorange sequences [11].Practical applications of position spoofing are shown in experiments with Unmanned Aerial Vehicles (UAVs) [12] or vehicles in motion [13], [14].In the work of [15], experiments with UAVs entail the hacker injecting counterfeit signals wirelessly and misleading the computed location 50 meters off from the true position.
Research has demonstrated that spoofing can simultaneously impact position and timing, even if only one of the two is targeted.For instance, in [10], the generation of spoofing signals to mislead the z-coordinate position estimate leads to deviations in other coordinates along with clock bias estimates for a stationary receiver scenario.Likewise, to mount a successful TSA, some deviation in the position estimates may be allowed [16].Therefore, there is a need for novel anti-spoofing techniques that can withstand combined attacks against position and timing.
According to [3], most of GNSS spoofing attacks can be mathematically modeled as the sum of an authentic transmitted signal y authentic , the attacking signal from the spoofer y spoof , and white Gaussian noise n, which can be equivalently written as y total = y authentic + y spoof + n where y total is the received signal at the antenna.Under the influence of y spoof , the victim receiver starts to produce inaccurate user location or timing solutions.In such environment, spoofing countermeasures can be categorized according to whether they 1) recognize and classify the presented inflicting signal y spoof (detection); or additionally 2) retrieve the authentic signal y authentic (mitigation).Once both procedures are successfully performed, accurate estimates of true PVT states are produced.

A. PRIOR ART
A variety of detection and mitigation methods have been developed in literature; see e.g., [3] and [17] for reviews.Generally, the reported approaches can be categorized into (1) system-level protection, and (2) signal-level countermeasure techniques.Majorly dominated by encryption methods, system-level protection furnishes symmetric or asymmetric encrypted code-phase keys on GNSS transmitter and receiver [18], [19].For instance, the work of [20] showcases encryption-based protection against replay-type spoofing with application to civilian and military receivers.
Signal-level techniques aim to detect anomalies in the receiver observables on the GPS baseband or navigation domains [21].The baseband domain countermeasures, which commonly exploit real-time correlation peaks, rely on power and automatic gain control [22], [23], vector tracking loop (VTL) [24], and Least Absolute Shrinkage and Selection Operator (LASSO) minimization [25].Anti-spoofing protection techniques operating in the baseband domain typically require physical modification of the receiver circuitry; thus harder to integrate into COTS receivers.
It is desirable to develop techniques that rely on readily available data at the navigation domain in the receiver, such as pseudoranges, pseudorange rates, and satelliterelated information.Navigation domain countermeasures exploit statistics in pseudorange residuals or carrier-to-noise (CNR) ratio [26], [27], [28], drift-monitoring [29], sparsity characteristics of realistic spoofing profiles [21], [30], and robust estimation based on pseudoranges and pseudorange rates [31].
The works in [27] and [29] focus mostly on developing spoofing detection schemes and do not readily yield corrected PVT estimates.The works in [21], [30], and [31] focus on spoofing detection and mitigation of TSAs only.However, when spoofing attacks target both time and position domains, the aforementioned algorithms [21], [30], [31] are not generally capable of isolating and mitigating the TSA impact even for a static receiver.A need therefore arises to develop spoofing mitigation methods addressing both position and timing related domains.

B. CONTRIBUTIONS
The present paper develops a detection and mitigation technique against timing and position attacks in static and low-dynamic receivers.The static receiver is defined as having constant position and zero velocity.Meanwhile, the low-dynamic receiver model entails that its velocity follows a random walk [32].
The method operates in the navigation domain and relies on characterizing spoofing attack profiles according to the order of derivative in which they appear to be sparse.Most realistic spoofing profiles can therefore be classified as first, second, on third order attacks.The characterization of spoofing according to sparsity in the derivative was first introduced in [21] and [30] for TSAs only and is extended here to include attacks against position in the x, y, and z-coordinates and associated velocities.To the best of our knowledge, this is the first paper that analyzes a model of joint spoofing against position and timing domains.
The estimation of PVT states involves a set of nonlinear measurement equations.To deal with this challenge, a new linearized model is developed that is different than traditional ones used in the GPS literature [32].Specifically, the new linearized model relates spoofing attacks against position, velocity, clock bias, and clock drift with pseudoranges and pseudorange rates.Building on the aforementioned linearized model, an optimization formulation is developed to jointly estimate the authentic PVT states and the multi-domain spoofing attacks.Appropriate sparsity promoting regularizations are included to this end.The overall formulation amounts to a convex quadratic program that can be easily solved.
The contributions of this paper are summarized as follows: 1) A novel GPS dynamical model is designed to capture position, velocity, and timing related states of stationary and low-dynamic receivers featuring a linearization method incorporating spoofing effects.2) An efficient optimization model is developed to estimate PVT states, and joint TSAs and position spoofing attacks.With proper choices of tuning parameters, the optimization technique can capture sophisticated and smooth shapes of spoofing profiles, potentially higher than third order attacks, regardless of measurement integrity.
3) The technique is tested on two receiver platforms representing the commercial and research environments, namely, Huawei Pro30 and an in-house software defined radio (SDR) receiver [33].The remainder of this paper is organized as follows.Section II introduces the GPS dynamic model.Section III models the sparse characteristics of spoofing profiles and integrated joint attacks.Section IV develops the anti-spoofing technique and lists the tuning parameters of the optimization model.Section V presents the testbeds and methodology for assessing the novel anti-spoofing technique.Section VI showcases the performance of the optimization model on realistic spoofing scenarios.Limitations of the present work are articulated in Section VII.Section VIII concludes the paper and provides pointers to future work.

II. GPS PVT DYNAMIC MODEL
This section provides an overview of the GPS dynamics, which comprise the PVT states dynamic equation and the linearized measurement equations.

A. MEASUREMENT AND STATE DYNAMIC MODEL
During satellite-to-receiver signal transmission, the receiver collects binary-coded navigation data including satellite clock information and orbital parameters to determine satellite position and velocity [34].Then, the receiver estimates the user location at transmission epochs indexed by k = 1, 2, . . ., K .Let n = 1, 2, . . ., N denotes the number of visible satellites at k.The n th satellite position is expressed as Earth Fixed (ECEF) coordinates.In the same manner, the receiver position is denoted by The true distance d n between the receiver and satellite n can be calculated using the Euclidean distance formula as d n = ||p n − p u || 2 .An alternative way of expressing d n is based on the difference between absolute GPS times of signal transmission t GPS n and reception t GPS R , respectively recorded by the satellites and receiver, as follows where c represents the speed of light.Imperfect clocks at the transmitter and receivers are modeled by including appropriate biases.Specifically, upon introducing biases b n and b u respectively affecting the satellite and receiver clocks, the perceived signal transmission and reception times are written as The pseudorange is computed based on the biased clock times as ρ n = c(t R − t n ).Combining the previously stated equations together, the pseudorange is written as while ϵ ρ n refers to measurement noise.
In addition to the pseudorange, the receivers measure the pseudorange rate, which depends on the Doppler shift.As satellites and the receiver maneuver, the relative velocity difference created by both ends is manifested on carrier frequency discrepancy.The Doppler rate, or equivalently pseudorange rate, can be expressed by the relative velocity difference projected on the satellite line-of-sight vector and adjusted by frequency offsets, as follows where 1 n represents aforementioned line-of-sight vector defined as In practice, the receiver measures the pseudorange and range rate using separate circuitry.If the pseudoranges and pseudoarange rates are computed from authentic GPS signals, then the two sequences of measurements must pass the measurement integrity check as follows where t represents the discretized time interval between consecutive measurements.In equations ( 2) and (3), navigation data provide the satellite positions p n and satellite clock biases b n , and the user seeks to solve for the receiver position, velocity, and time.The unknown state variable is denoted by x ≡ p u b u v u ḃu ⊤ and are typically computed with traditional techniques such as Weighted Least Squares (WLS) [35] and Kalman filter (KF) [36].
The low-dynamic receiver model entails that the velocity follows a random walk, and the position follows accordingly as the integral of the velocity [32].Specifically, the state 138988 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.dynamical equation is given as follows: where F k is known state transition matrix and w k refers to the state transition noise.
The ensuing section introduces a linearization of the measurement equations ( 2) and (3) that is different than the traditional linearization used in iterative estimation algorithms such as WLS and KF.

B. LINEARIZATION OF MEASUREMENT EQUATIONS
Linearizations of the pseudorange and range rate equations are performed with respect to a reference receiver position (p u,ref [k]) and velocity (v u,ref [k]).Theoretically, for stationary receiver applications such as PMUs, the user can program the position to be the known location of device.In this study, however, the reference is extracted from the initial WLS PVT estimation, whose horizontal dilution-of-precision (HDOP) is reported to be less than 1 meter.While stated in IV-A, spoofing-free initial position and velocity estimates guarantee trustworthy data to initiate linearizing process.This way, the proposed algorithm does not rely on the GPS receiver being static.
Additionally, if the receiver is considered to have a pre-defined route or tends to maneuver in a trajectory that the user can readily compute [37], the linearized model also can be applied via updating position and velocity references in every appropriate interval.Thus, by developing the following effective yet precise linearized measurement model, the formulated problem is considered to be applicable to GPS receivers with stationary or pre-determined dynamics.
The developed linearization approach is somewhat different from conventional iterative techniques used in GPS receivers, in which consecutive iterations of the receiver are linearized with respect to the previous and closest PVT solution.The reference in the current spoofing context indicates some known spoofing-free estimate which is either known from previous receiver iterations or alternative sources.In particular, the receiver can iterate with a frozen coarse reference estimate for several cycles to detect and mitigate spoofing if it appeared during these cycles.
The linearized form of the pseudorange measurement equation ( 2) is stated as where we define is a known quantity as it will be explained shortly: Details on the derivation of the linearization are provided in Appendix.The position vector p u [k] and the receiver clock bias b u [k] are the only unknown variables that need to be computed from (6) (apart from the noise).The right-hand side of ( 7) only includes known quantities as follows: • ρ n [k]: real-time observed pseudorange value given by the receiver at every time k.
• p u,ref [k]: reference point for position estimate; it can be constantly updated at every time epoch k if provided by the end user, otherwise set to p u,ref [1].
• ρ n.ref [k]: reference range corresponding to the linearization point calculated as In the same manner, the pseudorange rate equation ( 3) is linearized as where we define || 2 and the left-hand side of ( 8) is given by The unknown quantities in (8) are and ḃu [k].The following list specifies known values in (9): • ρn [k]: real-time observed range rate sequence provided by the receiver; the user must ascertain that the observed range rate has measurement integrity in regards to ρ n .
• v u,ref [k]: a reference user velocity point where linearization on range rate occurred; this also can be updated on every time epoch k if already provided, otherwise set to v u [1].
The introduced linearization method is mainly conducted upon the assumption that the receiver is characterized by a low-dynamic movement, in which the user position does not fluctuate severely.While the static receiver scenario has practical relevance, this work is not limited to the scenario in which the reference position remains constant throughout the processing time.However, if the user has prior knowledge of where the receiver is highly likely to be located at, the reference pair (p u,ref [k], v u,ref [k]) can be continuously updated.
The linearized equations ( 6) and ( 8) differ from measurement equations commonly used in the literature-see e.g., [32]-in which the models are driven by an interval of PVT estimation.Interestingly, the virtue of the developed linearized measurement sequences hinges upon the fact that ψ n and ψn maintain integrity when computed based on authentic GPS signals: Modeling of joint attacks against timing and position is the theme of the next section.

III. INTEGRATED ATTACK
Signal-level spoofing can generally be modeled as injecting modifying signals to the measurement observables: while s ρ n and s ρn represent spoofing signals on pseudorange and range rate, respectively.When the spoofed pseudoranges and range rates, ρ n,s [k] and ρn,s [k], are used as inputs to conventional PVT estimation algorithms, such as WLS or the KF, the effect is that the receiver will report erroneous PVT solutions to the user.This section analyzes certain characteristics of the injected signals s ρ n and s ρn and the resulting spoofed PVT solutions.Upon application of the linearization developed in Section II-B on ρ n [k] and ρn,s [k], the following set of spoofed linearized measurement equations are obtained where ψ n [k] and ψn [k] are given by ( 6) and ( 8), respectively.The spoofing signals s ρ n and s ρn induce spoofed PVT estimates.The relationship between the attacks on the pseudoranges and range rates on one hand and the resulting unwanted alterations of the PVT estimates on the other can be modeled after the linear equations ( 6) and ( 8) relating the pseudoranges and range rates with the PVT estimates . Specifically, introducing the notations s p , s b , s v , and s ḃ for the respective spoofing-induced deviations of the receiver PVT estimates, the following relationships are modeled after ( 6) and ( 8): The proposed algorithm in this work aims at precisely estimating The following subsections discuss the characteristics in s ρ n and s ρn that enable the formulation of a suitable optimization problem to this end.

A. SPOOFING PROFILES
In this subsection, characteristics of spoofing profiles that pertain to TSAs introduced in our previous work [21], [30] are reviewed and extended to spoofing attacks against positioning.
The work of Time Synchronization Attack Rejection and Mitigation via Sparse Technique (TSARM-S) [30] introduces a novel way of defining shapes of TSA spoofing, or concisely, spoofing profiles.Specifically, it is observed that even smooth spoofing signals exhibit identifiable sparse-like behavior in higher-order discrete-time derivative domains.Then accordingly, the spoofing profiles can be categorized on the basis of the smallest order derivative where the attack exhibits evident spikes.For instance, an attack exhibiting step-like behavior is referred to as first order attack since sparsity occurs in the first derivative of position, that is, velocity.Likewise, second order and third order attacks are identified from the sparse behavior occurring in acceleration and jerk domains, respectively.
The proposed anti-spoofing technique in the present paper attempts to detect sparse spikes not only in derivative domains of s b -which was the theme of our previous work [30]but also those in s p .The following set of equations describes the construction of derivative sequences up to third order on the basis of position z-domain alteration, namely, s z in s p : Equations ( 13) can also be applied to position alteration in s x and s y .For instance, Figure 1 depicts the attack profile s ρ n on pseudoranges and its manifestation as a second order attack that alters the position z-domain alongside its derivatives constructed by (13).The developed spoofing countermeasure in this paper searches for sparsity in the acceleration domain and mitigates the attack shown in s z .The spoofing attacks may satisfy the following integrity equations: If the inflicting signals on velocity or clock drift, s v or s ḃ, satisfy (14), they are referred to as consistent attacks.Otherwise, they are called non-consistent attacks.
In the works of [21] and [31], TSA countermeasure techniques are introduced which utilized the consistency characteristic in clock bias and drift.In contrast, the present work can detect and mitigate consistent and non-consistent attacks against both time and positioning.Also, it is worth noting that non-consistent attacks render pseudorange and pseudorange rates that could fail the integrity check of (4) and may therefore be detected by the receiver's rudimentary FIGURE 1. Second order attack applied on position z-domain.Fig. 1a shows the injected pseudorange spoofing signals s ρ n to create 600m of second order attack on z-domain.Fig. 1b shows the attack s z on z-domain and the first and second derivatives of s z .Note that sz is a sparse signal.
defense mechanism, but without the attack being rejected.On the other hand, the proposed technique offers an adequate mitigation functionality for the receiver by recovering the true PVT states, as will be described in Section IV.In the remainder of the present section, the characteristics of attacks against time and positioning and their relevance in developing the new spoofing mitigation technique are further elaborated.

B. TIME SYNCHRONIZATION ATTACK (TSA)
TSAs strive to solely mislead the estimation of clock bias and drift, leaving spatial-related information intact.TSAs can be modeled by substituting zero sequences for s p and s v in ( 12), yet allowing non-zero values for s b and s ḃ.It follows that TSAs are manifested when the injected pseudorange and range rate spoofing signals on all satellites are identical to s b and s ḃ, respectively:

C. POSITION SPOOFING
Different implementations of position spoofing are developed in the literature, for example, by hijacking the correlation peaks to deceive the receiver's location [26].Independently of the spoofing mechanism, the interpretation of position spoofing in the navigation domain simply translates into deviations appearing on position estimates in x,y, and z coordinates.In order to achieve intermediate spoofing signals [38], the inflicting signals must be applied to all open GPS satellite channels.
Experiments have manifested that subtly fabricated spoofing attacks have the ability to only affect user position estimates, while leaving velocity solution intact [39], which constitutes a non-consistent attack.Fig. 1 depicts the pseudorange deviations on 7 visible satellites for 450 seconds to fabricate 600m deflection on z domain.The relationship between the attack on pseudoranges and the deviation in the z domain is captured by (12), where line-of-sight unit vectors (1 n,ref ) pertaining to different satellites correspond to different sequences s ρ n .This work considers more general position attacks where the carry-off magnitudes and times, as well as the corresponding derivative order where sparsity shows up, may all be different for s x , s y , and s z .

D. JOINT ATTACK
The practicability of joint spoofing, that is, simultaneous TSA and position attack, is attested in [40], in which the spoofer takes advantage of disparate PVT solution sources in a commercial Android mobile receiver.While the receiver accumulates GPS raw data that is transmitted from satellites, separate hardware reads the calculated user position that is threaded in Wi-Fi signals.The spoofer purposefully inserts inaccurate information in both signals, and the joint spoofing is successfully manifested in the commercial receiver.
A generic attack on pseudoranges and pseudorange rates s ρ n and s ρn can generally lead to alterations in all PVT states (x, y, z, b, ẋ, ẏ, ż, ḃ).Even in cases where a single state is targeted, deviations in non-targeted domains can occur.For instance, Fig. 2 depicts the PVT solution for TEXBAT DS4 scenario, in which the spoofer intends a 600-meter push only in position z-domain.However, such attack also creates maximum deviations of 66 meters in position x, 114 meters in position y, and 64 meters in clock bias domains.The pure-TSA countermeasures that we previously investigated in [21], [30], and [31] are not capable of detecting and mitigating such deviations presented in all PVT domains.However, the developed algorithm in this work is anticipated to catch not only the intended attack on z-domain, but also subtle deflections occurred in other domains.

IV. SPOOFING DETECTION AND MITIGATION
This section introduces the state-space dynamic model augmented with joint attacks.It also depicts the development of l 1 -norm minimization technique which detects and mitigates abnormal behaviors presented onto PVT solutions.

A. MODEL ASSUMPTIONS
To achieve successful spoofing detection and mitigation, it is meaningful to clearly indicate the assumptions the optimization model requires.These are summarized next: 1) The set of measurements must be collected from stationary or low-dynamic GPS receivers.
2) The position and velocity references must be reasonably accurate.
3) The set of measurements must include a segment of spoofing-free pseudoranges and range rates at the beginning of the batch window.4) The inflicting signal should present sparsity in some derivative domain (e.g., first, second, or third order).

B. DYNAMICAL MODEL WITH INTEGRATED ATTACK
The GPS receiver system consists of two models: the state dynamic model and the linearized measurement model which includes spoofing.The state update model is given by eq. ( 5) in vector form as x k = F k x k−1 + w k where F k represents the state transition matrix and w k is the state noise vector modeled as zero-mean white Gaussian noise with covariance matrix Q k which depends on the type of oscillator in the receiver clock.Upon combining ( 6), ( 8), (11), and ( 12) and defining , the spoofed measurement equations are written in matrix-vector form as follows where and 1 N ×1 is an N × 1 vector of all ones.
The vector s k in ( 15) includes all domains of the receiver PVT solution where TSAs and position attacks can be manifested: Additionally, ϵ k represents zero-mean Gaussian measurement noise with covariance matrix R k = diag(σ 2 ρ 1 , . . ., σ 2 ρ N , σ 2 ρ1 , . . ., σ 2 ρN ).The combined dynamic models based on ( 5) and ( 15) can be written as The ensuing subsection develops an optimization problem to estimate spoofing attacks in s k , and simultaneously furnish the authentic PVT solution x k .

C. OPTIMIZATION PROBLEM FORMULATION
The optimization problem that jointly solves for PVT estimates x = [x 1 , . . ., xK ] ⊤ and reconstructed spoofing attacks ŝ = [ŝ 1 , . . ., ŝK ] ⊤ is stated as where ||x|| 2 M = x ⊤ Mx is a quadratic norm.Note that in (18), s p and s b are subvectors of the optimization variable s defined as s p = [s p [1], . . ., s p [K ]] and 138992 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
The minimization in (18) consists of four objective functions to facilitate the joint estimation of authentic PVT solutions and spoofing attacks against position and time.The first summation term is derived from the measurement equation ( 15) and aims to fit the PVT solutions and attack estimations to the observed vector.The second term aims at producing PVT estimates that abide by the random walk model in (5).The third and fourth terms are penalty terms, each representing the higher-order domain where sparse spikes are likely to be displayed in position and clock sequences, respectively.Specifically, matrices It is worth emphasizing that the derivatives are taken from the receiver's position and clock bias, rather than velocity and drift.For such penalizing formulation intends to cover both consistent and non-consistent spoofing cases since the attacks guarantee the presence of inflicting signals on position or clock bias, while only the former type of attack manifests a corresponding attack in velocity or drift.Overall, the technique is anticipated to protect against TSAs, position attacks, and joint spoofing.
The minimization in ( 18) is a quadratic program (QP).The problem can be written in MATLAB and solved by off-theshelf convex programming solvers such as CVX [41] efficiently and to global optimality.However, custom numerical optimization methods can also run on the receiver's lowcomputing power CPU and memory [42].The next subsection describes a detailed selection of tuning parameters.
The algorithm provides protection in traditional GPS receivers against any spoofing behavior defined in Section III.Thus, regardless of spoofing detection techniques on the baseband layer, the algorithm solely relies on measurements at the navigation layer.With an overpowered attack in which the tracking loop is locked onto a counterfeit peak, the algorithm will still read defective measurements, and alarm the receiver.While signals are cleansed within the current window of time, it is recommended the receiver initiate signal re-acquisition for the upcoming window.It can easily be seen that matrix D formulates successive differences of the entries in vectors s p and s b yielding sequences of first-order derivatives of the respective vectors (scaled by t).Such matrix is called the total variation matrix of the discrete-time signal [43].Furthermore, D m formulates the m-th order derivative.Note that once a signal is sparse in the m-th order derivative, then it is sparse in any order greater than m.Matrices D p and D b can correspond to any derivative order, though the third order appears to be adequate in practice.
Moving on to the weights λ p and λ b , these values determine the sensitivities of minimization performance.The proper choice of λ p and λ b affects the level of sparsity in the designated derivative order, hence λ p and λ b must be tuned for each receiver.Furthermore, a balanced selection of the pair of weights renders an emphasis on position or clock states.
The state noise covariance matrix Q k is generally formulated as where every component in Q k is determined by the process uncertainties in position (S p ) and clock dynamics given by S g = 2π 2 h 2 and S f = h 0 2 .Position related covariance matrices are formed as follows: while covariances related to clock dynamics are set to The user has the freedom to select the values for S p , h 0 , and h 2 ; the latter two depend on the accuracy of clock.For the numerical simulations performed in Section VI, the values h 0 = 2 × 10 −21 and h 2 = 3 × 10 −24 are selected [44].Different values of S p may be appropriate depending on whether the receiver is static or exhibits lowdynamic movement.The particular selections for S p are given in Section VI.Furthermore, since receivers acquire GPS raw measurements in 1-sec intervals, t by default is set to 1 second.

V. EVALUATION METHODOLOGY
This section describes the experiments aiming at assessing the performance of joint PVT estimation and spoofing rejection presented in Section IV-C.The receiver testbeds for real-time GPS data collection along with navigation data post-processing methods are discussed in Sections V-A and V-B.Then, Sections V-C and V-D present the spoofing data generating mechanism and the evaluation metric to assess the performance of the antispoofing technique.

A. GNSS RECORDING TESTBEDS
This work relies on a set of testbeds for capturing raw GNSS measurements at the University of Texas at San Antonio (UTSA).The testbeds consist of two receivers representative of research and commercial hardware.
The first testbed exploits the real-time LabView based L1 single-frequency GPS SDR receiver [33] in which the end-user has the capability of performing controlled environment testing by tuning receiver parameters.Due to the complete set of internal processing modules, once TEXBAT binary files are replayed on the SDR, the raw measurements are automatically provided.Furthermore, given that the SDR receiver operates in noiseless offline mode, once the recorded signal replays over such test bed, similar signal quality is ensured, thus providing precise GPS measurements i.e.C/N0, pseudorange, and range rate.
The second testbed utilizes a Huawei Android-based mobile receiver with an embedded GPS chipset.The Huawei receiver requires the Android application GNSSLogger [45] to extract raw measurements.TEXBAT signals are transmitted over-the-air (OTA) by an NI PXIe-1075 Chassis with a PXIe 5673 RF Signal Generator via PCIe interface and a Vert 900 antenna and captured by the Huawei receiver.As opposed to the SDR test bed previously described, commercial receiver signal collection involves OTA transmission in the controlled yet noisy lab setting.Therefore, the processed measurements are degraded.Fig. 3 shows the trends in the average from the top four C/N0 values in TEXBAT CS data.As noted in this section, SDR measurements maintained higher and more stable signal power than those of Huawei.

B. MATLAB POST-PROCESSING
With GPS measurements replayed and captured by the two testbeds, the PVT estimation is solved in MATLAB as a postprocessing task.Both the SDR receiver and the Huawei P30 Pro output a sequence of raw measurements over K epochs.Then, the required observables such as pseudorange, Doppler rate, GPS seconds, and signal uncertainty are given as inputs to WLS and EKF to produce PVT estimates in MATLAB.In addition to the traditional techniques of WLS and EKF, the same inputs are given to the developed optimization problem, solved with the convex optimization modeler CVX [41], in which the algorithm runs in a single snapshot.The tuning parameters listed in Section IV-D are set based on the testing receivers.

C. EVALUATION ENVIRONMENTS
To evaluate the performance of the developed optimization model against joint attacks, we perform experiments in two environments, namely, synthetic controlled simulations and practical experiments with real spoofing data.
The synthetic simulations entail the generation of various types of joint attacks.To assess the performance of the novel technique, we consider an unimpaired stationary receiver scenario of the TEXBAT data set-clean static (CS) scenario-as the ground truth.Under Scenarios 1 through 7, the CS signals are initially replayed over the SDR receiver, and the resulting raw measurements are processed via WLS to acquire PVT solutions.Scenario 8 entails the mobile Huawei receiver under low-dynamic movement recorded at a UTSA parking lot.Then, the designed synthetic joint attacks are appended on the raw clean measurement pair as per (12).Specifics of the eight position-targeting and joint synthetic attack scenarios are listed in Table 1.Each scenario is designed to test the resilience of the technique against various shapes and magnitudes of the attack, as well as the presence or absence of consistency.
The practical experiment, on the other hand, uses TEXBAT provided authentically spoofed data, specifically, the static position push scenario DS4.The z-domain PVT estimates deviate about 600m off from the receiver's initial stance in that scenario, hence the device is fully captured by the spoofer.Because at least one of the satellites was still locked to the authentic signal, the x-coordinate and clock bias estimates exhibit minor offsets, yet not significant enough to be considered as spoofing attack.Our SDR receiver replayed the DS4 binary-coded file and observed analogous position and bias estimates to the TEXBAT report [10] shown in Fig. 2. Scenario 9 features the DS4 signal replayed over the SDR receiver, and Scenario 10 entails the DS4 signal transmitted over-the-air and received by the Huawei receiver.
138994 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.The simulations and experiments induce spoofing attacks and contain smaller deviations on more than one position and time domain.The results presented in Section VI indicate that our proposed technique is effective in not only capturing and rejecting intended malicious attacks but also in reducing unwanted modifications.

D. EVALUATION METRIC
The performance of the developed optimization model is assessed through the estimation error formulated as the difference between corrected estimates and the ground truth.The root mean square error (RMSE) for the estimated sequence of each domain is calculated for all scenarios.With K defining the total length of observation period, and i = 1, 2, . . ., 8 indexing the 8 states in the PVT vector x k , the RMSE is defined as

VI. SIMULATION RESULTS
Insights in the tuning parameter selection for the experiments in this paper are presented in Section VI-A.The results for synthetic attacks are presented in Section VI-B, followed by the ones with replayed attacks from TEXBAT in Section VI-C.The resulting RMSE values from the output of the optimization model for each domain of the PVT solution and each scenario are listed in Table 2.The position uncertainty parameter S p is set to 0.5 meters for static scenarios, which include synthetic cases 1 through 7 and real spoofing scenarios 9 and 10.For the low-dynamic scenario 8, S p is set to 10 meters.

A. TUNING PARAMETER SELECTION
For all test scenarios, the tuning parameters are carefully selected to promote the best performance.As mentioned in Section V-D, the chief analysis metric is RMSE value.For any given parameter choices, the RMSE value in the PVT domains of interest must be less than those of spoofed data.Though, the effort is made to reduce the RMSE value to be close to zero.A complementary metric utilized in the performance analysis process was the numerical value of penalty term in Eq. (18).Based on the λ weights and total variation metric selection, the optimization performance is reflected on the sum of arguments making up the objective function.However, due to the intrinsically large pseudorange values, the sum of four terms inadvertently becomes a larger number.For that matter, the sum of the third and fourth terms, that is, the penalization terms, can offer a clue to the performance of the spoofing mitigation process.With successful spoofing signal capture, the penalization terms (l-1 norms) are close to zero as well.Thus, with the aforementioned metrics, tuning parameters are determined through systematic experiment procedures.

1) SCENARIO 1: FIRST ORDER POSITION NON-CONSISTENT ATTACKS (D
The first synthetic scenario fabricates the altered spatial estimates on x, y, and z domains with the shape of first order attack in which sparse spikes are observed in velocity.While nominal conditions are maintained for the initial 100 seconds, the spoofer dominates the receiver's behavior in the period between 100 and 300 seconds.Note that the amplitude modifications are different for each domain.As reported in [46], achieving at least 600 m deviation from the initial location in the z domain meets the requirement of a fully deceived receiver position.As this simulation features a nonconsistent attack, the receiver velocities remain the same as in the clean static scenario.Fig. 4a demonstrates the clean and attacked outputs of WLS as well as the output of the novel optimization model, resulting in RMSE of 4.52 m on z coordinate.Additionally, zoomed-in box plots show that estimation errors are smaller than 10 m.Fig. 4b shows the estimated attack s z and its first derivative ṡz as we expect to observe sparsity on the velocity domain.2) SCENARIO 2: SECOND ORDER POSITION NON-CONSISTENT ATTACKS (D P = D 2 , λ P = 0.01) The second synthetic scenario also alters the position coordinate values in a non-consistent manner, but with second order shapes, in which sparsity occurs in the acceleration domain.Accordingly, we select second degree total variation matrix enhancing the sensitivity in searching for spikes in second derivative of the position.After 100 seconds of nominal conditions, the spoofer carries off the receiver from 100 to 300 seconds at a constant rate.Then, once the intended amplitudes are reached on every domain, spatial estimations maintain the steady state until the recording terminates.Fig. 5a depicts the clean, attacked, and corrected position solution.The algorithm significantly lowers the RMSE values as opposed to those of WLS.Specifically, the RMSE value in z-domain, which has the highest deviation magnitude, is 2.36 m.The applied attack on z-domain, s z , is depicted in Fig. 5b.The third synthetic scenario targets position solutions with amplitudes of 100m, 300m, and 600m in x, y, and z coordinates, respectively.The spoofer deceives the receiver through the end of the recording.Compared to Scenarios 1 and 2, the third order attack applied here displays smoother and even more gradual change, in which sparse spikes are most likely to be observed in the third derivative of position, that is, the jerk domain.The total deviation matrix is D 3 .Furthermore, with respect to tuning parameters, since the attacking signal is smoother than those of the previous two simulations, λ p and λ b are set to higher values.Fig. 6a depicts that corrected estimations are close to the CS scenario with a fairly small RMSE value of 3.8 m on z-domain.Fig. 6b depicts the estimated attack s z .

4) SCENARIO 4 SECOND ORDER SINGLE POSITION AND TSA JOINT CONSISTENT ATTACK (D
Unlike the previous three scenarios, the fourth synthetic spoofing simulation deals with a joint attack in which the zcoordinate and bias estimates deviate from the initial value simultaneously with respective magnitudes of 600 m.This scenario assumes that the spoofer not only fully captures both spatial and timing solutions on the receiver, but also impacts the velocity and drift accordingly.Fig. 7 depicts the outcome of synthetic spoofing along with mitigated solutions 138996 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply. in four domains, namely, position and velocity with respect to z-coordinate, clock bias, and drift.The total variation matrices are selected as D 3 since a second order attack exhibits sparsity in the jerk domain.With RMSE values for the z-coordinate and bias equal to 6.03 m and 12.9 m, respectively, this experiment demonstrates that the optimization model successfully mitigates joint consistent attacks.The fifth and sixth scenarios simulate joint attacks on all three spatial coordinates and a TSA with various orders of spoofing.A smooth third order attack with 600 m carry-off is applied on the z-domain, which guarantees spoofing deception of the receiver.The bias is altered by a 50-m step attack.Scenario 5 introduces non-consistent attacks.Scenario 6 introduces the same attacks as Scenario 5, but the attacks maintain consistency.The optimization outputs for each scenario are shown in Fig. 8a and 8b respectively.Both of D p and D b are selected as the third-order total variation matrix, or equivalently D 3 , along with suitable weights.The resulting in RMSE values for z and bias domains are 3.77 m and 12.90 m, respectively.It is concluded from the observed estimation outputs in Scenarios 5 and 6 that the proposed algorithm successfully recognizes and mitigates the anomalies in any PVT state regardless of spoofing profile and consistency.The reported RMSE values for scenarios 5 and 6 in Table 2 support the argument by indicating similar numbers across domains.

6) SCENARIO 7: INTEGRATED SECOND ORDER ALL POSITION CONSISTENT ATTACK WITH MULTIPATH EFFECT
Scenario 7 tests the performance of the optimization model when the receiver experiences spoofing and multipath simultaneously [47].Multipath is modeled here similarly to existing literature.Specifically, a 50-m first-order shape deviation is inserted in 4 out of 7 visible satellites between 120 to 140 seconds of recording.The spoofing attack amounts to 100, 300, and 600 meters on x, y, and z domains, respectively.Unlike the spoofing attack, multipath manifests as unintentional discrepancies in only certain channels of measurement data.The multipath results in anomalous behaviors of PVT states; thus, the proposed algorithm can detect and mitigate the combined effect owing to the sparse estimation.Fig. 9 depicts the aforementioned effects in position and bias domains.
Because the spoofing follows a second order profile, the position related TV matrix is set to D 2 with identical λ p value as in Scenario 2. Fig. 9 reveals that the optimization mitigates the presence of multipath and accurately estimates the true position and clock offset of the receiver.Further, as shown in Table 2, the RMSE values corresponding to all position coordinates are noticeably reduced.The measurements are processed via WLS, KF, and our optimization problem.Fig. 10a depicts the resulting position and bias solutions.The plot indicates approximately 40 meters of user movement in position x-domain.Without spoofing, the optimization results accurately follow the WLS estimates, are considered as the ground truth.Fig. 10b) depicts a 100-meter second order synthetic attack applied to position x in a non-consistent manner.Considering the shape of the attack, the position total variation matrix is set to D 3 with sensitivity parameter λ p = 0.35.The computed RMSE value of spoofed x domain features 47 m as opposed to 8.17 m with mitigation.As the spoofer overtakes the transmitted signals from 150 to approximately 280 seconds, the z-coordinate experiences a purposeful carry-off, and the remaining three states undergo perturbations as well.Subsequently, every coordinate reaches a steady state regardless of the magnitude.The optimization model successfully detects and rejects the attack and perturbations attaining minimal estimation errors (cf.Table 2).The total variation matrices are set to D 3 with comparatively higher weights than synthetic scenarios.receiver.As opposed to the previous scenario with the SDR receiver, the carry-off stages contain higher noise which results in visibly high peaks in the attacked position estimates.The total variation matrix is selected as D 3 .The fourth panel in Fig. 12 depicts two graphs, namely, the difference between the spoofed clock bias and the clean version (ground truth), as well as the difference between the corrected clock bias and the clean version.Due to the internal clock offset in the Huawei receiver, the reported raw clock bias sequences had decreasing trends reaching 10 6 meters magnitude.To illustrate the effect of spoofing in the clock domain more clearly, we manually reconstructed the bias, bu,recon [k], by removing the existing trend due to the average Overall, the results in Table 2 and Fig. 12 indicate that the optimization model detects and mitigates the attacks with estimation errors ranging from very small to acceptable, depending on the domain.With regard to the two scenarios featuring the authentic DS4 signals, Fig. 11 and 12 manifest that over-the-air captured signals by the Huawei receiver contain more noise or uncertainty than the ones received by the SDR.As a consequence, even the CS signal contains relatively large transients leading to higher RMSE values after spoofing rejection.With proper tuning of λ p , λ b , and TV matrix D n , the uncertainty was substantially overcome resulting in comparatively small RMSE values.

VII. LIMITATIONS ON THE PRESENT WORK
It is worth mentioning that the proposed spoofing countermeasure must operate in controlled environments.
Introduced spoofing detection and mitigation technique exploits the specific characteristics discovered in abnormal signal behaviors; thus, the algorithm must take initiative prior the spoofing attack to fully capture and recognize unexpected changes during the receiver run-time.Further, as mentioned earlier in Section II-B, receiver position and velocity references are computed via WLS and provided to a receiver.It must be noted that the present work fixates the reference estimates throughout the processing.Thus, the technique may induce high estimation errors for a longer duration of GPS recording.
Additionally, we stress that the applied doesn't necessarily cure the impaired signal nor reacquire the hijacked correlation peak in the baseband domain.Instead, by leveraging the signal patterns that manifest spoofing in the measurements for a certain period of time (tens of seconds), the proposed work can detect and mitigate the abnormality, and provide reasonably accurate PVT estimations.

VIII. CONCLUSION AND FUTURE WORK
This work presents feasible signal-level joint spoofing that purposefully inflicts one or more PVT domains with different combinations of attack profiles, carry-out magnitudes, and duration.Consistent and non-consistent attacks based on measurement integrity are defined.A suitable PVT dynamical model that includes joint attacks against position and time is developed.An optimization problem that leverages the aforementioned dynamical model and the sparsity characteristics of position and time attacks in higher-order derivatives is formulated.The optimization can be easily solved to yield corrected PVT estimates as well as reconstruction of the attack.To assess the new technique's robustness against both synthetic and authentic spoofing attacks, the TEXBAT dataset is utilized in two testbeds.Eventually, the PVT RMSE values are noticeably reduced.
In light of Section VII, future work will address the limitations of the proposed technique's applicability.While the approach works well with the considered spoofing models, available datasets, spoofing intensity, and durations, there are other scenarios and factors that may be limiting.In particular, additional comprehensive study on reference position is required.While described PV reference choice method in Section II-B worked well with given data sets, coarse reference based linearization may not be sensitive enough to capture low intensity attacks.Moreover, future work intends to include algorithm response against existing initial clock estimation error, which leads to imprecise position acquisition.Also if the spoofer is persistent for longer durations and applies variable dynamic spoofing patterns not related to actual measurements, the proposed technique 139000 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
will be applicable only for a limited time frame when measurement consistency can be still exploited.Another aspect future work will investigate is the signal degradation at the occurrence of spoofing injection such as sudden C/N0 drops or spikes in urban canyon.Though the proposed algorithm generally takes into consideration the signal quality through the noise covariance matrix, extensive study in how the proposed algorithm responds against tampered C/N0 must be followed.
Furthermore, since the proposed work processes a window of K measurements in a batch manner, real-time implementation is worth pursuing.Finally, we plan to expand the developed technique to highly dynamic receivers that involve larger movements and incorporate readings from additional sensors.

APPENDIX MEASUREMENT EQUATIONS LINEARIZATION PROCEDURE
The pseudorange and pseudorange rate measurement equations are written as follows: where we define the nonlinear functions As the nonlinearities in (21) complicate the formulation of a convex optimization problem, equations (22) are linearized with respect to user position and velocity p u [k] and v u [k] around reference points given by p u,ref [k] and v u,ref [k] as follows: The linearization is achieved via Taylor series expansion while acknowledging only up to the first degree power.Application of multivariable differentiation rules yields the following expressions  Equations ( 24) are rearranged so that the left-hand sides include only reference, satellite-provided, or measured quantities, while the right-hand sides depend only on the unknown user PVT quantities.Specifically, the left-hand sides are defined as With the aforementioned definitions, the linearized measurement equations take the following form: It is worth emphasizing that the right-hand sides of ( 26

FIGURE 2 .
FIGURE 2. TEXBAT scenarios CS and DS4 recorded by SDR receiver.The position and clock bias estimates are computed with recursive WLS.
D p and D b are constructed such that sequences of derivatives of the entries in s p and s b is formulated by D p s p and D b s b , respectively.The penalization of their ℓ 1 norms has the effect of recovering sparse vectors D p s p and D b s b , such as the sparse sequence in Fig. 1.The construction of matrices D p and D b is elaborated in the next subsection.

D
. TUNING PARAMETERS This section discusses the selection of three sets of parameters for the optimization in (18): 1) Matrices D p and D b 2) Weights λ p and λ b of the penalty terms 3) State noise covariance matrix Q k The matrices D p and D b are chosen to have the following form D m = D × D × . . .× D

FIGURE 3 .
FIGURE 3. Average of top 4 C/N0 of TEXBAT CS data recorded on two experimental test beds: SDR (blue) and Huawei (red).The tow on x-axis indicates GPS Time of Week.

FIGURE 4 .
FIGURE 4. Response to synthetic attack scenario 1: (a) comparison of position estimates; and (b) z-domain estimated attack and its first derivative.

FIGURE 5 .
FIGURE 5. Response to synthetic attack scenario 2: (a) comparison of position estimates; and (b) z-domain estimated attack and its second derivative.

FIGURE 6 .
FIGURE 6. Response to synthetic attack scenario 3: (a) comparison of position estimates; and (b) z-domain estimated attack and its third derivative.

FIGURE 7 .
FIGURE 7. Response to synthetic attack scenario 4. Plots indicate estimates of z-coordinate position and velocity, clock bias, and drift.

7) SCENARIO 8 :
SECOND ORDER SINGLE POSITION NON-CONSISTENT ATTACK IN LOW-DYNAMIC RECEIVER (D P = D 3 , λ P = 0.35, λ B = 10) Scenario 8 features the Huawei P30 mobile receiver and real-time data acquisition.The aim is to assess the performance of the spoofing countermeasure technique under low-dynamic behavior.For the data collected in an open-sky environment with the handheld device, a total of 220 seconds of GPS recording indicated about 37 dB-Hz C/N0, in which the position uncertainty was fairly low.

FIGURE 9 .
FIGURE 9. Response to synthetic attack scenario 7. Plots depict position and bias estimates.
C. SCENARIOS WITH REAL SPOOFING 1) SCENARIO 9: TEXBAT DS4 ON SDR (D P = D 3 , D B = D 3 , λ P = 0.001, λ B = 10) Scenario 9 amounts to an experiment with authentically spoofed signals, where TEXBAT DS4 signals are replayed and captured via the SDR receiver.Fig. 11 depicts the clean, attacked, and corrected estimates on each domain.

2 )
SCENARIO 10: TEXBAT DS4 ON HUAWEI PRO30 (D P = D 3 , D B = D 3 , λ P = 0.001, λ B = 10) This experiment features the TEXBAT DS4 signal transmitted over-the-air and captured by the Huawei Pro30 138998 VOLUME 11, 2023Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

FIGURE 10 .
FIGURE 10.Response to synthetic attack scenario 8.The real-time recording is acquired with the Huawei P30 receiver under low-dynamic movement.(a) Position and clock bias without spoofing; and (b) position and clock bias solutions in cases of spoofed and mitigated estimation.

FIGURE 11 .
FIGURE 11.Response to SDR replayed TEXBAT DS4 signal.Plots depict the position and clock bias estimates.

FIGURE 12 .
FIGURE 12. Response to Huawei captured TEXBAT DS4 signal.The top plots depict the position estimates.The bottom plot depicts the difference between the spoofed clock bias and its clean version, as well as the difference between the corrected bias and the clean version.

TABLE 1 .
Applied synthetic joint attack scenarios.

Table 2
lists the RMSE values resulting from the optimization model and from WLS for all 10 scenarios.Although WLS is not intended to mitigate spoofing, it is meaningful to display the results side-by-side with the outputs of the optimization model.The WLS results show the extent of unmitigated spoofing effects on the PVT solutions.Overall, the table indicates that the optimization model succeeds in mitigating the spoofing attacks and returns relatively accurate PVT states.For example, the average RMSE across all scenarios where the z-coordinate is attacked by 600 m (all scenarios except scenario 8) is 11.57m for the optimization model as opposed to 353.35 m for WLS.