#### Abstract

In this paper a ground safety assessment model is introduced based on the probability estimation of possible impact positions when unmanned aerial vehicle (UAV) crashes on the ground. By incorporating the random uncertainties during the descending process, risks associated with UAV’s ground crash are estimated accurately. The number of victims on the ground per flight hour is selected as the indicative index to evaluate the risk levels of the corresponding ground area. We mainly focus on the analysis of uncertainties that usually appear in drag coefficient which would generate a great amount of effects on the travelled horizontal distance from the failure point to the impact point on the ground, which further influences the possible impact positions. The drag force in the air, failure velocity of a UAV, and wind effects in the local area are all considered in the proposed model, as well as ground features, including sheltering effects on the ground, UAV parameter settings, and distribution of local population. Uncertainties in drag force when a UAV descends, UAV’s initial horizontal and vertical speeds at failure point, and local wind patterns are all considered as the indispensable factors in the proposed model. Especially the probability of fatality once hit by the UAV’s debris is explored to make the safety assessment more reliable and valuable. In the end, the actual UAV parameters and official historical weather data are used to estimate the risks in a real operation environment when a failure event happens at a legal flying height. Experimental results are given based on different types of UAVs and random effects in the descent. The results show that all the operations of all kinds of UAVs selected in the validation are so dangerous that the safety of people on the ground cannot be guaranteed, whose value is much bigger than the manned aircraft safety criterion 10^{−7}.

#### 1. Introduction

Unmanned aerial vehicle (UAV) is the fastest growing sector of the aviation industry, which has become the focus of many commercial and civilian applications including delivery of goods, search and rescue, disaster relief, infrastructure monitoring, precision agriculture, public safety, and weather monitoring [1–3].

Especially because of the UAV’s strong line-of-sight connection links, fast and flexible deployment and mobility, they have been widely used in various civil Internet of Things (IoT) applications. For example, in a smart city environment that consists of a set of targets, UAVs can be applied in multitarget tracking and sensing for the data collection from critical areas to report them to a central entity with high efficiency [4].

However, there are a lot of risks and hazards that are associated with a large amount of operations of UAVs in the airspace. Because of the uncertainties in the airspace and limited flying resources of UAVs themselves, a failure event often happens for UAVs at some specified altitude. It will introduce a lot of safety hazards to the people and properties on the ground, which further results in the imposition of a significant amount of operational restrictions on them [5–7]. It has been recognized and accepted widely by the Federal Aviation Administration (FAA) and European Aviation Safety Authority (EASA) that the risk-based approach should be adopted to the development of a regulatory framework for UAV’s operation [8, 9]. In this way how to assess and estimate the ground risks and hazards that are caused by a UAV’s failure in the airspace and crash on the ground becomes a key challenge and hot research field.

Until recently many methods and algorithms have been proposed to realize the safety assessment on the ground and estimation of ground risks caused by a UAV’s crash. Two UAV failure models were proposed by Burke et al. [10] to determine the nature of the risk imposed to people on the ground directly. Similarly, Barr et al. [11] explored multidependent failures that may result in ground impacts. These models however only assume a single system-level failure and associated failure mode to the ground. Hayhurst et al. [12] provide examples of functional and structural decompositions of unmanned aircraft system (UAS) that can be used to develop a failure and impact model. A similar approach is applied in the space launch industry [13, 14], which also shares the problem of little data and high uncertainties. Many other ground impact models [15, 16] were proposed, each of which used a geometric approach to evaluate the impact positions. This in turn requires some input parameters in the form of position, velocity, attitude, and so on [17]. A number of the existing models [18–20] do not take the size of the lethal area into consideration, assuming the impact occurs at a finite point in space and time. This assumption limits the risk evaluation to that of individual risk. Further assumptions must be made in order to evaluate measures of group or collective risks. Other models [21–23] provided calculations of the lethal area, which differs from the impact area by including overlapping dimensions of different UAVs. But there are no existing models accounting for impact locations given the activation of a contingency mechanism.

In this paper a ground safety assessment model is introduced based on the probability estimation of possible impact positions on the ground. By incorporating the random uncertainties during the UAV’s descending process, risks associated with UAV’s ground crash are estimated accurately. We mainly focus on the analysis of uncertainties that usually appear in drag coefficient which would generate a great amount of effects on the travelled horizontal distance from the failure event point to the impact point on the ground. The drag force in the air, failure velocity of a UAV, and wind effects in the local area are all considered, as well as ground features, including sheltering effects, UAV settings, and distribution of local population. Uncertainties in drag force, UAV’s initial horizontal and vertical speeds at failure point, and local wind pattern are all included as the indispensable factors in the proposed model. The probability of fatality that people once hit by the UAV’s debris is also explored. In the end, the actual UAV parameters, scenario settings, and official historical weather data are used to estimate the risks in an actual operation environment at a legal flying height. Experimental results are given and analyzed based on different types of UAVs and random effects during the descent.

The remainder of this paper is organized as follows. Problem description of the ground safety assessment and risk estimation is given in Section 2. Section 3 describes the modeling and realization for assessment and estimation in the proposed algorithm in detail. Simulations and verifications in actual operation scenarios are listed and analyzed in Section 4, especially including the analysis of possible uncertainties in drag coefficient during the UAV’s descent. Finally, the paper is concluded in Section 5.

#### 2. Problem Description

The management of risks on the ground that are imposed by the UAV’s crash is a significant foundation and traceability for the aviation regulations. The risk-based approach becomes an effective way for ground risk estimations when the UAVs are operating in the low altitude airspace. Once a UAV fails in the air, how to accurately estimate the positions where it may crash plays a key role in the assessments and estimations of ground impact hazard. The problem description is given in Figure 1.

During the UAV’s descent when a failure event happens, two kinds of primary forces would act on the body of a UAV, namely gravity and drag force. The gravity is a constant once the mass of a UAV is determined at the moment when a failure event happens. Because in this paper an actual scenario is adopted, many random uncertainties should be taken into account in the drag force, such as the randomness in initial horizontal and vertical speed of UAV at the failure point and even the wind effects in local areas. A very standard model can be used to describe the drag force with different velocities and shapes that is realistic for a crippled descending UAV, which is correspondingly defined as a second order model that is proportional to the square of the operation velocity. In order to estimate the possible impact positions, it is crucial to calculate the travelled horizontal distance from the failure event point to the impact point on the ground first. By the random effects of drag force and characteristics of UAVs as well as wind effects, the covered area on the ground where the descending UAV may crash becomes a geographical area instead of a line segment. How to accurately calculate the probabilities of the possible impact positions on the ground determines the accuracy of ground risk estimation.

After that the features of the ground should be then analyzed in detail, which could determine the final safety assessment and risk estimation. The affected area by the UAV’s debris is the foundation for the risk calculation. That is, because the risk level is highly related with the operations of different types of UAVs in the airspace, different UAVs will lead to different levels of risks. The UAV parameters that cannot be neglected include wingspan, length, operation velocity, and maximum takeoff mass (MTOM).

Another feature is the probability of fatal injuries to people. It is known that human body is able to sustain a certain level of external forces or injuries caused by some kinds of impacts. Not all the impacts will introduce fatality. So how to precisely evaluate this fatality level would generate a great amount of effects on the results of the ground risk estimation. The density of people and their distributions on the ground are both necessary when estimating the risks in a specified area. More dense distribution means more people would be affected at the same time. How to define this parameter would influence the safety level directly.

The energy possessed by a UAV on the ground impact point cannot be neglected either, which is a direct factor that will damage people and properties on the ground. Generally, the kinetic energy would be used which can be determined by the current mass and impact velocity. The obstacles on the ground may provide extra sheltering effects to people when a UAV hit happens. Many types of obstacles could keep people away from being attacked by the UAV’s crash, including reinforced concrete buildings, trees, and sparse trees. The worst situation is when no sheltering effects exist because of bare surface. It makes people get exposed to the crash one hundred percent.

The casualties per flight hour caused by UAV’s operation is selected as the index to represent the risk levels on the ground. The number of victims on the ground per flight hour is calculated under different scenarios. For manned aircrafts, the widely accepted value given by FAA is 10^{−7}, which can be listed as a reference for UAVs’ operations.

#### 3. Ground Safety Assessment and Risk Estimation

In this part based on the analysis in the previous section, safety assessment of ground impact caused by the crashed UAV is realized and the corresponding values of ground risks that are imposed on the people can be then calculated accurately. The mentioned models and algorithms proposed in this paper are analyzed in detail.

##### 3.1. Calculations of Impact Positions on the Ground

As mentioned above, it is so significant to deduce the travelled horizontal distance from the failure even point to the impact point on the ground. Once the travelled horizontal distance is obtained, the possible impact positions can be determined combining with the flight direction. In this part, an estimation model of the possible impact positions on the ground is introduced [24]. The desired horizontal distance can be divided into three portions, which is given in Figure 2. The first is the horizontal distance from failure point to the highest altitude (defined as top point), which can be denoted as *x*_{1}. The second is from the highest altitude to the point where the vertical speed is equal to the horizontal speed, namely the crossing point of vertical speed and horizontal speed in Figure 2, which is denoted as *x*_{2}. The third is from the crossing point to the impact point on the ground that can be denoted as *x*_{3}. For simplicity in this part the last two portions will be analyzed together as the horizontal distance from the top point to the impact point, namely (*x*_{2} + *x*_{3}).

###### 3.1.1. Horizontal Distance from Failure Point to Top Point

When a failure event happens, it is supposed that all the thrust is lost, which makes sense in the actual UAV operations. In this way the velocity of a UAV at the failure point could be divided into horizontal and vertical speeds respectively. The former one will exist and its value is set as positive all the time during the descent. But for the latter one, the downward is set as positive. In this situation the time (*t*_{top}) at the highest altitude shown in Figure 2 since the failure event happens can be obtained using following [24]:

In (1), is the initial vertical speed at the failure event point, which is set as a random value in the proposed model that is subject to a stochastic process. Due to the uncertainties in airspeed measurement and the effects of the onboard controller attempting to maintain not just speed but also altitude and heading, the actual velocity may easily vary somewhat from the expected value. Γ is denoted as which is the reciprocal of *γ*. *c* is a variable that comprehensively captures the air density (*ρ*), drag area (*A*), and drag coefficient (*C*_{D}), which in turn can be denoted as *c* = 0.5·*ρ*·*A*·*C*_{D}. Drag coefficient here expressed as *C*_{D} is also subject to a stochastic process to better incorporate the uncertainties usually in the drag force during the UAV’s descending process.

The horizontal distance (*x*_{1}) from failure point to the highest altitude can then be determined by (2). In this equation below is the initial horizontal speed of UAV which is also a stochastic variable that is similar to the initial vertical speed. It is used to capture the randomness in the proposed model. *m* is the mass of UAV when the failure event happens. *c* is as the same as above and the value can be obtained using *c* = 0.5·*ρ*·*A*·*C*_{D}.

The altitude gained from flight level to the highest altitude is denoted as *y*_{top} in Figure 2, which can be calculated using the following:

###### 3.1.2. Horizontal Distance from Top Point to Impact Point

Furthermore, the time (*t*_{drop}) given in Figure 2 from the highest altitude to ground impact can be obtained using (4). In (4) *y* is the flying height of a UAV when a failure event happens. And *H*_{d} can be determined by *H*_{d} = arctanh() and so as *G*_{d} = lncosh*H*_{d} correspondingly.

In this way combining (1) with (4), the total time (*t*_{im}) for a UAV from the failure event to the ground impact can be calculated with (5) as follows:

The horizontal speed () of the UAV at the highest altitude is then expressed with (6). In (6), is the initial horizontal speed of UAV which has been defined as a stochastic variable previously. Another important parameter is the moment when the horizontal speed is equal to the vertical speed, which is denoted as *t*_{c} in Figure 2. It is also called the crossing point of vertical speed and horizontal speed. Based on the analysis above, the required time *t*_{c} can be obtained by solving (7):

The horizontal distance (*x*_{2}) travelled from the highest altitude to the point where the vertical speed is equal to the horizontal speed can be then determined by

In this situation, the horizontal and vertical speeds at the time of crossing point can be calculated using (9) and (10) respectively. Theoretically the values of the two speeds should be the same, namely . But there is a fact that the two are not one hundred percent the same because the crossing time *t*_{c} is approximated as given in (7):

In order to calculate the distance (*x*_{3}) from the crossing point of vertical speed and horizontal speed to the impact position on the ground, (11) can be finally obtained:

###### 3.1.3. Total Distance from Failure Point to Impact Point

Based on the derivation of (2), (8), and (11), the total travelled horizontal distance as described in Figure 2 from a failure event point projected on the ground to the impact point is then denoted as

In the end the horizontal and vertical speeds ( and respectively) at the impact point on the ground can be calculated using (13) and (14), which can then be used to calculate the kinetic energy of crashed UAVs at the impact point.

###### 3.1.4. Wind Effects on Impact Positions on the Ground

The wind pattern of the local areas in the UAV’s descending process will generate a great amount of effects which should also be considered. To simplify the model, it is reasonably assumed that the air volume only moves at a constant horizontal speed with a constant direction, which are both independent of the altitude. Of course, there are uncertainties in the wind effects. In the proposed model, both of the wind’s moving speed and direction are set as stochastic processes with some specified distributions.

To make the results more reliable, historical weather data from local official government website can be used to define the randomness in the wind effects based on different application scenarios. Taking the hypothesis above into account, only the horizontal movement of the descending UAV is affected by the wind, which is further determined by the wind’s moving speed and direction, and the total drop time from the failure point of a UAV to the ground impact point.

In this way because of the unavoidable uncertainties existing in the local wind, the impact point will be moved away from the line coincident with the UAV’s flight direction, which results in the affected area on the ground becoming a geographical area instead of a line segment. Finally, the impact point relative to the projection of the UAV’s failure point on the ground can be obtained with (15), which inevitably becomes a vector **P**(*y*).

In (15), *y* is the altitude when a failure event happens which should follow the local regulations of UAV operations, *θ* is the flight direction of a UAV, and *x*(*y*) is the total travelled horizontal distance which is obtained using (12). and *λ* are the wind speed and wind direction respectively. These two parameters can both be set as random variables or can be determined by the local years’ historical weather data. *t*(*y*) is the total time that costs from the failure event to the impact on the ground, which can be obtained by (5) in the proposed model above.

##### 3.2. Derivation of Ground Risk Caused by UAV’s Crash

As described in Section 2, the features on the ground could offer protections for human beings to avoid being attacked by a crashed UAV or reduce the damages imposed on the ground people. To some extent, ground features and distributions of the population could affect the risk estimations. As a result, many factors should be extracted and modelled precisely in the risk assessment algorithm. Casualty areas and fatalities of affected people on the ground are both considered in the proposed algorithm.

###### 3.2.1. Casualty Area on the Ground

The first is the casualty area that could be affected by the UAV’s debris. Normally two types of impacts happen when a UAV crashes on the ground, namely vertical impact and horizontal impact. For the former one, the casualty area is seen as a circle whose radius is the sum of the radius of a human being and the radius of a circle with area equal to the largest cross-sectional part of the UAV piece [25]. For the latter one, the traversed distance since the UAV contacts the ground should be calculated, which is highly related with horizontal and vertical speeds at the ground impact point that are already obtained by (13) and (14). The dimensions of failed UAVs should also be taken into account. In this way the covered area by a UAV’s debris can be calculated with (16). In this equation *r*_{P} is the average radius of a human body and *h*_{P} is the average height of standing human bodies on the ground. These two parameters are determined based on different operation scenarios. *R*_{UAV} is the maximum radius of a UAV dimension.

###### 3.2.2. Fatalities of Affected People on the Ground

It is known that not all the impacts on human bodies are fatal because our bodies could hold a certain level of external force or energy caused by some impacts. In order to estimate the probability of fatality by a UAV’s crash, a comprehensive model [26] is used to define this vital index, which is given in (17). *P*_{f} is the corresponding probability of fatality for ground people and *P*_{S} describes the protection effects of shelters on the ground. The parameter *a* gives the impact energy required for a fatal probability of 50% for ground people when *P*_{s} = 6. *β* is the impact energy threshold that makes the crash fatal with the value of 34 J [27]. *E*_{C} is the kinetic energy that a crashed UAV generates at the impact point. *E*_{C} can be calculated by using (13) and (14), combining with the maximum takeoff mass (MTOM). *k* is the correction factor to adjust the intensity of the fatality, whose value is always between zero and one.

It is definitely that the shelters on the ground could protect people from being injured by the UAV’s attack. Different shelters could provide different covers for the ground population in the actual situations. There is no doubt that bare ground offers no protections which leads to the highest risk. In the previous research, four main kinds of sheltering effects are considered [28], namely the reinforced concrete buildings, trees, sparse trees, and areas without obstacles. The sheltering factor of each mentioned type is an absolute real number as given in Table 1.

###### 3.2.3. Estimation of Ground Risk Caused by Crashed UAVs

In this paper the number of victims per UAV flight hour is selected as a standard criterion to evaluate the risk levels on the ground, which also defines the value of ground risk. So far based on all the models described above, the total number of victims on the ground per UAV flight hour which is caused by the possible UAVs’ crash accidents can be finally determined, which is shown in

In (18), *N* is the number of victims on the ground per flight hour of a UAV’s operation. *P*_{A} is the probability of the possible impact positions on the ground. It stands for the ground area as described in (15) in which a UAV may crash when a failure event happens in the air. *A*_{C} is the covered area by the UAV’s debris as given in (16). *D*_{P} is the density of ground population where a UAV crashes, which can be derived accurately from the local official website. *P*_{f} is the probability of fatality caused by the UAV’s ground impact given in (17).

###### 3.2.4. Analysis of Computational Complexity

It is known that the computational complexity is highly concerned with the computable problems. It means that the problems are considered as computable if they can be solved in principle by a computing machine, which takes as input any precise mathematical statement and could decide whether the statement is true or false after executing a finite number of steps [30]. It is also concerned with the computational resources that are required to solve the corresponding computational problems.

In the calculation of horizontal distance from the failure point to the top point which is also one of the difficult and complicated processes in the proposed framework, once the inputs are given in advance the corresponding outputs can be obtained directly by using (1)–(3) without executing an infinite number of steps. Next as in the calculation of horizontal distance from the top point to the impact point, (7) is a little more complicated compared with other equations used in the model. It is obvious that even though (7) belongs to a quadratic equation, this core parameter can also be solved using existing mathematical formula with very little and limited computational cost. It means the time *t*_{c} can be calculated with high efficiency and low computational complexity. Once the time *t*_{c} is obtained, the other parameters will be calculated step by step easily. When considering the wind effects on the UAV’s impact positions on the ground, a matrix given by (15) is solved once the total distance and time cost from the failure point to the impact point are given by (12) and (5) respectively. Other models used in the derivation of ground risk caused by UAV’s crash are all the first order equations including the determination of casualty areas and fatalities of affected people on the ground, which can all be solved immediately once the inputs are given with one-time calculation. The output is only determined by one input in each equation which further proves that less implementation costs are needed in the proposed framework. All the models used in the algorithm are simple and easy to solve with low complexities and little time cost, which are quite suitable in the risk assessment for UAV operations.

#### 4. Experiments and Case Study

In this work different experiments under different settings are done with actual scenarios to realize the assessment of ground risk and estimation of the hazard levels when a UAV crash happens. The modeling and algorithm are coded in MATLAB and ran on a server with a 2.8 GHz CPU and 16.0 GB of RAM.

##### 4.1. Simulation Parameters

An actual environment is used to make the assessment results be more valuable. Changqing Campus of Shandong Jiaotong University is selected as the UAV’s operation environment [29], as well as the second divisions on the ground to improve the accuracy. The maximum flying height of all UAVs is set as 120 m, which satisfies with the current Chinese UAV regulation that was published in January, 2018. The percentage of different types of ground areas and central angles are the same as in [29], as well as the distribution of population in different divided zones on the ground.

Table 2 shows the percentage of different types of areas and central angles in detail. From the table, we can find that there are no concrete buildings in Zone 2, Zone 4, and Zone 6 because the percentages are all zeros. The type of trees is mainly distributed in Zone 2, Zone 3, and Zone 5. However, there are few sparse trees in Zone 3 and Zone 5. Most areas of Zone 1 are bare that nothing exists in it. Based on the data of central angle, the area of Zone 3 is the biggest with 88.7° and Zone 5 is similar to the central angle 87.0°. Zone 2 and Zone 4 are the two smallest areas of the six. Zone 1 and Zone 6 are in the middle position.

In the experiments, both of fixed-wing and rotary-wing UAVs are taken into account, including all four kinds of UAVs from three different manufacturers. The detailed UAV data can be found in Table 3 below. Wingspan, length, and maximum takeoff mass (MTOM) are used.

From the data we can see that the wingspans of the fixed wing UAVs are much bigger compared with those of rotary wing UAVs, all of which are wider than 1000 mm. But the differences of aircraft length are not so obvious which are all less than 900 mm. The maximum takeoff mass (MTOM) of Zenith ATX8 is much heavier than the other three with 9.65 kg. Firebird can only load the fewest with only 1.2 kg that becomes the lightest one. Based on the models proposed previously, different UAV parameter settings could generate a lot of different effects on the people on the ground.

In the following, the effects of uncertainties existing in the drag coefficient on the horizontal travelled distance are analyzed first. In order to describe the randomness precisely in drag coefficient, normal distribution is assigned to it. And then the safety assessment and estimations of ground risks are given in the end of this part. The number of victims on the ground caused by the UAV’s crash is set as the index to evaluate the risk levels.

##### 4.2. Effects of Uncertainties in Drag Coefficient

In this part, we mainly focus on how the uncertainties in the drag coefficient in the UAV’s descending process would influence the probability dense function (PDF) of the travelled horizontal distance when a UAV fails at some altitude in the air. In the experiments, without loss of generality, it is supposed that the drag coefficient denoted as *C*_{D} which appears in the equation *c* = 0.5·*ρ*·*A*·*C*_{D} is subject to normal distribution *N*(*μ*, *σ*) to reflect the randomness. Here *μ* is the mean value and *σ* is the standard deviation respectively. The former stands for the intensity of the drag coefficient and the latter reflects how the real value deviates from the means. There are two main kinds of drag coefficients that are considered in the experiments. One is the weak randomness with *N*(0.5, 0.1) and the other is a strong randomness with *N*(0.9, 0.3). All the four UAVs operate in the proposed scenario above and the flying heights are set as 30°m, 60°m, 90°m, and 120°m. It is illegal to fly the UAVs above 120°m without authorization in China right now. The results are given in Figures 3 and 4 respectively.

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

In Figure 3 it is known that the drag coefficient is set as a normal distribution which is subject to *N*(0.5, 0.1). From the curves in the figure, we can find that no matter at which altitude the UAVs fail, Zenith ATX8 always travels the longest horizontal distance from the event point to the impact point on the ground with the highest mean values. Then is Typhoon H, whose travelled horizontal distance is shorter than that of Zenith ATX8. The distance that Firebird travels is the shortest compared with the other three. The distance of X8 is a little bigger than that of Firebird.

It proves that the fixed-wing UAVs are more easily affected by the random drag coefficient than rotary-wing UAVs. Also, the probability distribution of fixed-wing UAVs is more concentrated with very small fluctuations. Oppositely the probability distributions of rotary-wing UAVs are more dispersed with bigger standard deviations. Also, from the four subfigures in Figure 3, when the UAVs operate at a higher altitude, they would travel longer distances since they fail. Especially when the flying height is set at 120°m, the travelled horizontal distance is longer than 100°m. It is known that longer travelled horizontal distance means more areas on the ground will be covered by the UAVs’ debris and more people would be affected under the same situations.

For the UAV of Firebird that operates at any height, there is a high probability that the shortest horizontal distance can be smaller than 5°m, which means it will fall from the flying height with a fast speed when a failure event happens. It results in smaller affected area and fewer people attacked on the ground.

What is more, for the four UAVs operating in the low altitude airspace, there are no doubts that Firebird is the safest of the four with the lowest ground risk values. Oppositely Zenith ATX8 is the most dangerous that will introduce the most damage on the ground of the four. Uncertainties in drag coefficient definitely generates a lot of effects on the travelled horizontal distance for all kinds of UAVs. It describes the randomness during the descending process of a UAV. Many factors play different roles in the proposed model. From the results we can find that by incorporating the stochastic process in the descent, the covered area on the ground will also follow some similar distributions. By using kernel density estimation (KDE), the precise distribution can also be obtained further.

When some strong uncertainties are added into the drag coefficient, which is defined previously, the results of the experiments are shown in Figure 4. The parameters of the scenarios in Figure 4 are the same as those in Figure 3 except the distribution of drag coefficient. It also follows a normal distribution which is subject to *N*(0.9, 0.3). Compared with the data in Figure 3, all the curves in Figure 4 are getting flatter especially for Typhoon H and Zenith ATX8, which means that the fluctuations in them are more drastic than those in Figure 3. For Firebird and X8, the differences between the two are getting smaller and smaller. The uncertainties in drag coefficient makes the two kinds of UAVs get closer and closer in the probability dense function. Stronger randomness could generate a lot of effects on the travelled horizontal distance.

Further there are more overlaps of the four curves in Figure 4. It is because more fluctuations can be introduced by stronger uncertainties in drag coefficient. The mean values of the four curves at any flying altitude are smaller than those in Figure 3, which is because the intensity of drag coefficient becomes stronger. Strong drag coefficient increases the drag force that makes all the UAVs crash into the ground with more time, which further slows down the UAVs’ crash. In this way, more time is needed from the failure point to the crash point on the ground. But from the probability distribution described in Figure 4, for Zenith ATX8 especially, there is a probability the travelled horizontal distance could achieve as far as 120°m and even longer. All in all, the randomness existing in the travelled horizontal distance becomes stronger because of more randomness in drag coefficient settings.

##### 4.3. Safety Assessment and Estimations of Ground Risks

In this part, ground safety is assessed and the risk is estimated with random uncertainties. The radius and the height of human beings on the ground are set as 0.25 m and 1.8 m respectively. Local official weather data of Changqing district in Jinan is used which can be captured from the government website.

In the experiment it is supposed that all the UAVs will fly with different directions, and their flight angles are 0, *π*/2, *π*, and 3*π*/2. The wind direction and speed are both random variables satisfying with the normal distributions. They are subject to *N*(5*π*/4, *π*/8) and *N*(3.4, 0.5) respectively. The experimental results are given in Figure 5.

**(a)**

**(b)**

**(c)**

**(d)**

As shown in Figures 5(a) and 5(b), all the six zones on the ground are affected by the UAVs’ crash. In Figure 5(a), Zone 1 is attacked by all the four UAVs and the biggest number of victims is over 1.25 × 10^{−2} that is caused by Zenith ATX8. X8 could generate the strongest damage in Zone 6, whose values are bigger than 2.50 × 10^{−2} and 2.10 × 10^{−2} in Figures 5(a) and 5(b) respectively. They are all above the accepted safety baseline 10^{−7} that is given to manned aircrafts by FAA. The operations of UAVs in the corresponding airspace introduce a great amount of risks to the ground in the situations mentioned above.

In Figure 5(c), Zone 5 is also affected by all the four kinds of UAVs with flight direction *π*. The risks of the four UAVs in Zone 5 are all above 2.0 × 10^{−3}. But Zone 2, Zone 3, and Zone 4 are safe enough that no UAVs crash in these areas, which means no people or properties are damaged by them. The most dangerous is still Zone 6 which is caused by X8 whose value is close to 2.4 × 10^{−2}. In Figure 5(d), four UAVs all crash into Zone 1. The corresponding risks are much bigger than 1.5 × 10^{−4} which is of a high level compared with 10^{−7}. Still no UAVs could affect Zone 2, Zone 3, and Zone 4. Zone 5 is attacked by Firebird, X8, and Typhoon H. The highest risk in Zone 6 is 2.7 × 10^{−2}.

All in all, from Figure 5 we can find that fixed-wing UAVs could introduce more risks on the ground compared with rotary-wing UAVs. The risks generated by any kinds of UAVs are all higher than the safety baseline 10^{−7} compared with the requirement given to manned aircrafts, which further proves the higher requirements for UAV operations are quite needed to ensure the safety of people and properties on the ground.

#### 5. Conclusions

We have presented a method to assess the ground risks caused by the crashed UAVs, in which the possible impact positions are estimated first by deducing the PDF of travelled horizontal distance from the UAV’s failure point to the ground impact point. Then by considering the uncertainties that usually appear in the UAV’s descending process and local wind patterns, the obtained risk values would be more practical and valuable, especially when the real UAV parameters are used as well as the ground features. Finally, by means of incorporating an actual flying environment, the proposed models and algorithms are tested and validated. From the case study we can see that the uncertainties in random drag coefficient would generate more effects on fixed-wing UAVs (Firebird and X8) than rotary-wing UAVs (Typhoon H and Zenith ATX8). Further the probability distribution of travelled horizontal distance given by fixed-wing UAVs (Firebird and X8) is more concentrated with very small fluctuations, compared with that given by rotary-wing UAVs (Typhoon H and Zenith ATX8) whose probability distributions are more dispersed with bigger standard deviations. What is more, because strong drag coefficient increases the drag force acting on the UAVs themselves, it would make all the UAVs crash into the ground with more time which further slows down the UAVs’ crash. Further, it could be found that rotary-wing UAVs (Typhoon H and Zenith ATX8) introduce less risk comparing with fixed-wing UAVs (Firebird and X8). But the experimental results show that all the four UAVs operating in the specified airspace are so dangerous that the safety of people on the ground cannot be guaranteed whose value is much bigger than the manned aircraft safety criterion 10^{−7}.

#### Symbols

x_{1}: | Horizontal distance from failure point to the highest altitude |

x_{2}: | Horizontal distance from top point to crossing point |

x_{3}: | Horizontal distance from crossing point to impact point |

t_{top}: | Time at the highest altitude since failure event happens |

: | Initial vertical speed at failure point |

ρ: | Air density |

A: | Drag area |

C_{D}: | Drag coefficient |

: | Initial horizontal speed at failure point |

m: | Mass of UAV when a failure event happens |

y_{top}: | Altitude gained from flight level to the highest altitude |

t_{drop}: | Time from top point to impact point |

y: | Flying height of a UAV when a failure event happens |

t_{im}: | Total time from failure event to ground impact |

: | Horizontal speed of UAV at the highest altitude |

t_{c}: | Time when horizontal speed is equal to vertical speed |

x: | Horizontal distance from failure point to impact point |

: | Horizontal speed at impact point |

: | Vertical speed at impact point |

P(y): | Position vector of impact point on the ground |

θ: | Flight direction of a UAV |

x(y): | Total travelled horizontal distance determined by y |

: | Wind speed |

λ: | Wind direction |

t(y): | Total time from failure event to impact determined by y |

r_{P}: | Average radius of a human body |

h_{P}: | Average height of a standing human body |

R_{UAV}: | Maximum radius of a UAV dimension |

A_{C}: | Covered area by the UAV’s debris |

P_{f}: | Probability of fatality for ground people |

P_{S}: | Protection effects of shelters on the ground |

E_{C}: | Kinetic energy crashed UAV generated at impact point |

N: | Number of victims on the ground per UAV flight hour |

P_{A}: | Probability of possible impact positions on the ground |

D_{P}: | Density of ground population where a UAV crashes. |

#### Data Availability

The datasets generated and/or analyzed during the current study are available from the first author on reasonable request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by National Natural Science Foundation of China (grant number 71731001), Shandong Provincial Natural Science Foundation (grant number ZR2020MF151), ZheJiang Key Laboratory of General Aviation Operation Technology (General Aviation Research Institute of Zhejiang JianDe) (grant number JDGA2020-6), National Natural Science Foundation of China (grant numbers U1933130, U1533119), and Research Project of Chinese Academy of Sciences (grant number ZDRW-KT-2020-21-2).