Nanosecond-precision Time-of-Arrival Estimation for Aircraft Signals with low-cost SDR Receivers

Nanosecond-precision Time-of- Arrival Estimation for Aircra Signals with lo w-cost SDR Receivers Roberto Calvo-Palomino IMDEA Networks Institute, Spain University Carlos III of Madrid roberto.calvo@imdea.org Fabio Ricciato University of Ljubljana, Slovenia fabio.ricciato@fri.uni- lj.si Blaz Repas University of Ljubljana, Slovenia br9404@student.uni- lj.si Domenico Giustiniano IMDEA Networks Institute, Spain domenico.giustiniano@imdea.org Vincent Lenders armasuisse, Switzerland vincent.lenders@armasuisse.ch ABSTRA CT Precise Time-of- Arrival (TOA) estimations of aircraft and drone signals are important for a wide set of applications including air- craft/drone tracking, air trac data verication, or self-localization. Our focus in this work is on TO A estimation metho ds that can run on low-cost software-dened radio (SDR) receivers, as widely deployed in Mode S / ADS-B crowdsourced sensor networks such as the OpenSky Network. W e evaluate experimentally classical TO A estimation methods which are based on a cross-correlation with a reconstructed message template and nd that these methods are not optimal for such signals. W e propose two alternative methods that pro vide superior results for real-world Mode S / ADS-B signals captured with low-cost SDR receivers. The best method achieves a standard deviation error of 1.5 ns. 1 IN TRODUCTION Aircraft and unmanned aerial vehicles continuously transmit wire- less signals for air trac control and collision avoidance purposes. These signals are either sent as responses to interrogations by sec- ondary surveillance radars (SSR) or automatically on a p eriodic basis (ADS-B). Both types of signals are transmitted over the so- called Mode S data link [12] on the 1090 MHz radio frequency . Over the last few years, sensor network pr ojects have emerged which collect those signals using a crowd of low-cost software- dened radio (SDR) receivers such as e.g. the OpenSky Network [ 20 ], Flightaware [ 5 ], Flightradar24 [ 6 ] and many others. These sensor networks can leverage the time-of-arrival (TOA ) of Mo de S sig- nals for various kinds of applications, including aircraft localiza- tion [ 20 , 22 ], air trac data verication [ 13 , 16 , 17 , 19 , 21 ], and self-localization [ 15 ]. In those applications, a set of cooperating receivers measure locally the TOA of the arriving signals and then send these data to a central computation server . By joint process- ing the TO A of the same signal arriving at dierent receivers, the central server is able to estimate the location of the transmitter , the location of the receivers, or the exact time when the signal was transmitted. The accuracy of these applications heavily depends on the pr e ci- sion of the TO A estimation, and in order to estimate positions up to a few meters it is necessar y to estimate the TO A with nanose cond precision. The goal of this work is to provide a method for the TO A estimation of Mode S signals that delivers nanosecond-level precision even with low-cost SDR receivers, such as the widespread RTL-SDR dongle [ 3 ]. W e show that existing TOA estimation ap- proaches based on a cross-correlation with a r econstructed signal template are sub-optimal in the particular conte xt of Mode S signals. In fact, the lo ose tolerance margins allowed by the spe cications on the shape and position of each individual symb ol within the packet (up to ± 50 ns) adds uncertainty to the reconstruction of the whole packet waveform at the receiver . W e propose two alternative methods that improve the precision and at the same time reduce the computational load. W e test dierent variants of TOA estimation on real-world signal traces captured with RTL-SDR, which is currently the cheapest SDR device on the market and widely use d by crowdsourced sensor networks. Our results show that the best propose d method delivers TOA estimates with a standard deviation error of 1.5 ns. W e further identify the limited dynamic range of the RTL-SDR device (less than 50 dB with 8-bit analog-to-digital converter ( ADC) and xed automatic-gain controller ( AGC)) as the main performance bottleneck, and show that sub-nanosecond precision is achievable for signals that are not clipped due the limited dynamic range of the device. 2 BA CKGROUND This se ction provides background on aircraft signals which we rely on to estimate the T O A, and the limitations of classical TO A estimation methods. 2.1 Mode S signal format Hereafter we briey review the physical-layer format of SSR Mode S [ 18 ] reply and ADS-B messages transmitted by aircrafts on the 1090 MHz channel. Both packet formats consist of a preamble of 8 µ s plus a payload of 112 or 56 bits (only for other SSR Mode S replies) sent at 1 Mbps rate, for a total duration of 120 µ s or 64 µ s , respectively . The information bits are modulated with a simple Binary Pulse Position Mo dulation (BPPM) scheme as illustrated in Fig. 1: the symbol p eriod of 1 µ s is divided into two “ chips" of 0.5 µ s, and the high-to-low and low-to-high transitions encode bits “1" and “0" , respectively . It is clear from Fig. 1 that the BPPM modulation produces two types of pulses of dierent duration, denoted hereafter as “T ype-I" and “T ype-II". T ype-I pulses have a nominal duration of one chip period and are produced by the bit sequences “00" , “11" and “10" . The preamble consists of four T yp e-I pulses. On the other hand, T ype-II pulses have a nominal duration of two chip periods and are produced exclusively by the “01" sequence. 1 µsec “1” symbol 1 µsec “0” symbol preamble 8 µsec payload 11 2 µsec Ty p e -I peaks Ty p e -II peak … … time Figure 1: Mode S packet structure with a binar y PPM modu- lation. On average, we expe ct approximately 112 / 2 = 56 T ype-I and 112 / 4 = 28 T ype-II pulses for a payload of 112 bits. The real-valued baseband signal is then modulated on the 1090 MHz carrier frequency and transmitted over the air . On the r e ceiver side, the deco ding process relies exclusively on the signal amplitude , since in BPPM the signal phase carries no information. 2.2 Limitation of standard TOA methods The standard “course bo ok" approach to TO A estimation in the Additive White Gaussian Noise ( A W GN) channel is a correlation lter [ 14 ]: the received signal is cross-correlated with a known template corresponding to the source signal, and the point in time maximizing the cross-correlation module is taken as TOA estimate. The correlation method relies on the assumption that the source signal can b e reconstructed very precisely at the receiver , based on the signal specications and knowledge of the payload bits p m . Under this assumption, the correlation method represents the Max- imum Likelihood Estimator (MLE) [ 14 ]. Howev er , this assumption is problematic in the particular case of r eal-world Mode S signals. In fact, the standard specications tolerate up to ± 50 ns jitter in the position of each individual pulse within the packet: such high tol- erance value is practically negligible for the decoding process, but not for the task of determining the TOA with nanosecond precision. As to the shape of each pulse, tolerance values of 50 ns are allowed for the pulse duration and rise time and up to 150 ns for the decay time , while pulse amplitude may var y up to 2 dB (approximately 60%). Such lo ose tolerance margins introduce uncertainty in the prediction of the shape and position of the pulses in the source signal. Considering that Mode S signals are typically received with high SNR, such an uncertainty might well prevail over the eect of additive noise. Consequently , the correlation-based approach with a known packet template is no longer guaranteed to b e optimal, motivating the quest for alternative, more pr e cise methods. legacy' receiver ' (e.g.'dump1090)' decoded bits p m high- precision timestamp t m High7Precision' TOA' es?ma?on' modified receiver with HPTOA function IQ samples s m Figure 2: Block diagram of improved receiver with high- precision TOA estimation. 3 OUR TOA ESTIMA TION METHODS In this se ction we describe the general approach to TOA estimation based on the de coded payload and received signal samples, and then present the dierent TOA estimation algorithms that were tested. 3.1 Signal acquisition architecture In the software domain, the high-precision TO A estimation process can be seen as an additional function that is optionally called within the receiver and remains independent from the main decoding pro- cess. As such, it can be implemented on top of any legacy receiver , including but not limited to the widely adopted open-source tool dump1090 [ 1 ]. The overall block diagram of the proposed scheme is exemplied in Fig. 2. The legacy receiver takes as input a stream of complex in-phase and quadrature (IQ) samples collected at sam- pling rate f s (for RTL-SDR hardware f s = 2 . 4 MHz). The legacy receiver seeks to detect and decode the incoming packet and, if successful, provides as output the decoded bit sequence p m along with the indication of the leading IQ sample of the detected packet. Denote by s m the se quence of complex IQ samples corresponding to the whole packet. The sequence includes approximately 300 samples since we also pick a few samples immediately before and after the packet in or der to mitigate edge eects. The sample vector s m and the decoded bit vector p m represent the input to our TO A estimation block. 3.2 Proposed methods: CorrPulse and PeakPulse Hereafter w e describe two novel TO A estimation algorithms specif- ically developed for Mode S signals. For a generic packet m we shall denote by K m the total number of pulses in the whole packet (preamble and payload). The input v e ctor of complex samples s m is preliminarily upsampled by a factor N and transformed into vector s ′ m (for a review of upsampling pr ocess see e.g. [ 10 ]). T o illustrate, Fig. 3 plots an excerpt of the amplitude of both vectors, namely | s m | (top plot) and   s ′ m   (bottom plot), for a generic packet found in a real-world trace. The key aspect of the propose d algorithms is that the actual temporal position ˆ τ k of the generic k th pulse within the packet is estimated independently from other pulses, with no need to re- construct a template for the whole packet. For each pulse k ≥ 2 , 2 Figure 3: Example of received signal amplitude correspond- ing to the preamble and initial payload of a real ADS-B packet. Original samples at f s = 2 . 4 MHz (top, red circles) and corresponding upsample d version (b ottom, blue line). we compute the individual shift ∆ τ k def = ˆ τ k − τ k , i.e., the dierence between the estimated and nominal pulse position relative to the (estimated) position of the rst pulse. Finally , the pulse shifts are averaged in order to obtain the nal TO A estimate: ˆ t = ˆ τ 1 + 1 K m − 1 K m Õ k = 2 ∆ τ k (1) The two propose d variants dier in the way individual pulse position estimates are obtained, and which type(s) of pulses are con- sidered. In the rst variant, labeled CorrPulse , each pulse position is determined through pulse-level cross-correlation of the upsampled vector s ′ m with the corresponding nominal pulse shape . Both T ype-I and T ype-II pulses are considered in the nal averaging. In the second variant, labeled PeakPulse , individual pulse posi- tions are determined by simply picking the local maximum point value within the pulse interval, with no cross-corr elation operation. In this variant only T yp e-I pulses are considered, while Type-II pulses are ignor e d. This is motivated by the fact that T ype-II pulses have lower curvature, hence their local peaks cannot be identied as reliably as for T ype-I pulses. 4 EV ALU A TION METHODOLOGY This section describes how we e valuate our new methods. First, we introduce the other competing methods taken as reference for the comparison. Then, we present the testbed setup with commercial low-cost hardware. Finally , we provide details on the procedure adopted to empirically assess the precision of the TO A measurement methods in the given setup. 4.1 Other metho ds for comparison 4.1.1 Correlation with whole-packet template: CorrPacket . This is the canonical cross-correlation method with a known signal template. For every packet m , the whole packet template is reconstructed from the decoded bits p m and then cross-correlated with the am- plitude of the incoming signal. Her e also, upsampling by a factor N is adopted to achieve sub-sample precision. Within the template, the k th pulse is positione d at the nominal time τ k . As to the pulse shapes, we have tested two dierent variants: “Rectangular" (R), and “Smoothed" (S). The two versions will be denoted by CorrPacket/R and CorrPacket/S . The rectangular pulses have a nominal duration of 0.5 µ s and 1 µ s for T ype-I and T ype-II pulses, respectively , and zero rise/decay times. The rectangular pulse mask is represented ex- clusively by “0" and “1" values, hence multiplications with another vector reduce to element selection, which saves on computation load. The “Smo othed" shap e corresponds to the output of a low- pass lter with passband of 2.4 MHz—matched to the bandwidth of the RTL-SDR receiver—when the input signal is a nominal T ype- I/T ype-II pulse with the minimum decay/rise time of 50 ns as per specications [11]. 4.1.2 Existing dump1090 base d implementations. W e also evalu- ate the precision of the timestamp reported by the mutability fork of the open-source tool dump1090 [ 1 ]. Furthermore, we test on our traces also the method adopte d by Eichelberger et al. in a recent A CM SenSys’17 pap er [ 15 ] which is also based on dump1090. Code inspection rev ealed that this method is based on a cross-corr elation (implemented in frequency domain) with a partial packet template consisting of the preamble plus one quarter of the payload, with rectangular pulses and upsampling factor N = 25 . 4.2 T estbe d setup The experimental setup consists of two identical sensors connected to a single antenna through a power splitter and cables of identical length. The sensors are located on the roof of a building as Figure 4 shows. Every sensor consists of one RTL-SDRv3 “Silver" model [ 4 ] attached to a Raspb erry Pi-3 [ 2 ]. The A GC gain is set to a xed value, manually tuned to maximize the packet decoding rate. The sampling rate was set to f s = 2 . 4 MHz, the maximum value that our setup is able to acquire with sample losses. Ev er y I and Q sample is represented with 8 bit. The full stream of IQ samples are r ecorded one and processed multiple times oine. Our results are based on a sample trace of 5 minutes collected in Thun (Switzerland) on 02- Aug-2017 at time 09:41. The number of ADS-B packets that are correctly de coded at b oth sensors by the dump1090 open-source to ol [1] amount to 26445 from 59 dierent aircraft. 4.3 Evaluation Metrics In this se ction, we briey describe the methodology adopted to assess the precision of the dierent TO A estimation methods. The problem is not trivial, since our receivers are not synchronized and the “true " TO A is unknown. Therefore, we developed an evaluation method which allows us to quantify the TOA precision without a ground truth. Denote by t m , i the true absolute arrival time of packet m to receiver i and by ˆ t m , i the corresponding measured TO A 3 Figure 4: Experimental setup. T wo identical receivers con- nected to the same antenna via a splitter are collecting Mode S messages sent by aircraft. (by the method under test). In general, the measured value ˆ t m , i is aected by two distinct sources of error , namely clock error and measurement noise: ˆ t m , i = t m , i + ξ i ( t ) | t = t m , i + ϵ m , i . (2) The term ξ i ( t ) models the clock error between the receiver clock and the absolute time refer ence, and can be modeled by a slowly- varying function of time. Its magnitude depends on the hardware characteristics of the device, and specically on the stability of the local oscillator . The term ϵ m , i represents the measurement noise in the TOA estimation process and is modeled by a random variable with zero mean and variance σ 2 TOA . The precision of the TOA estimate, dened as the reciprocal of the noise variance, is independent of the clock error . The goal of the present study is to reduce σ 2 TOA . The prob- lem of counteracting the clock error component remains outside the scope of the present contribution. Here it suces to mention that the clock error can be mitigated by adopting receivers with GPS Disciplined Oscillators (GPSDO), or it can b e estimated and compensated in post-processing [7–9]. Hereafter w e illustrate the methodology to experimentally quan- tify the empirical TO A standard deviation ˆ σ TOA notwithstanding the presence of a non-zer o clock error component. First, we need to get rid of the unknown true absolute arrival time t m , i in Equa- tion (2). Since we use two identical receivers attached to the same antenna, we can set t m , 1 = t m , 2 = t m and subtract the TO A mea- surements at the two sensors to obtain the corresponding time dierence: ˆ ∆ t m def = ˆ t m , 2 − ˆ t m , 1 = ∆ ξ ( t m ) + ∆ ϵ m (3) wherein ∆ ξ ( t ) def = ξ 2 ( t ) − ξ 1 ( t ) denotes the compound clock error , and ∆ ϵ m def = ϵ m , 2 − ϵ m , 1 the compound measurement error with variance σ 2 ∆ ϵ = 2 σ 2 TOA . At short time-scales, within the coherence time of the process ∆ ξ ( t ) , the clock error r epresents a systematic (a) Low upsampling factor (b) High upsampling factor Figure 5: ECDF of ∆ ϵ residuals. error , i.e. a bias term that can be estimated and removed in order to estimate the error variance σ 2 ∆ ϵ . W e do so by modeling the slowly- varying function ∆ ξ i ( t ) by a p olynomial whose coecients are estimated by standard order-recursiv e Least Squares (refer to [ 14 , Chapter 8] for details). After removing the estimated clock error component, we obtain a set of residuals { ∆ ϵ } . Their Mean Square Error (MSE) represents an empirical estimate of twice the TOA variance M S E ∆ ϵ = 2 · σ 2 TOA . Accordingly , their Root Mean Square Error (RMSE) pro vides a direct empirical estimate of the TO A error standard deviation, formally: ˆ σ TOA = 1 √ 2 R M S E ∆ ϵ ≈ 0 . 7 · R M S E ∆ ϵ . 5 N UMERICAL RESULTS W e now present the r esults on the precision of the dier ent TOA estimation methods as evaluated in our testbed. 4 Figure 6: TOA standard dev . error vs. upsampling factor N 5.1 Error distribution In Fig. 5 we plot the Empirical Cumulative Distribution Function (ECDF) of the residuals ∆ ϵ ’s obtained with dierent TOA estimation methods for all the packets in the test trace. The corresponding values of the TOA error standar d deviation ˆ σ TOA are reported in the leftmost column of T able 1. For those applications where the computation load is of con- cern, it is relevant to investigate the performance of the dier ent methods with moderate value of the upsampling factor ( N = 25 ). For CorrPacket and CorrPulse , we consider the rectangular pulse shape with binar y 0/1 values, due to lower computation load. Refer- ring to Fig. 5(a), we observe that the proposed PeakPulse algorithm achieves a R M S E ∆ ϵ = 3 . 15 ns, less than half the value of the canoni- cal CorrPacket/R method. It is remarkable that such go od r esult was obtained with no cross-correlation operation. Fig. 6 shows ˆ σ TOA for dierent values of the upsampling factor N . W e observe that the precision of the proposed methods PeakPulse and CorrPulse/R improves faster than CorrPacket/R with increasing N . These results indicate that PeakPulse should be preferred when computation load is at premium. Next we consider applications that enjoy abundant computation power , for which the main goal is to maximize pr e cision and compu- tation load is not of concern. For these, it is convenient to consider higher upsampling factors ( N = 83 in our case) and, for the cross- correlation methods, the more elaborated “Smoothe d" pulse shape. The latter matches more closely the pulse shape passed through the RTL-SDR front-end, leading to slightly higher precision than the simpler “Rectangular" shape, as can b e veried from T able 1. The ECDF of the residuals ∆ ϵ ’s for these methods are plotted in Fig. 5(b). It can be se en that the proposed CorrPulse/S method is more precise than the classical CorrPacket/S method, and achieves R M S E ∆ ϵ = 2 . 16 ns corresponding to ˆ σ TOA = 1 . 51 ns. 5.2 Error vs. signal strength In the following, we investigate the impact of signal strength on the TO A error obtained with the most precise method, namely Cor- rPulse/S with N = 83 . For a generic packet m and sensor i , we denote by γ m , i the av erage of the squared pulse height ov er all pulses — an indicator of the arriving packet strength. Furthermore , we denote by β m , i the number of pulses that are clipp ed in the receiver due to one or more of the corresponding IQ samples saturating the ADC. estimation metho d ˆ σ TOA [nanoseconds] all packets L M H legacy dump1090 45 . 20 44.94 45.19 45.43 SenSys’17, N = 25 5 . 90 6.11 5.88 5.78 CorrPacket/R , N = 25 4 . 98 5.48 4.85 4.94 CorrPacket/R , N = 83 2 . 14 3.04 1.78 2.35 CorrPacket/S , N = 83 2 . 07 3.00 1.68 2.275 CorrPulse/R , N = 25 1 . 89 2.75 1.56 1.86 CorrPulse/R , N = 83 1 . 63 2.72 1.04 1.77 CorrPulse/S , N = 83 1 . 51 2.60 0.79 1.77 PeakPulse , N = 25 2 . 20 3.36 1.70 2.23 PeakPulse , N = 83 2 . 12 3.44 1.62 2.17 T able 1: Empirical estimates of TO A error standard devia- tion ˆ σ TOA . Figure 7: Absolute error | ∆ ϵ m | vs. packet strength γ m . In Fig. 7, we plot for each individual packet m the absolute value of the r esidual error | ∆ ϵ m | obtained with CorrPulse/S ( N = 83 ) against the mean signal strength between the two sensors γ m def = γ m , 1 + γ m , 2 2 . Each packet is classied into one of three classes: packets with γ m ≤ 0 . 04 are labele d by “L" , packets with min i = 1 , 2 β m , i ≥ 10 are lab eled with “H" , and all remaining packets are labeled with “M" . The three classes are marked respectively with black, red and blue markers in Fig. 7. The estimated TO A error standard de viation obtained by each metho d for each class are reported in T able 1. On one extreme, timing estimates for “L" packets with lower str ength are impaired by quantization noise. On the other extr eme, packets received with high strength are subject to ADC clipping, a form of distortion that clearly degrades timing precision. As expected, these two classes yield higher error with all methods. Between the two extremes, the strength of “M" packets ts well the dynamic range: for these, the proposed method achieves ˆ σ TOA = 0 . 79 ns. In our traces, less than 60% of all packets fall into class “M". With better hardware, and specically with more ADC bits and larger dynamic range, it would be possible to tune the AGC gain so as to increase the fraction of packets falling in this class, thus improving the overall precision. The above results indicate that the received packet metrics γ m , 1 and β m , i can be used to provide, for each individual TOA mea- surement ˆ t m , i , also an indication of the expected precision, i.e., of 5 -4 -2 0 2 4 Stand ar d No r mal Qu antiles -20 -10 0 10 20 Quantiles of I nput Sample All pack ets -4 -2 0 2 4 Stand ar d No r mal Qu antiles -20 -10 0 10 20 Quantiles of I nput Sample L -cl ass -4 -2 0 2 4 Stand ar d No r mal Qu antiles -20 -10 0 10 20 Quantiles of I nput Sample M-cl ass -4 -2 0 2 4 Stand ar d No r mal Qu antiles -20 -10 0 10 20 Quantiles of I nput Sample H-c lass Figure 8: Quantile-quantile plot of empirical errors ∆ ϵ vs. normal distribution. the error variance ˆ σ 2 m , i aecting each individual measurement. In this way , algorithms that take TOA measurements as input (e.g., for position estimation) have the p ossibility to weight optimally each individual input measurement, as done e.g. in W eighted Least Squares methods [23]. Finally , we nd that within each class the empirical error distribu- tion is very well approximated by the Gaussian distribution, as seen from the normal Q-Q plots in Fig. 8. This justies the adoption of Least Squares (LS) methods for position estimation problems based on input TO A measurements [ 7 ], since for normally distributed in- put errors the LS solution coincides with the Maximum Likelihood estimate. 6 CONCLUSIONS AND OU TLOOK W e have presented two variants of a novel TO A estimation method for Mode S signals that does not rely on long cross-correlations with the template of a full packet. The most precise variant, namely CorrPulse/S , involves only short cross-correlation operations on individual pulses. The other variant, namely PeakPulse , is lighter to compute, involv es no cross-correlation operation and w orks well also with moderate upsampling factors. W e have shown that such algorithms can achieve TO A estimates with nanosecond-level precision even with r eal-world signals cap- tured with the cheapest SDR hardware that is curr ently available, namely RTL-SDR. A closer look at the test results rev eals that the main limiting factor for the achievable TO A precision with RTL- SDR is the limited dynamic range — less than 50 dB with 8-bit ADC and xe d A GC — resulting in a large fraction of packets being clipped or drowned into quantization noise. For packets that are received with signal strength well within the dynamic range of the receiver , the CorrPulse/S achieves sub-nanosecond precision. It can be expected that precision can b e further improved with better hard- ware. The PeakPulse method has be en implemented in C, integrated in the dump1090 receiver and is released as open-source 1 . 1 http://github.com/openskynetwork/dump1090-hptoa A CKNO WLEDGMEN TS This work has been funded in part by the Madrid Regional Gov- ernment through the TIGRE5-CM program (S2013/ICE-2919). W e would like to thank Manuel Eichelberger from ETH Zurich for shar- ing the code we used in our evaluation for comparison purposes. REFERENCES [1] 2016. dump1090, https://github.com/mutability/dump1090. https://github .com/ mutability/dump1090. [2] 2016. Raspberry Pi 3 Model B, https://ww w .raspberr ypi.org/products/raspberry- pi-3-model-b/. [3] 2016. Silver v3 specications . http://w ww .rtl-sdr .com/buy-rtl-sdr-dvb-t-dongles/. [4] 2017. RTL-SDR dongle v3, http://ww w .rtl-sdr .com/buy-rtl-sdr-dvb-t-dongles/. [5] 2018. FlightA ware . http://ww w .ightaware.com/. [6] 2018. FlightRadar24 . https://www.ightradar24.com/. [7] F. Ricciato, S. Sciancalepore , F. Gringoli, N. Facchi and G. Boggia. 2018. Position and V elocity Estimation of a Non-cooperative Source From Asynchronous Packet Arrival Time Measurements. IEEE Trans. on Mobile Computing (2018). [8] G. de Miguel, J. Portas, J. Herrer o. 2005. Correction of propagation errors in Wide Area Multilateration systems. In Proc. of 6th European Radar Conference . [9] G. Galati, M. Leonardi, I.A. Mantilla-Gaviria and M. T osti. 2012. Lower bounds of accuracy for enhanced mode-S distributed sensor networks. IET Radar , Sonar and Navigation (2012). [10] Fred J. Harris. 2004. Multirate Signal Processing for Communication Systems . Prentice Hall. [11] ICAO . 2007. A eronautical Communications, V olume IV , Sur veillance and Collision A voidance Systems. (July 2007). [12] ICAO 2008. T echnical Provisions for Mode S Services and Extended Squitter . ICAO . Doc 9871, First Edition. [13] K. Jansen, M. Schäfer , D. Moser, V . Lenders, C. Pöpp er , , and J. Schmitt. 2018. Crowd-GPS-Sec: Leveraging Crowdsourcing to Detect and Localize GPS Spoong Attacks. In IEEE Symposium on Security and Privacy (S&P) . [14] S. M. Kay . 1993. Fundamentals of Statistical Signal Processing, vol. 1, Estimation Theory . Prentice-Hall. [15] S. Tanner M. Eichelberger , K. Luchsinger and R. Wattenhofer . 2017. Indoor Localization with Aircraft Signals. In ACM SenSys . [16] V . Lenders M. Strohmeier and I. Martinovic. 2015. Lightweight Location V erica- tion in Air Trac Surveillance Networks. In ACM Cyber-Physical System Security W orkshop (CPSS) . [17] D. Moser , P. Leu, V . Lenders, A. Ranganathan, F. Ricciato, and S. Capkun. 2016. Investigation of Multi-device Location Spoong Attacks on Air T rac Control and Possible Countermeasures. In ACM Conference on Mobile Computing and Networking (Mobicom) . [18] RTCA 2011. Minimum Operational Performance Standards for Air Trac Control Radar Beacon System / Mode Select (A TCRBS / Mode S) Airborne Equipment . RTCA. DO-181E. [19] M. Schäfer , M. Strohmeier , V . Lenders, and J. Schmitt. 2015. Secure Track V eri- cation. In Security and Privacy (S&P), 2015 IEEE Symposium on . [20] M. Schäfer , M. Strohmeier , I. Martinovic V . Lenders, and M. Wilhelm. 2014. Bring- ing Up OpenSky: A Large-scale ADS-B Sensor Network for Research. In 13th ACM/IEEE Int. Conf. on Information Processing in Sensor Networks (IPSN’14) . [21] M. Strohmeier , V . Lenders, and I. Martinovic. 2015. On the Security of the Automatic Dependent Surveillance-Broadcast Protocol. IEEE Communications Surveys and Tutorials Journal (2015). [22] M. Strohmeier , V . Lenders, and I. Martinovic. 2018. A k-NN-based Localiza- tion Approach for Crowdsourced Air Trac Communication Networks. IEEE Transactions on Aerospace and Electronic Systems (TAES) (2018). [23] Tilo Strutz. 2010. Data tting and uncertainty: A practical introduction to weighted least squares and beyond . Vieweg and T eubner. 6