Theory of Automatic Dependent Surveillance–Broadcast Position Verification Using Time Difference of Arrival

Automatic dependent surveillance–broadcast (ADS-B) is an emerging means of aeronautical surveillance. Because ADS-B is susceptible to false aircraft information, position verification is essential. This article proposes a theoretical performance model for ADS-B position verification using the time difference of arrival (TDOA). The model was derived by considering the effects of receiver–aircraft geometry, TDOA accuracy and outliers, latency, self-localization accuracy, and detection threshold. The model was verified through a numerical simulation, flight experiment, and spoofing experiment.


I. INTRODUCTION
Automatic dependent surveillance-broadcast (ADS-B) is an emerging means of aeronautical surveillance for air traffic control (ATC) where aircraft periodically report their positions. Compared with conventional means of aeronautical surveillance such as secondary surveillance radar (SSR), ADS-B can provide a better positional accuracy and update rate; this is essential to support advanced operations in ATC to increase safety and capacity [1].
However, a crucial drawback of ADS-B is that anomalous position reports, due to avionics failure or spoofing, may disrupt its operation and, consequently, ATC services [2]- [7]. Various countermeasures have been proposed against these security issues [6], [7]. One solution is enhancing the signal format and protocol for authentication, encryption, and position verification [8]- [14]. However, this solution requires updating the avionics, which places a large cost burden on aircraft operators. Another solution is combining ADS-B with SSR or multilateration [15]- [21]. However, this solution does not allow air navigation service providers to fully enjoy the cost benefits of ADS-B for ground infrastructure.
On the other hand, position verification techniques that use radiofrequency (RF) properties [21]- [33] are a promising approach to achieve cost benefits for both aircraft operators and air navigation service providers. Position verification techniques are also important as a countermeasure by ground systems to protect the subsequent information processing, where additional countermeasures, e.g., use of flight plan [30], can be additionally implemented.
RF properties that are widely being considered for position verification include the time difference of arrival (TDOA) [21], [30]- [33], received signal strength [21], [24], angle of arrival [21], [22], [27] and its derivative [23], and features extracted from base-band or intermediatefrequency signals [26], [29]. TDOA is especially promising because techniques for accurately measuring the TDOA of an ADS-B signal are readily available from an existing system called wide-area multilateration [34]- [37]. However, the performance of a TDOA-based method has not been fully investigated in the open literature.
Graziano et al. [31] proposed an improved geometric calculation that they tested through numerical simulations and measurements. Strohmeier et al. [32] proposed a lightweight statistical method that they tested with measurement data. SESAR Joint Undertaking [21] summarized the prototype implementations of various methods. Neufeldt [30] discussed the general idea for ADS-B security from the manufacturer's point of view. Leonardi [33] used TDOA with a Kalman filter and clock-offset model to estimate the clock offset, which was used to detect spoofing. However, these works did not consider or verify the complete theoretical performance in terms of both probabilities of detection and a false alarm. In addition, the null hypothesis adopted by Leonardi [33], whose work is the closest to the present study, can be improved by incorporating the latency and correlation of the GPS position error, and it must be verified through measurements. In the field of navigation, countermeasures against GPS spoofing are studied [38], but the purpose is different from the ADS-B position verification. The former considers a signal used by the aircraft to obtain its position, whereas the latter considers a signal used by the aircraft to report its position.
The lack of a complete and verified performance model has hindered system design, risk assessment, and policymaking. The performance model is also important for clarifying the underlying characteristics and limitations of TDOA-based methods for future development. Therefore, this article proposes a theoretical performance model for TDOA-based ADS-B position verification.
The threats considered in this article are spoofing, which are transmissions of ADS-B signals by an adversarial user, and avionics failures that transmit wrong position reports. In the both cases, the reported position is likely to differ from the true emitter position, which is detected by the TDOA. The statistical distribution of the test statistic, probability of detection, and probability of a false alarm were derived. The model explicitly describes the effects of the receiver-aircraft geometry, TDOA accuracy and outliers, latency, self-localization accuracy, and detection threshold on the performance. A threshold that guarantees a predefined probability of a false alarm was also derived. The derived model was verified through a numerical simulation, flight experiment, and spoofing experiment. In the flight experiment, the verification was based on comparing the test statistic and false alarms with the true target position. To the best of the authors' knowledge, the results of such a flight experiment are not available in the literature.
This article is an extension of previous work [39]- [41] with more realistic assumptions for the derivation, new measurement results (flight experiment), and other updates.
The rest of the article is organized as follows. Section II derives a basic model under simplified assumptions. Section III derives an advanced model under more realistic assumptions. Section IV derives further advanced model considering outliers. Section V presents the numerical simulation. Section VI presents the flight experiment. Section VII presents the spoofing experiment. Section VIII concludes the article and discusses possible applications of the derived model.

A. Notation
H 0 and H 1 denote a null hypothesis and an alternative hypothesis, respectively. x ∼ N (μ, σ 2 ) means that the variable x follows a normal distribution with the mean μ and variance σ 2 . (x) and¯ (x) denote the cumulative distribution function (cdf) and the complementary cdf (CCDF) of x, respectively. By definition,¯ (x) = 1 − (x). (x; μ, σ 2 ) and¯ (x; μ, σ 2 ) denote the cdf and CCDF with the parameters μ and σ , respectively. (x; H) and¯ (x; H) denote the cdf and CCDF under the hypothesis H, respectively. Pr(·) denotes the probability when an event specified by the argument occurs. Q(x) and Q −1 (x) denote the CCDF of the standard normal distribution and its inverse function, respectively. The Q-function enables numerical calculation of the model by¯ The following equation will be often used: x means the measurement or estimation of x. Max[·] denotes an operation to find the largest value.

II. BASIC MODEL
In this section, a basic model was derived under simplified assumptions: 1) ADS-B from a legitimate aircraft has no error and 2) receivers are perfectly synchronized such that the measured TDOA is nonbiased. The decision rule and the assumptions are introduced in Section II-A. Then, the distribution of the test statistic is derived in Section II-B. This allows deriving probability of detection, probability of a false alarm, and detection threshold for a constant false alarm, which is presented in Section II-C.

A. Decision Rule and Assumptions
Let H 0 and H 1 are the hypotheses where the position in an ADS-B message is valid and an anomaly, respectively. Let l = [x l , y l , z l ] T be the position of the emitter, which is an aircraft under H 0 or spoofing transmitter under H 1 . Two receivers identified by the index i = {1, 2} receive an ADS-B signal and measure its time of arrival (TOA) t i,o . The two TOAs produce a TDOA by t o = t 2,o − t 1,o , which is the measured value. The position is verified by comparing the measured TDOA t o with the predicted TDOA t e , which is available from the ADS-B position and receiver positions. This test is written as follows: where γ is the detection threshold. t e is obtained as follows. Let r i = [x i , y i , z i ] T be the ith receiver position and l = [x l , y l , z l ] T be the ADS-B position. The distance from the ith receiver is defined by the following function: The difference between the distances from the first receiver (i = 1) and from the second receiver (i = 2) is defined by the following function: Substituting the ADS-B position into (5) gives t e as follows: where c is the speed of light. In the following, assumptions and model parameters are introduced. Most of the parameters can be measured by the method described in Section VI. Different ways for obtaining a parameter will be additionally mentioned if available. ASSUMPTION 1 The TOA measurement error at the ith receiver, denoted by i , follows an uncorrelated zero-mean Gaussian distribution with the variance σ 2 t : i ∼ N (0, σ 2 t ).
The Gaussian distribution was selected based on the previous study, where majority of the TDOA was well characterized by a Gaussian distribution [37]. Outliers will be considered in Section IV. The zero-mean assumption means that the receivers are perfectly synchronized. σ t is mainly determined by the TOA measurement performance of the receiver. Also, the literature and product catalogs can be referred for a nominal value. Assumption 3 was made because an anomalous ADS-B position is subject to the intent of the spoofing attacker or avionics failure and the statistical distribution of l is difficult to be expressed in the parameterized closed form. Instead, the model can be calculated for various values of l . For example, in Sections V and VII, l anom is discretely uniform in a specified area, and the probability of detection is calculated for each position.

B. Test Statistic
Distribution of the test statistic T in (3) is derived here. Let t i be the true TOA at the ith receiver. The true TDOA can be written with t 21 = t 2 − t 1 and as follows: The measured TOA is given by the sum of the true TOA and the noise i Then, (8) and (7) can be used to obtain Substituting (9) and (6) into (3) obtains The term ( 2 − 1 ) means a TDOA measurement error. From Assumption 1, it is characterized as follows: Under H 0 , the reported position and emitter position perfectly agree (see Assumption 2), which removes the first term in (10). Thus, substituting (11) into (10) obtains the following distribution: Under H 1 , substituting l anom in Assumption 3 and (11) into (10) obtains the following distribution:

C. Position Verification Performance
The threshold γ can be applied to the test statistic distribution in (12) and (13) to derive the position verification performance. First, the probability of a false alarm is given as follows: The detection threshold that achieves a constant false alarm is given as follows: The probability of detection is given as follows:

III. ADVANCED MODEL
This section extends the basic model in Section II to present an advanced model derived under more realistic assumptions: 1) each receiver has a bias in time synchronization and 2) a legitimate ADS-B position contains a small error caused by latency and self-localization error. The assumptions are introduced in Section III-A. Then, the distribution of the test statistic, probability of detection, and probability of a false alarm are derived in Sections III-B and III-C. By leveraging the model, two approaches to determining the threshold are presented in Sections III-D and III-E.

A. Assumptions
To consider imperfect synchronization among the receivers, Assumption 1 is replaced by the following assumption. ASSUMPTION 4 The TOA measurement error i follows an uncorrelated Gaussian distribution with the mean μ t,i and the variance σ 2 t : i ∼ N (μ t,i , σ 2 t ). μ t,i and its range are determined by the synchronization state and performance of the receiver, respectively. μ t,i is a constant, which is an approximation that the temporal change is slow enough compared to the time-duration of interest, e.g., one measurement or one flight. In the case where the temporal variation is not negligible, e.g., very long measurement, the model should be applied to segmented time durations.
To consider a small error in a legitimate ADS-B position, let = l − l be the position error comprising the elements = [ x , y , z ] T . Then, Assumption 2 is replaced by the following Assumptions 5-8.

ASSUMPTION 5 The error in an ADS-B position is modeled as the sum of the latency effect late and self-localization error loc
The latency is the elapsed time from when an aircraft measures its position to when it transmits an ADS-B message. The self-localization error is the difference between the true position and position estimated by the onboard navigation system. ASSUMPTION 6 The latency l follows a normal distribution with the mean μ l and variance σ 2 l : Although a more sophisticated model of latency is known [42], the normal distribution was adopted here to achieve analyticity. μ l and σ l are determined by the avionics configuration. σ l can be measured by the method by Kakubari et al. [46].

ASSUMPTION 7 The velocity vector of the aircraft
This constant assumption was made because velocity does not have random nature like the other stochastic terms, e.g., TOA. A velocity usually changes slowly according to the state of the flight. v can be set by prior knowledge of typical flight conditions, the velocity squitter, and stateestimation techniques.
ASSUMPTION 8 The self-localization error loc is modeled by a correlated Gaussian distribution with the mean μ loc and the covariance matrix loc loc ∼ N (μ loc , loc ) (19) μ loc and loc are determined by the state of the avionics, mainly the navigation equipment. The ranges of μ loc and loc are determined by the performance of the avionics. loc can be set by referring to accuracy measures such as navigation uncertainty category (NUC) [47] or performance of the navigation system, e.g. [48].
The next assumption is for relaxing the requirement on the knowledge of the model parameters. Obtaining exact values for the model parameters are often difficult. Even in such case, the proposed model can be used for calculating a threshold that satisfies a false alarm requirement. For this purpose, the upper bounds (worse case value) of the parameters are introduced instead of exact values as follows. ASSUMPTION 9 The following upper bounds of model

B. Test Statistic
The test statistic is derived under the assumptions. First, the effect of is considered by employing the linearization technique in TDOA-based localization [43]. When is small, the first term in (10) is approximated by a linear term to yield where

is the coefficient given by
In (30), A is calculated for a known emitter position. Alternatively, the ADS-B position can also be used for calculation The choice of (30) or (31) depends on the situation. For example, (30) can be used during the system design phase where a specific emitter position is assumed, whereas (31) can be used in a running system. From Assumption 5, (29) is written as follows: The latency term can be expressed by the velocity vector v and latency l as follows: Eq. (33) and Assumption 6 yield For the self-localization error, Assumption 8 yields For the TDOA measurement error, Assumption 4 yields Substituting (35)-(37) into (32) obtains the following distribution of the test statistic under H 0 : Comparing (38) with (12) shows that an erroneous ADS-B position produces a nonzero bias and larger spread in the test statistic. In addition, the distribution varies depending on the aircraft position because of A.
Next, substituting (37) and l anom into (10) obtains the following distribution of the test statistic under H 1 :

C. Position Verification Performance
Applying the threshold γ to the test statistic distribution in (38) and (39) obtains the position verification performance. First, the probability of a false alarm is given as follows: Next, the probability of detection is given as follows:

D. Threshold I: Empirical Method
Unlike the basic model in the previous section, solving (40) for γ is difficult. One alternative to determining the threshold is an empirical method. The test statistic T is measured in a real environment; the results are then used to empirically determine γ such that a constant false alarm is achieved. One advantage of this method is simplicity. Only one parameter is needed, which is the required probability of a false alarm. However, a disadvantage is that a physical interpretation of the obtained threshold cannot be given.
The threshold determined by this method is a fixed value that is applied to any ADS-B message. Therefore, it is referred to as a fixed threshold hereafter to distinguish it from the next threshold.

E. Threshold II: Upper-Bound-Based Method
Instead of trying to achieve a constant false alarm, the other approach to determining γ is to satisfy a predetermined upper bound using the proposed model. As derived in (38), the distribution of the test statistic varies depending on the aircraft position. Consequently, the threshold determined by this method is a dynamic value depending on the aircraft position. Therefore, this threshold is referred to as a dynamic threshold hereafter to distinguish it from the fixed threshold described in Section III-D. Also, this method relaxes the requirement on the knowledge of the model parameters; only the upper bounds are sufficient. The derivation is described ahead.
First, an upper bound of μ 0 is given by the following lemma.
where μ 0,u is an upper bound of μ 0 given by PROOF The proof is provided in Appendix A An upper bound of σ 0 is given by the following lemma.
LEMMA 2 where σ 0,u is an upper bound of σ 0 given by PROOF The proof is provided in Appendix B.
LEMMA 3 Assume P FA is given by (40) with the distribution parameters μ 0 and σ 0 that are upper bounded by |μ 0 | ≤ μ 0,u and σ 2 0 ≤ σ 2 0,u . The threshold γ that satisfies P FA ≤ P FA,u is given by PROOF The proof is provided in Appendix C.

IV. ADVANCED MODEL II
This section extends the advanced model in Section III to another advanced model that considers outliers in the TDOA. To do so, a mixture of Gaussian distributions is employed. The assumptions are introduced in Section IV-A. Then, the distribution of the test statistic, probability of detection, and probability of a false alarm are derived in Section IV-B. Approaches to determining the threshold are discussed in Section IV-C.
A. Assumptions ASSUMPTION 10 TDOA outliers are caused by TOA outliers. The TOA outlier follows a mixture of M Gaussian distributions with the mean μ t,m,i , variance σ t,m , and weight w m for the mth distribution: i ∼ N (μ t,m,i , σ 2 t,m ). Possible sources of the outliers are signal collision, multipath, and receiver noise. Therefore, μ t,m,i , σ t,m , and w m are mainly determined by the surrounding environment (signal rate, antenna location, and scatterers) and the receiver performance. μ t,m,i is affected also by the synchronization state of the receiver.
For the same purpose as Assumption 9, the following assumption is introduced to relax the requirement on the knowledge of the outlier parameters. ASSUMPTION 11 The following upper bounds of model parameters-μ t,m,u and σ t,m,u -are available: To treat the nonoutlier and the outlier in a unified manner, the following assumption is introduced. ASSUMPTION 12 m = 0 denotes the nonoutlier, which was modeled in Section III. That is ASSUMPTION 13 The weight w m satisfies M m=0 w m = 1.

B. Test Statistic and Position Verification Performance
The derivation can be achieved by leveraging the results in Section III. To do so, Assumption 4 is replaced by Assumption 10. Consequently, μ t,m,i and σ t,m replace μ t,i and σ t in (38), which yields the test statistic for the mth nonoutlier/outlier component. Let denote the resulting μ 0 , σ 0 , and cdf by μ 0,m , σ 0,m , and (x; μ 0,m , σ 2 0,m ), respectively. Then, the test statistic distribution under H 0 is given in the form of a cdf as follows: Similarly, μ t,m,i and σ t,m replace μ t,i and σ t in (39). Let denote the resulting μ 1 , σ 1 , and cdf by μ 1,m , σ 1,m , and (x; μ 1,m , σ 2 1,m ), respectively. This yields the distribution under H 1 as follows: Applying the threshold γ to the test statistic distributions in (53) obtains the probability of a false alarm as follows: (56)

C. Threshold
In Section III, the two approaches for determining the threshold were presented. The empirical method (the fixed threshold) is directly applicable with the model presented in this section. On the other hand, the upper-bound-based method (the dynamic threshold) needs to be modified accordingly.
To do so, the results in Section III-E are applied to each of the nonoutlier and M outliers separately. Eqs. (25) and (28) in Assumption 9 are replaced by Assumption 11. P FA of (40) is replaced by P FA,m of (55). Consequently, μ t,m,u and σ t,m,u replace μ t,u and σ t,u in (43), (45), and (46). The resulting threshold is denoted by γ u,m . From Lemmas 1-3, γ u,m satisfies Then, the most conservative value is selected based on the following lemma. LEMMA 4 P FA ≤ P FA,u is satisfied by PROOF The proof is provided in Appendix D.

V. NUMERICAL SIMULATION
A numerical simulation was conducted to verify the derived model. Two types of simulation were conducted: for a false alarm and for detection, which are presented in Sections V-A and V-B, respectively.
A. False Alarm Simulation 1) Setup: For the false alarm simulation, two receivers and a legitimate aircraft were assumed. A Monte Carlo simulation was conducted; an ADS-B position and TOA measurements were randomly generated, which were used for the TDOA-based position verification. The resulting probability of a false alarm was evaluated according to where N H 0 is the number of signals generated and N H 1 ;H 0 is the number of signals that were determined as anomalous among the generated signals. The parameters were substituted into (53) and (55) to obtain the theoretical distribution of the test statistic and the probability of a false alarm. The parameters are given in Table I. They were determined to represent the following practical scenario for ADS-B.
1) M = 1 was selected as the simplest model that can consider outliers. m = 0 and m = 1 were associated with the nonoutliers and outliers, respectively. 2) w m , σ t,m , μ t,m,i , and μ t,m,u were determined based on the flight experiment in Section VI. 3) μ l was determined according to Federal Aviation Authority (FAA) regulations [44], [45], which require the ADS-B latency to be below 0.6 s. σ l was determined from the measurement results of Kakubari et al. [46]. 4) μ loc was determined based on the flight experiment in Section VI. 5) σ x and σ y were determined from the navigation uncertainty category (NUC) of 11, which is a typical NUC value. σ z was determined from σ x , σ y , and a horizontal dilution of precision of 0.83 and a vertical dilution of precision of 1.35, which were reported by Renfro et al. [48] to be the average values in 2017. Correlation was assumed among the x, y, z components. 6) v was determined based on a level flight at cruising speed traveling in a direction perpendicular to the hyperbola passing through the true position. 7) The parameters needed for the dynamic threshold have the same values as the actual values (e.g., σ t,m,u = σ t,m ).
The simulation was conducted using MATLAB. 2) Results: Fig. 1(a)-(d) shows the test statistic and probability of a false alarm for different values of the fixed threshold γ . Good agreement was observed between the simulation results (blue dot-dash lines) and theoretical values (solid black lines), which verifies the derived model. Fig. 1(d) also plots γ u calculated by (58) as a vertical red dotted line. At the threshold, the simulation and theoretical values of P FA do not exceed the predetermined value P FA,u = 0.05, which verifies (58). Thus, the derived dynamic threshold can guarantee the desired probability of a false alarm.

1) Setup:
In the detection simulation, the emitter was a spoofing or malfunctioning aircraft. The ADS-B position was uniformly distributed, which created the situation of a denial-of-service attack. This is appropriate for a malfunctioned aircraft because it is difficult to predetermine the ADS-B position in such a situation. A Monte Carlo simulation was conducted for each ADS-B position. TOA measurements were randomly generated for TDOA-based position verification. The resulting probability of detection was evaluated at each ADS-B position as follows: where N H 1 is the number of signals generated and N H 1 ;H 1 is the number of generated signals that were determined to be anomalous. Substituting the parameters into (54) and (56) obtained the theoretical value for the probability of detection.
The parameters are given in Table II. The parameters common to the false alarm simulation had the same values except that the emitter was on the ground (z l = 0). The parameters specific to this detection simulation were determined as follows.
1) The ADS-B position was uniformly generated as follows: y l = y 0 , y 0 + y , y 0 + 2 y , . . . , y 1 (62) Two configurations were considered: a wide area and narrow area. 2) γ u was calculated with the parameters from Table I. 2) Results: Because the probability of detection was 1 at most ADS-B positions, the focus here was on positions with a low probability of detection. Fig. 2 shows the positions at which the probability of detection was below 0.95. A good agreement between the simulation results (red circles) and derived model (blue circle) is verified in (56). Positions with a low probability of detection were distributed in a hyperbola because they occurred when the TDOA was the same for the ADS-B position and emitter. This illustrates a fundamental limitation of the TDOA-based position verification. The average probability of detection for the generated ADS-B positions was also evaluated with where L is the number of the ADS-B positions. The results are given in Table III. A good agreement for the average probability of detection is verified in (56).

VI. FLIGHT EXPERIMENT
A flight experiment was conducted to verify the derived model in terms of the test statistic and false alarm. An experimental aircraft was used as a verified target with a known true position. The TDOA-based position verification was then applied. The setup of the experiment is presented  in Section VI-A. The evaluation of the test statistic and probability of a false alarm is presented in Sections VI-B-VI-D.
A. Setup Fig. 3 shows a schematic of the experiment. The experimental aircraft was a Beechcraft B300 equipped with an onboard global navigation satellite system (GNSS) receiver, and its measured position was assumed to be true. Receivers  were used to measure the ADS-B signals transmitted by the aircraft, and the measured signals were stored by the target processor. After the experiment, the measured signals and onboard GNSS data were processed offline for analysis. Fig. 4 shows the flight path and receiver positions. Fig. 5 shows photos of the receiver and antennas used in the experiment. The receivers were developed by NEC and the Electronic Navigation Research Institute (ENRI) as prototypes for investigating future operational systems in Japan. The receivers were synchronized with the global positioning system (GPS) using the common-view technique. The antennas were developed by Antenna Giken and ENRI. The antenna for Rx A had a horizontal beamwidth of 30 • , whereas that for Rx B had a horizontal beamwidth of 120 • . The position of Rx A was 35. After the measurement, the signals received by the two receivers were paired. Because TDOA measurements are known to contain outliers presumably caused by the multipath and cochannel interference [37], such outliers were removed. Two consecutive TDOAs were compared, and the latter was identified as an outlier when it differed by more than 5000 ns from the former. After the pairing of results and removal of outliers, the following measurements were obtained: the index of the paired signal k, ADS-B positioñ l k , aircraft true positionl k , and TOAt i,o,k .

B. Test Statistic
First, the test statistics were compared for the derived model and measurement. The test statistic was measured as follows:T where the TDOA was measured with Next, the proposed model was calculated according to the following parameters.
1) The term associated with the ADS-B position error in (29) was calculated with Here,˜ k is the measured and is obtained as the difference between the true and ADS-B positions 2) The term associated with the TDOA measurement error, ( 2 − 1 ) in (29), was calculated by subtracting T k,theo from the measurement 21 =T k,mea −T k,theo .
4) M = 1 was selected as the simplest model that can consider the outliers. m = 0 and m = 1 were associated with the nonoutliers and outliers, respectively. Table IV lists estimated parameters. Fig. 6(a) shows the temporal variation of the test statistic. The red points labeled as "theory" indicatesT k,theo . The variation was steep at the beginning before 400 s and at the end after 4100 s, where the aircraft flew near the receivers. The direction of the variation was symmetrical: negative at the beginning but positive at the end. This can be associated with the aircraft traveling direction. The derived model reproduced this trend and showed good agreement with the measurement (i.e., blue circles). Some scattered measurements that did not agree with the theoretical values are outliers. Fig. 6(b) and (c) shows the distribution of the TDOA measurement error term in the test statistic. Because the measured test statistics contained the outliers, the distribution had long tails. The derived model (i.e., dashed blue line) reproduced this trend and showed good agreement with the measurement (i.e., black solid line).

C. Probability of a False Alarm With the Fixed Threshold
The probabilities of a false alarm were compared for the derived model and measurement with a fixed threshold. The theoretical value from the derived model was calculated according to the following parameters. 1) v was estimated from the GPS position data for each k byṽ 2) μ l was estimated by minimizing the residual after subtracting the latency error from the position error as follows:˜ 3) μ loc and loc were estimated as the mean and covariance matrix of˜ k,r , respectively. loc contained the fluctuating components of both the latency and self-localization error. 4) σ l = 0 because the fluctuating component of the latency was already included in loc . 5) M, w m , σ t,m , and μ t,m,i were set to the estimated values in Section VI-B.
Steps 2-4 were compromises made because of the limitations of the flight experiment. Table IV lists estimated parameters. The previous parameters were substituted into the model, and the theoretical probability of a false alarm P FA,theo,k was calculated for the kth signal. Then, the average was taken fromP where K is the number of analyzed signals.
The measured probability of a false alarm P FA,mea was calculated withP where N H 1 ;H 0 ,mea is the number of signals that satisfies |T k,mea | > γ . Fig. 7 comparesP FA,theo (black solid line) and P FA,mea (blue dotted line). The derived model and measurement showed acceptable agreement, which verifies the model.

D. Probability of a False Alarm With the Dynamic Threshold
The probability of a false alarm with a dynamic threshold was examined. The threshold γ u,k was calculated for each k such that the desired probability of a false alarm P FA,u was achieved according to (58). Then, the measured probability of a false alarm (i.e., the achieved value) was obtained byP where N H 1 ;H 0 ,mea was the number of signals that satisfied |T k,mea | > γ u,k . γ u,k was calculated according to the following parameters.
x ,σ 2 y ,σ 2 z }], whereσ 2 x ,σ 2 y , andσ 2 z were the diagonal components of˜ loc . 5) μ x,u = Max[{μ 2 x ,μ 2 y ,μ 2 z }], whereμ x ,μ y , andμ z were the components ofμ loc . The results are plotted as a black solid line in Fig. 8, where the vertical axis is the measurement and the horizontal axis is the predetermined (desired) values. P FA < P FA,u was achieved. However, a large difference between P FA and P FA,u was observed, suggesting that γ u was pessimistic. To identify the reason, the variations inT k,mea and γ u were plotted, as shown in Fig. 9. γ u , γ u,0 , and γ u,1 are plotted by the green, magenta, and black lines, respectively. γ u was dominated by γ u,1 , which represents the contribution of the outliers. The outlier has a larger TOA variance, which resulted in the pessimistic threshold. However, the amount of the outliers is smaller than the nonoutliers, and there is room for lowering γ u,1 . This improvement was suggested as future work.

E. Summary of the Flight Experiment
In summary, the flight experiment led to the following results.
1) The measurement and derived model agreed in terms of the temporal variation and distribution of the test statistic.
2) The measurement and derived model with a fixed threshold showed agreement for the probability of a false alarm.
3) The dynamic threshold was effective, but it was pessimistic and improvement is suggested as future work.

VII. SPOOFING EXPERIMENT
A spoofing experiment was conducted to verify the derived model in terms of detection. In the experiment, false messages that contain a fake aircraft position were generated and mixed with real messages broadcasted from target of opportunity. The TDOA-based position verification was then applied. The setup of the experiment is presented in Section VII-A. The evaluation of the test statistic and probability of detection is presented in Sections VII-B-VII-D.
A. Setup Fig. 10 shows a schematic of the experiment. The assumed scenario had a spoofing emitter, targets of opportunity (all considered legitimate), and two receivers. This scenario was simulated by using two receivers at a local site (denoted by Rx C and Rx D'), a receiver at a remote site (denoted by Rx D), and a signal generator as the spoofing emitter. Rx D received messages from the opportunistic targets. Rx D' received the generated false messages. The messages from Rx D' were treated as if being from Rx D. Thus, Rx D and Rx D' virtually received a combination of opportunistic messages and generated false messages. For Rx C, the opportunistic messages and false messages were mixed by an RF combiner, and the mixed signals were input to the receiver. The received messages were sent to a target processor and stored, and the stored messages were processed offline for evaluation. Fig. 11 shows photos of the receivers and antenna used in the experiment. The receiver hardware for Rx D was the same as that for Rx A and B in the flight experiment. The receiver hardware for Rx D' and Rx C was enhanced. Rx D used an AS-177B/UPX antenna. The signal generator for the false messages comprised the software-defined radio transceiver USRP 2901 and a host computer. The spoofing emitter was assumed to be located at 36.22 Fig. 12(a) shows the spatial distribution of the measured test static. Most of the false messages, which were distributed as a grid, had a test statistic above the scale maximum of 2000 ns. In contrast, the test statistic of the target of opportunity was mostly below 1000 ns. Fig. 12(b) also visualizes the test statistic as a histogram. The bin was highest at 10 5 -10 6 ns and 10 2 -10 3 ns for the false messages and target-of-opportunity messages, respectively. Thus, a significant difference was observed in the value of the test statistic for the false messages and targets of opportunity. This supports the effectiveness the TDOA-based position verification method.

C. Probability of Detection With the Fixed Threshold
First, a fixed threshold was applied. To empirically decide the threshold, a measurement was conducted for 800 s without false messages. Based on the results, a threshold of γ = 985.4 ns obtained a probability of a false alarm of 0.05. 1 The previous result [41] was updated according to this improvement.  where N H 1 is the number of false messages that was processed by the position verification and N H 1 ;H 1 is the number of false messages that were determined as anomalies. With N H 1 ;H 1 = 559 519 and N H 1 = 560 529,P D,mea = 0.9982 was obtained. This performance showed that the position verification is effective as a spoofing countermeasure. Although false messages remained, they should be easily removed by additional methods in the following information processing phase. For example, messages can be tracked or correlated with flight plans. The obtained performance is also sufficient for protecting a system against a denial-of-service attack, which transmits false signals at a high rate to overflow the surveillance system chain. Considering the length of an ADS-B signal, a maximum of 1/120 us = 8333 signals can be sent per second. If the position verification can reject 99.82% of the signals, then 15 signal/s would remain. This amount is ordinary for ADS-B signals and should be below the capacity of a typical system.
Next, the results of the measurement and derived model were compared. The following parameters were used to calculate the derived model. 1) M, w m , σ t,m , and μ t,m,i were set to the estimated values in Section VI. 2) l was set to the emitter position. 3) l anom was set to the position where the false messages were generated.
The previous parameters were substituted into the model, and the theoretical value for the probability of detection P D,theo was obtained for a false position l anom . Then, the average was taken with where L is the number of false positions. As a result, P D,theo = 0.9983 was obtained, which showed good agreement withP D,mea = 0.9982. P D,theo andP D,mea were compared for different values of the threshold. Fig. 14 shows the results; good agreement was observed between the theoretical values (black solid line) The threshold is computable for any three-dimensional position, but it was computed and plotted only for the observed ones here. By this selecting a specific altitude was avoided. and measurements (blue circles). Thus, the derived model was verified experimentally for the fixed threshold case.

D. Probability of Detection With the Dynamic Threshold
The previous section presented the results with a fixed threshold. Here, the results are presented with a dynamic threshold.
The dynamic threshold was calculated for each received message according to the following parameters. The parameters that were the same as in the case of the fixed threshold are omitted. The previous parameters were substituted to obtain the dynamic threshold γ u for each signal. Fig. 15 shows the spatial distribution of γ u . Fig. 16 shows the results when (3) was applied with γ u , where Fig. 16(a) shows the positions determined as anomalous and Fig. 16(b) shows the positions determined as valid. Most of the false messages were successfully determined as anomalous, whereas the target-of-opportunity positions were determined as valid.
In the same manner as the fixed threshold case, the measurements were compared with the derived model. This resulted inP D,theo = 0.9955 andP D,mea = 0.9954, which indicated good agreement. The comparison was repeated for different values of P FA,u . Fig. 17 shows the results; good agreement was observed between the derived model (black solid line) and measurement (blue circles). Thus, the derived model was verified experimentally for a dynamic threshold.  Finally, the probability of a false alarm was evaluated as follows: where N H 0 is the number of target-of-opportunity messages that were processed by the position verification and N H 1 ;H 0 is the number of target-of-opportunity messages that were determined as anomalies. When P FA,u = 0.05, P FA,mea = 8.5 × 10 −4 was obtained with N H 1 ;H 0 = 437 and N H 0 = 512 058. Thus, P FA,mea < P FA,u was satisfied.

E. Summary of the Spoofing Experiment
In summary, the spoofing experiment led to the following results.
1) The generated false messages and target-ofopportunity messages were successfully separated. 2) The measurements and derived model agreed in terms of the probability of detection with both the fixed and dynamic thresholds.

VIII. CONCLUSION
A theoretical performance model was proposed for TDOA-based ADS-B position verification and verified. The statistical distribution of the test statistic, probability of detection, and probability of a false alarm were derived to consider the TDOA accuracy, accuracy of a reported position, and detection threshold. A threshold that guarantees a predefined probability of a false alarm was also derived. Furthermore, outliers were included in the model by using a mixture of Gaussian distributions. The derived model was verified through a numerical simulation, flight experiment, and spoofing experiment.
The proposed model has several possible applications. First, the detection logic and threshold in this article can be employed in operational systems used by an air navigation service provider. The proposed logic and threshold have the advantages of a verified and predictable performance. Second, it can be applied to system and coverage design. In the design phase of an operational system, parameters such as the receiver position, coverage, and TDOA accuracy (which is related to the cost and performance of a receiver) can be determined by referring to the predicted performance with the proposed model. Because the proposed model is based on a relatively simple detection logic, it can be interpreted as the baseline performance of other advanced TDOA-based methods. Third, it can be applied to risk assessment and policymaking. The probability of detection can be interpreted as a measure of risk where anomalous positions are missed. This value can be referenced when developing an implementation plan for ADS-B.
Although the proposed model is useful for the aforementioned applications, it has some limitations. First, the theoretical threshold determined by the model tends to be pessimistic due to the outliers. Second, it is assumed that the model verifies each signal independently; a statistical test using multiple signals or in combination with a tracking technique will improve the performance.
In future work, the proposed model will be applied in a Japanese surveillance environment to derive insights into the operation implementation plan. In addition, improvements in the model or detection logic will be investigated.