GPS SPOOFING DETECTION USING MULTIPLE ANTENNAS AND INDIVIDUAL SPACE VEHICLE PSEUDORANGES

Spoofing is the common term used for describing the intentional broadcasting of false radio frequency signals intended to disrupt and mislead systems that depend on accurate position, navigation, and timing information provided by the Global Positioning System (GPS). Spoofing is an increasingly recognized threat which is garnering increased interest from researchers and users, both military and civilian. This thesis presents a novel GPS spoof detection algorithm that exploits the geometric distribution of a horizontal array of GPS antenna-receivers and the geometric configuration of visible navigation satellites. Using a Neyman-Pearson hypothesis testing formulation, a spatial correlation test is developed that can accurately and dependably detect a GPS spoofing scenario. Analysis is conducted showing the performance effects of the number of receivers used, internal receiver clock bias estimation, and temporal and spatial locations of the detector. Simulations were conducted using theoretical definitions of false alarm and detection probabilities, a GPS simulator and receiver combination, and a live-sky experimental set-up. Experimental and theoretical performance results are presented. ACKNOWLEDGMENTS I would like to first acknowledge my advisor, Dr. Peter Swaszek for his guidance and mentorship during the course of this work and my studies at University of Rhode Island. His insight and intuition was invaluable in solving many of the challenges in the implementation of this research. The theory behind this detection solution is the product of Dr. Swaszek’s expertise and ingenuity; I appreciate the opportunity to contribute to the theory and demonstrate its effectiveness. Many thanks are due to Dr. Richard Hartnett, Capt USCG (Ret.) and Commander Kelly Seals, USCG, from the U.S. Coast Guard Academy Electrical Engineering section; because of their support, encouragement, and knowledge I have had the opportunity to pursue this research and an advanced degree in Electrical Engineering. For inviting me to participate in the team with themselves and Dr. Swaszek researching this problem, I am grateful. I would also like to acknowledge Senior Chief Electronics Technician Ken McKinley, USCG, and Electronics Technician Second Class Rob Gutzeit, USCG, for their support and technical knowledge; their assistance was instrumental in completing the simulations and experiments described in this study. It is due to the love and support of my family, Michael, Sharon, and Andrew Radin, that I am in a position to contribute this research to the pursuit of scientific knowledge. I couldn’t have asked for a more awesome family. And to my beautiful and loving wife Brenna, whose understanding, patience, and compassion is infinite, thank you for everything, I love you, and aloha!


LIST OF TABLES
The goal of this thesis is to develop and analyze an algorithm for detecting GPS spoofing attacks using pseudoranges outputted from an array of commercial-offthe-shelf GPS antenna-receiver combinations.

Motivation
This research is motivated by the need to develop noninvasive, inexpensive, and inter-operable software GPS spoof detectors.
A software detector is a desirable solution because of its relative ease of implementation and interoperability with legacy equipment. A laptop pre-loaded with required software could be connected to any commercial-off-the-shelf GPS receiver that outputs the requisite data and an instant spoof detector would be deployed.
It is easily recognized that this form of detection solution would save considerable costs in development and installation versus a hardware based detector.
Prior work on this subject provides many promising detection solutions at all stages of the receiver data processing operation. Many of these detection methods involve extracting data from receivers in an intrusive manner that would compromise the integrity of legacy or certified systems in active use. Figure 1 presents a simplified depiction of information transformation as it flows through a GPS receiver. A spoof detector could work with RF data directly [4], [5]. At the other extreme, the detector could use position solution data [6] [7]. An alternate and middle-ground possibility is a detector based on pseudorange data. Pseudorange measurements are readily available from most mid-upper range receivers through standard ASCII or NMEA code outputs, and thus satisfy the need for non-invasive data extraction.
The choice of pseudorange measurements is not an arbitrary one. Ease of access has its drawbacks; it is expected that performance is inversely correlated to the complexity of data extraction. A hypothetical performance comparison is shown graphically in Figure 2. This figure plots probability of detection (P D ) versus probability of false alarm (P F A ). While described in more detail later in this thesis, at this moment we note that the closer a performance curve is to the to the upper left corner the better. A detector based on the received RF, while harder to implement, should have the best performance because the RF data represents the complete data set. The position solution is a much reduced version of the received signal; it presents the worst performing detector in that much of the detail in the RF is 'averaged' away. Thus, it is expected that pseudorange detection methods will provide better performance than methods using position solutions; developing a pseudorange detector is the next logical step in both complexity and performance.
The algorithm presented in this thesis is designed to satisfy all of the above motivating factors.

Thesis Organization
This thesis conducts background discussions of the Global Navigation Systems (GPS) constellations, GPS spoofing and spoofing detection, navigation characteristics of the Global Positioning System, COTS receiver data processing, and random processes to support the techniques used for development of this test.
We will start by framing the problem and discussing some background information in Chapter 2. Chapter 3 provides the notation and conceptual development of the test. The testing and simulation methodologies are presented in Chapter 5 and followed by the performance evaluations of all the tests in Chapter 6. Finally, a discussion of conclusions and future work is presented in Chapter 7. Appendices A and B provide detailed derivations of the formulas used and detailed hardware and software configuration notes. others (India, Japan). The focus of this study is on GPS, but the conclusions and methodology could be applied to any of these systems with some modification.

Integrity of GNSS
Threats to the integrity and functionality of these systems can be classified into the broad categories of jamming and spoofing. Jamming is essentially the overpowering of the authentic signal with noise or clutter to prevent the receiver from obtaining a correct navigation decision. These attacks are well researched and understood with plentiful counter-measures available. The threat area of growing interest is spoofing, which effectively amounts to providing receivers with counterfeit GNSS data to confuse, disrupt, and/or mislead the receiver's position or timing solution. Some potential spoofing methods are described in detail in References [8] and [9]. Current research and literature have presented many methods of detection and, potentially, defense against a remote GPS spoofing attack. These methods consist of various carrier to noise ratio (C/N o ) detectors in the acquisition plane, [4], [5], spatial energy detection at the input [10], post-position solution detection [6], [7], and integrated inertial measurement unit data / GPS receiver systems [7] amongst others [11], [12]. Some methods use single antennas, others multiple antennas contained in both single platforms and platform networks [13]. These methods provide many promising detection techniques and attack the detection problem at almost all phases of the receiver data flow.
Much literature is available on the general theory and background behind GNSS and GPS in particular, as well as the tools required to analyze and evaluate the experimental data. For example, [1] and [14] provide in depth discussions on the topic, from the construction of the receiver components to the decision algorithms and code generation. We will commence on a brief discussion of GPS in the following paragraphs.

Global Positioning System
GPS navigation is a time of arrival (TOA) navigation system fundamentally based Control Segment consists of a series of terrestrial monitoring and control stations that maintain the health of the system. User Segment consists of GPS receiver equipment which receives, processes, and analyzes the data received from the Space Segment to provide position, navigation and timing (PNT) information to the user [15].

Figure 3: GPS Constellation
The GPS constellation consists of 32 active satellites (as of the drafting of this thesis) [16] in six equally spaced orbital planes. Figure 3 is a depiction of the satellites in the orbital planes [17]. These satellites broadcasts two primary signals, the legacy civil signal, L1 C/A (Coarse Acquisition) code and the military signal, P(Y) or Precise code. The L1 signal is broadcast at 1575 MHz. There are also new signals: L2 at 1227 MHz, and L5 at 1176 Mhz. These are the carrier frequencies, where the navigation data is embedded. The P(Y) code is broadcast on both L1 and L2 frequencies for military users [15]. An example of the spectrum of GPS broadcast codes is shown in Figure 4 [18].   Figure 5 shows an example of code-delay correlation (this figure is reproduced from [1]). The delay between the peak in the correlation function and the receiver's clock is the code delay; this measurement (in units of meters) is the pseudorange. The frequency of the received signal is shifted from the designed carrier frequency; this error from the designed carrier frequency is Doppler frequency shift. A search is conducted over frequencies around the designed frequency until a maximum correlation is found. Figure 6 graphically shows this 2-D grid search across code-delay and frequency (this figure is reproduced from [1]). This search is conducted as described in detail in Chapter 7.7 of [19]; the smaller the "frequency bin" the more precise the correlation and better signal detection, but longer processing time is required. The correlation peaks are found in the Doppler frequency vs. code delay plane and are exemplified by Figure 7. It is in this realm that a spoofer could create a falsified correlation peak and take a receiver hostage [8]. We discuss this point in further detail in the next section. Once a signal is identified as matching a known GPS satellite, the receiver compares its internal clock to the noted time of transmission, and this ∆t is used to calculate a range. Due to the many errors involved in this solution, (see  Once four pseudoranges have been measured by the receiver, a position solution is calculated based on the relationship where ∆x is the displacement from an previously calculated approximate position, H is a matrix of direction cosines pointing from approximated position to navigation satellites, and ∆ρ is the difference between pre-approximated pseudoranges and measured pseudoranges.
The next chapter picks up with a more in depth look at the pseudoranges and lays the groundwork for our test.

Spoofing
As mentioned in the previous section the commonly accepted method of spoofing begins with an attacker providing a counterfeit signal that matches the authentic one in time and frequency. This signal shall match the composition of the authentic signal in every aspect (with authentic information) so that one could not distinguish between the two in any manner. As shown in Figure 9 the attacker would then increase the power of his signal and slowly move the delay of his counterfeit signal away from the authentic one. Because the power of the spoofed signal will surpass that of the authentic one, the receiver's tracking mechanism would follow the new signal. The spoofer would have successfully hijacked the receiver at this point. This process is described in further detail in references [8] and [20].

Additional Theory
The tools used to conduct the analysis were coded throughout the length of this study, and use skeletons of functions provided in the GPS MATLAB toolbox from L3 Communications. Many of the rudimentary GPS calculations were completed using these functions [21]. Tactics for statistically evaluating the behavior of spoofed signals and the detector were drawn from reference texts [19], [22], [23], [24], [25] and [26].

CHAPTER 3
Notation and Concept Development

Notation
We begin by imagining an array of m antennas spaced around a circular horizontal platform of radius r. In local east-north-up (ENU) coordinate frame (discussed in detail in Appendix [C], the position of an antenna k is We assume equal angular spacing let θ represent the geographic orientation of the platform giving the angular position of the antennas, θ k , as: This horizontal array is implemented in the configuration example of Figure 10. It is important to note here that elevation is a relative term while azimuth is measured in degrees True (000 • T points towards geographic North). This is an important distinction because this is the way the Novatel GPS receiver we use in simulation reports its data and it allows for an orientation independent statistic to be developed. The Spirent GSS8000 simulator used for testing uses an axis as defined in Figure 11. In the local ENU format the position coordinates of each SV The range from SV n to the n th antenna is: d k,n = (d 0,n cos ψ n sin φ n − e k ) 2 + (d 0,n cos ψ n cos φ n − n k ) 2 + (d 0,n sin ψ n − u k ) 2 This equation is simplified to: (please see Section A.1 for full derivation) in which δ k,n is δ k,n = r cos ψ n [sin φ n sin θ k + cos φ n cos θ k ] (3.

1.4)
This relationship is shown in Figure 12 for k=1. Figure 12: Relationship between δ k,n , d 0,n , and d k,n The distance from any specific satellite to any antenna relative to its distance to the center of the array is approximately equal to −δ k,n . This starts to form the basis for our test. Figure 13 provides a graphical representation of these measurements. In this graphical example, the δ 1,n term is positive; if we had instead used Antenna 2, we would have seen δ 2,n appear longer than d 0,n and the resulting δ 2,n would have been negative. This alignment provides the spatial correlating effect of our test, as shown later, and plays a large part in our analysis. This is also the fundamental basis behind Equation 3.1.5.
Note: the following analysis depends on two summations of the δ k,n terms. These summations are valid for m ≥ 3. Section 3.3 explores these summations for m = 2.
where ρ k,n is the observed receiver pseudorange, b k is receiver clock bias, and w k,n is white Gaussian noise. As part of its solution method a receiver estimates b k , hence the receiver is estimating the actual range. We use the term d k,n to refer to the receiver's measured distance between an antenna k and a SV n. Assuming a perfect estimate of b k , this is We return to inaccurate estimates of b k later. Note: The COTS GPS receiver used in this study executes the ρ k,n − b k calculation before reporting the pseudorange [27]. Hereafter, d k,n will be synonymous with receiver observed pseudoranges.

Assumptions
Satellite range is much larger than the spacing of antennas in array (d 0,n r) The GPS constellation orbits Earth at a range of approximately 20000km above the Earth's surface. Our test is being developed for antenna arrays less than 20m in diameter. Thus, with 6 orders of magnitude difference, we ignore the contribution of the antenna spacing in this calculation.

Clock bias is estimated perfectly
Until discussed in Section 4.4, we will assume the receiver clock bias is perfectly estimated. This estimation process is discussed in detail in Reference [1] and is typically accurate to 10 −9 seconds. Removing the clock bias estimation simplifies the primary analysis and allows for a more concise formulation. The effect of the assumption is discussed in detail in Section 4.4.
Spoofing attacker is using a single transmitting source This assumption is accurate exempting only the most complex of spoofing attacks [8]. A multi-point transmitter broadcasting independent and unique satellite RF signals creates a very challenging detection environment and is beyond the scope of this research. What this assumption allows us to do is to declare that under spoofing every antenna in our array measures identical pseudoranges.
This point is discussed further in Section 4.1.

Special Case: Two-Antenna Test
The primary focus of this study is on antenna arrays of greater than two antennas, but there is also use in developing the theory for a 2-antenna array, for example, we could imagine an array consisting of an antenna on the front and rear of a shipping container, truck, or ship. This example lends itself to a two-antennae array more-so than a multi-antennae circular array.
The derivation of this special case is the equivalent to the multi-antenna case until we derive the distribution of the Test Statistic in Section 4.3.

CHAPTER 4
Development of Hypothesis Test

The Hypotheses
Our hypothesis test has two scenarios: the null hypothesis, H 0 , where no spoofing is present and the alternate hypothesis, H 1 , where spoofing is present.
H 0 : With no spoofer present each individual range measurement is an accurate estimate of the actual range for that antenna and satellite pair.
H 1 : With the spoofer present scenario, we assume the spoofer uses a simple oneradiator spoofing environment, and thus we assume that each individual antenna will receive the same identical RF GPS signals [8] and independent noise. A more complicated multi-radiator spoofer would be more unpredictable and this case is not addressed in this paper. The model of the receiver pseudorange is n is the spoofer's generated pseudorange to the n th satellite. All antennas will receive these identical signals at different times due to propagation range, but this difference will be estimated out through the receiver clock bias (effects of receiver clock bias are presented later in Section 4.4).
We model each noise term, w k,n , as independent Gaussian statistics with zero means and equivalent variances σ 2 under both hypotheses. We can thus model our psuedorange distributions for H 0 and H 1 in the standard Gaussian shorthand distribution format (x ∼ N (µ, σ 2 ), where µ is the mean and σ 2 is the variance) as: respectively.

The Hypothesis Test
We construct our test viewing the null hypothesis, H 0 , as a simple hypothesis of Neyman-Pearson formulation [25]. Using this method we will derive a test statistic as a function of the psuedorange data T d 1,1 , ..., d m,N and compare this statistic to a threshold value, λ, to declare spoofing or no spoofing.
With this method of testing we have two fundamental metrics for our test: Probability of False Alarm (P F A ), and Probability of Detection (P D ).

P F A A False Alarm is a type 1 error in statistical hypothesis testing and occurs
when the H 0 hypothesis is rejected when H 0 is in fact true. The Probability of False Alarm is known as the size of the test or level of significance (α).
This level of significance is the amount of risk of erroneous detection built into the test [24]. If this level is exceeded, then an event of significance has occurred and it can be concluded that the outcome was not the result of sampling error. Mathematically P F A is defined as where R 1 is the critical region, or set of values where H 1 is decided (H 0 rejected) [26].
P D A Detection is a successful rejecting of the H 0 hypothesis when H 1 is true, in other words, deciding H 1 when H 1 is correct. The Probability of Detection is known as the power of the test [26] and is defined as where R 1 is the region of possible values where H 1 is declared.
When using a Neyman-Pearson criterion, the optimum test is a likelihood ratio test (LRT) [23] which takes the form: where the notation d k,n signifies the set of all m · N range measurements.
The LRT is a ratio of the conditional probability density functions (PDFs) of the measurements under the two hypotheses. In Equation 4.1.1 we characterized both hypotheses as Gaussian Normal Distributions and so we use the Gaussian PDF to express our LRT as:  The GLRT takes the form [26]: In Section A.4 we show how the MLE of the unknowns is derived and discover that the best estimation of d

Because the test statistic Equation (4.2.3) is a linear combination of Gaussian
variables, the test statistic itself is also Gaussian distributed.
Under hypotheses H 0 and H 1 the distributions are: and Details of the distribution derivations are included in Section A.5.
A test with Gaussian statistics has a false alarm probability (size of test) in the form [24] Here we reference the Gaussian right tail probability Q function which is defined as Q(x) = 1−Φ(x), Φ being the cumulative distribution function (CDF) [22] The power of the test, or probability of detection (P D ) takes the form: Neyman-Pearson hypothesis testing typically fixes P F A and allows us to solve for the desired threshold: Plugging in expressions (4.3.4), (4.3.1), and square root of (4.3.2) into the P D equation yields a final performance expression: The terms present in the performance equation highlight a few key characteristics for us. Namely the performance of the detector increases as the radius of the array, r, is enlarged or the number of antennas in the array, m is increased. The performance likewise decreases as the user estimated range error (UERE), σ, or signal to noise ratio (SNR) is increased. These are obvious and expected.
One term remains in this equation that will be explored in more detail in the next section, and this is the satellite constellation/ satellite elevation dependent term, N n=1 cos 2 ψ n .

Sky Term
Let's define the Sky Term as: This expression is a function of the number of satellites in view, N , and their individual elevations,ψ n . Since 0 • ≤ ψ n ≤ 90 • , cos 2 (ψ n ) ranges from 0 to 1 and |cos 2 (ψ n )| is inversely related to the elevation of the satellite. Furthermore, we can state that It is a satellite-dependent predictor; we can look at the relevant GPS Almanacs and predict how well the test will perform. Figure 15 demonstrates the temporal variation in this detector, which nearly repeats every 24 hours. In this example plot we see a maximum Sky Term of just over 10 around hour 17, and a minimum Sky Term just under 4 around hour 6. The average Sky Term throughout the day is near 6; we see this result replicated in Figure 16. Performance spoiler alert: we use this average value of 6 to build Table 1 which is a quick, easy reference for what radius the antennas should be placed at to achieve desired performance with a guess-timated UERE and α set for .001. Using almanac data, Figure 16 shows the average Sky Term over the course of a day (August 20, 2014) in North America. This highlights the spatial variation in our test. Less obvious here is the implications on temporal variation; the GPS constellation has a period of 11 hours and 58 minutes thus shifting the bands of

Two-Antenna Arrays
Similar to the multi-antenna case, the two-antenna test statistic will be Gaussian distributed.

Clock Bias Effects
Let's re-state the range measurement equation, Equation 3.1.7, here: Up until this point, we assumed that the clock bias, b k , estimate was perfect, i.e.
b k = b k . If we remove that assumption and include this estimate, b k in the equation, we can analyze the effect of this imperfect estimation in the test. Restating the definition of the range measurement: The full analysis of this subject is presented in [28]; we will present the conclusions here: The mean under the imperfect clock bias estimation remains the same at µ 0 , but the variance is modified to: matrix generated by relevant satellites' azimuths and elevations.
The result is that the detector performance is reduced by a small amount (not particularly significant), independent of the actual clock bias.

Methodology of Performance Testing
In order to effectively present and substantiate the performance of the detection method, testing was completed in three phases of increasing complexity.
Theoretical The theoretical performance metrics (P F A and P D ) are used to generate expected performance curves.
Simulation Using a Spirent GSS 8000 GNSS simulator, simulated constellation RF was created and fed into a NovAtel receiver for processing. Extracted pseudoranges were fed into the testing algorithm, and a realized performance curve was generated.
Experimentation An array of three antenna-receiver combinations were deployed to measure real-time GPS signals. Extracted pseudoranges were fed into the testing algorithm for a spoof/no-spoof declaration.
A more in-depth discussion of methodology is undertaken in the following sections.
All data processing and function calls are completed in MATLAB®.

Almanac-based Simulations
The first level of testing uses GPS Almanac data provided by the United States Coast Guard Navigation Center [16] to determine visible SV position information for a fixed reference location and then calculate ranges to the array center and offset antennas. These range calculations are combined with some additive white Gaussian noise (AWGN) to produce pseudoranges. These pseudoranges combined with SV position information provide the inputs to the hypothesis test and a performance curve is generated.  Figure 17 shows these positions.
2. Using the initial position's ECEF coordinates, the almanac data can be used to determine which space vehicles are visible above the 5 • mask angle and their ECEF coordinates, relative elevations, and relative azimuths. Functions from [21] are used for this step.
3. Using the equations given in Chapter 4 the δ n,k 's are calculated.
(5.1.1) Figure 18 shows an example of various thresholds plotted against some T0 and T1 data. It is easy to visualize the above relationships on this graphic for each threshold.
6. Finally, the ROC curve is generated by plotting P D vs. P F A . This is the standard procedure for generating performance curves.
Results Shown in Chapter 6     The hardware configuration for this test is described in Appendix B.
Given inputs Desired P F A ; number of antennas m; radius of antenna array r; Procedure 1. The receiver/antennas were all spaced according to their H 0 radius, turned on at roughly the same time, and started collecting data which was extracted in the same ASCII codes as the simulation test.
2. When "spoofing" occurred the antennas were moved into a cluster position as close as possible. Figure 22 shows the transition graphically.  Hour of day Figure 24: Performance of test utilizing various Sky-terms Next, we present the performance of the two-antenna test. In Figure 25 ROC curves were generated for a two-antenna system at various θ orientations and compared to a three-antenna system. These performance curves are generated using a single snapshot of SVs. The constellation is shown in Figure 26. As expected, no orientation of a two-antenna array can compare to a three-antenna array, but performance is affected and reduced for poor geometric alignment. The more broadside the SVs are to the array orientation, the lower the performance.

Simulator Testing
The GSS8000 Simulator tests provided excellent results as expected. Shown below are the performance plots generated using 24 hours of data from the simulator. Figure 28 shows the observed simulation performance plotted against theoretical performance curves generated using maximum and minimum observed Sky Term bounds. For these tests, the parameters , σ = .1269 and r = .10m. were used. The performance in the plot is not all that impressive -this is due to the incredibly small radius (.1m) used to generate the data. A radius this small was required because the only noise experienced in this exercise was receiver noise at less than a tenth of what atmospheric noise would normally provide [1]. Any increased radius results in unity performance, which is awesome, but does not present well graphically. The curve is also a little jagged, and that is due to the snapshot method of performance detected described in the last chapter.
It is interesting to see how the observed curve fits cleanly between the maximum and minimum theoretical curves -and this makes sense. Because the test is time variant the performance curves are also time variant; if a ROC curve was drawn at each snapshot, we would see the curves fill the space between these bounds perfectly. As such, the observed simulation curve is in effect an average of all potential curves.
In Figure 29 we see how well the estimated P F A matches the desired exact P F A .
Ideally, this would be a straight line. The curvatures (errors) are due to an imprecise estimate of observed σ. This is a graphical realization of the noise estimation problem discussed in Section 6.3.1.
As the estimation of noise improves, the performance of the test improves. This dilemma manifests itself again in the next test. For this test, the following parameters were used: σ = 2.61, r = 3.307, andm = 3.

Real-Time Multi-Antenna Testing
The results of this experiment match predictions and demonstrate some key points of the test. Performance was set for 1% false alarm rate for all experiments. The Sky View for each test is shown in Figure 30 and theoretical performance curves (using the sample mean Sky Terms) for each case are compared in 31.

Case 1: Obstructed View Experiment 1
In Case #1 the antennas were exposed to an obstructed western sky-view that blocked all SV's to the West. This greatly reduced the average Sky Term to a quite minimal 2.805. This is far below the average predicated Sky Term of that position (given an assumed unobstructed sky-view) of 6.0898 as shown in Figure 32. The performance is obviously very poor. Figure 33 shows    The sky view was unobstructed (the few nearby buildings didn't rise above the 5• mask angle.) and so our experienced SkyT erm was much improved. Performance was also much improved, as you can see in Figure 35. This is expected.

Noise Measurement
As discussed in Section 4.3, the threshold is a function of the UERE. The UERE is difficult to measure directly, and must be estimated in some manner. The NovAtel receiver provides us with a root-mean squared (RMS) average of the standard deviation (1σ) of range measurement errors. [29] This value is provided in field three of the GPGST ASCII string referenced in Section B.3.
Since our test is inherently weighted towards lower-elevation satellites, we want to weight our noise estimate equivalently. A mini-experiment to determine this scale factor was conducted by taking a four hour time sample and looking at the range residuals as a function of the elevation of the satellite. The resultant residuals grouped by low-elevation satellites (< 45 • ) and high-elevation satellites (> 45 • ) is shown in Figure 37.  The equipment used were off-the-shelf pieces found in an undergraduate EE communications laboratory; these were not expensive and custom items. Furthermore, they were units the U.S. Coast Guard actually has installed on various assets: we know the test works with this equipment! Given a wideenough array, this test would provide very dependable results at a fraction of the cost of an integrated hardware system. We can thus declare the objectives met.

Future Work
There is still work to be done on this detection method.
The analysis can (and should be) extended to include a 3-dimensional array.
If this test works for low-elevation satellites using a horizontal antenna plane, then it should work for high-elevation satellites using a vertical antenna plane.
What is the optimal configuration under this scenario? Furthermore, a more in-depth study of pseudorange noise estimation needs to be conducted to nail down a procedure for procurring a real-time UERE value for detection calculations. The last step of analysis is an obvious one -how does this test function with moving platforms? My suspicion is that it will take the form of the noncoherent unknown orientation problem discussed in [28] and will experience a small drop in performance (will likely exacerbate clock bias estimation performance loss).
Better testing methodologies can be employed to solidify performance promises; we need to develop real-time testing software and experiment with it. An obvious testing improvement would be to implement a broadcasting spoofer to test our array with. When someone markets the iSpoof we can realize this goal. Next, we need to test with a moving platform, and identify methods for precisely monitoring spoof vs. no-spoof test conditions so delays in detection identification can be measured.

A.1 Multi-Antenna Array Distance
The range from SV n to Ant k is: Expanded and simplified: d k,n = d 2 0,n − 2rd 0,n cos ψ n (cos θ k cos φ n + sin θ k sin φ n ) + r 2 With d 2 0,n >> r 2 we can neglect the r 2 term . d k,n ≈ d 2 0,n − 2rd 0,n cos ψ n (cos θ k cos φ n + sin θ k sin φ n ) Defining δ as: δ k,n = r cos ψ n (cos θ k cos φ n + sin θ k sin φ n ) Note: the maximum value of each δ k,n is r.
The expression is further simplified to: And finally simplified to: An approximation of this expression is conducted using Taylor's Theorem [30] letting x = 2δ k,n ρ 0,n * about x = 0: * We can truncate this approximation after 2 terms because the size of x is quite small relative to the range. max (δ k,n ) = r, r d 0,n ≈ 1 10 6

A.1.1 Multi-Antenna Array Facts
Mirroring [28], we develop some initial facts for simplifying the analysis: Calling again on the cos(α + β) = cos(α) cos(β) − sin(α) sin(β) identity and some further simplification: Now let's show the derivation for m k=1 sin θ k = 0, starting with substituting Equation 3.1.2 for θ k : Using the trigonometric identity for sin(α + β) from above: Pulling out the scalars: It's clear that this case now mirrors the Z c case derived in detail above.
And thus we conclude derivation of Fact 1. We can see this has the components Z s and Z c from Fact 1. Applying these facts show that this Fact is also 0.
It is clear to see that the 2πk m trig terms are equivalent to the analysis from Fact 1 and will simplify to zero. Thus, we can conclude at this point with our solution: Deriving Z cc begins with the power reduction identity for cosine, cos 2 θ k = 1+cos 2θ k It is clear to see that the 2πk m trig terms are equivalent to the analysis from Fact 1 and will simplify to zero. Thus, we can conclude at this point with our solution: Deriving Z sc begins with the product-to-sum identity for sin/cosine, sin α · cos β = Removing the zero term and consolidating: With some foresight, we can convince ourselves that this derivation will follow the same as the previous one. Thus we declare the conclusion:  Using cos(α + β) identity twice: (r · cos φ n [sin φ n sin θ k + cos φ n cos θ k ]) = r · cos φ n sin φ n 2 k=1 sin θ k + cos φ n 2 k=1 cos θ k Here we can apply Fact 1 of the two-antennae scenario and conclude: sin θ k cos θ k = sin 2θ Mirroring our work from the Multiple-Antenna Fact 3 derivation above, we begin with Z ss : cos 2 ψ n $ $ $ $ $ $ $ $ $ X 1 sin 2 φ n + cos 2 φ n − sin 2 φ n cos 2θ + 2 sin φ n cos φ n sin 2θ + cos 2 φ n cos 2θ = r 2 2 k=1 cos 2 ψ n 1 + − sin 2 φ n +$ $ $ $ $ X 1 − sin 2 φ n cos 2 φ n cos 2θ + 2 sin φ n cos φ n sin 2θ = r 2 2 k=1 cos 2 ψ n 1 + 1 − 2 sin 2 φ n cos 2θ + 2 sin φ n cos φ n sin 2θ

A.4 Generalized Likelihood Estimator and Maximum Likelihood Estimators
Given a conditional probability density function, p r|a (R|A) the goal of the MLE is to find "the value of A that most likely caused a value of R to occur" [23].
The well known process of evaluating an MLE is to take the log-likelihood function of the distribution, differentiating the expression, and setting it equal to zero. This Because the test statistic is Gaussian, we characterize it by the means and variance under both hypothesis.

Multi-Antenna Array
Starting with the null hypothesis, H 0 ,

Two-Antenna Array
Starting with the null hypothesis, H 0 , cos 2 ψ n 1 + 1 − 2 sin 2 φ n cos(2θ) + 2 sin φ n cos φ n sin(2θ) alternate hypothesis, H 1 : cos 2 ψ n 1 + 1 − 2 sin 2 φ n cos(2θ) + 2 sin φ n cos φ n sin(2θ) We previously determined that the variance of d k,n is equivalently σ 2 for both hypotheses; likewise the test statistic variance is equivalently σ 2 T under both hypotheses. Real-time data is not flawless; the receivers used in simulations and experiments were configured (as described in an upcoming section) to report data in one second intervals. Occasionally, a message code was skipped and the requisite data missing, or the codes are printed in the wrong order. Both error types cause the data to be incorrectly classified and invalidates the entire collection. These errors must be recognized and corrected prior to data analysis. As the functions stand, the corrections must be made by hand by either duplicating a neighboring time-set's data or re-ordering the codes in the original data file.
This grooming process is quite intuitive, and not altogether surprising, but the author would like the reader to be aware of any and all pitfalls in the creation of these experiments.

Data Alignment
The importance of aligning the data correctly for the computations described in this text that involve actual receiver-generated data cannot be overstated and thus I wanted to highlight a few key points here. I write this section for the purpose of easing the process of recreating the experiment.

Timing
That the data needs to be organized by time-sample is intuitive but I want to provide some support for this; I will briefly describe my process. The key to this requirement is that the Test Statistic and Threshold values both change with time, related to the GEO-Term as previously discussed.
All the experiments done here used post-processed data and the time-alignment was a two step process. First the data from each "receiver" was parsed from its ASCII strings into a .mat file organized by time-step. Each section contained all the observation data from each satellite vehicle at each time.
Secondly, we had to ensure that all of the receivers were reporting the same timesample at the same iteration during the data analysis. This was slightly more challenging because at times the receiver output could provide corrupted data (in various forms) that were thrown out during the parsing process leaving some time-

B.2 Hardware
The hardware used for the experiment was: • Three (3)    The antennas were spaced using a measuring tape set for 5m and oriented facing due North (000 • T ). These measurements were done at the best accuracy of the author and assumed to be accurate; the effect of any measurement errors would be absorbed by the test and were likely to be immeasurably small.    NovAtel Connect was configured to write ASCII strings with standard receiver processing data at intervals of 1 second. These strings were recorded into a predetermined text file for post processing. A screen shot of the data collection window with typical settings is shown in Figure B.7. The third value in the GPGST string (2.41) is the RMS pseudorange noise term referred to in Section 6.3.1.

APPENDIX C Geographic Coordinate Transformations
In Chapter 3 we presented Equation 3.1.1 which is the ENU coordinate frame version of our individual antennae locations. In Chapter 5 we talk about converting a latitude/longitude reference point into ENU coordinates for antenna array generation, and then converting these ENU coordinates to the ECEF coordinate frame for calculations involving space vehicles. This section is to present the theory behind the coordinate frames and mathematically how these transformations are conducted.

C.1 Geodetic
The Latitude-Longitude-Height coordinate system is one of the most common methods of describing a geographic position on Earth and is known as the

C.2 Local Earth-North-Up
The local Earth-North-Up coordinate frame is a cartesian coordinate system that consists of a horizontal plane placed tangent to Earth's surface at any reference position. This plane does not follow the curved surface of Earth, and so the further a point is placed from the reference position in a local ENU system the greater the "error" or difference from the actual geographic location of the point. Figure C.2 demonstrates this coordinate frame along with the relationship between ENU, ECEF, and Lat-Long.
When generating a Local ENU frame, we start with a given Latitude/Longitude to serve as a reference position and is represented as (0, 0, 0) local. Other ENU points can be given as unit vectors from this reference point. In other words, the process is to convert geodetic to ECEF, and then ECEF to ENU. Likewise, to go back to geodetic coordinates we convert ENU to ECEF and ECEF to geodetic.
Converting a point from ECEF to ENU creates a unit vector from a reference point (h, θ, φ) to the new point. The conversion looks like this:

C.3 Earth-Centered, Earth-Fixed
Earth-Centered, Earth-Fixed (ECEF) coordinate frame is one which the center of the earth is at position (0, 0, 0) and coordinates are presented in a (X, Y, Z) format. The X ECEF axis points out from the center of Earth through the point Converting from geodetic coordinates to ECEF utilizes the following equations: (N (θ) + h) · cos θ · cos φ (N (θ) + h) · cos θ · sin φ We are going to stop the formulation here and refer the reader to the references for a more technical discussion of Earth's geodetic coordinate system characteristics , [31]. The transformation from ECEF back to geodetic coordinates was not used in this paper but if the reader is interested in the topic I refer her to References [32] and [1].