Performance Analysis of Direct GPS Spoofing Detection Method with AHRS/Accelerometer

The global positioning system (GPS) is an essential technology that provides positioning capabilities and is used in various applications such as navigation, surveying, mapping, robot simultaneous localization and mapping (SLAM), location-based service (LBS), etc. However, the GPS is known to be vulnerable to intentional attacks such as spoofing because of its simple signal structure. In this study, a direct method is proposed for GPS spoofing detection, using Attitude and Heading Reference System (AHRS) accelerometer and analyzing the detection performance with corresponding probability density functions (PDFs). The difference in the acceleration between the GPS receiver and the accelerometer is used to detect spoofing. The magnitude of the acceleration error may be used as a decision variable. Additionally, using the magnitude of the north (or east) component of the acceleration error as another decision variable is proposed, which shows better performance in some conditions. The performance of the two decision variables is compared by calculating the probability of spoofing detection and the detectable minimum spoofing acceleration (DMSA), given a pre-defined false alarm probability and a pre-defined detection probability. It turns out that both decision variables need to be used together to obtain the best spoofing detection performance.


Introduction
The global navigation satellite system (GNSS) is an essential technology for positioning and timing, and its application covers various areas such as navigation, surveying, mapping, robot simultaneous localization and mapping (SLAM), location-based service (LBS), etc. The global positioning system (GPS) is the original GNSS and its full operational capability (FOC) was declared in 1995 in the United States of America. The legacy L1 C/A code signal of GPS is very weak at the Earth's surface and has a simple structure to implement [1,2]. Thus, the GPS signal is vulnerable to intentional interferences such as jamming and spoofing. While jamming attacks make the GPS receivers malfunction, spoofing attacks make the target receiver unaware of it being attacked by the spoofer. Spoofing threats have garnered attention since the initial finding of the 2001 Volpe Report [3]. GPS spoofers may cause significant damage to the target receiver by transmitting counterfeit navigation data which can result in erroneous navigation. Thus, spoofing attacks are a significant problem to users and many studies on spoofing attacks and anti-spoofing methods have been carried out since 2001.

Block Diagram to Obtain the Acceleration Error from GPS Receiver and Accelerometer
In this study, GPS spoofing detection is done by using the comparison of accelerations obtained from the GPS receiver and accelerometers. The block diagram of Figure 1 shows the procedure to obtain the difference of the two acceleration measurements. The accelerometers measure the specific forcef b acc and then, it is changed intof n acc through the transformation matrix C n b . Inf b acc andf n acc , the superscript b denotes the body frame and n denotes the navigation frame. The navigation frame uses the north(x)-east(y)-down(z) (NED) frame in this study. AHRS produces the transformation matrix C n b by using the sensor outputs and Kalman filter. C n b denotes the direction cosine matrix from the body frame to the navigation frame. The hat (ˆ) denotes measured or calculated values containing errors and Ψ denotes the skew symmetric matrix of the attitude error.
It is assumed that IMU calibration and initial alignment is performed in advance depending on the characteristics of various inertial sensors since there are many kinds of gyroscopes, such as ring laser gyro, fiber optic gyro, hemispherical resonator gyro, and low cost micro-electro-mechanical system (MEMS) gyro, and so on, and different gyroscope has different error sources, and accelerometer has also many types, such as pendulous type, vibrating type, silicon type, and MEMS type. Figure 1 shows the IMU calibration and initial alignment with the dotted block, which will not be considered in this paper. Thus misalignment, bias, scale factor, and others are assumed to be compensated in advance. The IMU calibration and initial alignment is an essential process in the inertial navigation system and thus there are much research results which have been already performed [19,20]. Only GPS spoofing detection will be considered in this paper.
Matrices use capital letters, and vectors, small letters. Matrices and vectors will use bold letters and scalars use plain letters.
Sensors 2020, 20, 954 3 of 22 and 5 using the two performance measures, the probability of detection and the DMSA. The conclusion is presented in Section 6.

Block Diagram to Obtain the Acceleration Error from GPS Receiver and Accelerometer
In this study, GPS spoofing detection is done by using the comparison of accelerations obtained from the GPS receiver and accelerometers. The block diagram of Figure 1 shows the procedure to obtain the difference of the two acceleration measurements. The accelerometers measure the specific force ̂ and then, it is changed into ̂ through the transformation matrix . In ̂ and ̂, the superscript b denotes the body frame and n denotes the navigation frame. The navigation frame uses the north(x)-east(y)-down(z) (NED) frame in this study. AHRS produces the transformation matrix by using the sensor outputs and Kalman filter. denotes the direction cosine matrix from the body frame to the navigation frame. The hat (^) denotes measured or calculated values containing errors and denotes the skew symmetric matrix of the attitude error. It is assumed that IMU calibration and initial alignment is performed in advance depending on the characteristics of various inertial sensors since there are many kinds of gyroscopes, such as ring laser gyro, fiber optic gyro, hemispherical resonator gyro, and low cost micro-electro-mechanical system (MEMS) gyro, and so on, and different gyroscope has different error sources, and accelerometer has also many types, such as pendulous type, vibrating type, silicon type, and MEMS type. Figure 1 shows the IMU calibration and initial alignment with the dotted block, which will not be considered in this paper. Thus misalignment, bias, scale factor, and others are assumed to be compensated in advance. The IMU calibration and initial alignment is an essential process in the inertial navigation system and thus there are much research results which have been already performed [19,20]. Only GPS spoofing detection will be considered in this paper.
Matrices use capital letters, and vectors, small letters. Matrices and vectors will use bold letters and scalars use plain letters.

GPS Kalman Filter
The GPS receiver usually provides position and velocity information. To obtain acceleration from the GPS receiver, the Kalman filter is used by including acceleration as a state variable. The dynamics of the GPS receiver can be described by the state-space model as in Equation (1), which has 11 state variables such as 3-dimensional position, velocity, acceleration, GPS receiver's clock bias , and drift .

GPS Kalman Filter
The GPS receiver usually provides position and velocity information. To obtain acceleration from the GPS receiver, the Kalman filter is used by including acceleration as a state variable. The dynamics of the GPS receiver can be described by the state-space model as in Equation (1), which has 11 state variables such as 3-dimensional position, velocity, acceleration, GPS receiver's clock bias c b , and drift c d . .
x, y, and w x , w y , w z , w b , w d are independent white noises. Pseudo range is the distance between the GPS satellite and the receiver. The difference between the measured pseudo range ρ i and the estimated pseudo rangeρ i is used as the measurement in the Kalman filter for i-th GPS satellite is the estimated user position, and v i is the white noise. The whole measurement equation is given as follows: From the dynamic Equation (1) and the measurement (2)

AHRS
AHRS provides the attitude and heading and thus, the direction cosine matrix C n b can be calculated uniquely if the rotation sequence of roll, pitch, and heading is pre-defined. Many approaches have been proposed for AHRS [22][23][24][25] using accelerometer, gyroscope, and magnetometer. Accelerometers provide roll and pitch, and magnetometers provide heading. Hence the roll, pitch, and heading obtained from accelerometers and magnetometers can be compared with those from the gyroscope, and thus the Kalman filter can be used to estimate attitude and heading.
In many cases, quaternion is used to avoid deadlock and to save time. Quaternion q is defined as one scalar and a three-or four-dimensional vector as follows: The direction cosine matrix C n b is related with the quaternion q as in (3).
The quaternion q is updated as the following differential equation Sensors 2020, 20, 954 The direction cosine matrix C n b , which is the output of AHRS, can be described aŝ where C n b is the true direction cosine matrix, and Ψ denotes the orientation error of AHRS and is a skew symmetric matrix as follows: where φ is roll, θ is pitch, and ψ is heading. As discussed in Section 2.1, it is assumed that IMU calibration and initial alignment is performed in advance before AHRS block as in Figure 1. Thus deterministic and some random errors are compensated in the IMU calibration and initial alignment block. Then the remaining orientation error δφ, δθ, and δψ can be assumed to have Gaussian distribution as in (4).

Definition of the Decision Variables and the Structure of the Second-half of the Proposed Direct GPS Spoofing Detection
This section describes the acceleration error equation for the direct GPS spoofing detection and defines two decision variables to decide whether a spoofing signal exists or not. One decision variable is the magnitude of the horizontal acceleration error and the other is the magnitude of the north (or east) direction acceleration error. The probability density functions and the thresholds for spoofing detection are given for the two decision variables.

Acceleration Error Equation
The acceleration estimated from the GPS receiver is described as follows, f n gps = f n gps + n gps (5) where f n gps is the true moving acceleration (plus spoofing acceleration if any) in the navigation frame and n gps is the white noise. The specific force measured from the accelerometersf b acc in the body frame is transformed into the navigation framef n acc by using the direction cosine matrixĈ n b obtained from the AHRS as follows, where f n acc is the true specific force, Ψ is defined in (4), and n acc is the white noise of accelerometers. The Coriolis effect is assumed to be negligible with the assumption of low moving velocity for brevity. Notice that the z-component of f n acc contains the gravity. The acceleration error equation is described from the difference off n gps andf n acc as follows: Sensors 2020, 20, 954 6 of 22

Decision Variable z mag as the Magnitude of the Horizontal Acceleration Error
Suppose that the hypothesis H 0 denotes the case of the absence of the spoofing signal, and H 1 denotes the case of the presence of the spoofing signal. For hypothesis H 0 , the acceleration error is denoted as z 0 , and for hypothesis H 1 , the acceleration error is denoted as z 1 . Then, z 0 and z 1 can be described from Equations (5) and (6) as follows: where f n s,n , f n s,e , f n s,d are north, east, down components of spoofing acceleration. Only the x, y components of the acceleration error equation (8), (9) are used in this study. The acceleration errors z 0 and z 1 are expressed in the navigation frame and the superscript n will be omitted henceforth, for brevity.
In Equation (9), the random variables z 1n and z 1e , which are north and east components of z 1 , have Gaussian distribution and the probability density function (PDF) is as follows, 3.2.1. Probability Density Function of the Magnitude of the Horizontal Acceleration Error z mag z 1mag is defined as the magnitude of the horizontal component of z 1 as follows: , where z 1n ∼ N f s,n , σ 2 n and z 1e ∼ N f s,e , σ 2 e with σ 2 n and σ 2 e being defined in (10) and (11).
The variances σ 2 n and σ 2 e have different values and depend on the AHRS attitude accuracy times moving acceleration. Thus σ 2 n and σ 2 e are time-varying if the moving acceleration varies with time. The PDF of the z 1mag could not be found in the literature and thus it was derived in this paper and called as Shim PDF. Lemma 1 shows the PDF of z 1mag .

Lemma 1. (Shim PDF)
Consider the independent Gaussian random variables X and Y with X ∼ N m 1 , σ 2 1 and Y ∼ N m 2 , σ 2 2 . Then, the magnitude Z = √ X 2 + Y 2 has the following PDF: Sensors 2020, 20, 954 7 of 22 Proof. The variable change for the independent Gaussian random variables X and Y is as The idea of using s and b above comes from Hoyt's paper [26]. Consider the joint PDF f UV (u, v).
can be obtained as follows: From the relation between random variables Z and R as Z = sR, the PDF of Z can be obtained from f Z (z) = 1 s f R z s as: From algebraic manipulation of Equation (13) z 2 + α 2 is obtained, and the integral term in Equation (13) becomes the integral term of I 1 (z), which results in Equation (12).
Rayleigh PDF and Rice PDF are well-known PDFs as the relation to Gaussian, and those two PDFs are special cases of Equation (12), which becomes Rayleigh PDF with the condition of m 1 = m 2 = 0 and σ 1 = σ 2 , and becomes Rice PDF with the condition of σ 1 = σ 2 .
Defining z 0mag as the magnitude of the horizontal component of z 0 as follows, z 0mag = z 2 0n + z 2 0e , where z 1n ∼ N 0, σ 2 n , z 1e ∼ N 0, σ 2 e . The PDF of z 0mag can be obtained from Lemma 1 with m 1 = m 2 = 0 and the result is shown in Corollary 2.

Corollary 2.
Consider the independent Gaussian random variables X and Y with X ∼ N 0, σ 2 1 and Y ∼ N 0, σ 2 2 . Then the magnitude Z = √ X 2 + Y 2 has the following PDF: where I 0 (z) = 1 π π 0 e z cos θ dθ, Proof. Equation (14) can be obtained easily from the PDF in Equation (12) with m 1 = m 2 = 0 and manipulation in the I 0 (z) part.
Equation (14) can be found in Reference [27] without proof.

Threshold to Detect GPS Spoofing for the Decision Variable z mag
In this study, the probability of false alarm is used to obtain the threshold for the detection of spoofing. The threshold γ mag to detect GPS spoofing is defined according to the pre-defined probability of false alarm P f a as follows : where and the probability is calculated from the integral of equation (14) from γ mag to infinity. Whether a spoofing signal exists or not is decided by the following decision rule: The variable z mag above is said to be a decision variable since it is used to decide whether there is a spoofing signal or not. Defining z 1absN and z 0absN as the magnitude of the north component of z 1 and z 0 , respectively, z 1absN = |z 1n | and z 0absN = |z 0n | (similarly, z 1absE = |z 1e | and z 0absE = |z 0e |).
The PDF of z 1absN (or z 1absE ) and z 0absN (or z 0absE ) can be obtained as Equations (17) and (18), which are called folded Gaussian [28], since z 1n and z 0n have Gaussian density functions as z 1n ∼ N f s,n , σ 2 n and z 0n ∼ N 0, σ 2 n . The threshold γ absN to detect GPS spoofing is defined according to the pre-defined probability of false alarm P f a as follows : here and the probability is calculated from the integral of Equation (18) from γ absN to infinity. Whether a spoofing signal exists or not is decided by the following decision rule:

The Structure of the Second-Half of the Proposed Direct GPS Spoofing Detection
This subsection shows the structure of the proposed second-half of direct GPS spoofing detection method in Figure 2, which is drawn after the rightmost signal in Figure 1. The analysis of the proposed structure shown in Figure 2 will be given in Sections 4 and 5 in detail.
In Section 4, it will be observed that the decision variable z absN (or z absE ) shows a higher detection probability than z mag in the condition that both moving acceleration and spoofing acceleration head within roughly 25 • from the north (or east). Section 5 shows that when DMSA is used for performance measure, the decision variable z absN (or z absE ) shows a smaller DMSA than z mag in the condition that both moving acceleration and spoofing acceleration head within roughly 25 • from the north-south direction (or east-west direction). From these results, a direct GPS spoofing detection method is proposed as follows: If any of the three decision variables z mag , z absN, and z absE are larger than or equal to the corresponding thresholds, then a spoofing signal is declared to exist.
Note that the threshold γ mag (t), γ absN (t), and γ absE (t) in Figure 2 are time-varying, not constant. The threshold γ mag (t) is obtained from Equation (15) given the probability of false alarm P f a , where the PDF is Equation (14) with σ 1 = σ n and σ 2 = σ e . The north and east variances σ 2 n and σ 2 e given in Equations (10) and (11) contain the moving acceleration and thus σ 2 n and σ 2 e are time-varying, which is why γ mag (t) is time-varying. The threshold γ absN (t) is obtained from Equation (18) and the PDF contains σ n , which is time-varying. Thus γ absN (t) depends on the moving acceleration and becomes time-varying. Similarly γ absE (t) is time-varying. The red line and arrow in Figure 2 means that the threshold γ mag (t), γ absN (t), and γ absE (t) depend on the moving accelerationf n acc .

The Structure of the Second-Half of the Proposed Direct GPS Spoofing Detection
This subsection shows the structure of the proposed second-half of direct GPS spoofing detection method in Figure 2, which is drawn after the rightmost signal in Figure 1. The analysis of the proposed structure shown in Figure 2 will be given in Section 4 and 5 in detail.
In Section 4, it will be observed that the decision variable (or ) shows a higher detection probability than in the condition that both moving acceleration and spoofing acceleration head within roughly 25° from the north (or east). Section 5 shows that when DMSA is used for performance measure, the decision variable (or ) shows a smaller DMSA than in the condition that both moving acceleration and spoofing acceleration head within roughly 25° from the north-south direction (or east-west direction). From these results, a direct GPS spoofing detection method is proposed as follows: If any of the three decision variables , , and are larger than or equal to the corresponding thresholds, then a spoofing signal is declared to exist.
Note that the threshold ( ), ( ), and ( ) in Figure 2 are time-varying, not constant. The threshold ( ) is obtained from Equation (15) given the probability of false alarm , where the PDF is Equation (14) with 1 = and 2 = . The north and east variances 2 and 2 given in Equations (10) and (11)

Performance Analysis of the Decision Variables using the Probability of Detection
This section shows the performance of the proposed direct GPS spoofing detection method for the two decision variables and (or ) which are defined in Section 3.

Performance Analysis of the Decision Variables using the Probability of Detection
This section shows the performance of the proposed direct GPS spoofing detection method for the two decision variables z mag and z absN (or z absE ) which are defined in Section 3.

Detection Threshold According to Moving Acceleration
Suppose that the probability of false alarm P f a is pre-defined. Then, spoofing detection thresholds γ mag and γ absN (or γ absE ) are determined according to P f a as in Equations (15) and (19). Taking the PDF Equations (14) and (18) into account, Equations (21) and (22) are obtained from Equations (15) and (19) to further obtain γ mag and γ absN .
where σ 2 n = 2 gps,n + σ 2 ψ f 2 acc,e + σ 2 θ f 2 acc,d + 2 acc,n and σ 2 e = 2 gps,e + σ 2 ψ f 2 acc,n + σ 2 φ f 2 acc,d + 2 acc,e . To see the detection performance result clearly, the vertical moving acceleration is supposed to be zero and the gravity is compensated before the acceleration error is obtained. Thus the following variances in Equation (23) are used in the simulation from now on. σ 2 n = 2 gps,n + σ 2 ψ f 2 acc,e + 2 acc,n and σ 2 e = 2 gps,e + σ 2 ψ f 2 acc,n + 2 acc,e The detection thresholds γ mag and γ absN depend on the variances σ 2 n and σ 2 e which are functions of moving acceleration f acc as in (23). Thus, the detection threshold γ mag and γ absN are not constant but vary according to the magnitude and direction of the moving acceleration f acc . Figure 3 shows detection thresholds γ mag and γ absN according to the direction of f acc with two cases of magnitude, (a) f acc = 0.2 m/s 2 and (b) f acc = 0.4 m/s 2 . The threshold is the distance from the origin for the corresponding direction of f acc . It is observed that for the same magnitude of moving acceleration, γ mag has maximum values in the north and east directions and γ absN has the minimum value in the north direction. Similarly, γ absN has the minimum value in the east direction.
To see the detection performance result clearly, the vertical moving acceleration is supposed to be zero and the gravity is compensated before the acceleration error is obtained. Thus the following variances in Equation (23) The detection thresholds γ mag and γ absN depend on the variances σ n 2 and σ e 2 which are functions of moving acceleration as in (23). Thus, the detection threshold γ mag and γ absN are not constant but vary according to the magnitude and direction of the moving acceleration . Figure 3 shows detection thresholds γ mag and γ absN according to the direction of with two cases of magnitude, (a) | | = 0.2 / 2 and (b) | | = 0.4 / 2 . The threshold is the distance from the origin for the corresponding direction of . It is observed that for the same magnitude of moving acceleration, γ mag has maximum values in the north and east directions and γ absN has the minimum value in the north direction. Similarly, γ absE has the minimum value in the east direction.

Effects of Moving Acceleration on the Performance of Spoofing Detection
This subsection analyzes the effects of moving acceleration on the performance of spoofing detection. The effects of moving acceleration, magnitude and direction are separately examined, for both decision variables and which are defined in Section 3. The probability of detection is used for the performance of spoofing detection with the pre-

Effects of Moving Acceleration on the Performance of Spoofing Detection
This subsection analyzes the effects of moving acceleration on the performance of spoofing detection. The effects of moving acceleration, magnitude and direction are separately examined, for both decision variables z mag and z absN which are defined in Section 3.
The probability of detection P d is used for the performance of spoofing detection with the pre-defined probability of false alarm P f a . When the detection threshold γ mag and γ absN are obtained from P f a , the corresponding detection probabilities P d,mag and P d,absN are defined as follows: P d,mag = prob z mag ≥ γ mag H 1 and P d,absN = prob z absN ≥ γ absN H 1 (24) where the probability density functions (12) and (17) are integrated from the detection threshold to infinity.

Performance of Spoofing Detection According to the Magnitude of Moving Acceleration
The probability of detection P d depends on both moving acceleration f acc and spoofing acceleration f s . The effect of the magnitude of moving acceleration is focused on in this subsection.
Suppose that the probability of the false alarm is pre-defined as P f a = 0.001, the AHRS attitude accuracies of roll, pitch, and heading are 2 • , 2 • , and 4 • , respectively, and moving acceleration is heading north. Figure 4 plots P d according to the magnitude of f s and shows that P d increases as f s increases. The big arrow in cyan color in the upper left corner of Figure 4 denotes the moving acceleration f acc and the narrow arrow in red color denotes the spoofing acceleration f s . Figure 4a shows the case of both f acc and f s heading north, and plots four curves, two pink in color and two black in color, where the two pink curves show P d,mag and the two black curves show P d,absN . The two black curves are same as one, implying that two cases of f acc = 0.1 m/s 2 and f acc = 0.2 m/s 2 do not cause any effect on P d,absN , while the two pink curves show different results. As f acc changes from 0.1 m/s 2 to 0.2 m/s 2 , the performance of using z mag , which is P d,mag , deteriorates. Figure 4a shows that the performance of using z absN is always better than that of using z mag when both f acc and f s head north. Figure 4b shows similar results as Figure 4a when f s heads 20 • east from north.
The reason the two black curves are the same in Figure 4 despite two different f acc s is because the variance σ 2 n in (10) does not contain f acc,n , but f acc,e . When f acc heads north, f acc,e = 0 and thus for a different north speed, the threshold γ absN is same, which provides the same P d,absN . Figure 4a shows that the performance of using is always better than that of using when both and head north. Figure 4b shows similar results as Figure 4a when heads 20° east from north.
The reason the two black curves are the same in Figure 4 despite two different | |s is because the variance 2 in (10) does not contain , , but , . When heads north, , = 0 and thus for a different north speed, the threshold is same, which provides the same , .

Performance of Spoofing Detection According to the Direction of Moving Acceleration
In this subsection, the effect of the direction of moving acceleration on the spoofing detection performance is presented. When the direction of f acc changes from north to northeast of 30 • , Figure 4a changes into Figure 5, where two black curves are distinct. Since the direction of f acc changes from north to northeast of 30 • , the east component f acc,e exists and thus σ 2 n is different for different speeds, which results in a different threshold γ absN and thus different P d,absN . For both z mag and z absN the performance P d deteriorates as f acc increases from 0.1 m/s 2 to 0.2 m/s 2 . Looking into the black curve of f acc = 0.1 m/s 2 carefully in Figures 4a and 5, it is observed that P d,absN with f acc heading north is bigger than P d,absN with f acc heading northeast of 30 • . For both Figures 4 and 5, the performance of P d,absN is better than that of P d,mag . from north to northeast of 30°, the east component , exists and thus 2 is different for different speeds, which results in a different threshold and thus different , . For both and the performance deteriorates as | | increases from 0.1 / 2 to 0.2 / 2 . Looking into the black curve of | | = 0.1 / 2 carefully in Figure 4a and Figure 5, it is observed that , with heading north is bigger than , with heading northeast of 30°. For both Figure  4 and Figure 5, the performance of , is better than that of , .

Effects of Spoofing Acceleration on the Performance of Spoofing Detection
This subsection analyzes the effects of spoofing acceleration on the performance of spoofing detection. The effects of the direction of the spoofing acceleration are examined for both decision variables and . Figure 6 shows the spoofing detection probability when spoofing direction changes from 0° to 40° from north in the case of heading north as (a) and heading 30° east from north as (b). The horizontal axis is the magnitude of spoofing acceleration. It shows that , does not depend on the direction of spoofing acceleration while , decreases as the spoofing direction changes from 0° to 40° from north. This is because the north component of spoofing acceleration decreases as the spoofing direction gets far away from the north. It is observed that , is greater than , when the direction of is less than 20° for both Figure 6a and 6b.

Effects of Spoofing Acceleration on the Performance of Spoofing Detection
This subsection analyzes the effects of spoofing acceleration on the performance of spoofing detection. The effects of the direction of the spoofing acceleration are examined for both decision variables z mag and z absN . Figure 6 shows the spoofing detection probability when spoofing direction changes from 0 • to 40 • from north in the case of f acc heading north as (a) and heading 30 • east from north as (b). The horizontal axis is the magnitude of spoofing acceleration. It shows that P d,mag does not depend on the direction of spoofing acceleration while P d,absN decreases as the spoofing direction changes from 0 • to 40 • from north. This is because the north component of spoofing acceleration decreases as the spoofing direction gets far away from the north. It is observed that P d,absN is greater than P d,mag when the direction of f s is less than 20 • for both Figure 6a

Effects of Sensor Accuracy on the Performance of Spoofing Detection
This subsection analyzes the effects of sensor accuracy on the performance of spoofing detection. Figure 7a considers AHRS accuracies of 2°, 2°, and 4° for the roll, pitch, and heading errors while Figure 7b considers AHRS accuracies of 1°, 1°, and 2°. It is observed that , (pink color) decreases as the magnitude of increases or AHRS accuracy deteriorates. When both and head north, , is greater than , in Figure 7. Better the AHRS accuracy, the better is the spoofing detection performance and this can be explained in Figure 8 which shows that with better AHRS

Effects of Sensor Accuracy on the Performance of Spoofing Detection
This subsection analyzes the effects of sensor accuracy on the performance of spoofing detection. Figure 7a considers AHRS accuracies of 2 • , 2 • , and 4 • for the roll, pitch, and heading errors while Figure 7b considers AHRS accuracies of 1 • , 1 • , and 2 • . It is observed that P d,mag (pink color) decreases as the magnitude of f acc increases or AHRS accuracy deteriorates. When both f acc and f s head north, P d,absN is greater than P d,mag in Figure 7. Better the AHRS accuracy, the better is the spoofing detection performance and this can be explained in Figure 8 which shows that with better AHRS accuracy, the threshold is smaller, which results in higher detection performance.
(a) (b) Figure 6. Probability of spoofing detection , , and , according to ∠ with | | = 0.6 / 2 , where is red color and is cyan color; (a) is heading north, (b) is heading 30° east from north.

Effects of Sensor Accuracy on the Performance of Spoofing Detection
This subsection analyzes the effects of sensor accuracy on the performance of spoofing detection. Figure 7a considers AHRS accuracies of 2°, 2°, and 4° for the roll, pitch, and heading errors while Figure 7b considers AHRS accuracies of 1°, 1°, and 2°. It is observed that , (pink color) decreases as the magnitude of increases or AHRS accuracy deteriorates. When both and head north, , is greater than , in Figure 7. Better the AHRS accuracy, the better is the spoofing detection performance and this can be explained in Figure 8 which shows that with better AHRS accuracy, the threshold is smaller, which results in higher detection performance.

Performance Analysis of the Decision Variables Using the Detectable Minimum Spoofing Acceleration (DMSA)
This section compares the performance of the proposed direct GPS spoofing detection method for the two decision variables and (or ) using the minimum threshold of spoofing detection with pre-defined false alarm probability and detection probability.

Spoofing Detection Threshold According to Moving Acceleration
The detection threshold γ mag and γ absN according to the direction of for the two cases of | | = 0.2 / 2 and | | = 0.4 / 2 are shown in Figure 3, where two decision variables and are used. Figure 9 shows Figure 3a and 3b together again upon adding the case of | | = 0.6 / 2 under the condition of P =0.001 and the AHRS attitude accuracy of 2/2/4°. Figure 9 shows

Performance Analysis of the Decision Variables Using the Detectable Minimum Spoofing Acceleration (DMSA)
This section compares the performance of the proposed direct GPS spoofing detection method for the two decision variables z mag and z absN (or z absE ) using the minimum threshold of spoofing detection with pre-defined false alarm probability and detection probability.

Spoofing Detection Threshold According to Moving Acceleration
The detection threshold γ mag and γ absN according to the direction of f acc for the two cases of f acc = 0.2 m/s 2 and f acc = 0.4 m/s 2 are shown in Figure 3, where two decision variables z mag and z absN are used. Figure 9 shows Figure 3a,b together again upon adding the case of f acc = 0.6 m/s 2 under the condition of P f a = 0.001 and the AHRS attitude accuracy of 2/2/4 • . Figure 9 shows the exact threshold γ mag and γ absN and is obtained using Equations (21) and (22) given the pre-defined p f a .

Acceleration (DMSA)
This section compares the performance of the proposed direct GPS spoofing detection method for the two decision variables and (or ) using the minimum threshold of spoofing detection with pre-defined false alarm probability and detection probability.

Spoofing Detection Threshold According to Moving Acceleration
The detection threshold γ mag and γ absN according to the direction of for the two cases of | | = 0.2 / 2 and | | = 0.4 / 2 are shown in Figure 3, where two decision variables and are used. Figure 9 shows Figure 3a and 3b together again upon adding the case of | | = 0.6 / 2 under the condition of P fa =0.001 and the AHRS attitude accuracy of 2/2/4°. Figure 9 shows the exact threshold γ mag and γ absN and is obtained using Equations (21) and (22) given the predefined . Figure 9. The threshold of spoofing detection γ mag and γ absN according to , the angle denotes the direction of , i.e., ∠ .

Definition of DMSA
This section defines the detectable minimum spoofing acceleration (DMSA) and compares the DMSA for the decision variables z mag and z absN .
DMSA is the magnitude of the minimum spoofing acceleration to obtain a pre-defined detection probability P d given a pre-defined false alarm probability P f a . P f a = 0.001, P d = 0.99, and AHRS attitude accuracy of 2/2/4 • are used for DMSA in the simulations. Figure 10 shows an example of the computation results of DMSA, where the big arrow in cyan color is the moving acceleration f acc and the angles of 0 through 360 • denote the angle of spoofing acceleration f s . The contour of red '+' is the set of DMSA for all directions of f s . For example, the 'x' point means that when f s comes from 330 • direction, the DMSA is the distance from the origin to 'x' point, which guarantees P d = 0.99.

Contour of DMSA Using The Decision Variable z mag Depending on The Moving Acceleration
For a given DMSA in Figure 10, the values of P f a , P d , and AHRS attitude accuracy are fixed in advance, and thus the moving acceleration is the only remaining parameter that can affect the DMSA.  Figure 11c appears like a circle since the north and east component of f acc are same and thus σ n = σ e . To check the effect of the AHRS accuracy, the DMSA is calculated for two sets of AHRS attitude accuracy of 1/1/2 • and 2/2/4 • in Figure 11d which shows that the better the accuracy, the smaller the contour of DMSA.
DMSA is the magnitude of the minimum spoofing acceleration to obtain a pre-defined detection probability P d given a pre-defined false alarm probability P fa . P fa = 0.001, P d = 0.99, and AHRS attitude accuracy of 2/2/4° are used for DMSA in the simulations. Figure 10 shows an example of the computation results of DMSA, where the big arrow in cyan color is the moving acceleration and the angles of 0 through 360° denote the angle of spoofing acceleration . The contour of red '+' is the set of DMSA for all directions of . For example, the 'x' point means that when comes from 330° direction, the DMSA is the distance from the origin to 'x' point, which guarantees P d = 0.99.

Contour of DMSA Using The Decision Variable
Depending on The Moving Acceleration For a given DMSA in Figure 10, the values of P fa , P d , and AHRS attitude accuracy are fixed in advance, and thus the moving acceleration is the only remaining parameter that can affect the DMSA. To check the effect of the AHRS accuracy, the DMSA is calculated for two sets of AHRS attitude accuracy of 1/1/2° and 2/2/4° in Figure 11d which shows that the better the accuracy, the smaller the contour of DMSA.  Figure 12 shows the contour of DMSA for ∠ = 0°. When heads north, the variance does not depend on | | as shown in Figure 9, so the threshold and DMSA are the same for different

Contour of DMSA Using The Decision Variable Z absN Depending on the Moving Acceleration
This subsection shows the contour of DMSA for the decision variable Z absN = |z n |. Figure 12 shows the contour of DMSA for ∠ f acc = 0 • . When f acc heads north, the variance σ n does not depend on | f acc as shown in Figure 9, so the threshold and DMSA are the same for different magnitudes of f acc as in Figure 12a,b. The contour of DMSA is a line passing through the minimum point of the north direction. Figure 13 shows the contour of DMSA for ∠ f acc = 30 • . In this case, the east component of f acc has an effect on the σ n , which results in a different threshold and DMSA according to different | f acc as in Figure 13. The bigger the acceleration, the bigger is the DMSA.

Optimal Combined Contour of DMSA Using Both and (or )
When both decision variables and (or z absE ) are used, the optimal combined contour can be obtained by combining Figure 11 and Figure 12. Here, the optimal combined contour, colored in pink, is the inner most combined contour from the two contours of DMSA using and . Figure 14a shows the optimal combined contour in the case of | | = 0.2m/s 2 and ∠ = 0°. When the magnitude is increased to | | = 0.4m/s 2 while maintaining the direction, Figure 14b shows that the threshold γ mag is increased, but the threshold γ absN does not change. In Figure 14a

Optimal Combined Contour of DMSA Using Both and (or )
When both decision variables and (or z absE ) are used, the optimal combined contour can be obtained by combining Figure 11 and Figure 12. Here, the optimal combined contour, colored in pink, is the inner most combined contour from the two contours of DMSA using and . Figure 14a shows the optimal combined contour in the case of | | = 0.2m/s 2 and ∠ = 0°. When the magnitude is increased to | | = 0.4m/s 2 while maintaining the direction, Figure 14b shows that the threshold γ mag is increased, but the threshold γ absN does not change. In Figure 14a

Optimal Combined
Contour of DMSA Using Both z mag and z absN (or z absE ) When both decision variables z mag and z absN (or z absN ) are used, the optimal combined contour can be obtained by combining Figures 11 and 12. Here, the optimal combined contour, colored in pink, is the inner most combined contour from the two contours of DMSA using z mag and z absN . Figure 14a shows the optimal combined contour in the case of | f acc = 0.2 m/s 2 and ∠ f acc = 0 • . When the magnitude is increased to | f acc = 0.4 m/s 2 while maintaining the direction, Figure 14b shows that the threshold γ mag is increased, but the threshold γ absN does not change. In Figure 14a,b, it is observed that when f acc heads north, as | f acc increases, the difference γ mag -γ absN becomes larger and the range of angles where DMSA absN < DMSA mag holds, becomes larger. Figure 14c,d show the case of ∠ f acc = 30 • and ∠ f acc = 90 • , respectively. When f acc heads east, i.e., ∠ f acc = 90 • as in Figure 14d, z absE and γ absN should be used instead of z absN and γ absN .  Figure 15 shows the combined contour of DMSA using both and according to the AHRS attitude accuracy. It shows that as the attitude accuracy enhances, the combined contour of DMSA shrinks.  Figure 15 shows the combined contour of DMSA using both z mag and z absN according to the AHRS attitude accuracy. It shows that as the attitude accuracy enhances, the combined contour of DMSA shrinks.  Figure 16a and Figure 16b show the case of ∠ = 30° and ∠ = 120°, respectively, with | | = 0.2, 0.4. ,0.6 m/ 2 . The decision variable is used for ∠ = 30° and is used for ∠ = 120°. It is observed that as | | increases, the range of angles using (or ) becomes larger.  Figure 17 shows the case of ∠ = 45° and the decision variables , and are all necessary to obtain the optimal combined contour.  Figure 16a,b show the case of ∠ f acc = 30 • and ∠ f acc = 120 • , respectively, with | f acc = 0.2, 0.4, 0.6 m/s 2 . The decision variable z absN is used for ∠ f acc = 30 • and z absE is used for ∠ f acc = 120 • . It is observed that as | f acc increases, the range of angles using z absN (or z absE ) becomes larger.  Figure 16a and Figure 16b show the case of ∠ = 30° and ∠ = 120°, respectively, with | | = 0.2, 0.4. ,0.6 m/ 2 . The decision variable is used for ∠ = 30° and is used for ∠ = 120°. It is observed that as | | increases, the range of angles using (or ) becomes larger.  Figure 17 shows the case of ∠ = 45° and the decision variables , and are all necessary to obtain the optimal combined contour.  Figure 17 shows the case of ∠ f acc = 45 • and the decision variables z mag , z absN and z absE are all necessary to obtain the optimal combined contour.

Conclusions
In this study, a direct GPS spoofing detection method is proposed, with AHRS and accelerometers via the difference of the acceleration estimated from GPS receiver and the acceleration measured from IMU. From the acceleration error expressed in the navigation frame, two decision variables are defined for spoofing detection. One decision variable , which may be commonly used, is defined as the magnitude of the horizontal acceleration error. The other decision variable (or ) is defined as the magnitude of the north (or east) component of the acceleration error. The spoofing detection performance can be evaluated using the detection probability, which can be calculated from the probability density function of both decision variables. The decision variable shows higher detection probability than in the condition that both moving acceleration and spoofing acceleration are heading within roughly 25° from the north or south. Similarly, the decision variable shows higher detection probability than in the condition that both moving acceleration and spoofing acceleration are heading within roughly 25° from the east or west.
When detectable minimum spoofing acceleration (DMSA) is used, the decision variable (or ) shows smaller DMSA than in the condition that both moving acceleration and spoofing acceleration head are within roughly 25° from the north-south direction (or east-west direction).
The spoofing acceleration can happen to be any direction. Thus, given a pre-defined false alarm probability, the best algorithm to detect GPS spoofing is that the three decision variables , , and are calculated and compared with the corresponding threshold, and declare the existence of the GPS spoofing if any of the three decision variables exceed the corresponding threshold.
The proposed GPS spoofing detection method in this paper depends on the acceleration error. If a ground vehicle runs across road irregularities such as potholes, bumps, and rubble, etc., then accelerometers may show large changes and deteriorate the GPS spoofing detection performance. Therefore, the flying or driving environment may have an effect on the GPS detection performance.

Conclusions
In this study, a direct GPS spoofing detection method is proposed, with AHRS and accelerometers via the difference of the acceleration estimated from GPS receiver and the acceleration measured from IMU. From the acceleration error expressed in the navigation frame, two decision variables are defined for spoofing detection. One decision variable z mag , which may be commonly used, is defined as the magnitude of the horizontal acceleration error. The other decision variable z absN (or z absE ) is defined as the magnitude of the north (or east) component of the acceleration error.
The spoofing detection performance can be evaluated using the detection probability, which can be calculated from the probability density function of both decision variables. The decision variable z absN shows higher detection probability than z mag in the condition that both moving acceleration and spoofing acceleration are heading within roughly 25 • from the north or south. Similarly, the decision variable z absE shows higher detection probability than z mag in the condition that both moving acceleration and spoofing acceleration are heading within roughly 25 • from the east or west.
When detectable minimum spoofing acceleration (DMSA) is used, the decision variable z absN (or z absE ) shows smaller DMSA than z mag in the condition that both moving acceleration and spoofing acceleration head are within roughly 25 • from the north-south direction (or east-west direction).
The spoofing acceleration can happen to be any direction. Thus, given a pre-defined false alarm probability, the best algorithm to detect GPS spoofing is that the three decision variables z mag , z absN , and z absE are calculated and compared with the corresponding threshold, and declare the existence of the GPS spoofing if any of the three decision variables exceed the corresponding threshold.
The proposed GPS spoofing detection method in this paper depends on the acceleration error. If a ground vehicle runs across road irregularities such as potholes, bumps, and rubble, etc., then accelerometers may show large changes and deteriorate the GPS spoofing detection performance. Therefore, the flying or driving environment may have an effect on the GPS detection performance.