{'fraction_non_alphanumeric': 0.029655619908197352, 'fraction_numerical': 0.008002691595972219, 'mean_word_length': 4.850133558273583, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Probabilistic Determination of Maximum Safe Altitudes for Unmanned Traffic Management

This paper considers the problem of automatically determining maximum safe altitudes for unmanned aircraft operating within the aerodrome environment. First, the spatial-temporal behaviour of manned aircraft is modelled at each airspace location defined by a planer polar grid with centre at the aerodrome reference point. Second, an arbitrarily defined vertical separation threshold, unmanned altitude uncertainty and target level of safety (defined as the encounter likelihood) is used to back derive the corresponding maximum safe altitude at each location. Example results are presented for a major Australian aerodrome, including a comparative study with an existing quantitative approach currently considered by regulators. Results highlight the utility of the proposed approach for practical applications that support unmanned traffic management such as automated facility map generation and flight approval.


I. INTRODUCTION
Unmanned aircraft have recently been proposed as a feasible approach to many tasks (emergency response, infrastructure inspection, media etc.) and as an alternative transportation solution for people (Urban Air Mobility) and packages within large urban areas. To reach a scale where the business case becomes viable, unmanned aircraft will need frequent access to more airspace and receive a range of safety related services provided by unmanned traffic management (UTM) systems [1]. Due to collision risk concerns and technology gaps, airspace access is currently restricted through the establishment of static exclusion zones or maximum safe altitude limitations cast over a grid-based airspace partitioning (LAANC 1 and LAANC-like facility maps). Flight authorisation and approval is then granted through direct negotiation with the regulatory authority and Air Navigation Service Provider (ANSP) or via an authorised UTM service provider. How to determine these maximum safe altitude is non-uniform across regulators with both qualitative and quantitative methods being proposed. Additionally, derivation of the maximum safe altitude can be time-consuming and costly in the absence of an automated or semi-automated approach.
In manned aviation, probabilistic collision risk models have been used to examine the implications of changes to airspace designs, route locations, aircraft performance characteristics and separation standards [2]. By adopting manned aviation processes for airspace design and separation standard development [3]- [5], it is reasonable to suggest similar probabilistic collision risk models could offer a feasible solution to maximum safe altitude determination in the unmanned case.
The problem is essentially two-fold. First, an appropriately parametrised model is required to estimate the encounter likelihood, encounter rate or collision probability. Second, the same model needs to be adapted or restructured to determine maximum altitude limits. To avoid over-simplification, relevant uncertainty must be captured in the model. To provide utility, parametrisation should only require key variables of interest such as target level of safety and separation threshold. This paper provides a solution to these challenges by answering the following two questions: A. How can the encounter rate E and associated collision likelihood L between manned and unmanned aircraft be estimated given uncertainty on the spatial position of all aircraft, frequency of operation and desired separation?
B. How can the maximum unmanned operating altitude H * be found such that a given target level of safety, defined as an encounter rate E * or likelihood L * , be achieved?
To evaluate the proposed approach and provide a performance analysis, a case study is provided for a major Australia aerodrome. To this end, a final question is posed and answered: C. What is the difference in maximum unmanned operating altitude found using the proposed method and alternate quantitative methods and do the results align to expectations?
To the best of the authors knowledge, the proposed approach is currently the most comprehensive and aviation-compatible approach to determining maximum safe unmanned altitudes. Subject to appropriate scrutiny and similar to their manned counterparts, the derived models can provide a uniform and repeatable approach across international jurisdictions with the added benefit of transparency and relatively straightforward implementation using modern technology and digital platforms. It is hoped that this work further compliments existing risk assessment tasks posed and conducted by operators, regulators and air navigation service providers (JARUS, ICAO, CASA, Airservice Australia etc.) [6], [7]. This paper is structured as follows. Section II provides the problem background and literature review. Section III details the data processing and spatial-temporal traffic modelling approach used for manned and unmanned aircraft. Section IV and V present the methodology used to determine the encounter likelihood (or rate) and maximum unmanned operating altitude respectively. Section VI provides a comparative study of the proposed approach for maximum unnamed altitude determination with an alternate quantitative method. Further notes regarding the utility of the proposed approach and conclusions are given in section VII and VIII respectively.

II. BACKGROUND
Determining the maximum operating altitudes for unmanned aircraft in urban airspace is a very new and important problem. Indications are that methodologies should be based on encounter likelihood due to mid-air collision risk concerns [6], [7], yet limited information is available on how to accomplish this. Perhaps the most relevant quantitative approach to determining maximum safe altitudes was recently commissioned by the FAA and presented in [8]. The method leverages manned traffic surveillance data and attempts to include uncertainty on both the manned and unmanned vertical position. The approach is simplistic and easy to implement, but does not explicitly consider temporal information or provide estimates of encounter likelihood. Essentially, the encounter likelihood is unknown and does not directly underpin the maximum altitude such that a variable level of safety is achieved between adjacent airspace volumes.
Determining the maximum operating altitude (i.e. airspace class limits) for manned aircraft is a complex airspace design problem. Embedded in this problem are a number of complimentary tasks, whose solutions could be mapped across or augmented to address the related unmanned problem. The most relevant tasks relate to route design and separation standard development. To solve these problems, a set of appropriate analytical models backed by operational data are used to derive suitable solutions. For example, separation standards have been derived using purely analytical probabilistic encounter models [9]- [10], before data-driven techniques were used to refine and better parametrise the models and verify suitability. The models are comprehensive, with the ability to capture occupancy and position uncertainty through data analysis and flight performance characteristics. By combining these features, the models can derive required separation standards that meet a given target level of safety defined by an acceptable encounter rate or likelihood.
Despite the significance of such methods in manned aviation, limited work has attempted to explore or extend the such approaches for unmanned applications. A reasonable suggestion would be to consider similar underlying encounter models to derive separation standards for unmanned aircraft operations. This was suggested for unmanned-unmanned separation standard development in support of unmanned traffic mangement [11], before the encounter models were augmented and analysed in [12]. Although useful, the latter approaches cannot be used without modification to investigate manned-unmanned separation. Initial work in this direction was presented in [13] to characterise encounter liklehood in rectangular grids, then refined in [14] to investigate the impact of unmanned navigation error on encounter likelihood when operating beneath runways splays. In both cases, the encounter rate estimation methods required the unmanned altitude to be defined. Finding both a separation standard and safe altitude would then require numerous parametric studies.
In this paper, we consider the above findings and propose a new method to automatically determine maximum unmanned operating altitudes in dense air traffic regions like aerodromes. The approach is based on findings from [14] and can overcome some of the drawbacks of [8]. Importantly, the proposed approach is quantitative and strongly aligned to well-established and accepted techniques used in manned aircraft separation standards development. Example results are provided using real data and a comparison study to [8] is given. Further insights are provided as to the suitability of the approach and some common drawbacks.
Further Information Other notable collision risk modelling approaches include gas models and thier derivates [15], [16] but they do not scale well to dense airspace and small operational volumes. Simulation based methods [17], [18] or those that require simiulating specific unmanned tracks [19], [20] can require an infeasible number of simulations (or unmanned tracks) and multiple parameter sets to accurately quantify encounter likelihood. Importantly, it is not immediately evident how any of the aforementioned approaches can be augmented to derive maximum safe altitude directly (i.e. without further simulations and/or parameter sweeps).

A. Data Processing
Air traffic data recorded by The Australian Advanced Air Traffic System (TAAATS) is considered in this paper. The data consists of position reports, flight plans and other relevant information to support Air Traffic Management in Australia. Data processing occurs in a number of stages. First the data is partitioned to include only position reports (surveillance data) around Brisbane International Airport. This data is filtered to isolate all traffic movements in December 2015 within five and a half nautical miles and below 3000' AMSL of all runway thresholds. Aircraft identification tags and report timing are then used to condense the position reports into an ordered set of flight tracks F j where j ∈ {1, ..., J}. The data is further refined to remove all flights with less than two data points such that J = 15, 804 total trajectories.
Each trajectory is then checked for intersection with each radial band or boundary r * ∈ {0.0, 0.5, 1.0, 1.5, ..., 5.0} nautical miles (i.e. Δr = 0.5 nm) centred at the aerodrome reference point. The intersection points are then used as query points to re-sampled the trajectories in position, time and velocity using a linear scattered interpolant without extrapolation. The result is a set of trajectories represented by a sequence of range r, azimuth θ, altitude h and speed v values with a corresponding time stamp. Data points at each radial are then merged before being partitioned into a set of azimuth bins θ * ∈ {0, 20, ..., 320, 340} degrees (i.e. Δθ = 20deg).
The subset of data D r,θ denotes all data for radial r and azimuth bin θ such that aircraft behaviour within circular segments measured from the airport reference point can be analysed. Control of the segment volume is via Δr and Δθ, but the segments will not be uniform like a rectangular grid. An example result of the data processing for a subset of radials is shown in Fig. 1.

B. Manned Traffic Modelling
Consider the altitude of manned aircraft at each radial and azimuth to be a random variable H(r, θ). Given the potential for error on the surveillance data measurements, manned altitude estimates can be expressed as whereĤ(r, θ) is the altitude from the resampled position reports and v ∼ N (0, σ 2 m ) captures the associated uncertainty on the surveillance data. Assuming the measurement error to be equal at each radial and azimuth, random samples from the error distribution can be added to the resampled data. The probability distribution of H(r, θ) denoted by f r,θ can then be approximated using non-parametric approaches by fitting a kernel density estimate [14], [21] to the augmented data such that H r,θ ∼ K(n, w) and where n is the data sample size and Selecting the kernel function K in K w as a Gaussian distribution with w as the bandwidth or smoothing parameter then For each radial and azimuth segment, the bandwidth w of the kernel is optimised by minimising the mean integrated squared error (MISE) (or L 2 risk function) between the kernel estimate and an unknown underlying density function.
To ensure the data is appropriately bounded, each distribution must be truncated to avoid negative altitude values with a probability greater than zero. Additionally, the measurement and resampling error on the surveillance data can be included in the truncation to enable a variable level of conservatism to be explicitly included in the manned traffic model. The truncated manned altitude distribution f r,θ (h) at each radial and azimuth can be defined such that with support a and b such that and D r,θ (h) is the set of resampled data points at each radial and azimuth. A depiction of the truncation process is given in Fig. 2(a) and the resulting normalised histogram of samples taken from each f r,θ (h) for dataset D is shown in Fig. 2 As expected, we see the data disperse over a wider altitude range for aircraft located further from the aerodrome reference point. Additionally, the shape of the underlying distributions resemble data verified analytical distributions used in manned separation analysis [3], [4]. The crossing rate n r,θ defines the number of crossing per hour at each radial and azimuth. The crossing rate is more difficult to approximate due to subjectively in how to define the variable and implications when data is missing. Taking a conservative approach defines the maximum crossing rate recorded over any hourly interval and any day where t ∈ {1, ..., 23} and d ∈ {1, ..., M } defines the set of hourly intervals and days respectively. Taking a more refined yet suitably conservative approach, n r,θ is given by the rate parameter of a Poisson process such that where N t ∼ Po(λ r,θ,t ) is the number of crossings at hour t ∈ {0, ..., 23} on any day. In this case, the crossing rate is essentially defined by the maximum average crossing rate recorded over any hourly interval. The remaining option is to simply define a fixed crossing rate based on domain knowledge or operational experience such that This can enable a non-zero but negligible vertical overlap probability and simplify the calculation of p r,θ , L r,θ and H r,θ .

C. Unmanned Traffic Modelling
Consider the unmanned aircraft altitude at any location to be a random variable H subject to altimetry system errors and weather variation. The actual altitude of the unmanned aircraft can then be modelled with an arbitrary probability distribution f u (h ), which may consist of a Nominal, Gaussian, Double Exponential (Laplacian), Double Double Exponential or Mixture model (see Remark 2). Assuming a Double Exponential model, then H ∼ L(μ u , σ u ) and the unmanned aircraft altitude probability density function f u is given by where μ u is the mean, σ u is the standard deviation and d is the diversity parameter. Note, the model assumes the unmanned altitude distribution and thus altitude keeping ability is the same for all locations and all unmanned aircraft. If the unmanned altitude keeping ability or uncertainty is assumed to be dependent of location then H → H (r, θ). The derivation presented here still applies, however each radial and azimuth region would contain a unique manned f r,θ and unmanned f u → f u,r,θ altitude distribution. Using f u,r,θ would allow the unmanned aircraft altitude keeping performance to vary which may be required subject to the navigation equipment and operating region (GNSS degradation at low altitude in urban environments etc.). Remark 2 Typically double exponential and Gaussian models are used to model individual aircraft Altimetry System Errors (ASE) in manned aviation [3], [4]. Aircraft fleets and the inclusion of additional error terms then take on more complex models (i.e. Double Double Exponential, Mixture models etc.).
In the absence of unmanned data, it is reasonable to assume a similar set of distributions for unmanned altitude with different parameterisations. This is not a drawback but instead offers flexibility in the modelling regarding parameter selection [14]. The models derived in later sections can be used to inform unmanned altitude keeping error requirements to meet a target level of safety and maximum altitude level.

IV. ENCOUNTER RATE AND LIKLIHOOD DETERMINATION
This section derives an algorithm that can compute the encounter rate E r,θ and associated encounter likelihood L r,θ between manned and unmanned aircraft at location x(r, θ). The customisable (user defined) input parameters p include vertical separation distance s z , manned f r,θ (h) and unmanned f u (h ) position distributions. If traffic data is available, measurement uncertainty and processing errors (i.e. σ m ) can be used to shape f r,θ (h). Otherwise, predefined distributions can be used for both f r,θ (h) and f u (see Remark 1). As an example, a vertical separation distance s z , manned altitude data error σ m , unmanned altitude distribution f u , mean μ u and error σ u can be defined to return a corresponding E r,θ and L r,θ .

A. Methodology
An encounter occurs when the manned and unmanned aircraft are separated by less than a defined separation standard in all dimensions simultaneously. A Cartesian coordinate frame is typically used such that lateral s y , vertical s z and longitudinal s x separation thresholds are defined. Separation thresholds could also be specified in polar coordinates such that radial s r , azimuthal s θ and vertical s z separation thresholds apply (see Remark 3). Given the analysis being pursued, it may be appropriate to consider some of these dimensions having assumed overlap in the remaining dimensions.
Consider a manned and unmanned aircraft co-located at a given radial and azimuth. In this case, the aircraft do overlap in the radial and azimuth dimensions so an encounter is determined by whether the aircraft overlap in the vertical dimension. The likelihood or probability of an encounter then depends on the altitude distribution of the manned f r,θ and unmanned f u aircraft and the vertical separation threshold s z .
To determine vertical overlap probability, the approach presented in [14] is adopted. The altitude for manned and unmanned aircraft are independent continuously distributed random variables with probability distributions defined earlier as f r,θ (h) and f u respectively. The vertical separation distance z between the manned traffic and unmanned traffic is also a random variable Z(r, θ) with probability distribution g r,θ (z) such that where (14) defines the convolution between f r,θ and f u and is of the same form defiend in [4], [3]. Depending on how the manned and unmanned distributions are defined, it may be simpler to sample each distribution then fit an appropriate (non-parametric) distribution to the data set Z i (r, θ), i ∈ {1, ..., K}. The resulting approximationĝ r,θ (z) for g r,θ (z) is well aligned provided the sample size K is sufficiently large. The probability p r,θ that the separation is less than the threshold s z can then be found using the cumulative distribution G r,θ (z) such that (15) where the lower limit at −∞ is retained as separation distances less than zero infer unmanned aircraft are above the manned aircraft by (13). This is undesirable and can be considered equivalent to a loss of separation. Alternately, if an expression for g r,θ (z) is not available, G r,θ (s z ) can be approximated by dividing the number of samples K z whereẐ(r, θ) < s z by the total number of samples K (see Remark 5).
In either approach, the resulting probability is not an encounter rate, and instead an instantaneous encounter likelihood or probability p r,θ that assumes a single manned and unmanned exist at each radial and azimuth for all time. To find the encounter rate E r,θ , the number of manned aircraft crossings n r,θ at each radial and azimuth per unit time is required. The expected encounter rate can be approximated as E r,θ = n r,θ p r,θ (16) where each crossing event is independent and the probability of overlap is constant for each crossing. Furthermore, encounters between multiple unmanned aircraft are neglected. The collision likelihood L r,θ at each radial and azimuth can then be adjusted to account for the crossing rate such that the number of collisions per hour M follows a Poisson approximation where and the collision likelihood L r,θ is given by The probability of no encounters P (M ≤ 0) can be found by applying the Chernoff bound argument such that when substituted into (18) gives Using (16) and (20) the encounter rate and likelihood can be found respectively, where each expression has embedded uncertainty on the manned and unmanned locations and the separation threshold. The aforementioned parameters are essentially inputs to the model such that the approach is generic.

Remark 3 Separation distances can be mapped between
Cartesian and polar coordinate frames. The lateral separation distance s y maps to arcs along the radials associated with an angle s θ from the origins such that s θ = 2sin −1 (s y /2r). The longitudinal separation distance s x is analogous to the radial separation s r . The vertical separation distance s z is identical in both frames.

B. Example Results
The first study considers the quantification of encounter likelihood L r,θ between manned and unmanned aircraft at each radial and azimuth. The manned altitude distribution f r,θ is found using the non-parametric density estimation method presented earlier where measurement error is included such that σ m = 30 feet. The unmanned altitude distribution f u is defined such that f u ∼ L(μ u , σ u ) where L represents a truncated double exponential distribution with μ u = 400, σ u = 30 feet and support vector (μ u − 2σ u μ u + 2σ u ). The separation threshold s z = 100 feet is fixed. Results using data set D are depicted in Fig. 3 Table VIII. Summarising the results and their implications with respect to Fig. 3:

(a)-(c) with parameters defined by Analysis A in
No Crossing Rate (a): The instantaneous encounter likelihood p r,θ generally increases toward the aerodrome reference point. Higher probabilities are observed along azimuth segments that correspond to the approach departure paths of the two runways located above (32/14) and to the right (01/19) of the aerodrome reference point. Given the instantaneous encounter likelihood assumes that a manned and unmanned aircraft exist at each radial and azimuth, these results are intuitive and align to expectations as aircraft tend to be lower near the runways (see Fig. 1).

Nominal Crossing Rate (b):
The encounter likelihood L r,θ when crossing rate has been considered provides a more realistic representation of operations. For a nominal crossing rate, segments aligned with the major runway, L r<3,180−200 have doubled in some cases whilst segments further from the main flight path L r<3,220 have reduced by an order of magnitude. Very near the aerodrome an unexpected reduction in L r<1.5,θ is observed, but can be explained by the reduced data at the inner radial and the discrete nature of the segments. The incomplete (not missing) data essential lowers n r<1.5,θ and thus L r<1.5,θ , but has lesser impact on p r,θ . The manned altitude distribution f r,θ is unlikely to have changed significantly with the reduced data at these locations, provided the sample is a sufficient representation of expected altitude (see further notes as well).
Fixed Crossing Rate (c): The encounter likelihood L r,θ using a fixed crossing n r,θ = 25 improves likelihood estimates for inner radial segments and segments covering major flight paths. Essentially the incomplete data issues using the nominal crossing rate are avoided. However, the fixed rate also produces overly conservative estimates for segments further from the runways and main flight paths. Importantly, using a fixed rate could provide a maximum or worst case encounter likelihood estimate by simply fixing n r,θ to be the maximum expected rate within the aerodrome region.

V. MAXIMUM ALTITUDE DETERMINATION
This section derives an algorithm that can compute an upper altitude that when enforced would guarantee the encounter rate E r,θ or collision likelihood L r,θ is below a target level of safety E * or L * for locations x(r, θ). Two cases can be considered with respect to the target level of safety. The encounter likelihood may be explicitly defined to some acceptable level (1 × 10 −6 , 1 × 10 −4 etc.). Alternatively, the encounter likelihood may be unknown but implicitly defined by current exclusion or unmanned fly zones. For example, many aerodromes prohibit unmanned aircraft flying within a nominal radial distance from the aerodrome reference point or runway thresholds. A nominal altitude (i.e. 400 feet) may be allowed outside and on the boundary of the excluded region. The encounter rate associated with operating at this boundary and maximum allowable operating altitude could be considered as the target level of safety. As an example, a desired encounter rate E * or likelihood L r,θ would be specified alongside the vertical separation distance s z and standard deviation on unmanned σ u and manned σ m altitude to return H * r,θ .

A. Methodology
Estimating the maximum safe unmanned aircraft altitude H * r,θ for each radial and azimuth requires two main steps. No further assumptions are required on the input parameters and the process can be simplified by acknowledging the statistical properties associated with the summation of random variables.
The first step is to find the required instantaneous encounter likelihood given by a target level of safety L * and subject to the observed crossing rate. To do this, consider substituting the target level of safety defined as an acceptable encounter likelihood L * into (20) and rearranging such that As the crossing rate is fixed or observed, the expression can be evaluated to give the required instantaneous vertical overlap probability that would result in meeting the target level of safety. Recall, that p r,θ was the output from evaluating the cumulative distribution G r,θ (s z ), so we could use the inverse cumulative distribution or quantile function G −1 r,θ (p r,θ ) to find the unmanned altitude H * r,θ that corresponds to p r,θ such that where the quantile function is defined as and gives the separation distance where p r,θ is less than the target level of safety. It should be evident that this approach requires the nominal unmanned altitude to be defined upfront (i.e. μ u ) and then adjusted as per (24). This is not necessarily a drawback as selection of another unmanned nominal altitude would have just resulted in a different adjustment. However, we can use some intuition to remove this requirement. Re-consider the derivation of the distribution for vertical separation. Instead of defining the unmanned altitude distribution f u (h ), a separation distribution f s (z) of the same type and standard deviation could have been defined but with a mean equal to the separation threshold s z . If the same convolution is conducted as in (14) where Z is a random variable drawn from distribution f s (z) then leaving an expression for the unmanned altitude distribution at each radial and azimuth location f u,r,θ (h ). Essentially, the uncertainty on the unmanned altitude is shifted to uncertainty on the separation threshold s z such that Remark 2 translates accordingly. Following the same process as above, the required unmanned altitude H * r,θ to meet L * is thus

B. Example Results
The second study considers the determination of H * r,θ for unmanned aircraft at each radial and azimuth. Again, the manned altitude distribution f r,θ is found using the nonparametric density estimation method including added measurement error σ m = 30 feet. A separation distribution f s is defined such that f s ∼ L(μ s , σ u ) where L represents a truncated double exponential distribution with μ s = 100, σ u = 30 feet and support vector (μ u − 2σ u μ u + 2σ u ). The separation threshold s z = 100 feet is fixed. Results using data set D are depicted in Fig. 4(a)-(c) with parameters defined by Analysis B in Table VIII. Summarising the results and their implications with respect to Fig. 4: Equivalent TLS (a): Unmanned altitudes are derived by considering the TLS to be equivalent to the maximum observed encounter likelihood at a nominal radius of 3nm from the aerodrome (L * = max{L 3,θ } ≈ 5 × 10 −2 ). With such a liberal TLS, H * r<1,θ = 0 north and south of the reference point as expected as this constitutes the major approach departure paths. To the east of the aerodrome reference point H * r<1,θ remains low whilst the west contains higher altitudes. This can be explained by the runway configuration and the limited traffic in the area between the runways and tower (i.e. inner west). From approximately 1-3.5 nm the segments are quite variable and beyond 3.5nm H * r>1,θ ≥ 400 with the exception of one segment surrounding a western heliport (medical).
Specified TLS (b)-(c): Unmanned altitudes are derived by specifying a desired TLS L * = 1 × 10 −4 which corresponds to 1 encounter in 10,000 hours. In (b) a nominal crossing rate is used whilst in (c) a fixed crossing rate is used. As expected the general trend is a reduction in H * r,θ , with the application of a fixed crossing rate providing more conservative maximum altitude estimates. Zero or very-low altitude segments span the majority of the runway splays, with greater variability surrounding the minor runway (32/14). Even with the maximum crossing rate applied and the resulting overestimate of encounter likelihood, many segments west and south east of the aerodrome from 3.5 nm exist where H * r,θ ≥ 350 providing useful airspace for unmanned operations.
Further Observations (a)-(c): The results suggest the proposed methodology is a feasible approach to automatically determine safe altitudes and clearly show the ease at which customisation w.r.t risk tolerance can be achieved. The method is consistent across the airspace segments so can at least provide a relative measure of encounter likelihood between segments. For example, there are a number of segments at 5nm where H * 5,θ < 300. This may be surprising but is consistent with the encounter likelihood estimate as the corresponding segments in Fig. 3(a)-(c) have a probability in the order of 10 −2 . Lastly, the results quantitatively support observations that could have been derived from inspection of the raw or re-sampled data shown in Fig. 1. The data suggest there may be more segments at which unmanned aircraft could operate to higher altitudes with low risk west of the aerodrome.

VI. COMPARATIVE ANALYSIS
This section provides a comparative analysis between the proposed method for maximum safe altitude determination and an alternate quantitative method considered by regulators [8].
The objective is to further identify the merits and drawbacks of the two approaches and help determine their suitability for use and important considerations in their application.

A. Alternate Methodology
The alternate approach employs a simply data search and summation to estimate the maximum safe altitude in any grid cell or radial segment such that whereD r,θ denotes a reduced data set with the lowest 5% of altitudes removed from D r,θ and s * z is a scalar variable that encapsulates the separation threshold and the uncertainty on the manned and unnamed altitude. This approach is implemented using data set D with two changes to ensure the comparison is meaningful. Firstly, the data extrapolation method and ground elevation adjustment are not employed (see Remark 6) and secondly a radial grid is used instead of a rectangular grid. With these alterations, the two approaches can be compared as the same data source, processing and representation is employed.

Remark 6
The extrapolation scheme in [8] produces linear trajectory segments directly toward the aerodrome reference point from positions within the 2nm radial. This approach will introduce erroneous and unrealistic data. A similar extrapolation scheme that aligns the missing data to the appropriate runway would be more suitable and improve issues with missing data outlined earlier. This is not pursued in this paper and left as further work.

B. Example Results
The third study considers finding the difference ΔH * r,θ in maximum altitude at each radial and azimuth, derived by subtracting H * for the proposed approach from H * for the alternate approach. Again, f r,θ is found using the non-parametric density estimation method with added measurement error σ m = 30 feet. A separation distribution f s is defined such that f s ∼ L(μ s , σ u ) with μ s = 100, σ u = 30 feet and support vector (μ s − 2σ u μ s + 2σ u ). The separation threshold s z = 100 and separation-uncertainty threshold s * z = 200 are fixed. Essentially the manned and unmanned uncertainty is considered to account for 100 feet such that a separation buffer equal to s z remains to keep comparisons consistent. Results using data set D are depicted in Fig. 5(a)-(f) with parameters defined by Analysis C in Table VIII. Summarising the results and their implications with respect to Fig. 5: Complete Data Set (a)-(c): Using D for both the alternate method and the proposed method provides some interesting results. Using (a) and (c) suggests H * from the alternate method aligns well with H * with L * = 1 × 10 −4 using the proposed approach. The majority of segments have a zero altitude difference which suggests the alternate method implicitly assumes a target level of safety of approximately 1 × 10 −4 . Encounter likelihoods were not explicitly calculated or stated in the alternate method, and the proposed provides this valuable insight. A small number of segments in the main runway (01/19) path show a slightly increased maximum altitude when using the proposed approach yet L * = 1× 10 −4 is retained. This is because crossing rate is not considered in the alternate approach. A single aircraft crossing at a low altitude dominates the calculation in the alternate approach, whilst the proposed approach essentially has knowledge of the relatively low frequency of such occurrences. This is the question of possibility versus probability. Using (b) where L * = max{L 3,θ } ≈ 5 × 10 −2 provides a much higher risk tolerance and as expected using the proposed approach suggests many locations with up to 500' difference in altitude.

Reduced Data Set (d)-(e):
In this comparison, only the alternate approach uses the reduced data setD. Using (d) and (f) suggests H * from the alternate method aligns well with H * with L * = 1 × 10 −4 using the proposed approach for inner radials. This suggests that the removed altitude data was very similar to the data that remains. As the radius increases, this is no longer the case. Focusing on the flight path associated with the main runway, through the arbitrary removal of data the alternate approach has essentially increased the risk tolerance. One explanation is that the removed data contained multiple unusually low altitudes that are appropriately modelled via f r,θ . Comparison of maximum unmanned operating altitudes determined using proposed approach and alternative approach presented in [8]. (a) and (d) depict the maximum altitude defined by the alternative approach using 100% and 95% of the data respectively ( -). (b) and (e) provide the altitude difference between the alternate approach using 100% and 95% of the data respectively, and proposed approach when L * = max{L 3,θ }. (c) and (f) provide the altitude difference between the alternate approach using 100% and 95% of the data respectively, and proposed approach when L * = 1 × 10 −4 . Zero differences ( ) are visible along with numerous segments where the alternate approach provides higher ( -) and lower ( -) maximum altitudes.
in the proposed approach. Essentially, when data is removed the alternate approach does not provide a uniform adjustment of encounter likelihood across the segments. Further evidence for this can be seen in (e). The same segments in the runway path now have the equivalent or higher H * , but the comparison is now made with a risk tolerance two orders of magnitude greater (i.e. L * = max{L 3,θ } ≈ 5 × 10 −2 ). Of note, it would be useful to compare the alternate approach with the proposed approach usingD to further verify this observation.

VII. FURTHER NOTES AND CONSIDERATIONS
The methodologies, example results and comparative studies provide an overview of the merits of the proposed approach. Some further considerations are given below, some of which apply to general data and grid-based methods.

A. Overlap Dimensions
The encounter likelihood and rates assume the positional overlap in the radial and azimuth dimension to be unity. The assumption in the radial dimension is justified as a single radial band is used. The assumption is not well justified in azimuth as the aircraft pair could be anywhere within the segment such that separation could be achieved in this dimension. The situation is more likely the further away from the aerodrome. If the probability of overlap in azimuth were included in the model, the encounter likelihood and rate would decrease. However, care should be taken to consider the change in size of the data at each radial. Additionally, this may not be desired as only assuming the vertical dimension will provide a more conservative estimate.

B. Data and Grids
The size and completeness of the raw data needs to be considered in applying the proposed method, particularly with respect to crossing rates. Care should be taken to either provide a suitable extrapolation scheme or individually fix the crossing rates per segment based on operational experience. The data size should be sufficient to cover all aerodrome runway configurations or partitioned into specific time windows to produce altitude maps relevant to specific conditions. The grid size, type and reference point will also alter the results. Consider a grid boundary under a main traffic flow. It would be reasonable to assume that half of the data points exist in neighbouring segments, so the resulting instantaneous encounter likelihood in each segment may be representative, but the corresponding crossing rates will be halved. Additionally, the non-uniform and discrete nature of a polar representation may be practical, but means the arc length used to count crossings and fit altitude data varies at each radial and can skew results. Lastly, there is inherent ambiguity on the appropriate altitude on grid boundaries if adjacent grids are dissimilar. Such discritisation issues are observed with any grid based method.

C. Parameters
The tunable parameters include the uncertainty on manned and unmanned altitude, target level of safety and separation threshold. As a result, the approach can be used to conduct parametric studies on any of the parameters. The unmanned altitude distribution was fixed in this analsyis but could have been varied as in [14]. Conducting such a parametric study would help uncover the required navigation performance for unmanned aircraft (in altitude) should a specific target level of safety be specified. The altitude maps derived in this paper suggest that all unmanned aircraft with altitude keeping ability to within less than 30' can operate at derived altitudes.

D. Applications
The proposed approach offers a uniform method to aid development, implementation and substantiation of low-level airspace characterisations or Airspace Encounter Categories (AEC) within risk assessment frameworks such as those proposed by JARUS [7]. This directly supports regulators and Air Navigation Service Providers (ANSPs) in the determination of where unmanned aircraft can safely operate, including altitude restrictions, low-level airspace designs (facility maps) and automated flight approval/authorisation. VIII. CONCLUSIONS This paper presented a probabilistic encounter modelling technique that can automatically determine maximum safe operating altitude for unmanned aircraft. The techniques were evaluated using real surveillance data and compared to a recently proposed alternative approach supported by some regulators. The results demonstrated the utility of the techniques which could be used to develop new procedures, products and services for unmanned traffic integration and management.
This work is aimed at moving towards an automated, unified and reliable approach for determining maximum safe operating altitudes for unmanned aircraft, such that the process is efficient and repeatable whilst balancing safety and flexibility. Additionally, the work provides evidence to support quantitative encounter modelling, which can then help verify qualitative methods and provide a common language for air navigation service providers, operators and regulators to discuss air risk.