Safety Assessment and Risk Estimation for Unmanned Aerial Vehicles Operating in National Airspace System

This paper proposes an effective approach for modelling and assessing the risks associated with unmanned aerial vehicles (UAVs) integrated into national airspace system (NAS). Two critical hazards with UAV operations are considered and analyzed, which are ground impacts and midair collisions. Threats to fatalities that result from the two hazards are the focus in the proposed method. In order to realize ground impact assessment, a multifactor risk model is designed by calculating system reliability required to meet a target level of safety for different UAV categories. Both fixed-wing and rotary-wing UAVs are taken into account under a real scenario that is further partitioned into different zones to make the evaluation more precise. Official territory and population data of the operation scenario are incorporated,as well as UAV self-properties. Casualty area of impacting debris can be obtained as well as the probability of fatal injuries on the ground. Sheltering factors are not neglected and defined as four types based on the real scenario. Whenmidair collision fatality risk is estimated, a model of aircraft collisions based on the density of civil flight in different are constructedtocalculateexpectedfrequencyof fatalities for eachprovincecorrespondingly.Truthfuldata with differentnumbers of UAVs is incorporated in the model with the expected number of fatalities after a collision is included. Experimental simulations are made to evaluate the ground impacts and midair collisions when UAVs operate in the NAS.


Introduction
Unmanned aerial vehicles (UAVs) are a kind of aircraft without pilots onboard but can be remotely controlled or can fly autonomously based on preprogrammed flight plans [1,2], which is a viable and operational technology in the future [3,4]. UAVs are one of the fastest growing aviation sectors. Due to their capabilities to work remotely and under extremely hazardous environments [5,6], the UAVs could accomplish a wide variety of missions ranging from military reconnaissance to environmental monitoring [7,8].
The integration of UAVs into the NAS presents a number of challenges currently addressed by Federal Aviation Administration (FAA) in the context of the operational safety implications [9]. Like all technologies, there are two primary risks associated with their applications [10]. One is the impact with people or structures situated on the ground when UAVs crashed. The other is the collisions with other aircraft in the air.
Until now UAVs have to be operated in a segregated airspace based on FAA regulation. That is because there are inherent safety concerns with UAVs due to the lack of onboard human pilots [11]. Since a final rulemaking has not yet been completed, FAA only approves UAVs' access to the national airspace on a case-by-case basis [12]. FAA provides this approval through three different means, public or civil certificates of waiver or authorization (COA), Section 333 exemptions, and special airworthiness certificates in the experimental category and the restricted category [13]. One of the important steps in obtaining the approval is a proof that the UAS operations can be conducted at an acceptable level of safety [14,15]. Risk assessments and estimations of UAVs operating in the NAS become the basement and the first key step. Because it will fundamentally transform existing aviation patterns and its public perceptions. In the past few years, many methods have been proposed to realize the assessments and estimations. Anno [16] investigated midair collision risk using a random collision theory and compared the results with historic collision data from 1969 to 1978. McGeer [17,18] did hazard estimation research on the Aerosonde UAS. In the existing research, uniform traffic densities were incorporated. Two different constant densities were used for UAVs themselves and the background traffic. A study on both midair collisions and ground impacts based on the collision model of gas molecules was proposed by Lum and Waggoner [19]. Also in Lum's another paper [20], actual UAVs' trajectories for a ground impact analysis and a realistic distribution of average glide angle were used to refine the expected values of ground fatalities. A comprehensive system-wide study presented by Weibel and Hansman [21] used a ratio of the volume swept by the background aircraft to the total airspace volume. In his other related work [22], conditional probabilities were included to develop a separation standard model using the uncorrelated encounter model [23]. Hak-Tae Lee [11] proposed a distributed traffic model with actual traffic data, which was collected over a one-year period to enable a probabilistic approach for risk assessments. Sheridan [24] proposed a model to estimate the relative collision probability between any two aircraft at the closest point of their approaches based on Gaussian density functions.
The main contribution of the work is that an effective approach for modelling and assessing the risks associated with UAVs integrated into NAS is proposed. Two critical hazards of UAV operations are taken into account, namely, ground impacts and midair collisions. The corresponding hazard analysis is conducted and we focus on threats to fatalities generated by the two hazards. A ground impact assessment model is proposed considering system reliability required to meet a target level of safety for different UAVs. Both fixed-wing and rotary-wing UAVs are taken into account under a real scenario, which is further partitioned into different zones to make the evaluation more precise. In the model, territory and population data, casualty areas, and sheltering factors are all indispensable. To estimate the midair collision fatality risk, a model of aircraft collisions based on density of civil flight in different regions over China is proposed. A relative collision area and operating speed between UAVs and manned aircraft are constructed to obtain expected frequency of fatalities for each province using official government data with different numbers of UAVs. Experimental simulations are made to evaluate the ground impacts and midair collisions when UAVs operate in the NAS. The models in this paper provide a generic framework that can be used to structure the development of safety cases for any UAV operation.
The remainder of this paper is organized as follows. Ground impact assessment, including problem description, model establishment, and result analysis, is provided in Section 2. Section 3 describes the midair collision risk estimation models and the final simulation verifications. The paper is concluded in Section 4.

Ground Impact Assessment
A model of ground impact assessment for UAVs was proposed to investigate the influence of different factors on the equivalent level of safety (ELOS) in terms of ground fatalities per hour of UAV system operations. The model incorporated total system reliability, UAV size and kinetic energy, and population characteristics and probabilities of fatality in different vicinity of operation.

Problem Description.
Once a failure occurs on the UAVs, there will be uncontrolled ground impacts. It is assumed that there will be a number of fatalities on the ground once ground impact events take place, which are further defined as hazardous events [25] by FAA criteria [26]. There are many factors that should be taken into account, which is shown in Figure 1. The first is the properties of UAVs themselves, including the probability of failure and the casualty area on the ground. The other are ground features where UAVs crashed. Population density and probability of fatality are also included.
Probability of failure means whether a hazardous event will happen. It could be used to evaluate the frequencies that a UAV malfunctions. When the casualty area is calculated, dimension of UAVs, maximum flying height and velocity, and maximum takeoff mass (MTOM) are incorporated based on the UAV themselves. Casualty area is also strongly affected by glide angle for different UAVs. Combining with the population density, the average height and radius of a human body on the ground could also influence the final assessments. The probability of fatality is highly related to the sheltering factors considering different types of crash areas. Higher sheltering factors will lead to lower probability of fatality.

Model Analysis
(1) Ground Impact Equivalent Level of Safety (ELOS). The safety level of UAVs that operate in national airspace system (NAS) is usually measured by fatality rate that is given as ground fatalities per hour of operations of the UAVs. For manned aircraft the fatality rates based on NTSB data for the period from 1983 to 2006 are in the order of 10 −8 , although it does not include fatalities after emergency landings, ditching, and other situations. The equivalent level of safety (ELOS) is closer to 10 −7 if the latter are included [27]. Usually the ELOS for UAVs used for the ground impact analysis is set an order of magnitude beyond that required for manned aircraft systems. UAVs will not be accepted by the public unless it is safer than manned aircraft. In this paper 10 −8 is set for ground impact assessment of UAVs, which exceeds the current level of safety of manned aircraft operations [21].
(2) Probability of Failure. Probability of failure for UAVs means whether a hazardous event will happen. It is highly related to the safe operating time in the NAS. Longer safe operating time means lower failure probability. Different UAVs will have different failure probabilities. Obviously failure probability is derived from the inverse of the time between two consecutive failures [21]. All general types that could result in ground impact are considered as failures, which includes failures of any components on the UAV system or human errors. Suppose the failure probability of UAV is denoted as and t F is the time between two consecutive failures, the following equation could then be obtained: (3) Casualty Area of UAV Debris. There are two kinds of impacts on the ground when a UAV crashes. One is vertical impact and the other is horizontal impact. We will illustrate the casualty areas for these two kinds of impacts, respectively. The casualty area of a vertically falling piece of debris is a circle whose radius is the sum of the radius of a circle with area equal to the largest cross sectional area of the piece and the radius of a human being [28]. As given in Figure 2, r P is the average radius of human body and R UAV is the maximum radius of UAV dimension. Then the casualty area of a vertically falling UAV debris could be calculated with the following equation: The basic horizontal casualty area is illustrated in Figure 3. Note that d is horizontal distance that UAV debris travels as it falls the height of a person and the impact angle is an angle that the velocity vector makes with the horizontal plane or surface impacted [28]. H P is the average height of standing human bodies. Based on the illustrations in Figure 3, horizontal casualty area can then be obtained by using In (3) d is the horizontal distance that is calculated using (4). In (4), is the approaching angle of UAVs, which is derived from their corresponding features. (4) Probability of Fatal Injuries. It is known that human body is able to sustain a certain level of force or injury. In this way, the presence of a person in an area affected by a ground impact cannot guarantee a fatality. Moreover, obstacles on the ground such as buildings and trees may provide shelters, which will correspondingly increase the chances of survival.
In order to describe the probability of fatality of a person exposed to a UAV crash on the ground more accurately, the human vulnerability and sheltering effects are indispensable in the proposed model. It is evident that both the two aforementioned factors should be taken into account. This model is provided by the following equation [27]: Probability of fatality P f can be calculated by using (5). In (5), P S stands for the sheltering factor, whose value is selected from zero to infinity. Because it is rescaled, so that its average value matches closer the Feinstein average lethality curve. The parameter is the impact energy required for a fatality probability of 50% for human beings when = 6. The parameter is the impact energy threshold required to cause a fatality for a person, which can be considered to be a constant with value 34 J [29]. E C is the kinetic energy of crashed UAVs at the impact point, which is determined by their own masses and velocities. Another k is a correction factor that is used to improve the estimations given for low kinetic energies, especially those close to, or below the threshold limit of 34 J. It can be further obtained by (6). The objective is to maximize the robustness while satisfying some desired mission requirements. Especially in (5), sensitivity and robustness analysis are performed by changing different values of a, k, and , which is associated with the task and environmental complexities on value and risk models demonstrated the robustness and reliability of these models.
(5) Kinetic Energy Estimation at Impact. In this paper two kinds of UAVs are taken into account, one is fixed-wing and the other is rotary-wing. Once the velocity at the crash point is calculated, the kinetic energy can then be obtained. The key problem is how to calculate the velocity accurately. Figure 4 gives the velocity determination for fixed-wing UAVs. h max is the maximum flying height. The glide angle or the approaching angle is set as 45 ∘ [30]. In this way, the velocity at the crashing point can be calculated with [30] Figure 5: Velocity determination for rotary-wing UAVs.
V falling from h MAX is set as equal as the freefall velocity from the maximum flying height [30]. V max op is the maximum operation velocity which is determined by UAV itself.
For rotary UAVs, the impact velocity is obtained with another way as shown in Figure 5.
is the maximum operation velocity and is the freefall velocity from the maximum flying height. In this way, (8) can be formulated.
There is no doubt that can be calculated using the equation V y = • . Combining Figure 5 and (8) the following equation can be used to calculate the approaching angle : When the velocities of all kinds of UAVs at the crashing point are obtained, the corresponding kinetic energy could be then determined as given in the equation below as well as the maximum takeoff mass (MTOM): Once all the factors are taken into account, the ground risk assessment model could be completed using N is the number of on ground victims per flight hour as well as ground fatalities, which should satisfy the safety level given by FAA. Based on the analysis above, its value in this paper is set as 10 −8 . D P is the population density where the UAVs operate. In this way, by using (11), failure probability for UAVs can be determined. Further, the time between two consecutive failures t F could be obtained with (1) for each UAV operating in the NAS.

Model Simulations
The scenario considered in this paper is a real case, which is Changqing Campus, Shandong Jiaotong University, Jinan, China. The flying area is shown in Figure 6, where three circles are drawn. The green circle is the operation area with radius 200 meters, the yellow one is the buffer area with radius 350 meters, and the adjacent area with radius 450 meters is in red color. All the flight operations took place in the operation area, namely, in the green circle.
The Changqing Campus is characterized by different population densities and offers different kinds of shelters Journal of Advanced Transportation 5 Figure 6: Simulation scenario. for people on the ground. In the proposed model, four types of shelters are defined, namely, reinforced concrete buildings, trees, sparse trees, and areas without obstacles. The reinforced concrete buildings could offer high population density but as well as high values of sheltering factor. The others are characterized by low population densities and sparse trees can only be able to offer poor sheltering effects. The sheltering factor is an absolute real number as mentioned before. It is evaluated according to a qualitative estimation of the operative scenario. Different sheltering factors [31] for different types of areas in the proposed model can be found in Table 1. Suppose the total area where UAVs operate is S; the separate areas of reinforced concrete buildings, trees, sparse trees, and without obstacles are denoted as S 1 , S 2 , S 3 , and S 4 , with their own sheltering factor P 1 S , 2 , P 3 S , and P 4 S . We can easily find that = 1 + 2 + 3 + 4 . In this way, overall sheltering factor of the operation area can be calculated using following equation: In order to evaluate the average population density and sheltering factor accurately, the area has been partitioned into six separate flying zones, which is shown in Figure 7. In each zone, the percentage of different types of areas and central angles of the six zones are different based on the real scenario case. The partition in detail can be found in Table 2.
Combining with Figure 7, we can find that the six zones are determined by their corresponding central angles, whose units are degree. From Table 2, there are no concrete  Figure 7: Partition of the flight area. buildings in Zone 2, Zone 4, and Zone 6, which means the sheltering effects are very weak. Correspondingly, the density of population in these flying zones will be low. Another scenario parameter that should be paid attention to is the density of population in each zone. Based on the official data from the school website, the total number of people in Changqing Campus, Shandong Jiaotong University, is 26335. The mathematical statistics is used to find out how many people are there in different zones at a special time during the workday and the distribution of all the population is analyzed. Table 3 gives the distribution of all the people in different operation zones. From the table, we find that most people are in Zone 1 during the daytime, namely, as high as 40%. The number of people in Zone 3 is smaller than that in Zone 1. That is because dormitories and canteens locate in Zone 3. There are outdoor field and school gymnasium in Firebird and X8 are two fixed-wing aircraft with similar length. But X8 is larger than Firebird with wider wingspan, as well as heavier maximum takeoff mass. X8 could fly at 110 km/h, which is much faster than Firebird at 83 km/h. Typhoon H and Zenith ATX8 belong to rotary-wing UAVs, which are much smaller than fixed-wing UAVs. Their wingspans are 457 mm and 600 mm, respectively. The two also have similar length as 520 mm and 600 mm. However, the maximum takeoff mass of Zenith ATX8 is much heavier than that of Typhoon H as 9.65 kg and 1.98 kg, respectively. The flying speed of Zenith ATX8 is 20 m/s, which is faster than that of Typhoon H with 13.5 m/s. Another parameter is the flying height of UAVs in the operation airspace. Considering the performances of the four UAVs in Table 4, it is illegal to fly above 120 m in China based on a new regulation published in January 2018. For the baseline case, the maximum height is set at 120 m in the proposed model, which means all the UAVs will fall from the height of 120 m.
The average height and radius of human beings on the ground that is needed in (3) and (4), namely, H P and R P , are set as 1.8 m and 0.25 m, respectively, during the experiments. The sliding angle for fixed-wing UAVs is to be seen constant as the same as 45 ∘ .

Result Analysis
(1) Fatality and Sheltering Factors in Different Zones. Figure 8 shows the probability of fatality when the four kinds of UAVs fly in six different operation airspace in Changqing Campus, Shandong Jiaotong University. There is no doubt that it is more fatal in Zone 6 and Zone 4 compared with other four zones. That is because in Zone 6 30.28% is bare and 62.03% is covered by sparse tress. The sheltering effect is very low in this zone. In Zone 4, 36.95% is bare and there is a hill in this field. Also 21.10% of this zone are sparse tress. Zone 3 is the safest   because there are a lot of trees and buildings, which make the sheltering factor high. The total percentage of the two types could arrive 76.96%. Zone 1 and Zone 2 are similar as given in Figure 8. From the histograms in the figure, Firebird is the most fatal compared with the other three no matter where it operates. Typhoon H is a little better than Firebird. Zenith ATX8 is the safest of the four. Actually, the differences between Firebird and Typhoon H are not huge at all which are given in Figure 8. Figure 8 also shows the sheltering factors in different flying zones. Zone 3 could provide the best protection for the people on the ground with the highest sheltering factor, which is because 76.96% of this zone are buildings and trees. Zone 6 is the most dangerous with the lowest sheltering factor, which is led by bare ground and a lake. Different simulation results as given in Figure 8 are presented to show the robustness of the proposed model.
(2) Ground Impacts for UAVs in Different Flying Zones. Based on the parameters given in the previous part, four different UAVs operate in the six operation airspace. Then combining with (11), the failure probability for each kind of UAV can be calculated. The results are given in Table 5. Table 5 shows the area proportions, density of population, and sheltering factors of different zones. Zone 3, Zone 5, and Zone 1 are the three biggest areas with the percentage of 24.64%, 24.17%, and 21.67%, which are student dormitories, teaching buildings, and sports field. Correspondingly, the density of population in Zone 1 is the highest and then is    Zone 3 with the values 0.0764 and 0.0504, respectively. The third is the sports field. It means most of the population will be in the classrooms or dormitories, which is a true situation in Chinese universities. Also in Zone 3 most of this area is concrete buildings as high as 27.27%, which makes its sheltering factor as high as 21.428 that is shown in Table 5. It is similar for Zone 5 and Zone 1. The two types of concrete buildings and trees are the main areas in the two zones.
In Table 6, Table 7, Table 8 and Table 9, the kinetic energy of Zenith ATX8 at the impact point is the highest, which could arrive at 13278.4 J. It is because the maximum takeoff mass of Zenith ATX8 is the heavies. X8 flies the fastest of the four and the maximum takeoff mass is 4.2 kg. In this way, its impact energy is the second large but much larger compared with Firebird and Typhoon H. Firebird has the smallest energy at the impact point on the ground because of the smallest maximum takeoff mass as well as slow moving speed. The covered area on the ground of X8 is the biggest, which is mainly determined by the dimensions of UAV itself as analyzed in (3) and Typhoon H influences the smallest area on the ground with only 2.63 m 2 . The other two are in the medium positions. In all rotary-wing UAVs would cover much smaller areas on the ground compared with fixed-wing UAVs.
When the UAV Firebird flies in Changqing Campus, Zone 6 is the most fatal with = 0.37 and Zone 3 is the safest with the smallest = 0.11. The sheltering factors in the two zones are so different that one is the lowest and the other is the highest. The probability of failure in Zone 1 is the lowest, which means the safety requirement on UAVs in this zone is the highest. For Firebird, the safe operation time between two consecutive accidents should arrive 1.72 × 10 7 hours at least. The lowest value is 1.72 × 10 6 hours, which is in Zone 2. The requirements in Zone 3, Zone 5, and Zone 6 are similar but much higher than that in Zone 4. Table 7 shows the simulation results for X8. When X8 flies, Zone 6 is the most fatal with the value 0.35. Oppositely, Zone 3 is the safest only with 0.08. The highest failure probability is in Zone 2 as high as 3.32 × 10 −7 . In this way, the required safe flying time would be 3.01 × 10 6 hours. If X8 flies in Zone 1, its failure probability should be as low as 3.26×10 −8 , which means the time between two consecutive accidents should be no less than 3.07 × 10 7 hours. But Zone 3, Zone 5, and Zone 6 are in the same order of magnitudes, which is lower than Zone 4. Zone 3 is relatively safer compared with the Zone 5 and Zone 6 when X8 operates in Changqing Campus.
The similar results can be found for Typhoon H in all the six zones. Zone 6 is still the most fatal and Zone 3 is safest of the six. The probabilities of failure in Zone 1, Zone 3, Zone 5, and Zone 6 are in the same order of magnitudes. Zone 1 needs Typhoon H to fly safely without any malfunctions up to 3.48×10 6 hours which is the longest time. Zone 2 and Zone 4 have the similar requirements with the order of magnitudes 10 5 , which is much lower than that in the other zones. Table 9 gives the results for Zenith ATX8 in the six different zones. Zone 6 and Zone 4 are the two most dangerous because of the low sheltering factors. Zone 1 and Zone 6 require ATX8 operate a long safe time without any malfunctions. Zone 2 and Zone 4 would tolerate poor performance of ATX8, which need it fly safely up to 3.32 × 10 5 and 6.47 × 10 5 hours, respectively, which are much shorter compared with other zones.
Based on Table 6, Table 7, Table 8 and Table 9, no matter which kind of UAVs operate in the six zones, Zone 6 is always the most fatal to the ground and Zone 3 is the safest with the lowest fatality probability. UAVs will produce fewest injuries on the ground in Zone 3. Every UAV should fly the longest time without any failures in Zone 1, because the failure probability in this zone is the lowest for all kinds of UAVs. Zone 2 and Zone 4 are the two medium zones in Changqing Campus, Shandong Jiaotong University.

Midair Collision Fatality Risk Estimation
In this part an estimation model of midair collisions between UAVs and other civil manned aircraft are proposed. Actual civil flight data of taking off and landing in 2015 over Chinese airspace are incorporated to develop the model. The expected frequency of fatalities in the midair is used as an evaluation standard.

Problem Description.
In this paper the midair collision risk estimation was based on the use of an existing general gas model [32] of aircraft collisions to estimate the expected frequency of fatalities per hour of flight. In the proposed method a model of midair collisions between UAVs and other aircraft was developed. In the model a UAV is equally likely to be located anywhere in the airspace. Additionally, its velocity is assumed to be small compared to the civil aircraft. Once civil aircraft fly within the airspace, a potential collision volume will be extruded. In this way the expected level of safety in terms of potential collisions per hour of UAV operation is then becoming the ratio of volume extruded by threatened civil aircraft per hour to the volume of airspace.
The model incorporates domestic air traffic density data in 2015 over Chinese airspace. UAVs are supposed to be deployed randomly in a specified airspace. Once there are threatened civil aircraft entering the space, the collision volumes of UAVs will be extruded. The number of UAVs and civil aircraft in the airspace determines the occupation rate of the airspace. Relative collision area and speed between UAVs and civil aircraft could generate effects on the final collision frequency.

Model Analysis
(1) Expected Number of Fatalities after a Collision. The number of people exposed to the midair accident as well as the probability of them sustaining fatal injuries depends on the aircraft that are involved in the accident and the passengers they carry. In this way, it is difficult to get a good estimate without a priori knowledge of all air traffic in the area of operations. Based on the NTSB accident data from 1983 to 2006, the value of the expected number of fatalities after a midair collision is closer to 0.58 [27]. Moreover, if the onboard fatalities after a collision with obstacles other than aircraft are ignored, the expected number of fatalities per accident drops to below 0.09. It should be noted that this estimate can be considered conservative because in contrast with the accident data it was derived from, the midair collisions of interest will always involve at least one aircraft that is unoccupied. In this paper, the expected number of fatalities after a collision is set as 0.58, which is widely accepted in the existing research.
(2) Number of UAVs and Density of Civil Aircraft. In the proposed model, for the midair collision problem all the UAVs in the airspace are assumed to be uniformly distributed from sea level to 50,000 ft to simplify the calculation of the midair collision risk, in which elevation is neglected. The number of UAVs in the airspace, N UAV , is a variable. Density of civil aircraft AC reflects the number of aircraft per cubic meter and per hour, which can be obtained using following equation. A total year's data of Chinese domestic civil flights are used.
In (13), N total is the total number of domestic civil flights taking off and landing over Chinese airspace in the whole year of 2015. N day and N hour are the number of the days in a year and number of hours in a day respectively. V as is the volume of the specified airspace, which could be calculated by 50,000 ft multiplying territory acreage.
(3) Collision Areas of UAVs and Manned Aircraft. In the gas model of aircraft collisions [32], each civil aircraft will fly a distance, d AC , through the airspace segment under consideration. There is an area of exposure for threatened civil aircraft A AC as well as for UAVs A UAV , which representing the contact area that in a collision. The collision area is important for the calculation of collision times. For the preliminary analysis in the estimation model, the area of exposure for civil aircraft A AC is estimated as the frontal area of a Boeing 747, approximately 180 m 2 . The area of exposure for UAVs A UAV varies based on different kinds of UAVs. Collison area could be calculated as follows: (4) Relative Speed between UAVs and Civil Flights. Another parameter is the relative speed between UAVs and civil flights. It is used to calculate a potential collision volume combining with collision areas given in (14). Here we suppose that V UAV and V AC denote the speeds of UAV and civil aircraft respectively. The relative speed can be obtained as in In (15) is the relative direction angle formed by UAVs and civil aircraft, where = 2 / in the specified airspace. Combining (14) and (15) the potential collision volume can be further obtained by (16), in which T is the relative time unit that UAVs and civil aircraft operate.
Combining with all the analysis above, the expected number of collisions per hour of UAV flights in the NAS with By using the equations above, (17) could be further rewritten as follows for simplification: Journal of Advanced  Finally based on analysis above, the midair collision risk estimation model can be generated. The expected frequency of fatalities per hour of UAV flights in the NAS can be obtained with In midair collision risk estimation, since severe uncertainty often occurs in uncertain and dynamic environments, the related factors should be incorporated in the model in order to reduce the influence of severe uncertainty. Based on the existing works, it turns out that the model is often quite sensitive to even minor errors in the transition probabilities. As given in (19), based on the gas model and extruded collision volume the model proposed becomes less sensitive to uncertainty and simulation results also proves the robustness.
Simulations over the whole Chinese airspace using real data will be given in the next part.

Model Simulations
(1) Different Parameter Settings. In this part territory areas of each province in China are used, which could be found on the government official websites publicly. The airspace is limited below 50,000 ft. In this way, the airspace volumes of provinces in every Chinese province can be calculated. In order to simulate the number of civil aircraft in the airspace of each province in China, actual civil flight data of taking off and landing in 2015 over Chinese domestic airspace is used, which could be easily found out from CAAC website. Another in the proposed model all the types of aircraft operating in the airspace are supposed to be the same as Boeing 747, which means its frontal area is exactly 180 m 2 given in technical manual from Boeing company. In the same way, the cruse speed will also be the same as 900 km/h. All the UAVs operating in the NAS are set as large UAVs. The core parameters needed in the proposed model, such as frontal exposer area, flying speed and their numbers, are given in Table 10, which are all obtained from the manufacturer manuals. Once all the parameters are determined, the collision risk can be accomplished by (19).
(2) Simulation Results. To develop a preliminary estimate of midair collision risk, the variation of frequency of fatalities spatially over China was investigated, assuming that the UAV was equally likely to be located from sea level to 50,000 ft, neglecting effects of elevation. The model of midair collisions given by (19) was applied to all air traffic from sea level to 50,000 ft over China. The resulting expected level of safety over several regions of the country is shown in Figure 9. Figure 9 shows the expected frequency of fatalities in all the Chinese provinces when different numbers of UAVs operate in the airspace. The numbers of UAVs are set as 2000, 4000, 6000, and 8000, respectively.
There is no doubt that Shanghai and Beijing are the two most dangerous in China no matter how many UAVs are there in the airspace. When 8000 UAVs operate in the airspace, the highest value could arrive at 5.98 × 10 −4 in Shanghai and 2.06 × 10 −4 in Beijing. The former one is two times higher than the latter one. If there are 2000 UAVs operating, the values will be 1.49 × 10 −4 and 5.15 × 10 −5 in Shanghai and Beijing.
Tibet and Qinghai are the safest of all with the lowest expected frequency of fatalities under all the situations. The reason is that the density of civil flights in the airspace is low over the whole year. Their expected frequency of fatalities are 9.85 × 10 −9 and 1.77 × 10 −8 respectively if 2000 UAVs exist in the airspace. It is obvious the expected frequency in Tibet and Shanghai with 2000 and 8000 UAVs are not in the same order of magnitudes.
Tianjin, Shandong, Jiangsu, Zhejiang, Fujian, Henan, Chongqing are seven provinces, which are much safer than Shanghai and Beijing. But these seven provinces are still dangerous compared with the other provinces, which are determined by the density of manned aircraft in the airspace. The expected fatal frequencies of the seven provinces are all above 10 −5 .
Neimenggu, Xinjiang, and Gansu are the three provinces that are more dangerous than Tibet and Qinghai, but much safer compared with the others under all kinds of UAVs settings. The safety levels of Heilongjiang, Jilin, Hebei, ShanxiT, Anhui, Jiangxi, Hubei, Hunan, Guangxi, Sichuan, Yunnan, Guizhou, ShanxiX, and Ningxia are in the medium range.
Based on the analysis above, we could conclude that the majority of the midair collision risk is concentrated over metropolitan areas with major airports by approximately an order of magnitude. The structure of air traffic in the NAS is clear, with large collision risk along several welltraveled routes. The expected level of safety calculated by using this method does not adequately capture the expected level of safety in low density regions. Civil flight density, and therefore collision risk, is expected to be highest on major flight levels and within the airway boundaries, reflecting the operation of the majority of air traffic along airways in the NAS. The structure of operations on flight levels and along airways is likely to create local regions of increased density in dimensions, which is not analyzed by this method in this paper.

Conclusions
This paper has introduced an effective approach for modelling and assessing the risks associated with UAVs integrated into NAS. Both ground impact hazard and midair collision fatality risk are estimated, in which threats to fatalities generated by the two hazards are the focus. Based on system reliability required to meet a target level of safety for different UAVs, a ground impact assessment model is proposed. Both fixed-wing and rotary-wing UAVs are considered under a real scenario. In the model territory and population data, casualty areas, and sheltering factors are all indispensable. Since the fatal injuries yields the probability that an impact may cause a fatality, a random number generation process using this probability, is used to determine if the fatality happens or not for each simulation. A model of aircraft collisions is designed to estimate the midair collision fatality risk based on density of civil flight in different regions over China. The relative collision area and operating speed between UAVs and manned aircraft are constructed to obtain expected frequency of fatalities for each province with official government data. Experimental simulations are made to evaluate the ground impacts and midair collisions when UAVs operate in the NAS. The models in this paper provides a generic framework that can be used to structure the development of safety cases for any UAV operation.

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

Conflicts of Interest
The authors declare that they have no conflicts of interest.