Statistical test for GNSS spoofing attack detection by using multiple receivers on a rigid body

Global navigation satellite systems (GNSS) are being the target of various jamming, spoofing, and meaconing attacks. This paper proposes a new statistical test for the presence of multiple spoofers based on range measurements observed by a plurality of receivers located on a rigid body platform. The relative positions of the receivers are known, but the location and orientation of the platform are unknown. The test is based on the generalized likelihood ratio test (GLRT) paradigm and essentially performs a consistency check between the set of observed range measurements and known information about the satellite topology and the geometry of the receiver constellation. Optimal spoofing locations and optimal artificial time delays (as induced by the spoofers) are also determined.Exact evaluation of the GLRT requires the maximum-likelihood estimates of all parameters, which proves difficult. Instead, approximations based on iterative algorithms and the squared-range least squares algorithm are derived. The accuracy of these approximations is benchmarked against Cramér-Rao lower bounds.Numerical examples demonstrate the effectiveness of the proposed algorithm and show that increasing the number of GNSS receivers makes the attack easier to detect. We also show that using multiple GNSS receivers limits the availability of optimal attack positions.


Introduction
Global navigation satellite system (GNSS) technology provides real-time positioning for various civil and military applications. Due to the low received power (about − 160 dBW), GNSS is highly susceptible to intentional and unintentional interference. In addition, the signal and modulation formats of civilian GNSS are publicly available. For these reasons, a wide range of attacks on GNSS are viable.
Attacks on GNSS are categorized into jamming, meaconing (replay), and spoofing [1]. In a jamming attack, the adversary overshadows the received GNSS signals with a higher power noise-like signal to make the receiver unable

Techniques using auto-correlation characteristics
The first group of techniques uses the received signal auto-correlation function characteristics to detect the presence of a spoofer. These methods detect the spoofer when it tries to take over the tracking loop of the GNSS receiver. There are in turn three main approaches: • Methods that use the signal quality monitoring mechanism to measure the distortion in the auto-correlation function of the received signal and detect the spoofer [3][4][5][6][7][8][9][10][11][12]. The paper [12] studies the applicability of multi-path mitigation techniques for anti-spoofing purposes. A different approach is taken in [13], where the cross-correlation among the GNSS signal of multiple GNSS receivers is used as input to a machine-learning classifier (a support vector machine specifically) to detect the spoofing attack. • Considering that the spoofer must use higher power than the original GNSS signals to make the receiver lock on the counterfeit signals, monitoring the received power level of the auto-correlation function can be used to detect the spoofing attack [14][15][16][17][18][19]. The dynamic process of taking over the victim tracking loop is used in [19] to detect a spoofing attack. As a new metric, [18] proposes to measure the variance of the signal quality to detect when the spoofer starts taking over the tracking loop of the victim. • One may combine both power and distortion monitoring of the auto-correlation function to detect an attack [20,21].
These approaches have a number of difficulties. First, multi-path fading can cause distortion and power level fluctuations in the auto-correlation function of the received GNSS signal similar to those caused by the spoofer. In addition, the spoofer can intelligently change its power level to mimic the power fluctuation of multipath fading. Second, these techniques require access to the received raw GNSS signals.
Considering a common source of spoofing signals, the work [29] uses the time difference of arrival of the GNSS signal to detect spoofing. The work in [30] uses array processing along with multi-path detection algorithms to estimate the direction of arrival to detect the presence of one spoofer. A statistical test is used in [31] to estimate the power and angle-ofarrival of the GNSS signal in order to detect spoofing.
For example, carrier phase differences can be used to estimate this direction of arrival [22,25,27]. • One may exploit that counterfeit signals from a single spoofing source are spatially correlated and measure the spatial correlation between the received signals from different satellites [32][33][34]. The work [34] uses correlation at multiple receivers to separately classify authentic and spoofing signals; then, double differences of carrier phase measurements are used to detect the spoofer. • Using rotating antennas to measure the spatial correlation of the received GNSS signals [35,36]. The paper [35] uses rotating antennas to measure the correlation between the carrier phases, while [36] uses a single rotating antenna to perform spatial power measurements.
The techniques in [22-28, 35, 36] require modification of the GNSS receiver hardware. In addition, the methods [22-28, 35, 36] are based on the property that the signals coming from a single-antenna spoofer are correlated. Attacks using multiple spoofers can result in spatially uncorrelated spoofing signals, in which case these methods can fail.

Methods that use multiple GNSS receivers
The last category of spoofing mitigation methods uses multiple GNSS receivers to detect the presence of the spoofer. These works can be further grouped as follows: • Using multiple GNSS receivers to detect the presence of a spoofer, exploiting the fact that all GNSS receivers show the same position while being spoofed [37][38][39][40][41][42]. • Range measurements may be used directly to detect the spoofing attack [43][44][45][46][47]. More specifically, [43] considers multiple vehicles where each of them is equipped with a GNSS receiver and a range measurement sensor. The range measurements and GNSS positions are fused to detect the presence of one spoofer while assuming that only one of the vehicles is subject to a spoofing attack. In [44], the range measurements from multiple GNSS receivers are used to detect a spoofer, assuming that all range measurements from a spoofer are the same. Considering the effect of clock knowledge, the authors in [45] develop a technique to detect a spoofer using range measurements of multiple GNSS receivers. The authors of [46] propose to use the time difference of arrival of the GNSS signals derived from pseudorange measurements to detect the presence of one spoofer. The technique is based on the fact that signals coming from a spoofer have similar time differences of arrival. The work [47] studies detection of one spoofer using differential pseudorange and (2020) 2020:8 Page 3 of 16 carrier phase measurements by a double antenna receiver. The effect of synchronization between the measurements on spoofing detection is investigated.
All above methods except [41] consider the detection of a single-antenna spoofer. The paper [41] develops a technique for multi-antenna spoofer detection and presents implementation results.

Contributions
In this paper, we propose a signal processing approach that uses range measurements from multiple satellites to multiple GNSS receivers to detect the presence of multiple single-antenna spoofers. Each spoofer emulates one specific satellite. The GNSS receivers are assumed to be fixed on a rigid body platform (with known relative positions), and the range measurements gathered by all receivers are processed jointly to detect the presence of a spoofing attack. The technique developed here can work as a second layer of security to further strengthen methods that rely on properties of the auto-correlation function of the GNSS signal to detect spoofing attack (Section 1.1.1). We make the following specific contributions: • We cast the spoofing detection problem as a statistical hypothesis test, specifically a generalized likelihood ratio test (GLRT), based on a set of observed range measurements. • We use squared range-least squares (SR-LS) approach [48] to approximately find maximum-likelihood estimates of the parameters under two hypotheses corresponding to spoofing and no spoofing, respectively. We also calculate the Cramér-Rao lower bound (CRLB) for the estimated parameters under both these hypotheses and compare the empirical results with these bounds. • We determine optimal locations of the spoofers (from the adversary's perspective) that best counteract the proposed defense mechanism.
We present and evaluate all methodology using a twodimensional model of the world, leaving the extension to three dimensions-which incurs some nontrivial technicalities-to future work.

Notation
Uppercase and lowercase bold-faced letters are used to denote matrices and column vectors, respectively. The superscripts (·) T , (·) * , (·) H , and (·) † denote the transpose, conjugate, Hermitian, and Moore-Penrose pseudo-inverse operators, respectively. I N×N denotes an N by N identity matrix, diag(a) denotes a diagonal whose diagonal elements are the elements of the vector a, 0 is the all-zero vector, · is the Frobenius norm, and | · | represents the absolute value of a scalar. For a symmetric matrix A n×n

System model and problem description
We consider a two-dimensional scenario with N GNSS receivers and I satellites. The GNSS receivers are mounted on a rigid platform with fixed (and known) mutual distances. We assume that the clocks of the GNSS receivers are synchronized with those of the satellites. Spoofers may be present, and if so, each spoofer emulates one specific satellite. The emulation is done by receiving the GNSS signal of a satellite, modifying the ephemeris data, adding an artificial time delay, and re-transmitting the signal [49]. We assume that the spoofers use higher transmit power than the satellites, so that the GNSS receivers lock on the spoofing signals instead of the legitimate satellite signals. When the receivers are synchronized, they can operate as a virtual array to perform angle-of-arrival estimation [50]. This can be used for spoofing detection where emulated GNSS signals arrive from a specific direction. In contrast, our approach can detect a spoofing attack where multiple spoofers emulate GNSS signals originating from geographically different locations. In addition, performing array processing requires accessing to the baseband signal of the GNSS receivers, which requires hardware modifications and precludes the use of commercial off-the-shelf (COTS) receivers. Let p n be the position of the nth GNSS receiver. Since the GNSS receivers are fixed on a rigid body, we can parameterize their positions in terms of their locations relative to the platform, p n 0 , a translation vector (that determines the position of the platform), b 0 , and a rotation matrix (that determines the orientation of the platform), T, according to The rotation matrix, T, is parameterized in terms of a rotation angle θ as and we write b 0 = b 01 b 02 . Hence, {p n 0 } are known a priori, while b 0 and T are not. Figure 1 illustrates the relation between the variables {p n }, b 0 , T, and {p n 0 } in (1). Each GNSS receiver takes M range measurements to each satellite. To detect the possible presence of spoofers, we apply a binary hypothesis test to the so-obtained range measurements. In regular operation (no spoofing), the delay of a signal from a satellite to a receiver is entirely due to the physical distance. In contrast, in the presence of spoofing, the spoofer is generally located at a different position than the satellite, and it can add an artificial, a priori unknown, time delay to the GNSS signal to simulate a different physical distance. Therefore, the two hypotheses, regular operation (H 0 ) and spoofing (H 1 ) cases, can be formulated: where 1 ≤ n ≤ N, 1 ≤ i ≤ I, and 1 ≤ m ≤ M is the number of GNSS signal samples captured by each GNSS receiver. We assume that the rigid body does not move during the sampling. Here, n n,i,m are measurement noise samples that we model as identically distributed zeromean Gaussian random variables with variance σ 2 , and we assume that these noise samples are independent among n, i, and m. Also, s i is the true position of the ith satellite, s i is the forged position of the ith satellite (when a spoofer is present), and τ i is the artificial time delay caused by the ith spoofer. Note that the artificial time delay introduced by the spoofer is perceived by the GNSS receivers as an additional propagation distance. We assume that the noise power, σ 2 is known by the GNSS receivers (Table 1).
Under H 1 , the translation vector b 0 is the same for all measurements, and hence we can absorb b 0 into the unknown satellite positions s i by rewriting (4) as where To further simplify (5), we can factor out the rotation matrix T, exploiting that T T T = I, to rewrite (5) as where We first average the M measurements taken by each GNSS receiver for each satellite, to obtain: Since under the Gaussian-noise assumption, {r n,i } are sufficient statistics for localization [51], we have the equivalent hypothesis test based on averaged measurements: where n n,i = 1 M M m=1 n n,i,m is averaged noise with zero mean and variance σ 2 M . Averaging over the samples decreases σ 2 M and improves the signal-to-noise ratio. Under H 0 , the rotation matrix T and the translation vector b 0 are unknown, while under hypothesis H 1 , the vectors s 1 , ..., s I , and τ 1 , ..., τ I are unknown. In the next section, we devise a statistical test that can discriminate between H 0 and H 1 .

Generalized likelihood ratio test for spoofer detection
We formulate a generalized likelihood ratio test (GLRT) to detect the presence of a spoofing attack. The GLRT compares the likelihoods of the data under H 0 and H 1 , with all unknown parameters replaced by their maximumlikelihood estimates. The GLRT approach is appropriate for cases when there are unknown parameters under each hypothesis, and no prior knowledge of these parameters is available.
Since the measurements in (7) and (8) follow a normal distribution, and since all averaged noise samples are independent, the PDFs of the measurements under H 0 and H 1 are and p r; s 1 , ..., s I , τ 1 , ..., τ I , H 1 = 1 The GLRT can be written as and Taking the logarithm of (11) and simplifying yields where As we see in (14), GLRT calculates the difference between the measurements and the best fit assuming that H 0 and H 1 , respectively, are true.

Approximate parameter estimation
Evaluation of the GLRT requires the maximum-likelihood estimates of the parameters in (12) and (13). However, finding these estimates is a non-tractable problem. As a remedy, we resort to using the SR-LS technique [48] as a building block in an iterative algorithm. SR-LS was originally devised as a source localization method that finds the best position estimate (in a least squares sense) based on squared range measurements. We are going to use estimates from SR-LS here as a substitute for maximumlikelihood estimates required in (12) and (13). In this context, it is important to stress that under a Gaussian assumption on the measurement errors, SR-LS is asymptotically equivalent (for large numbers of samples, which (2020) 2020:8 Page 6 of 16 in our setting is equivalent to small σ 2 /M) to maximumlikelihood. The range-squaring idea in [48] in ingenious and has also been used by other authors for joint position and orientation estimation of a rigid body from range measurements [52,53].

Estimation of parameters under H 0
Under H 0 , as an approximation to (12), we formulate the following problem, defined in terms of the squared range measurements: We then solve (15) using cyclic optimization, alternating between minimization with respect to b 0 and with respect to T.

Minimization of (15) with respect to b 0
Considering the rotation matrix T to be given, we first minimize (15) with respect to b 0 . Expanding the norm in (15) Using the fact that T T T = I and ordering (16) with respect to b 0 yieldŝ To write (17) as squared norm, we introduce an auxiliary variable α and a constraint: The optimization problem in (18) can be written in a compact form aŝ where and Note that A should have full column rank, i.e., A T A should be non-singular for the solution to exist. The optimization problem in (19) is non-convex. However, its global optimum solution can be found using the result in [54] as where λ is the unique solution of and φ is defined as The search interval for λ consists of the values for which the expression A T A + λD is positive definite. As a result, the search domain V will be Based on [54], the function φ (λ) is strictly decreasing over the domain V. Hence, we can use bisection to find the root of φ (λ).

Minimization of (15) with respect to T
To minimize (15) with respect to the rotation matrix T, for given b 0 , we again write, similarly to (17): Defining c n,i = p T n 0 Using the structure of the rotation matrix in (2) and cos θ − sin θ sin θ cos θ p n 01 p n 02 = 2 p n 01 a i1 + p n 02 a i2 cos θ + 2 p n 01 a i2 − p n 02 a i1 sin θ = α n,i cos θ + β n,i sin θ.
Using (29), estimating T becomes equivalent to estimating θ as Problem (30) is non-convex with respect to θ and may have multiple local minima. To find the global minimum, we calculate the derivative of (30) with respect to θ and find its roots. The derivative of the objective function in (30) is which can be simplified into n,i sin θ cos θ + α n,i β n,i 2cos 2 θ −1 − α n,i c n,i sin θ + β n,i c n,i cos θ.
Summation over the indexes n and i in (32) results in f (θ)=e 1 sin θ cos θ +e 2 2 cos 2 θ−1 +e 3 sin θ +e 4 cos θ, with The roots of f in (33) We square both sides and substitute sin 2 θ = 1 − cos 2 θ to obtain e 2 1 + 4e 2 2 cos 4 θ + 2 (e 1 e 3 + 2e 2 e 4 ) cos 3 θ + −e 2 1 − 4e 2 2 + e 2 4 + e 2 3 cos 2 θ which is now a trigonometric function with only cos θ terms (which are periodic with period 2π). Replacing x = cos θ in (36) results in The fourth order equation in (37) has at most four roots, which can be found using the Ferrari's method [55]. After solving (37), we can use θ = arccos x to find the values of θ corresponding to the local optimums of (30). Note that since cos θ in (36) is periodic with period 2π, we look for the solutions in the domain [0, 2π] when using θ = arccos x. The value of θ corresponding to the global optimum of (27), i.e.,θ , is then

Iterative algorithm for minimization of (15)
Minimizing (15) jointly with respect to b 0 and θ (or equivalently, T) is not tractable in closed form. Instead, we propose an iterative optimization approach in Algorithm 1. Since the problem is highly non-convex with multiple local optima, the point obtained as solution could depend on the initialization values (especially of θ). To tackle this, we repeat Algorithm 1 with different initial values of θ.
Algorithm 1 Iterative approach to minimization of (15) with respect to b 0 and T.

Go to 3; 8: end if
It is shown in [48] that the SR-LS results in the global optimum to the approximated problem based on squared range measurements. Hence, in the minimization with respect to b 0 , for given T, in each iteration inside of Algorithm 1, we find the global optimum. The same is true of θ, as its global solution is given by the root of (37). However, the overall problem is non-convex and the joint estimate of b 0 and T may be suboptimal.

Estimation of parameters under H 1
Next, to find approximate solutions to (13), we consider again a cyclic optimization algorithm that alternates between minimization with respect to s 1 , ..., s I and with respect to τ 1 , ..., τ I and where SR-LS is used as a building block in the first-mentioned step.

Minimization with respect to
where r 2 n,i = r n,i − τ i 2 . We expand the norm expression in (39) as ŝ 1 , ...,ŝ I = arg min and add the terms in (40) By completing the square, we can cast (43) aŝ where and

Minimization with respect to τ i
Next, we find the artificial time delays caused by each spoofer for given values of s 1 , ..., s I . The problem is We expand (47) The objective function in (48) is sum of positive and independent terms. Again, the problem decouples and we can find the minima with respect to τ 1 , ..., τ I separately for i = 1, ..., I. The corresponding problem for the ith variable, τ i , is and has the (s) solution

Iterative algorithm for approximate solution of (13)
The whole approach to estimate the parameters under H 1 is summarized in Algorithm 2. Since the estimation problem is non-linear and can have multiple local optimum points, we need to repeat Algorithm 2 for different initialization points of τ 1 , ..., τ I . Similarly to the case under H 1 , the minimization with respect to s 1 , ..., s I for given τ 1 , ..., τ I in Algorithm 2 returns the global optimum. Also, for given s 1 , ..., s I , (50) yields the globally optimal values of τ 1 , ..., τ I . However, (2020) 2020:8 Page 9 of 16 the joint estimates s 1 , ..., s I and τ 1 , ..., τ I delivered by Algorithm 2 may be suboptimal. A more granular initialization range for the parameter θ increases the chance of finding the global optimum rather than a local optimum.

Complexity analysis of the proposed algorithms
In this section, we provide the Big-O complexity for each step of Algorithms 1 and 2, respectively. Algorithm 1 is comprised of two parts where b 0 is first estimated using (23) and θ is estimated using (37), which is solved using Ferrari's method. Adding up the complexity of these two parts, the complexity of each iteration of Algorithm 1 is: where r = log 2 ε 0 ε , ε 0 is the search domain of the bisection algorithm in (26), and ε is the required accuracy of the bisection algorithm. The overall complexity of (51) is the number of tested initial values of θ 0 multiplied by the complexity of each iteration provided in (51).
Similarly, Algorithm 2 is comprised of two parts. First, the values of s 1 , ..., s I are estimated and then τ 1 , ..., τ I are computed. Adding the complexity of these two parts, the complexity of each iteration of Algorithm 1 is:

Optimal design of attack parameters
As we develop signal processing methodology to detect an attack by multiple spoofers, it is expected that the adversary devises the best strategy to counteract the spoofing detection technique. In this section, we develop analytical solutions to design the optimal spoofing parameters.
In this analysis, we assume that the adversary has perfect knowledge of the initial positions of the GNSS receiver nodes, p n 0 .
To minimize the chance of detection, the adversary needs to satisfy two types of constraints. First, the adversary needs to find spoofed locations, p sp n , such that the relative coordinates of the GNSS receivers is preserved after spoofing. To achieve this, the spoofer selects values of the spoofed translation vector, b sp 0 , and spoofed rotation matrix, T sp . Then, it uses the knowledge of p n 0 to calculate the values of the spoofed locations, p sp n , as Second, the locations of the spoofers and forged positions of the satellites need to chosen such that the range measurements at all GNSS receivers result in the target spoofed location derived in (53), and the spoofed formation of the GNSS receivers is preserved. To satisfy this, the simulated distance by a specific spoofer, caused by its physical distance to a GNSS receiver and artificial delay, needs to be equal to the distance of forged satellite, provided in the simulated GNSS signal, to the estimated position of that GNSS receiver. To implement the former in the best way, the adversary first picks values for the forged positions of the satellites, s sp i , and uses the values of p sp n derived in (53) to satisfy the following constrains: where p A i is the position of the ith spoofer, p sp n is the spoofed location of the nth GNSS receiver, s sp i is the forged position of the ith satellite, and τ i is the artificial delay produced by the ith spoofer. The better the spoofers satisfy the corresponding equations in (54), the closer the estimated parameters in (3) to that set by the adversary. Consequently, this leads to a better fit between the measurements and the data model. As we see in the developed detection test in (14), this helps reducing the left hand side of the detection test and not reaching the detection threshold.
Each spoofer chooses values for s sp i since these values need to be in a specific range, e.g., GEO satellites, and then solves the related equations in (54) to get the optimal values of p A i and τ i . To illustrate, the best position for the first attacker is illustrated in Fig. 2 for N = I = A = 2. As we see, satisfying the set of constraints in (54) results in the optimal position for the spoofer being the intersection of two circles where the center of each circle is p n for n = 1, 2. For N ≥ 3, there are three circles which may not have a common intersection point. Therefore, the adversary needs to find the best point for each spoofer which best satisfies (54) for all GNSS receivers. One way is to use a least squares fit, and calculate the derivatives of the objective in (55) with respect to p n x , p n y , and τ i . After algebraic simplifications, we get where c n,i = y i − p n y 2 − μ n,i − τ i 2 , Relations (58) imply an intertwined system of non-linear equations in (57). This system of equations can be solved using classic approaches such as the Newton method. To find the optimal values of p A i and τ i for each spoofer, we perform a coarse search over multiple initial guess points and choose the best one among them.
Nonetheless, in practice, it may not be possible for each spoofer to be in the optimal location at each moment. This is due to the fact that the GNSS receivers are installed on a moving platform and physical barriers as well as unknown velocity and speed can prevent the spoofers from being at the required positions. We quantify the spoofing detection performance with respect to the deviation from the optimal location of the spoofers in Section 6.

Simulation results
In this section, we present numerical examples to quantify the performance of the proposed spoofing detection mechanism. Unless otherwise stated, the parameters of the simulation setup are as in Table 2.
To initializeθ for Algorithm 1, we select a set of coarsely spaced values within the range [ 0, 2π]. We repeat Algorithm 2 for different initial values of τ 1 , ..., τ I and then   To demonstrate the accuracy of the parameter estimates used in lieu of the maximum-likelihood estimates required in (12)-(13), we first compare the empirical estimation variance to the Cramér-Rao lower bound (CRLB). The CRLBs under H 0 and H 1 are derived in Appendix A and Appendix B, respectively. Due to space constraints, we present these comparisons only for the parameters θ for H 0 and s 31 for H 1 . Figures 3 and 4 show the results. As we can see, parameter estimates get close to their respective CRLBs as the noise variances decrease.  To generate the simulation results and find the threshold of the GLRT, we proceed as follows. First, we run a case without spoofers and empirically calculate the probability of false alarm using the GLRT test in (14) for various values of γ th . Next, we simulate the presence of spoofers and calculate the probability of detection using the GLRT test in (14) for various values of γ th . Finally, we plot the so-obtained probabilities of detection and false alarm.
In the first simulation scenario, we investigate how the difference between the true and spoofed positions of the GNSS receivers affect the spoofing detection probability. We present the probability of detection versus false alarm in Figs. 5 and 6 for different values of the spoofed trans- The results show that as the adversary tries to mislead the victims more from their true locations, the chance to spot the spoofing attack increases.
In the next scenario, we quantify the performance of the proposed technique when the position of each spoofer deviates from the optimal designed positions. The probability of detection versus false alarm is shown in Fig. 7 when each spoofer is moved by 2p A i / p A i from the optimal position. As we see, the probability of detection increases as the spoofers fail to occupy their optimal locations.
Finally, we investigate the effect of the number of GNSS receivers on the performance of the proposed algorithm. We consider the same satellite formation with I = 3 and three spoofers, i.e., A = 3, while increasing the number of GNSS receivers from three to five. The performance of the proposed algorithm when increasing the number of GNSS receivers is shown in Fig. 8. Compared to Fig. 5, we see that adding two extra GNSS receivers increases the detection performance of the proposed scheme considerably.
As we see in Figs. 5-8, increasing the number of samples taken by the GNSS receivers consistently increases the detection probability for a given false alarm probability.

Conclusions
We proposed an anti-spoofing approach for GNSS based on a statistical test. The test exploits multiple GNSS receivers mounted on a rigid-body platform (with a priori unknown position and orientation) and essentially performs a consistency check between all pairs of measured receiver-satellite distances and available prior knowledge about the relative positions of the receivers on the platform. Numerical simulations proved the feasibility of our method and specifically showed that the more the spoofers try to manipulate the estimated GNSS receiver positions from their nominal locations, the higher is the probability of attack detection. Also, the more GNSS receivers on the platform, the higher the probability of detecting a spoofing attack.
We furthermore showed that using multiple GNSS receivers limits the feasible attacker position to few locations, and we designed a framework for finding the optimal attack positions as well as the artificial time delays of the spoofers. Simulations showed that when the spoofers do not occupy their optimal locations, it is easier to detect them. For analytical tractability and to achieve a first proofof-concept of our statistical test methodology, we considered a two-dimensional geometry. Future work may include extensions to a three-dimensional world or combining our approach with anomaly detection in the autocorrelation function of the received signals, in order to enhance overall detection performance. In addition, new analysis can be performed by considering synchronization errors among the clocks of the GNSS receivers and the satellites. Also, the noise variance may be treated as an unknown parameter in the GLR tests.
It would be interesting to test our proposed approach in practice. We hope that this paper will stimulate both further theoretical research and experimental work.

Proof of CRLB for H 0
In this appendix, we calculate the CRLB for estimated parameters under hypothesis H 0 . Since the range measurements in (5) follow a normal distribution, we can use the Slepian-Bang formula [56] to calculate the CRLB as, and C −1 (α) = σ −2 since the noises across the GNSS receivers and the measurements are assumed to be independent.
To calculate the CRLB for H 0 , first, we calculate the values of ∂μ(α) ∂α j for j = 1, 2, 3. After some algebraic calculations, these values are derived as  23 . Using the calculated elements of the Fisher matrix, we can derive closed-form expressions for the CRLB by inverting the Fisher information matrix. For the sake of conciseness, here, we avoid showing these expressions.

Proof of CRLB for H 1
In this appendix, we calculate the CRLB for the parameters under hypothesis H 1 . Since the range measurements in (5) follow a normal distribution, we can use the Slepian-Bang formula [56] to calculate the CRLB.
Based on Algorithm 2, the first step is to find initial values of the parameters τ i for i = 1, .., I. As we see in (45), estimation of s i depends on r 2 n,i for n = 1, ..., N where r 2 n,i depends on τ i for n = 1, ..., N. Hence, only τ i is used in (43) to estimate s i for i = 1, ...I.
In the next step of Algorithm 2, calculated values of s i are used to estimate τ i for i = 1, .., I. Based on (50), only s i is used to estimate τ i for i = 1, .., I. According to the previous explanations, each pair (s i , τ i ) is estimated independently for i = 1, ..., I. In the following, we provide the details to calculate the CRLB for H 1 .