Quantifying Risk of Ground Impact Fatalities for Small Unmanned Aircraft

One of the major challenges of conducting operations of unmanned aircraft, especially operations beyond visual line-of-sight (BVLOS), is to make a realistic and sufficiently detailed risk assessment. An important part of such an assessment is to identify the risk of fatalities, preferably in a quantitative way since this allows for comparison with manned aviation to determine whether an equivalent level of safety is achievable. This work presents a method for quantifying the probability of fatalities resulting from an uncontrolled descent of an unmanned aircraft conducting a BVLOS flight. The method is based on a standard stochastic model, and employs a parameterized high fidelity ground impact distribution model that accounts for both aircraft specifications, parameter uncertainties, and wind. The method also samples the flight path to create an almost continuous quantification of the risk as a function of mission flight time. The methodology is exemplified with a 180 km flight in Danish airspace with a Penguin C aircraft.


Background
As drones are becoming ubiquitous in the airspace and as more applications of drones are focused on longer flights the need for reliable, detailed, and quantitative risk assessment methods is growing. While there are several methods available (borrowing from manned aviation) for qualitatively identifying and categorizing hazards and mitigating risks there is comparatively little methodology available for actually determining the probability of fatality for a particular drone operation in particular methods that lend themselves to reverse engineering to allow for pinpointing how a given flight scenario could be altered to reduce the risk. This is in large part due to the seemingly endless list of variables that enters such a method were it to be completely comprehensive.
This particular work is prompted by a decision by the Danish Transport Construction and Housing Authority Anders la Cour-Harbo alc@es.aau.dk 1 Aalborg University, Fredrik Bajers Vej 7C, 9220 Aalborg East, Denmark in 2016 to conducts a number of trial BVLOS flight operations in Danish airspace with the expressed intent of using these to pave the way for routine BVLOS operations by companies and other non-state actors. At that time no permissions had been granted for BVLOS flight in Denmark except for flights in Greenland and individual flights confined to specific routes and dates. The trial flights would be based on thorough analyzes of the risks, and the idea was that subsequent flights would be conducted routinely using the same risk assessment methodology.
This work presents the developed methodology. It is based in part on a previous publication [21] applying the methodology to a power line inspection flight done by the Danish company Heliscope. The present work expands and details the methodology, and the example flight has a different and somewhat longer flight path as well as a more professional drone platform.

Previous Work
There are numerous works on how to conceptually approach the challenge of determining the risk of an unmanned aircraft flight. Much is borrowed from the world of manned aviation that has been conducting risk management for decades. A number of examples of risk assessments and quantifications for unmanned aircraft include the following.
In [4] a study for ground impact fatalities resulting from power failure and subsequent uncontrolled glide is presented. A study on the impact area for a general uncontrolled descent, including a buffer zone, is presented in [17]. A method for automatically finding a proper landing area for an aircraft in emergency descent is shown in [33,34], and the ability of a fixed wing aircraft to glide to a designated emergency landing area is presented in [11]. In [3] a method for determining a no-thrust flight trajectory to reach a particular landing spot is presented. The barrier bow tie model also used in manned aviation risk assessment is presented in [8]. A study of trajectory models for explosive debris [29] attempts to determine the impact point based on initial conditions. [10] addresses the lack of an accepted framework and provides some guidelines for how to apply existing models to manage the risk. In [9] a comprehensive description of how to manage the risk of unmanned aircraft operations, including 'the systematic application of management policies, procedures and practices to the activities of communicating, consulting, establishing the context, and assessing, evaluating, treating, monitoring and reviewing risk.' This work also presents a series of quantification of existing risks for various types of aviation. Metrics for safety, including hazard metrics and risks metrics are presented in [22], in [14] a software safety case is developed, and in [13] a generic safety case is presented based on experience with NASA unmanned aircraft missions. The uncontrolled descent of unmanned aircraft into populated areas have been the subject in a number of publications. This includes [18] that investigate larger aircraft through an equivalent level of safety analysis. [35] specifically looks at distribution of possible impact positions based on simulation, and [7] uses the standard statistical setup (which is also used in this work) and applies a normal distribution approach using aircraft glide parameters to model the impact position.

Current Work
The aim in this work is to go beyond a qualitative approach that merely provides a framework for risk assessment, and instead apply a quantitative approach to determine as accurately as possible, the level of safety for a given flight operation using the metric 'fatalities per flight hour', similar to what is used in manned aviation. The modeling of the probability of fatalities (POF) is done with a stochastic approach similar to many of the previous works listed above. However, the determination of the individual probabilities in the model is done using georeferenced probability density functions with high fidelity and (almost) continuously along the flight path to provide not just a single probability for the entire flight, but indeed a fatality rate along the flight path itself. This in turn allows for easy identification of what parameters adversely affect the flight, with the possibility for reconfiguring the flight path to reduce the risk of fatalities.

Methods
The basic approach used in this work is similar to numerous previous works, namely a stochastic model that joins probabilities in the causal chain from drone malfunction to a potential fatality. The specific design of this stochastic model varies from work to work. As described in Section 2.1 below, we use a fairly simple setup, where the focus is on the probabilities related to the ground impact. The model is used for four different types of flight terminating events, all described in Section 2.2, each with their own high resolution ground impact probability density function, see Sections 2.2.1 to 2.2.4, including the effect of wind, see Section 2.3, and associated event probability, see Section 2.4. Rather than a priori assuming an average population density for the entire flight (as seen in many previous works), we employ high resolution population density maps generated to fit the spatial extend of each event type, see Section 2.5. We then use the stochastic model on a (sufficiently densely) sampled flight path to determine the probability of impacting a person as a function of the flight path. Finally, mapping from person impact to POF (probability of fatality) is applied based on work in the field of forensic science, see Section 2.7.
The POF for the entire flight is determine by summing over the flight path relative to the flight time between each sample in the sampled flight path. This is done separately for the different types of flight terminating events to accommodate the varying lethality parameters associated with the manner in which the aircraft descents in each of the event scenarios. The is described in more detail along with the results in Section 3.

Overall Modeling Approach
For computing the probability of fatalities during the flight we will use this formulation similar to what is used in several of the previous works listed in Section 1.2.
where p event is the probability per time of a given event (of which we will use four), p impact person is the conditional probability that given an occurrence of one of the events that a person will be impacted as a result, and p fatal impact is the conditional probability, given one or more persons are impacted as a result of one of the events that this person suffers a fatal injury. The primary focus in this work is on the probability that a descending aircraft will impact a person. We also briefly treat the probability of fatalities as a result of impact. The probability for each event occurring during flight is measured in events per flight hour and will be assumed known, since this probability is not the focus of this work. Determining such probabilities is no simple task, and in Section 2.4 is a brief discussion of the numbers used for exemplifying this work in Section 3. For each type of event the probability of impacting persons on the ground is determined in 4 steps: 1. Based on a model of a given event type a probability density function (PDF) is computed that determines the probability of ground impact relative to the position where the event happens. 2. This PDF is subjected to probabilistic wind, resulting in a new, often bigger, two-dimensional PDF, which can be interpreted as a georeferenced probabilistic impact map. 3. This map is correlated with a population density map of sufficient resolution for the same geographical area. 4. The resulting map is integrated over the entire area and the result is modified to account for various factors pertaining to each event type.
This gives a probability for a given event to result in an impact of a person on the ground relative to the event position. As a consequence this impact probability depends on the type of the event, aircraft parameters, wind, and flight path of the aircraft. In the following sections the above steps are discussed in more detail.

Descent Event Types
We consider four rather different types of uncontrolled descents of an unmanned aircraft. Each is a result of varying types of malfunction, and each has it own model (as indicated in step 1 above). We only consider events where the operator has no control authority and the aircraft descent path therefore is governed mainly by aerodynamics or the autopilot. In either case the aircraft is not expected to return home nor to head for any designated safe impact zone.
1. Ballistic descent This is a situation where the aircraft has lost most of its lift, for instance by a wing breaking off or a motor physically separating from the aircraft. The aircraft will then enter a (close to) ballistic descent governed solely by the aerodynamics of the crippled aircraft. 2. Uncontrolled glide For a fixed wing aircraft this is loss of thrust as well as loss of power for the flight control surfaces. For a helicopter type aircraft this event could be loss of thrust on the main rotor with an autopiloted autorotation descent. In any case the airframe is structurally intact and the aircraft is assumed to enter a descent path governed by the glide ratio/autorotation descent angle as well as wind. For a fixed wing the deflection surfaces are assumed to be in a close to neutral position such as to give a straight or perhaps slightly curved glide. 3. Parachute descent This is a descent with a fully deployed parachute. It is assumed that the deployment is a result of a malfunction detection, and thus that motor(s) are turned off giving a descent path based solely on the aerodynamic properties of the parachute. 4. Flyaway This is complete loss of operator control authority of the aircraft while the autopilot continuous to operate in a mode that maintains the stability of the flight. The motion of the aircraft is controlled by the autopilot, and it may fly to its maximum range in any direction, including vertically up.
One can envision additional failure scenarios such as a spin due to a faulty actuator, loss of a vital sensor that deprives the aircraft of any useful navigation, and so on. We assume that each of the four above listed events will happen with a probability independent of where the aircraft is along the flight path. The impact point on the ground of a descending aircraft is determined by the descent models relative to the event point. The event point is defined as the spatial location of the aircraft at the time it suffers a malfunction that gives raise to one of the above scenarios. Depending on the context the event point may refer to the projection of the point onto the ground (i.e. zero AGL). The ground impact point position is modeled with a number of uncertainties that reasonably may affect the descent, and these are described below in more detail.

Ballistic Descent
The ballistic descent is modeled as described in [19]. The model has a second order dependence on speed and acts solely under influence of gravity and drag. The flight velocity prior to the event is based on the georeferenced flight path. The model for the ballistic descent is a multiple stage formulation that includes altitude, mass, drag coefficient, and frontal area. The descent path is roughly a second order polynomial modified to account for the aerodynamic properties of the aircraft. It also has probabilistic assumptions on the horizontal and vertical initial speeds as well as on the drag coefficient. The output of this model is a 1D PDF with the probability of the aircraft impacting the ground a given distance away from the event point and along a line coinciding with the traveling compass direction of the aircraft. The model is much too comprehensive to be reiterated here, so do see the reference for further details.

Uncontrolled Glide
The uncontrolled glide assumes the aircraft is descending as a glider. The aircraft is without thrust and in an aerodynamic equilibrium configuration with actuators in neutral or close to neutral positions (that is, no banking with ailerons and no turning with the rudder). The horizontally traveled distance per vertically descended distance is called the glide ratio γ , and in this work is given as a normally distributed variable with mean equal to the estimated glide ratio and a variance to accommodate for variations in elevator deflection angle and variation in drag due to possible modifications of the aircraft to accommodate missions sensors. The horizontal distance x traveled (in the wind frame) in an uncontrolled glide from altitude y is simply x(y) = γ y and the drop time is t drop (y) = x(y)/v g , where v g is the horizontal glide speed. To get the ground impact PDF the wind variations as described in Section 2.3 are applied.

Parachute Descent
Descent with a parachute is essentially a vertical drop in the wind frame. We assume that the parachute is deployed in case of a detected emergency situation with a short delay of d seconds, and that the parachute once triggered will deploy instantaneously and reduce the horizontal velocity to zero instantly (in the wind frame). While this obviously in not the case in real life the distance traveled during the deployment and deceleration phase is negligible in relation to the impact area for a parachuted descent. Assuming a standard second order drag model. The drop time from altitude y is where m is mass, g is gravitational constant, A p is parachute area, and C d,p is the parachute drag coefficient. In this work the latter is a probabilistic value. Note that if the drop speed v drop for a given parachute configuration is known the formula is very simple. As with the uncontrolled glide the ground impact probability density function (PDF) is achieved by applying the wind variations as described in Section 2.3 plus offsetting the result according to the flight direction and deployment delay.

Flyaway
A flyaway event will potentially take the aircraft to the limit of its range, which dependents largely on the fuel left and on the wind speed and direction. The model for a flyaway is in this work composed of two contributions.
First, the probability of ground impact is assumed to decreases linearly with distance from the event point, reaching zero at the maximum range. The simplicity of this model is primarily due to lack of knowledge of how a flyaway will progress. On the one hand if a flyaway most often will continue until there is no fuel left the probability of ground impact should be relatively higher close to the circumference of the flight range circle. On the other hand if the flyaway is not in a straight line, but rather more with random movements, and the ground impact point therefore can be assumed to be equally likely at any distance from the event point, the probability should decrease with the square of the distance. As a compromise a linear relation is chosen. This probability is then modified according to the wind speed and direction. Mathematically, this is modeled as where θ is the wind average direction, R max is the maximum flight range given the available fuel, v c is the aircraft cruise airspeed, and p = (p N , p E ) is a north-east position relative to the event point. Second, a probability is added that the aircraft will ascend more or less vertically, either as 'helicopter climb' or spiraling up with a fixed wing. This causes the ground impact point to be close to the event point, and this part is modeled as a normal distribution for distance from the event point with mean 0 and standard deviation σ va (vertical ascend) from the event point, and uniformly in direction. That is, the impact PDF becomes The two scenarios as linearly combined with a relation factor β, such that β = 0.5 makes each contributor with equal probability. The resulting ground impact PDF becomes An example of the resulting 2D PDF is shown in Fig. 1.

Wind
The influence of wind is significant for all four event types. Ballistic descent, uncontrolled glide, and parachute descent are all modeled directly in the wind frame, while the flyaway event is not directly assumed to be in the wind frame, but the wind speed does in Eq. 3 affect the ability of the aircraft to fly distances in certain directions. Since the knowledge about the wind for a given flight will vary somewhat depending on the circumstances we propose to use one of three wind models: 1. Direction and speed modeled with normal distributions. 2. Direction unknown (modeled uniformly), speed modeled with normal distribution. 3. Direction and speed based on historic data. The first option is useful for missions in a particular geographical area where some wind statistics is available, and for computation just prior to a flight where the actual wind might be available. The second options is useful for more generic scenarios where the geographical locations is yet to be determined. The third option applies to scenarios that take place at a known location at a known time, and where historic data for a the given time period (say a particular week or month) is available.
The models for the ballistic descent, uncontrolled glide, and parachute descent all provide a series of drop times (the time it takes the aircraft to reach the ground from the given flight altitude) resulting from the probabilistic nature of the models. For all three models the offset caused by the wind is dependent solely on these drop times, since all other aerodynamic properties are already accounted for in the models.
One way to practically compute the effect of the wind is as follows. The range of possible drop times is sampled as {t k } k=0,...,n and for each drop time t k a sampled georeferenced PDF is generated as a matrix M k representing the probability of the aircraft doing a purely vertical drop in the wind frame to impact the ground at the geolocation represented by each entry in the matrix. The ground grid of the PDF (sampling density) should match the event type and the population density matrix that appears later in the computation (see Section 2.5). The drop times comes from the individual models with a set of probabilities {p k } k=0,...,n that gives the probability of each drop time. Thus, to generate the wind effect on the individual model, we simply sum up the PDF matrices as while keeping score of the georeference of each matrix such that the resulting matrix M is also georeferenced. This approach is particularly useful for flights where the same drop times occurs often, since the M k matrices need not be recalculated as long as the drop times are within the same range.
To demonstrate how the wind effect is in practice applied to a descent an example using the parachute descent is given in Fig. 2. The procedure of applying wind to an uncontrolled glide and a ballistic descent is essentially the same.

Probability of Events
The events that render the aircraft uncontrollable and eventually lead to a uncontrolled descent each has a probability attached to them. These probabilities are difficult to estimate and difficult to measure (as that would require many flight hours on precisely the same setup). For the computations in this work we simply assume given values that per aviation tradition are measured in 'per flight hour'. The probabilities used here are estimates based on the works of others (see below) as well as the experience of persons at the institutions listed in Acknowledgements. The event probability appears in the probability computation in Eq. 1 as a scalar. As a consequence the effect of changing the event probability is simply a scaling of the resulting probability.
The probability of a flight terminating malfunction on an unmanned aircraft has been studied by a number of groups. A reliability assessment of an Ultra Stick 120 is made in [15] and [16] using failure mode effect and analysis (FMEA), with particular attention to the control surfaces and servos. No specific probabilities for an uncontrolled   7, 2). The fourth graph shows a linear combination as given by Eq. 6, and this would be the non-offset (see Section 2.2.3 on parachute descent) georeferenced ground impact PDF for a parachute descent relative to the event position at (0, 0) descent are provided, but are considered to be high. In [24] probabilities related to military unmanned aircraft are reported, and the probability of a flight terminating event is in the order of 10 −4 -10 −2 per flight hour, with the probability for smaller aircraft being somewhat higher that for the larger aircraft. A group of students showed in [26] using FMEA based on component failure rates that their Ultra Stick 120 has on average 2.17 catastrophic (flight terminating) failures per 100 flight hours. The types of failures considered relate to the ballistic descent and uncontrolled glide in the present work. In [27] the same group showed how a dedicated reconstruction of the aircraft based on a fault tree analysis could theoretically reduce the failure rate by a factor 20, and they were able to implement changes to the physical aircraft to achieve a catastrophic failure rate of 0.76 per 100 flight hours. Actuators and control surfaces are investigated in [30] and [31] where the probability of having an uncontrollable aircraft is modeled using a servo fault detection algorithm.
In [23] a method for estimating mechanical failure rates of small unmanned aircraft is presented, and an example is provided based on the 25 kg SPAARO aircraft. This example explicitly lists the used probabilities for failures of servos and deflection surfaces as well as failures of engine and battery. These probabilities were provided by two experienced RC pilots. In addition probabilities of failure for a wing bolt and main spar are theoretically derived.

People Density
Density of people on the ground is the main factor in the probability of impacting a person in the event of a crash, and for BVLOS flight operations that stretches over longer distances variation in people density can be significant during the flight. The example flight presented in Section 3 is specifically chosen to include this feature. As the impact areas for each event change over the course of the flight as well as in between flight as a result of changes in altitude, wind conditions, flight speeds etc. it is important to represent the actual people density with a reasonably fine resolution. For this reason it is insufficient to assume a fixed people density, even for smaller geographical areas, and if such assumption is made the resulting POF may be misleading and changes in the flight planning that might actually have an impact on the fatality rate may go unnoticed. The impact area for different types of descents varies hugely, i.e. a ballistic descent is typically close to the event position, whereas a flyaway can result in a descent hundreds of kilometers from the event position. Consequently, for some events the resolution of the people density map must be fairly high to give accurate estimates of the person impact probability, while for other events the resolution can be more coarse and still give accurate results.
A list of geographical coordinates of all addresses in Denmark is publicly available and this has been used to generate people density maps with varying resolution to fit the different types of descents. While a fine grained resolution will of course work for any type of event the computation time grows significantly, so maps are generated that suits the spatial extend of the impact area for each type of event. In Fig. 3 are shown three examples of such maps.
While these maps do show where people live they obviously do not show where people actually are. As this information is evidently very difficult to obtain we will make the assumptions that people are, with some probability, in the vicinity of their home, and with some probability are outside exposed to a small unmanned aircraft potentially descending. Inspired by [12] an appropriate probability of people being exposed is around 30%. This is also referred to as the shelter factor. We will also assume that the number of people associated with each address is equal to the average number of people in a Danish household. This number is 5.75 million people divided by 2.65 million households, equal to 2.17 people per household. The number of addresses is 3.3 million as some addresses are not households, but rather businesses and industry. The density map used in this work is not adjusted to account for this.

Probability of Impact Persons
For a given event type we now have a ground impact probability density function (PDF) measured in meters relative to (0, 0). By offsetting this PDF relative to the coordinates of the event point we obtain a georeferenced impact map in WGS84 coordinates (latitude and longitude). This PDF matrix is then entry-wise multiplied with the population density map (appropriately sampled matrix D) for the same area and the result is a map of the probability of impacting a 1 m 2 large person, since our population density is measured in people per m 2 . We assume that a person takes up a particular area A person that depends on the expected impact angle, and this value is multiplied onto the result. Additionally, it is multiplied with the shelter factor S, which accounts for the probability that a person is sheltered by being indoors, under a tree, in a car etc. In this work we assume a fairly high probability of being sheltered (typically indoors), since the population density map is based on addresses where people live. This will then provide the probability p impact person of impacting a person (see Eq. 1 above) given a particular event at the given event point. The computation can be formulated as where • is the Hadamard product and · is the scalar product.

Probability of Fatality When Impacted
When a drone impacts a human there is a probability that the impact will inflect injuries that will result in a fatality. Determining this probability for a given person and a given drone is not simple, partly because of the many different ways the impact can occur, and partly because the easily determined parameters, such a speed and mass, do not have a simple correlation to injury severity, because the human body reacts differently depending on the impacted body part, and the fact that injuries primarily relate to how fast and where the kinetic energy is dissipated in the body, not the kinetic energy of the impacting object itself. For a review of literature on drone-like human injuries, see [20] and [2]. A number of reasonably accurate and empirically verified The color scale is the same and goes from 0 (white) through 1 (dark blue) to 40.000 person/km 2 (dark red). Note how the density tends to grow with increased resolution due to the same number of people being registered in still smaller squares. Semi-transparent topographical information is overlaid; roads are brown, urban areas in dark yellow, forest in green, and municipality border in black. The town in the center of view is Thorsø models have been developed. One model that fits well to a drone (chest) impact scenario is [32] which uses a lumpedmass thoracic model to develop a VC parameter, where V is thorax compression velocity and C is compression relative to chest depth. The VC parameter for a given impact maps well to injury severity. For this work we have chosen to use the blunt criterion (BC) from [6] for impacts at relatively high speeds (covering ballistic descent, uncontrolled glide, and flyaway) and the area weighted kinetic energy methods (AWKE) from [2] (covering parachute descent). It would be relatively easy to substitute these for other methods, since the modeling approach used here provides all necessary parameters, such as mass, speed, impact angle, impact area, etc.

Blunt Criterion
The blunt criterion (BC) is useful because it does map kinetic energy to injury severity. It is defined as where E is kinetic energy, W is mass of impacted object, T is thickness of the body wall (in cm), and D is diameter of impacting object (in cm). According to [28] T = kW 1/3 with k = 0.6 for females and k = 0.7 for males, and according to [5] we have AIS = 1.33 · BC + 0.6. And by interpolating the fatality rates normally associated with the AIS scale [1] we can now map kinetic energy to fatality rate. Note that an adaption of BC to impacts of drones is done in [25], where a generic drone design is used to develop formulas specifically for thorax and head impacts. It does not map all the way to fatality rate, though.

Area Weight Kinetic Energy
This method is an adaption done by [2] of earlier work to better represent the posture of a person when impacted, and also maps kinetic energy to POF. Unlike the BC criterion it does not account for the size of the impacted area, and as such is more suitable for impacts of larger areas. The actual mapping used in the present work is a cubic spline with a derivative equal to zero at the end knots applied to the numbers in Table 1, which in turn is copied from [2]. Figure 4 shows the maps from impact speed to POF for the Penguin C aircraft. The impact speeds for the example flight ranges from 4.5 m/s for parachute descent to over 20 m/s for ballistic impacts. As expected the lethality of the Penguin C aircraft is close to 1 for any descent type.

Approximation of WGS84 Coordinates
While a population density map would typically be in geographical latitude and longitude coordinate system the impact PDFs are in a local north-east coordinate system, since the models operate with a relative distance measure in Euclidean metric. In order to multiply those two maps a conversion of either one is required. In this work we convert the impact PDF to latitude/longitude coordinates. This conversion does require a significant amount of computation as the location of each entry in the PDF matrix must be converted. A very fast and simple approximation is to simply linearly interpolate lat/lon coordinates between two diagonally opposite lat/lon corner coordinates of the PDF. This approximation is fairly accurate for a PDF spanning single digit kilometers, in the sense that the distance error between the true position and the approximated position (measured as Euclidean distance) is somewhat smaller than population density resolution. However, as the PDF size grows the error soon becomes significant. Figure 5 shows the difference between a full conversion and interpolation for two example PDF sizes.  Table 2

Results
The proposed method for quantifying POF lends itself to a wide range of unmanned aircraft flight scenarios. It does require reasonably good knowledge on a number of aspects on the flight, including aircraft specifications, a fairly finegrained population density map, specific flight path, and assumptions on the probability of the flight terminating events. In this section the method is demonstrated using an imaginary, albeit quite realistic flight scenario where all the above parameters are assumed available, either as available specifications or as reasonable estimates. It seems sapient to

Aircraft and Flight Path
The example flight is a transport scenario where a Penguin C aircraft, shown in Fig. 6, is operating a service between  Table 2.

Flight Path
The flight path is a route from Kolding to Aalborg specified with 68 waypoints in latitude and longitude in WGS84 coordinates. The altitude AGL is mostly 100 meters (being the maximum altitude for flights outside urban areas in Denmark), but varies in some places to demonstrate the consequence of flights at higher or lower altitudes. The path is over areas with very low (forest areas) as well as fairly high (city area) population density, also for demonstration purposes. The flight path is shown in Fig. 7, including two excerpts over areas with low and high population density. The altitude of the path is shown in Fig. 8. The flight path is also given in WGS84 coordinates in Table 5 as 'Original path'.

Flight Path Sampling
For the purpose of computing POF during the flight the entire path is sampled at equidistant points between the waypoints under the assumption that the flight path consists of straight lines between WPs. The sample density is chosen in relation to the geographical extend of the probable impact area for each of the events described in Section 2.2. For small impact areas a higher sampling density is chosen such as to capture any change in population density that occurs along the path. The size of the probable impact area for the ballistic event is in the order of 100 m by 100 m, so the flight path sampling density for this event is set to 25 m, and the population density map used is 25 m by 25 m. A more dense sampling  Table 3. Examples of georeferenced impact areas for each of the event types are shown in Fig. 9 for ballistic, uncontrolled glide, and parachute, and for flyaway in Fig. 10. All four examples use the same event point, namely WP 39 over the city of Silkeborg. Notice how the size of the areas vary significantly, and all show signs of the east-north-east wind direction.

Computing Probability of Fatality Along Flight Path
For each flight path sample point the probability of fatal injury p fatal impact can now be computed. First we use Eq. 7 to determine the probability of impacting a person. Based on the descent parameters for each event type, given by Table 4, the lethality estimation described in Section 2.7 is applied to determine the POF given that an event occurs. This probability is then used in Eq. 1 along with the event probability p event to give the unconditional POF p fatality . This probability is valid for that brief period of time where the aircraft is in the vicinity of the flight path sample point; for computational purposes we will assume that this probability is valid for the flight from one flight sample point to the next. Table 4 lists the time distance between path samples for each event type assuming cruise speed. The entire time-varying POF for each event type is shown in Fig. 11. It is immediately obvious that there are distinct differences as well as similarities between the four graphs. Each graph is discussed in the following four subsections, highlighting similarities as well as differences.

Ballistic Descent
The ballistic descent is targeting a relatively small impact area with a comparatively high probability, as is exemplified in Fig. 9. At the same time the resolution of the population density map is small, and because the flight is mostly over a thinly populated area many of the (small) cells in this map have zero density. Consequently, for many of the sampled event points along the flight path the impact area for the ballistic descent only covers zero density cells, resulting in zero POF. However, since the y axis is in logarithmic scale this is not evident from the graph. Whenever the aircraft passes over a farm or cluster of houses, the probability increases briefly, which causes 'spikes' in the POF graph. The transit over Silkeborg at WP 37 through 41 gives a more prominent increase in POF.

Uncontrolled Glide
The uncontrolled glide impact area is somewhat larger (roughly one to two order of magnitudes depending on flight altitude) than that of the ballistic descent, and the population density map used is comparatively more coarse, and as a consequences the POF graph is less spiky. Since the flight path is over mainland Denmark the population density may be low, but virtually never zero for the impact PDF covering tens of thousands of square meters. This results in the POF graph always being non-zero, albeit mostly with quite low

Parachute Descent
The parachute descent impacts a tear drop shaped area stretching from close to the event point and in the mean wind direction. This is clearly seen in Fig. 9. This impact area is in the vicinity of the ballistic descent impact area and we also do see some correlation between the POF graph for the two events. The parachute graph is a little less spiky because the impact area is somewhat bigger, and there are a few places where the parachute POF is indeed zero, most prominently between WP 31 and 33, which is over forest area.

Flyaway
The flyaway impact area is quite large, as can be seen in Fig. 10, covering hundreds of square kilometers. Therefore it changes relatively little during the flight. The model for the flyaway includes the possibility of a vertical ascent (and subsequently descent 'close' to the event point), and this noticeably affects the POF when passing over Silkeborg, where a slight increase in the POF is seen in the graph. The path sample density is 1000 m, significantly higher than the three previous graphs, and this shows up as a distinct stair case pattern in the graph. Some of the slightly larger steps in the graph occurs at waypoints, and is caused by the maximum flight distance being updated (the maximum flight distance is reduced a little as a little less fuel is available) only at every waypoint, not at every sample point along the flight path. While this is visually a bit confusing it has negligible effect on the resulting average POF. Note that this graph is in linear scale. At the bottom of this graph is also shown the time position of the 68 waypoints.

Joint Probability of Fatality
For each of the four events the average POF is computed by averaging over all flight path samples relative to the flight time for each sample. These four average values are summed to give the total probability for the entire flight. The figures are shown in Table 6 as 'POF original path'.
It is evident from Fig. 11 that the POF can be reduce somewhat by avoiding the obvious higher population density areas. To quantify the effect of this the flight path has been modified slightly to circumvent all such areas. The WGS84 coordinates are listed in Table 5 as 'Modified path' and the resulting POF are shown in Table 6. Essentially, all the larger 'spikes' in Fig. 11 have been removed, and as expected the three event types with relative small impact areas have significantly reduced POF, while the flyaway event is largely unaffected.
It is important to note that the probabilities given in Table 6 are rather approximative for two reasons; 1) they are based on a series of assumptions with varying degrees of accuracy and obviously the end result is no more accurate than these assumptions (see Section 3.5), and 2) the joint probability should have been conditional in the sense that the probability of an event at any given time is conditional on any other event has not yet occurred. For instance, a flyaway at a given time is conditional on the aircraft not having experienced a ballistic descent prior to this time. However, as all the probabilities are indeed relative small the error caused by summing the individual event POF is negligible compared to the inaccuracy caused by the previously mentioned assumptions.

Interpretation and Validity of Results
The list of assumptions enabling the computation of the probabilities shown in Fig. 11 and in Table 6 is fairly long and holds values that are estimates based on work done by others on other drones, some estimates are based on general knowledge about small unmanned aircraft, and some of the parameters are pure guesswork and remains unproven. It would be valuable to determine how these uncertainties affect the uncertainty of the results. For some parameters this is easy, such as the event probability that affects the POF merely by scaling (as evident from Eq. 1), but other   parameters enter the computations in a highly nonlinear fashion (such as wind direction and speed). Therefore the effect of varying such parameters is not easy to determine. A study of the sensitivity of the individual parameters remains as future work. The substantial uncertainty aside the estimated parameters in Table 3 and 4 have deliberately been chosen slightly conservative to reduce the risk that the derived probabilities are unrealistically low. Coupled with the above considerations on the lack of knowledge of actual values, it is indeed likely that POF estimates are perhaps as much as an order of magnitude wrong. Still, the precision is within the uncertainty range that applies to many similar considerations for manned aviation.

Conclusion
The application to the example Penguin C flight from Kolding to Aalborg demonstrates how it is possible with the proposed methodology to quantify an estimate of the POF for a specific flight. The computations are fully parameterized, and are thus easy to repeat for another flight path, for another aircraft, for other wind parameters, etc. It is important to note that the risk assessment here does not cover all possible flight terminating events, as does it not include midair collision and impacts with ground obstacles. Also, there might be other event types depending on the type of aircraft. However, it is relatively easy to include additional events (which is certainly necessary for other aircraft types, such as rotorcraft). For approval of BVLOS flights this method contributes in a tangible way to assist the authorities in determining the risk associated with a given type of flight operation, and as indicated above the Danish Transport, Construction, and Housing Authority accepts this method as a valid tool to anyone applying for permit to conduct BVLOS operations in Danish airspace. There remains substantial future work to improve and refine the method, as well as including more events, and more types of aircraft. Also, more accuracy on assumptions will be beneficial for the resulting probabilities.
The Penguin C aircraft is picked at random for this work. The author is not affiliated with UAV Factory. The results cannot be used 'as is' for a quantitative correct risk assessment of the aircraft since some of the aircraft parameters most likely are not correct.