Opportunistic-Target-Measurement-Based Narrowband Statistical Modeling of Civil Aviation Surveillance Signal at 1090 MHz

Automatic dependent surveillance—broadcast (ADS-B) is one of the next-generation aeronautical surveillance systems for air traffic control. ADS-B requires an aircraft to periodically broadcast its own position to other aircraft and ground stations, thereby enabling high-performance surveillance. In this study, the received signal strength (RSS) of the ADS-B signal was measured and characterized for opportunistic flights. The RSS distribution for a single aircraft was modeled as a sum of three components: nominal RSS, bias in the equivalent isotropically radiated power (EIRP), and fading. Then, the RSS distribution for multiple aircraft was defined, considering that the EIRP bias and fading parameter are different for different aircraft. To this end, a composite distribution was employed, the parameters of which were estimated from the bias and fading statistics. Furthermore, a practical approximation of the model was proposed. The proposed model is suitable for large-scale data, enabling the realization of aircraft-by-aircraft analysis. Also, it can provide a clear explanation of the mechanism by which the RSS distribution is formed.


I. INTRODUCTION
A ERONAUTICAL surveillance systems are used to provide aircraft information to air traffic controllers, and thus, they constitute an essential infrastructure for ensuring safe flights. While conventional radars have been a reliable means of surveillance, the need for advanced systems is increasing to meet the requirements of growing air traffic. One such system is the automatic dependent surveillancebroadcast (ADS-B), which requires aircraft to periodically broadcast their information to ground stations and other aircraft. The signal used for ADS-B at 1090 MHz is called the extended squitter [1]- [3]. The signal conveys 112 bits, which are pulse position modulated. Different types of extended squitters can be defined depending on the information con- tained in the signal, such as position, velocity, and identification squitters. The operation of ADS-B has been either planned or initiated worldwide.
A key enabler for reliable ADS-B operation is appropriate coverage design based on deep understanding of air-ground propagation channels [4]. In addition, security enhancement in ADS-B has been extensively investigated recently [5]- [16], to which contribution from the field of air-ground propagation is expected. However, sounding of air-ground channels requires flight experiments [17]- [24], which requires a high cost. Performing channel sounding via unmanned aerial vehicles (UAVs) can reduce the measurement cost, as reported in [25] and [26]; however, UAVs are limited in terms of the flight altitude. Therefore, Naganawa et al. [27], [28], Naganawa and Miyazaki [29]- [31] have proposed measurement and analysis of opportunistic target as a complementing approach. Although ADS-B equipage of aircraft has not been mandated in many countries yet, presently, many aircraft have started broadcasting ADS-B signals in preparation for future operation. The measurement of the ADS-B signals results in a large amount of aeronautical data with relatively low cost. To manage these data, architectures for ADS-B data collection and analysis were proposed in [32]- [34]. Applications of ADS-B data analysis for meteorology [35], [36], flight performance [37], pilot behavior [38], and link-level performance [33], [39], [40] were proposed. Research attention on ADS-B data analysis is, thus, increasing.
Radio propagation is an attractive application of the ADS-B measurement. Compared with flight experiments, the ADS-B approach has opposite characteristics; it is less accurate; however, it can be performed under several measurements conditions, at lower operational costs [27]. Therefore, the ADS-B approach is suitable to complement flight experiments. However, detail analysis of ADS-B data from the aspect of radio propagation has not been reported yet. For example, measurements from a 14-day period were analyzed in [33] and [34], but the resulting log-distance model was relatively simple because the uncertainty of the transmission antenna type (bottom or top), transmit antenna gain, and transmit power hindered the realization of a detailed analysis. Another example is [28], in which visual comparison with the free space model was presented for one aircraft. In this study, a statistical model of the received signal strength (RSS) of the ADS-B signal was achieved based on opportunistic measurement. To the best of the authors' knowledge, this is the first study that conducted a detailed analysis of ADS-B data from the aspect of radio propagation. It is noted that this article is an extended version of the authors' conference article [30] and technical report without peer review [31].
First, the probability distribution of the RSS was modeled for a single aircraft as a sum of three components: the nominal RSS, bias in equivalent isotropically radiated power (EIRP), and fading. Then, the RSS distribution for multiple aircraft was presented. Because the EIRP bias and the fading parameter differ for each aircraft, the distributions of the EIRP bias and the fading parameters were obtained from the measurement result. This statistical treatment of the aircraft-dependent parameters results in a compound distribution as the model for multiple aircraft. Furthermore, an approximation of the compound distribution by a normal distribution was proposed. The proposed model was verified against the measurement data. The proposed model is expected to be suitable for large-scale data, enabling the realization of aircraft-by-aircraft analysis. Also, it can provide a clear explanation of how the RSS distributions of individual aircraft relate to the aggregated distribution.
The remaining article is organized as follows. Section II describes the measurement setup and data processing technique. Section III defines the statistical model for a single aircraft. In Section IV, the proposed model is further extended for multiple aircraft. In Section VI, the results obtained by using additional measurement sets are presented. The conclusions are presented in Section VII. Fig. 1 shows an image of the antenna used and a map that indicates several measured tracks of aircraft and the receiver position. The map is intended to show an overview of the measured traffic flows, where 208 aircraft selected from the 1657 aircraft measured during 24 h are plotted. The measured aircraft corresponded mainly to arrivals and departures from the Tokyo International Airport and Narita International Airport.

A. Measurement Setup
The antenna is directional in the horizontal plane for reducing the effect of signal collisions and increasing the gain. The horizontal beamwidth is 28.7 • . The vertical pattern is cosecant-squared with a 7.6 • beamwidth. The peak gain is 19.0 dBi. The measured pattern of the antenna was available for analysis and modeling. Fig. 2 shows the view from the antenna. The antenna is located on the rooftop of the antenna test tower of the Electronic Navigation Research Institute with a good line-of-sight (LOS). A wire and a VHF antenna, which are potential scatterers, are present in the vicinity of the antenna.
The antenna was connected to a receiver via a 30-m cable with a loss of 3.1 dB. The receiver consists of a softwaredefined radio (SDR) transceiver (USRP 2901), a host computer, and a rubidium oscillator (SRS FS 725). The SDR transceiver was calibrated against a spectrum analyzer before the measurement. The SDR transceiver continuously captured baseband signals with a center frequency, sampling frequency, and bandwidth of 1090, 10, and 8 MHz, respectively. The host computer performed preliminary signal detection with the baseband signals and recorded potential ADS-B signals. The recorded signals were demodulated offline by using the enhanced reception technique defined in DO-260B [42], with a detection threshold of −80 dBm. The demodulated signals were decoded according to [1] and [3]. For each valid ADS-B signal, the RSS was measured by the technique defined in [29], which mitigates the effect of signal collisions.

B. Data Preprocessing
The measured signals were preprocessed to generate the subsets to be analyzed by employing the following steps.
1) Track construction: The aircraft addresses and aircraft positions (longitudes, latitudes, and barometric altitudes) were obtained from the airborne position squitters. The barometric altitudes were converted to the geometric altitudes by using the information available in the airborne velocity squitters. In this manner, the aircraft tracks were constructed. 2) Aircraft and signal selection: Arrivals to the Tokyo International Airport were selected for analysis. Further, the signals when the aircraft stayed within the main lobe of the antenna were obtained by clipping the track.

C. Measurement Set
Table I describes measurement sets. Seven measurement sets, each for a duration of 1 h, and a 24 h measurement set were obtained. The measurement sets 1-2 and 3-8 were measured under the south-wind configuration and north-wind configuration of the Tokyo International Airport, respectively. Because the proposed model did not work appropriately when mixing the two configurations, the results obtained using the measurement sets 3-8 are presented in Sections III-V. The results obtained using the measurement sets 1-2 are summarized in Section VI.
In Table I, the measured signals/aircraft are labeled as "valid." As a result of the preprocessing in Section II-B, the subsets of signals/aircraft for analysis were obtained, which are labeled as "analyzed" in Table I. The number of these signals and aircraft are 112 354 and 497, respectively, under the north-wind configuration. However, it was found that two aircraft were possibly uncompliant in Section IV-B. Therefore, having excluding these aircraft and signals, the remaining 495 aircraft and 111 896 signals were used in the modeling in Sections III-V except for Section IV-B. They are labeled as "modeled" in Table I. A set of RSS values of the signals selected for the analysis and modeling is called a sample in the rest of this article. The sample size is 226 per aircraft on average.

A. Model
The modeling of the distribution of the RSS pertaining to a single aircraft is described in this section. The RSS is denoted by P r . The proposed model is a sum of three components: nominal RSSP r , EIRP bias B, and fading X, as follows: (1) Fig. 4 illustrates the relationship among the RSS and the three components. The nominal RSSP r is deterministic and calculated according to [2], which employs a free-space path loss model with a nominal transmission power and a nominal transmit antenna gain. The EIRP bias B is a constant that depends on the aircraft; it mainly models the differences between the nominal and actual values of the transmission power and transmission antenna gain. The fading component X is a random variable, which models the instantaneous RSS variation caused by multipath fading and aircraft maneuvering. Candidates for the statistical distribution of X were the normal distribution and the Nakagami-Rice distribution. The normal distribution was selected owing to its mathematical simplicity. The Nakagami-Rice distribution was selected owing to the agreement between its physical interpretation and the measurement condition. The best fit distribution function was selected from the measurement result, as explained in Section III-D.
Because the center of the RSS variation is determined by the EIRP bias, only parameters that characterize spread of the RSS variation are employed, specifically, the standard deviation and K -factor K for the normal and Nakagami-Rice distribution, respectively. In this study, and K were termed fading parameters and considered to depend on the aircraft. For the convenience of analysis, the stochastic part of the model was defined as = B + X. Here, it is noted that statistical modeling the rest of this article mainly concerns .
does not include the effect of distance. Using this, (1) was written as follows: The three components are described as follows.

B. Nominal RSS
According to [2], the nominal RSS can be calculated usinĝ where P t is the transmission power, G t is the transmission antenna gain, L atm is the atmospheric attenuation loss, L fs is the free-space propagation loss, G r is the receiver antenna gain, and L c,r is the cable loss at the receiver. G t =0 dB was selected, assuming that the transmission antenna is isotropic. P t = 51 dBm, which is the minimum power for A1, A2, and B1 classes of aircraft [1], was selected, as it was presumed to be suitable for a majority of the measured aircraft. L atm was 0.0090 dB per 1 NM (1.852 km). L fs = 20 log 10 (4πd/λ), where d is the distance between the receiver and aircraft, and λ is the wavelength. For the other receiver parameters, specifically, G r and L c,r , the measured values were employed.

C. EIRP Bias
The EIRP bias was estimated from the measurements. The measured values of the RSS and the stochastic part were denoted byP r and˜ , respectively. From (2),˜ was calculated by subtracting the nominal RSS: This eliminates the effect of distance variation on the RSS. Then, the estimate of the EIRP bias was calculated as a median value of the differenceB The bias estimate in the presented case was 6.2 dB.

D. Fading
Measurements of the fading component were determined as the part remaining when EIRP bias was subtracted from the stochastic partX Then, the fading parameters and K were estimated by the maximum likelihood estimator provided by MATLAB. In the presented example, the estimates were˜ = 1.2 dB andK = 14.1 dB. Fig. 5 shows the empirical and fit cumulative distribution functions (CDFs), and a good agreement can be observed between them. To evaluate the goodness of fit, the Kolmogorov-Smirnov (KS) test [48] with a significance level of 0.05 was applied, and neither of the normal and Nakagami-Rice distributions were rejected. Furthermore, the Akaike Information Criteria (AIC) was applied to identify the best fit distribution [43], [44], and the Normal distribution was selected.
The above procedure was repeated for the 495 aircraft selected for modeling under the north-wind configuration. The normal distribution was noted to be the best fit distribution for 288 aircraft, whereas the Nakagami-Rice distribution was selected for 207 aircraft. The KS test results indicated that the normal distribution was not rejected for 438 aircraft, whereas the Nakagami-Rice distribution was not rejected for 406 aircraft. These results indicated that the normal distribution was a better fit compared to the Nakagami-Rice distribution. To make sure the applicability of the normal distribution, the Lilliefors test, which improves the KS test for a normal distribution with estimated parameters [48], was additionally applied. The normal distribution was not rejected for 278 aircraft. Therefore, the normal distribution was determined as a model for the fading component IV. CHARACTERIZATION FOR MULTIPLE AIRCRAFT As defined in Section III, the RSS distribution for a single aircraft was modeled using the EIRP bias B and the fading parameter . However, the actual values of these parameters are different for each aircraft. Therefore, the aggregated RSS distribution formed by multiple aircraft was modeled, as defined in this section. It should be noted that only the stochastic part was considered; the nominal RSSP r is deterministic and not within the scope of modeling. Also, it is assumed that the number of aircraft that constitute the distribution is sufficiently large.

A. Model
The compound probability distribution was employed. When a random variable Y follows a distribution, the parameter(s) of which are also random variable(s), the distribution of Y is known as the compound probability distribution. In the present case, the RSS distribution formed by multiple aircraft was a compound distribution, and the bias and fading parameter were randomly distributed parameters. It should be noted that the compound distribution has been used for combining different propagation phenomena, for example, the Suzuki distribution [45]; however, the modeling of opportunistic aircraft is a novel application.
The probability density functions (PDFs) for the bias and fading parameter were denoted using p(B) and p( ), respectively. The joint PDF for the parameters were denoted by p(B, ). The RSS distribution for multiple aircraft can then be defined as a compound distribution obtained by marginalizing the parameters as follows [46]: where p( ) and P( ) are the PDF and CDF of , respectively. The ranges of integration are (−∞, ∞) and (0, ∞) for B and , respectively. The two parameters were assumed to be statistically independent, which writes (8) as follows: The closed-form solution was difficult to be obtained. Instead, numerical integration can be computed, for which quad2d function in MATLAB was used in this study. In the following, estimation of the distribution of B and is described.

B. Distribution of EIRP Bias
The EIRP bias was estimated for the 497 aircraft that were initially selected for analysis under the north-wind configuration. Among them, 495 aircraft had positive EIRP biases, whereas two aircraft had negative EIRP biases. A positive EIRP bias suggests that the aircraft is compliant with the minimum transmission power requirement. In other words, the aircraft with the negative EIRP bias possibly did not comply with the requirement, and they were removed from the analysis as outliers.
The bias distribution formed by the remaining 495 aircraft was fit using the normal, exponential, gamma, Weibull, and log-normal distributions. Among these, the normal distribution was selected as it had the lowest AIC value. Fig. 6 shows the fitting result, in which the CDF is plotted. Agreement between the normal distribution and the measurement was observed by visual inspection. In the Lilliefors test, the normal distribution was not rejected. Thus, the bias was modeled using with the estimated parameters beingμ B = 5.0 andσ B = 1.1.

C. Distribution of Fading Parameter
The same procedure as that for the bias was repeated to determine the fading parameter . Among the five distributions, the log-normal distribution was selected as it had the lowest AIC value. Fig. 7 shows the CDF, in which agreement between the log-normal distribution and the sample was observed by visual inspection. In the Lilliefors test, the log-normal distribution was not rejected. Thus, the following model was obtained: ∼ Lognormal μ ln , σ 2 ln (11) where μ ln and σ ln are the mean and standard deviation of the logarithmic values. The estimates wereμ ln = 0.30 and σ ln = 0.20. These values can be converted into the mean and standard deviation of non-logarithmic values, as follows [47]:

D. Correlation Between Fading and EIRP Bias
Correlation between the EIRP bias B and the fading parameter was evaluated. Fig. 8 shows a scatter plot for B and , in which each point represents an analyzed aircraft. By visual inspection, any correlation was not observed. In addition, the correlation coefficient was −0.04. Therefore, it was concluded that B and are not correlated. This result verified the assumption of (9).

E. RSS Distribution and Verification
Substituting the distributions estimated as described in Section IV-B and IV-C into (9), the model of the RSS distribution for multiple aircraft was obtained. To verify the proposed model, it was compared with the empirical distribution obtained by aggregating sample over all the aircraft. The size of the aggregated sample was 111 896. Fig. 9(a) shows a comparison between the proposed model and aggregated sample, in which a good agreement was observed. The maximum error between the CDFs was 0.010, which is sufficiently small considering that the main applications of the model are link budget calculation and coverage design. These results indicate that (8) and (9) are valid.
Furthermore, in order to analyze the source of the error, the proposed model was computed with measurements of B and by substituting the following equations into (9): whereB i and˜ i denotes the measurements of B and , respectively, for i th aircraft, N denotes the number of analyzed aircraft, and δ is the Dirac delta function. The result is plotted in Fig. 9(b), in which vicinity of the maximum error is enlarged. The error of the proposed model was reduced by approximately half when the measured B and were used. The removed error corresponds to modeling error of the B and distributions, while the remaining error corresponds to the RSS distribution model of each aircraft.

V. APPROXIMATION
The proposed model in Section IV requires numerical integration with respect to B and . Alternatively, an approximation of the model was developed, as described in this section.

A. Proposed Approximation
The compound distribution of given by (9) was approximated using a normal distribution where μ and σ are the mean and standard deviation, respectively. According to [46], the mean and variance of a composite distribution can be calculated from parameter statistics. To do so, μ B , μ , σ B , and σ , which were all defined in Section IV, are used. In addition, let μ |B, and σ |B, denote the conditional mean and standard deviation for a given B and , respectively. Then, the mean of the compound distribution can be calculated as follows: where p( )d = 1. Next, the variance of the compound distribution can be calculated as follows: To calculate the first term in (18), the following relationship can be used: where E is an expectation. The second term in (18) is the definition of the standard deviation of B. Thus, (18) can be reduced to the following expression:

B. Verification
Substituting the parameters estimated as described in Section IV-B and IV-C into (17) and (20), the approximated model was obtained. Fig. 10(a) shows a comparison between the proposed approximation and measurement sample, in which a good agreement was observed. The maximum error between the CDFs was 0.013, which is acceptable considering the main applications of the model. However, it should be noted that the Lilliefors test rejected the normal distribution for the aggregated sample.
To analyze the source of the error, the proposed approximation was compared with direct fitting, i.e., the normal distribution with statistics estimated from the aggregated sample, and the proposed model without approximation. Fig. 10(b) shows a comparison among them, in which vicinity of the maximum error is enlarged. Comparison between the proposed models with and without the approximation revealed error caused by the approximation, which was 0.0042. This amount of error is acceptable considering that the convenience of using the normal distribution outweighs the accuracy obtained using numerical integration, which is computationally intensive. Furthermore, comparison between the approximation and the direct fitting revealed error caused by the parameter estimation, which was negligible. Agreement in the parameter estimates was excellent indeed; the means μ for the proposed model and direct fitting were 5.038 and 5.044, respectively. The standard deviations σ for the proposed model and the direct fitting were 1.760 and 1.755, respectively. From the above results, it was concluded that (16)- (20) are valid.

C. Discussion
The advantages of the proposed model are summarized here. First, large-scale analysis can be easily realized by decomposing a sample. In the presented case, the aggregated sample, the size of which was 111 896, was divided into the 495 aircraft. The size of the divided samples was 226 per aircraft on average. This aspect is valuable, especially for opportunistic measurement, which can result in an enormous amount of data.
Second, the compound distribution can provide a clear explanation of how the RSS distributions of individual aircraft relate to the aggregated distribution. Equation (17) indicates that the center of the RSS distribution depends only on the EIRP bias and is not related to fading. Equation (20) indicates that the spread of the RSS distribution depends both on the EIRP bias and fading. With σ = 0.28, μ = 1.38, and σ B = 1.1, all the terms in (20) have a significant contribution.

VI. RESULTS FOR SOUTH-WIND CONFIGURATION
Sections III-V showed the results obtained using the measurement sets 3-8, which were measured under the northwind configuration of the Tokyo International Airport. This section shows the results obtained using the measurement sets 1-2, which were measured under the south-wind configuration. The main purposes of showing these results are to examine the potential applicability of the proposed model and potential changes of the parameter estimates in different measurement scenarios. It should be noted that these results are not intended as a complete characterization of the south-wind configuration because the number of analyzed aircraft was 51, which is far less than that of the north-wind configuration. Fig. 11(a) and (b) shows representative tracks in the north-and south-wind configurations, respectively. The flying directions in the north-and south-wind configurations were northeastward and eastward, respectively, as a whole. The altitude in the south-wind configuration tended to be higher than those in the north-wind configuration.

A. Applicability of the Model
First, the same analysis as described in Sections III and IV were applied to the measurement sets 1-8, i.e., the northand south-configurations were mixed. However, the proposed model was not applicable in this case. For example, Fig. 12(a) shows the fitting of , in which a disagreement between the measurement and log-normal distribution is observed. The Lilliefors test rejected the log-normal distribution.
Therefore, the measurement sets 1-2 were used independently. In this case, the proposed model was applicable. For example, Fig. 12(b) shows the fitting of , the agreement between the measurement and log-normal distribution is acceptable, contrary to Fig. 12(a). The Lilliefors test did not reject the log-normal distribution. Other remarkable results are as follows.
1) The normal distribution was determined as a model for the fading component. 2) Correlation between B and was not observed with the correlation coefficient of 0.027. 3) Fig. 13(a)-(c) shows the results that correspond to Figs. 6, 9, and 10, respectively. Because the proposed model and the approximation show a good agreement with the measurement, it was verified against the southwind data too. 4) The normal distribution was the second-best fit as an EIRP bias distribution in terms of AIC, while the lognormal distribution was the best fit, which was a different result from the north-wind configuration. However, a good agreement of the normal distribution can be observed by visual inspection of Fig. 13(a). In addition, the Lilliefors test did not reject the normal distribution. Therefore, the normal distribution is still acceptable as an EIRP bias model. These results indicate that the proposed model can be applied to different measurement scenarios. Table II summarizes the estimated parameters for the north-and south-wind configurations. Based on these results, the effect of the measurement scenario on the parameters can be discussed as follows.

B. Effect on Parameter Estimate
1)μ B andσ B did not change. This result was reasonable because the EIRP bias is a constant associated with each aircraft and should not change depending on the flight route. This result also indicated that the EIRP bias was successfully estimated as intended. 2)μ ln increased, which means fluctuation due to the fading was severer. One possible explanation for this is reflection and scattering by primary wings of the aircraft. The receiver was located almost perpendicular to the flying direction in the south-wind configuration, which possibly resulted in more significant effect of the primary wings. Another possible explanation is that the propagation channel changed according to the flight route. In either way,μ ln can change depending on the measurement scenario, and extensive characterization is suggested as future work. On the other hand,σ ln did not change, which is an interesting result. Becauseσ ln is a parameter for the fading, it can change according to the flight route. In order to judge whether this result is a special case or not, characterization in various scenarios is suggested as future work. 3) A good agreement between the proposed model and direct fitting was observed in terms ofμ andσ . Therefore, it is expected that the proposed model can substitute the direct fitting in other measurement scenarios too. 4)μ did not change, which can be explained by (17) withμ B unchanged. On the other hand,σ increased, which can be explained by (20) withμ ln increased. As demonstrated in these results, the proposed model can explain the mechanism by which the RSS distribution is formed.

VII. CONCLUSION
The RSS distribution for a single aircraft was modeled using the sum of the nominal RSS, EIRP bias, and fading. The RSS distribution for multiple aircraft was modeled by a compound distribution because the EIRP bias and the fading parameter are aircraft dependent and were modeled as random variables. Furthermore, the compound distribution was approximated by a normal distribution. The proposed model and its approximation demonstrated a good agreement with measurement in terms of the CDF. The proposed model enabled statistical analysis after the division of the large sample into smaller samples. Also, it provided a clear explanation of how the RSS distributions of individual aircraft relate to the aggregated distribution. The obtained model can be employed to calculate a link budget and in the coverage design of surveillance systems. For example, adding intoP r which is calculated for a candidate receiver position and expected flight yields a prediction of the RSS distribution.
Suggested future work is investigation on the statistics of the proposed model in various measurement scenarios. Parameters affecting the statistics should be investigated together with the propagation mechanism behind. As a data source for the future work, voluntary measurement platform such as OpenSky [33] is promising. However, variations and uncertainty in measurement condition such as the surrounding environment and characteristic of the receiver and antenna need to be eliminated or considered in the modeling. Also, automatic grouping of flights needs to be developed so as to achieve meaningful statistical analysis. These are technical challenges to realize a truly large-scale analysis.