Urban unmanned aerial vehicle ground risk assessment method based on high-precision risk factors

Abstract

The invention provides an urban unmanned aerial vehicle ground risk assessment method based on high-precision risk factors, and belongs to the technical field of unmanned aerial vehicle control. The method comprises the following steps of: 1, classifying satellite images by adopting a K-means clustering method, fusing a classification result with risk factor data, and establishing a refined risk map layer by integrating multi-source data; 2, introducing two risk cost models, namely a personnel death risk cost model and a property loss risk cost model, and accurately quantifying the two types of risks; secondly, comprehensively considering a personnel death risk and a property loss risk, establishing an unmanned aerial vehicle operation comprehensive risk cost model, and realizing unmanned aerial vehicle operation risk assessment; and finally, based on an unmanned aerial vehicle operation risk assessment result, a multi-layer stacked unmanned aerial vehicle operation risk map construction method is adopted. According to the method, the spatial precision, interpretability and availability of the risk map are remarkably improved.

Description

Technical Field

[0001] This invention provides a ground risk assessment method for urban unmanned aerial vehicles (UAVs) based on high-precision risk factors, belonging to the field of UAV control technology. Background Technology

[0002] In recent years, unmanned aerial vehicle (UAV) technology has developed rapidly and has been widely used in various industries, such as logistics delivery, urban inspection, environmental monitoring, and emergency rescue. However, with the rapid increase in the number of UAVs, the risks associated with their operation have attracted widespread attention. In particular, ensuring the safe operation of large numbers of UAVs in urban low-altitude airspace has become crucial. Therefore, researching UAV operation risk assessment methods and their operational mechanisms has become a primary issue.

[0003] Regarding ground risk assessment for unmanned aerial vehicles (UAVs), existing technologies focus on ground injury probability and kinetic energy modeling. Scholars generally use calculations based on the range of kinetic energy parameters under UAV crash or out-of-control conditions to assess the potential degree of injury. Population density, as a key factor, is often based on static population distribution data, combined with shading coefficients to assess the impact of obstructions such as buildings on risk propagation. Some studies introduce dynamic exposure models to reflect changes in population activity. In terms of ground risk map construction, traditional methods rely on static population heat maps. In recent years, some scholars have introduced semantic maps and multi-source data fusion techniques, combining dynamic population flow information with Geographic Information System (GIS) data, and using rasterization to refine the spatial distribution of risk. Furthermore, for risk distribution at different altitude levels, some studies have used hierarchical raster models to characterize UAV crash sites and kinetic energy decay characteristics. In recent years, deep learning methods have been used to predict dynamic population flow trends, assisting in the spatiotemporal updates of dynamic risk assessments.

[0004] The existing technical solutions have the following main drawbacks: Disadvantage 1: Single risk source. Most solutions focus on static population density or simplified kinetic thresholds, without systematically incorporating key exposure factors such as shielding coefficients, facility vulnerability, vehicle / pedestrian flow, activity periods, and temporary gatherings. They also lack sufficient identification of abnormal scenarios such as construction, events, and festivals, resulting in incomplete coverage of risk triggers and an underestimation of actual risks during peak periods and special events.

[0005] Disadvantage 2: Insufficient accuracy of risk maps. They often rely on low-resolution statistical data or a single raster scale, lack semantic segmentation and cross-layer registration, and ignore height / three-dimensional occlusion effects and differences in building form. Insufficient parameter calibration and uncertainty propagation processing result in "overly smoothed" risk value space and blurred boundaries, which do not match the granularity of the management threshold, making it difficult to support refined permitting and control decisions.

[0006] Disadvantage 3: Lack of validation in complex real-world urban scenarios. Evaluations often remain at the level of simulation or small-sample field tests, lacking joint verification with real flight paths, weather, pedestrian flow, and accident / failure records; there is insufficient validation of generalization capabilities across urban areas, time periods, and seasons, and a lack of reproducible engineering evidence chains, making it difficult to prove stable and reliable implementation in high-density, rapidly changing urban environments. Summary of the Invention

[0007] This invention aims to provide a ground-based risk assessment method for urban unmanned aerial vehicles (UAVs) based on high-precision risk factors. Addressing the common shortcomings of existing technologies—such as single risk sources, insufficient risk map accuracy, and lack of verification in real-world complex scenarios—this invention constructs an integrated mechanism encompassing "refined element extraction—multi-layer modeling—cost assessment—3D mapping—resolution conversion." The specific objectives are: (1) A method for fine extraction of ground features based on high-resolution satellite imagery is proposed. K-Means clustering and similarity matching are used to generate a surface semantic classification layer. Multi-source layers such as population density layer, occlusion coefficient layer, and ground obstacle layer are integrated to form an scalable multi-layer risk data model. (2) Establish a risk cost model for death and a risk cost model for property loss, provide a comprehensive risk calculation framework based on weighted normalization, and output the classification risk and comprehensive risk value of each grid. (3) Form a three-dimensional (3D) ground risk map, realize the hierarchical risk distribution estimation at meter-level spatial resolution, and provide a resolution conversion method to adapt to different operational needs and assessment scales.

[0008] By achieving the above technical objectives, this invention can provide high-precision, quantifiable, and verifiable ground risk assessment results in complex urban environments, providing a reliable basis for determining urban drone operation permits, mission planning, and refined management.

[0009] This invention adopts a technical approach of "multi-source data layering → risk cost assessment → three-dimensional raster synthesis": first, semantic classification and layer construction are performed on satellite imagery and urban data; then, quantitative assessment of each raster is performed based on the risk cost of death and property loss; finally, three-dimensional layered synthesis is completed to generate a meter-level high-precision ground risk map.

[0010] The specific technical solution is: a ground-based risk assessment method for urban drones based on high-precision risk factors, including the following steps: Step 1: Extraction of Urban Environmental Risk Map Elements and Layer Modeling First, K-means clustering is used to classify satellite images. The classification results are then fused with risk factor data to create a refined risk map layer by integrating multi-source data.

[0011] (1) High-precision urban satellite image classification (1.1) Land feature classification mechanism based on K-Means clustering The distance from each pixel to each cluster center is calculated using Euclidean distance. u i The distance is calculated, and the optimal clustering result is obtained by minimizing the objective function, the mathematical expression of which is as follows.

[0012] (1)

[0013] J It is the objective function. K For the number of clusters, C i Indicates the first i The number of pixels in each cluster p ij Indicates the first i In the cluster, the first j Data points, u i Indicates the first i The centroids of each cluster. In the K-Means clustering method, the cluster center is the centroid of the cluster. u i By minimizing the objective function J Iterative updates will be performed.

[0014] (1.2) Land cover type identification mechanism based on similarity matching Similarity matching was used to obtain an accurate match between each cluster center and the corresponding land cover feature type. Based on this, 15 land cover types were defined from the selected satellite images. Using pixels u ’ i ( i = 1, 2, ¼, 15) represents land cover types, and the 15 defined land cover types constitute a set. U' = { u' 1, u' 2, ¼, u' 15 The set of cluster centers obtained by the K-Means clustering algorithm is calculated. U With the defined set of land cover types U' The Euclidean distance between them yields the similarity matrix. D It is used to quantify the similarity between cluster centers and land cover types.

[0015] In the similarity matrix D In, each element d ij Indicates the first i The cluster centers and the first jThe Euclidean distance between different land cover types is calculated using the following formula: d ij = || u i - u' j The spatial similarity between cluster centers and land cover types is quantified by calculating the Euclidean distance, thereby achieving accurate matching between cluster centers and corresponding land cover types. The objective function for similarity matching is shown in the formula.

[0016] (2)

[0017] In the above formula, Z It is the objective function. c ij = 1 indicates the first i The cluster centers and the first j Matching of land cover types, c ij =0 indicates the first i The cluster centers and the first j The land cover types do not match.

[0018] (2) Multi-layer modeling of urban risk maps By combining the urban satellite image land cover classification results in (1) with various risk factors, a multi-layer model of urban risk map is established, including a population density layer, a ground vehicle layer, a shading factor layer, and an obstacle height layer, with each layer corresponding to different urban geographical and social characteristics.

[0019] (2.1) Population density layer N p This represents the total population within the city's ground grid. It is the first grid cell i Population density corresponding to crop type A i It is the first i The total area of ​​crops within the grid; establish a linear equation based on the formula to describe the total population in the urban ground grid. N p Area of ​​each land cover type A i Population density corresponding to each land cover type The relationship between them is shown in the formula.

[0020] (6)

[0021] A n'k Represents ground grid n The Middle kArea of ​​crop type, It is the first k Population density corresponding to crop type N n It is the ground grid n The total population; solved using multiple linear regression, where the intercept term is zero, ultimately yielding the population density corresponding to each land cover type, i.e. .

[0022] The population density of the same land cover type in different areas of the city is corrected using WorldPop data, as shown in the formula.

[0023] (7)

[0024] To correct the population density of a certain land cover type in different locations within the city, Population density corresponding to land cover type This represents the population density at that location in the WorldPop data. The weight of population density corresponding to land cover type. The weights for WorldPop population density are denoted as , where . and .

[0025] (2.2) Ground vehicle level Ground vehicle density is simulated using a normal distribution, as shown in the formula. Vehicle density of the urban surface transportation network. The average vehicle density of the urban surface transportation network. The standard deviation of vehicle density within the urban surface transportation network. It follows a normal distribution The random value.

[0026] (8)

[0027] (2.3) Masking factor layer The fatality rate based on the kinetic energy of a drone impact is introduced, as shown in the formula. R fatality For the fatality rate, E k The kinetic energy of the drone when it hits a person on the ground. S f It is a masking factor. , yes S f = 0.5 is the impact energy required to achieve a 50% fatality rate. yesS f = 0 is the threshold of impact kinetic energy required to cause death to ground personnel.

[0028] (9)

[0029] Different shading factors were calculated based on the formula. S f The mortality rate is below, of which = 10 6 J , = 100 J .

[0030] (2.4) Obstacle Height Layer The 15 land features were further divided into three types: the first type is open and bare areas with a shading factor of zero; the second type is woody plants; and the third type is urban buildings. The first category of land cover consists of open, bare areas, with a height value of zero. For the second category, woody plants, a normal distribution is used to simulate the actual height of woody plants in the urban environment, randomly generating their height values. The height range of woody plants is set between 3 and 8 meters, with a normal distribution mean of 5.5 meters and a standard deviation of 1.25 meters, as shown in the formula. h tree The height of woody plants within the urban area. The average height of woody plants. The standard deviation of the height of woody plants. h random It follows a normal distribution The random value.

[0031] (10)

[0032] The third category of urban buildings directly uses the building height dataset and employs a random forest model to estimate building height, generating building height data.

[0033] Step 2: Risk Assessment and Risk Map Construction for Urban Operation of Unmanned Aerial Vehicles First, two risk cost models were introduced: a personal death risk cost model and a property loss risk cost model, which accurately quantified these two types of risks. Second, by comprehensively considering both personal death risk and property loss risk, a comprehensive risk cost model for UAV operation was established to achieve UAV operation risk assessment. Finally, based on the results of the UAV operation risk assessment, a multi-layered overlay method for constructing a UAV operation risk map was adopted.

[0034] (1) Mortality Risk Cost Model The formula provides a cost model for the risk of death, in which the risk of death is further divided into the risk of death of pedestrians on the ground and the risk of death of people inside vehicles.

[0035] (11)

[0036] In the above formula, C r_f For the cost of personnel death risk, C r_p The risk of death due to an unmanned vehicle colliding with a pedestrian on the ground. C r_v The risk of death due to drones colliding with ground vehicles; (2) Property loss risk cost model The formula provides a model for the cost of property loss risk, in which C r_b For the cost of property loss risk, This is the building type coefficient, with a value of [value missing]. , The function is a log-normal distribution, which can be obtained from the formula. h b Indicates the height of the building. and These are the mean and standard deviation of the logarithmic variable, respectively.

[0037] (25)

[0038] (26)

[0039] The coefficient reflects the differences in economic losses suffered by different types of buildings after being hit by drones. (3) Comprehensive risk cost model for UAV operation Based on the personnel mortality risk cost model and the property loss risk cost model derived from the formula, a comprehensive risk cost model for UAV operation is presented; a normalization coefficient is proposed. Its value range is (0, 1), which converts the two to the same order of magnitude.

[0040] (27)

[0041] In the above formula, C u The overall risk cost of drone operation; and These are the weighting coefficients for the risk of personal death and the risk of property loss, respectively. , + = 1; The normalization coefficient for risk cost is defined as follows: ,in C i_max It is the maximum value of all types of risk costs; C i These are the risks of personal death and property loss, respectively. C 1= C r_f , C 2= C r_b .

[0042] To construct a risk map of drone operations in urban areas, the urban airspace is divided into several voxel grids. Based on the comprehensive risk cost of drone operations obtained from the formula, the corresponding risk value is calculated for each voxel grid, and each voxel grid is defined as having a center point. C xyz The rectangular prism; the resolution of the risk map is converted, and the conversion method is shown in the formula. C u,n_n_m ( i , j , k ) represents the average risk value of the voxel grid in the new risk map after resolution conversion, where the size of the voxel grid is . n ´ n ´ m ; i , j and k The new map voxel grid is in x , y and z Index of direction; a , b and c This is used to traverse all voxel locations of the original high-resolution risk map within the new grid to calculate the average risk value of the new grid.

[0043] (28)

[0044] The technical effects of the present invention are as follows: 1. Construct a method for extracting elements and modeling layers in urban environmental risk maps, unifying ground risk elements such as population density, ground vehicle exposure, shading factors, and obstacle height into the same data and mathematical framework, and outputting meter-level raster under a unified coordinate and time benchmark; forming a complete link from data processing and parameter calibration to permit threshold docking and audit traceability, significantly improving the spatial accuracy, interpretability, and usability of risk maps.

[0045] 2. A weighted normalized assessment framework for the risk costs of personal death and property loss is proposed. This framework generates a 3D ground risk map based on multi-layer overlay and provides resolution conversion capabilities, reducing two types of biases: "overly conservative" and "underestimating risk." An engineering verification process (data processing → parameter calibration → result reproduction) is also provided, allowing for incremental updates and cross-regional reuse by replacing data sources such as population and transportation as needed. This provides a high-precision, quantifiable, and reusable risk map for operational permit determination, mission pre-planning, and site selection, reducing third-party ground risks and operational uncertainties. Attached Figure Description

[0046] Figure 1 This is a schematic diagram illustrating the composition of the urban risk map layers of the present invention; Figure 2 This is a schematic diagram of the population density calculation method of the present invention; Figure 3 The effect of impact kinetic energy on lethality under different shielding factors according to the present invention; Figure 4 The ground area affected by the crash of the drone of this invention; Figure 5 The present invention addresses the safety threat posed by unmanned aerial vehicles (UAVs) to road traffic. Figure 6 This is a schematic diagram of the risk map resolution conversion of the present invention. Detailed Implementation

[0047] A ground-based risk assessment method for urban drones based on high-precision risk factors includes the following steps: Step 1: Methods for Extracting Elements and Modeling Layers in Urban Environmental Risk Maps The risks associated with drone operations in cities are complex and diverse, encompassing features such as urban population, ground vehicles, and buildings. Accurately extracting risk factors is a fundamental prerequisite for drone risk assessment and risk map construction. First, K-means clustering is used to classify satellite images, providing a foundation for extracting detailed risk factor data, such as population density and building height. Second, the classification results are fused with the risk factor data to establish a refined risk map layer by integrating multi-source data.

[0048] (1) High-precision urban satellite image classification In order to effectively extract key elements of urban environmental risk maps, it is necessary to accurately identify and classify urban environmental features in order to accurately assess the risks of drone operation in cities.

[0049] (1.1) Land feature classification mechanism based on K-Means clustering First, satellite images are acquired via a map API. The images consist of three channels ( R ,G , B It consists of three features per pixel, with a brightness range of 0 to 255, thus eliminating the need for data normalization. It is assumed that the pixel... P ( x , y The eigenvectors of ) are P i =( p 1, p 2, p 3) The area of ​​each pixel in the actual map is 1. m × 1 m To divide the data into clusters, Euclidean distance is used to calculate the distance from each pixel to each cluster center. u i The distance is calculated, and the optimal clustering result is obtained by minimizing the objective function, the mathematical expression of which is as follows.

[0050] (1)

[0051] In the formula, J It is the objective function. K For the number of clusters, C i Indicates the first i The number of pixels in each cluster p ij Indicates the first i In the cluster, the first j Data points, u i Indicates the first i The centroids of each cluster. In the K-Means clustering method, the cluster center is the centroid of the cluster. u i By minimizing the objective function J The goal of iterative updates is to assign data points to the nearest cluster centers and optimize clustering results by adjusting the positions of the cluster centers.

[0052] The pseudocode for the K-Means clustering algorithm is shown in Table 1. First, input all prime points. P and initial random cluster centers U Next, the iteration begins. Each data point is assigned to the nearest cluster center. The position of the cluster center is updated based on all the data points assigned to that cluster. This iteration continues until the cluster center positions no longer change. Finally, the results include cluster partitioning. C and the corresponding set of cluster centers U .

[0053] Table 1. Pseudocode for the K-Means clustering algorithm

[0054] (1.2) Land cover type identification mechanism based on similarity matching Since the K-Means clustering algorithm cannot directly determine the urban land cover type corresponding to the cluster center, similarity matching must be used to obtain an accurate match between each cluster center and the corresponding land cover feature type. Based on this, 15 land cover types were defined from the selected satellite images, as shown in Table 2. The land cover types were divided into 5 categories: Category 1 includes grassland, lakes, concrete surfaces, loess, and blue construction sites; Category 2 is roads; Category 3 includes woody plants, blue roofs, and red roofs; Category 4 includes low-rise buildings (brown, gray, and red) and industrial areas (silver-gray and blue); and Category 5 is high-rise buildings.

[0055] Table 2 Classification of Urban Land Feature Types

[0056] Using pixels u ’ i ( i = 1, 2, ¼, 15) represents land cover types, and the 15 defined land cover types constitute a set. U' = { u' 1, u' 2, ¼, u' 15 The set of cluster centers obtained by the K-Means clustering algorithm is calculated. U With the defined set of land cover types U' The Euclidean distance between them yields the similarity matrix. D It is used to quantify the similarity between cluster centers and land cover types.

[0057] In the similarity matrix D In, each element d ij Indicates the first i The cluster centers and the first j The Euclidean distance between different land cover types is calculated using the following formula: d ij = || u i - u' j The spatial similarity between cluster centers and land cover types is quantified by calculating the Euclidean distance, thereby achieving accurate matching between cluster centers and corresponding land cover types. The objective function for similarity matching is shown in the formula.

[0058] (2)

[0059] In the above formula, Z It is the objective function. c ij = 1 indicates the first i The cluster centers and the first j Matching of land cover types, c ij =0 indicates the first i The cluster centers and the first j Since the land cover types do not match, the similarity matching method can effectively match the cluster center with the land cover type, thereby improving the accuracy of land cover type identification.

[0060] (2) Multi-layer modeling of urban risk maps To accurately describe various risk factors in the urban environment and construct an urban risk map based on them, the urban satellite image land cover classification results in (1) above are combined with various risk factors to establish a multi-layer model of the urban risk map, such as... Figure 1 As shown, the map includes layers for population density, ground vehicles, occlusion factors, and obstacle height. Each layer corresponds to different urban geographical and social characteristics, and each layer of the urban risk map has a significant impact on the safety of drone operations.

[0061] (2.1) Population density layer Population density layers quantify the ground population density and its distribution within urban areas, which is crucial for drone ground risk assessment and reflects the number of ground personnel who may be hit by drones. Therefore, how to construct population density layers with acceptable resolution has become a key issue.

[0062] (2.1.1) Data Acquisition and Preprocessing To effectively improve the accuracy and authenticity of urban population density, urban population data was obtained through two methods. First, WorldPop population density distribution data, which provides global population distribution and density, is sourced from various channels such as satellite remote sensing images, geographic statistics, and ground surveys, greatly improving the accuracy and reliability of the data. Second, the total population data from the Seventh National Population Census Bulletin. The Seventh National Population Census is an official comprehensive census conducted by the Chinese government every ten years, covering data on basic population information, economic status, education level, living conditions, and other aspects nationwide.

[0063] First, we used 2020 population density data for China provided by WorldPop. This population density data uses the WGS84 geographic coordinate system, dividing the urban area into several grids, with each grid cell being 100 pixels in size. m × 100 mWorldPop provides the population count within each grid cell; secondly, it uses the total population from the Seventh National Population Census to correct the population density data provided by WorldPop.

[0064] (2.1.2) Population density estimation Because the population density data provided by WorldPop is outdated and cannot describe the current urban population characteristics in real time, a correction coefficient for urban population density is proposed based on national census data to effectively solve this problem, as shown in the formula.

[0065] (3)

[0066] In this formula, C c This is the population density correction factor, which is the total population provided in the national census data. N census Compared to the total population given in WorldPop N WorldPop The ratio of this to the population density correction factor C c The population density of a city is corrected by multiplying it by the population density of each city's ground grid in WorldPop, as shown in the formula.

[0067] (4)

[0068] In the formula, This is the corrected population density of the urban ground grid. This is the initial population density of the city's ground grid in WorldPop.

[0069] On the other hand, the population size within each urban ground grid is related to the population density corresponding to each land cover type within that grid, such as... Figure 2 As shown in the figure, the top layer represents the total population within the city's ground grid, with a grid size of 100. m ×100 m The middle layer represents the area of ​​different land cover types within the ground grid, with each land cover type having an area of ​​1. m × 1 m The bottom layer represents the population density corresponding to each land cover type within the ground grid, and the ground grid size is 1. m × 1 m Therefore, a relationship model can be established between the total population, land cover type, and population density corresponding to land cover type within each urban ground grid, as shown in the formula.

[0070] (5)

[0071] In the formula, N p This represents the total population within the city's ground grid. It is the first grid cell i Population density corresponding to crop type A i It is the first i The total area of ​​crops within this grid, and since the urban ground grid size is 100 m ×100 m With a resolution of 50 × 50, a total of 2500 linear equations can be established according to the formula to describe the total population in the urban ground grid. N p Area of ​​each land cover type A i Population density corresponding to each land cover type The relationship between them is shown in the formula.

[0072] (6)

[0073] In the formula, A n'k Represents ground grid n The Middle k Area of ​​crop type, It is the first k Population density corresponding to crop type N n It is the ground grid n The total population can be calculated using a multiple linear regression model, since the formula is based on this model. It's important to note that the intercept term in this regression model is zero. Ultimately, the population density for each land cover type can be obtained. .

[0074] Since the same land cover type has different population densities in different locations in the urban environment, for example, the sidewalks near large shopping malls have a large difference in population density compared to the sidewalks in residential areas, and the population density of the same land cover type may vary in different locations, this invention uses WorldPop data to correct the population density of the same land cover type in different areas of the city, as shown in the formula.

[0075] (7)

[0076] In the formula, To correct the population density of a certain land cover type in different locations within the city, Population density corresponding to land cover type This represents the population density at that location in the WorldPop data. The weight of population density corresponding to land cover type. The weights for WorldPop population density are denoted as , where . and .

[0077] (2.2) Ground vehicle level To calculate the risk of death of occupants in vehicles caused by a drone crashing into a ground vehicle, it is necessary to calculate the vehicle density of urban roads to further construct a ground vehicle layer model. This model describes the distribution of vehicles in the urban area and quantifies the risk effect of a drone colliding with a ground vehicle. Specifically, based on the similarity matching-based land cover type identification mechanism proposed in (1.1), an accurate urban ground traffic network can be obtained. How to calculate the vehicle density within the urban ground traffic network becomes a key issue.

[0078] Due to the lack of real-world urban vehicle density data, a normal distribution is used to simulate ground vehicle density, as shown in the formula. Vehicle density of the urban surface transportation network. The average vehicle density of the urban surface transportation network. The standard deviation of vehicle density within the urban surface transportation network. It follows a normal distribution The random value.

[0079] (8)

[0080] (2.3) Masking factor layer The shielding factor refers to the degree of protection that buildings or trees can provide to people on the ground when a drone crashes. The shielding factor is a positive value that reflects the strength of the shielding effect. The higher the value, the higher the degree of protection provided. The shielding factor is crucial in assessing the risk of human death because it can reduce the impact kinetic energy when a drone crashes, thereby reducing the risk of human death. If the shielding factor is ignored, the risk of drone operation in cities cannot be accurately assessed.

[0081] To accurately quantify the shielding factor, a fatality rate based on the impact kinetic energy of the drone is first introduced, as shown in the formula. R fatality For the fatality rate, E k The kinetic energy of the drone when it hits a person on the ground. S f It is a masking factor. , yes S f = 0.5 is the impact energy required to achieve a 50% fatality rate. yes S f = 0 is the threshold of impact kinetic energy required to cause death to ground personnel.

[0082] (9)

[0083] To further analyze the relationship between ground-based personnel fatality rate and shading factor, this invention calculates different shading factors based on the formula. S f The mortality rate is below, of which = 10 6 J , = 100 J The result is as follows Figure 3 As shown in the figure, the occlusion factor... S f The larger the size, the higher the mortality rate. R fatality The smaller the value, the higher the mortality rate. R fatality The larger the occlusion factor, the better. S f When the kinetic energy of the impact remains constant, E k The increase in mortality rate R fatality It grows exponentially.

[0084] This invention categorizes urban landforms into 5 major categories and 15 subcategories, and defines shading factors for each landform type, as shown in Table 3. Fixed values ​​of 0, 0.25, 0.50, 0.75, and 1 are selected to evaluate the degree of protection provided by different landform types to ground personnel during drone crashes. Specifically, grasslands, lakes, concrete surfaces, loess soil, and blue construction sites have a shading factor of 0, meaning no protection is provided when a drone crashes in these areas; roads have a shading factor of 0.25, providing some protection; woody plants, blue roofs, and red roofs have a shading factor of 0.50; low-rise buildings (brown, gray, red) and industrial areas (silver-gray, white-blue) have a shading factor of 0.75; and high-rise buildings have a maximum shading factor of 1, providing the best protection.

[0085] Table 3. Shading factors corresponding to ground features

[0086] (2.4) Obstacle Height Layer To generate accurate obstacle height layers, the 15 terrain features in Table 2 are further divided into three types: the first type is open and bare areas, such as grasslands and lakes, with a shading factor of zero; the second type is woody plants; and the third type is urban buildings, which are usually more than 4 meters high and will significantly affect the safety and efficiency of UAV operation.

[0087] Since the first type of land cover consists of open, bare areas, no height value needs to be assigned to it; that is, the height value is zero. For the second type of woody vegetation, due to the lack of uniform and accurate data, this invention uses a normal distribution to simulate the actual height of woody vegetation in the urban environment, randomly generating its height value. The height range of woody vegetation is set between 3 meters and 8 meters, with a normal distribution mean of 5.5 meters and a standard deviation of 1.25 meters, as shown in the formula. h tree The height of woody plants within the urban area. The average height of woody plants. The standard deviation of the height of woody plants. h random It follows a normal distribution The random value.

[0088] (10)

[0089] The third category of urban buildings directly uses the China Building Height Dataset (CNBH-10 m). This dataset integrates multi-source, all-weather ground observation data such as radar, optical, and nighttime light imagery, and uses a random forest model to estimate building heights, generating building height data for the whole country in 2020. It can effectively characterize the refined height features of buildings in urban areas.

[0090] Step 2: Methodology for Risk Assessment and Risk Map Construction of Urban Unmanned Aerial Vehicles To accurately assess the risks of drone operation in cities, this invention first introduces two risk cost models: a personal death risk cost model and a property loss risk cost model, which precisely quantify these two types of risks. Second, by comprehensively considering both personal death and property loss risks, a comprehensive risk cost model for drone operation is established to achieve drone operation risk assessment. Finally, based on the results of the drone operation risk assessment, a multi-layered overlay method for constructing a drone operation risk map is proposed.

[0091] (1) Mortality Risk Cost Model The risk of human death quantifies the number of deaths caused by a drone colliding with a pedestrian or vehicle on the ground, expressed as the number of deaths per hour. The formula provides a cost model for the risk of human death, in which the risk of human death is further divided into the risk of death of pedestrians on the ground and the risk of death of people inside vehicles.

[0092] (11)

[0093] In the above formula, C r_f For the cost of personnel death risk, C r_p The risk of death due to an unmanned vehicle colliding with a pedestrian on the ground. C r_v The following section will address the risk of death from drone collisions with ground vehicles. C r_p and C r_v The modeling methods will be discussed in detail.

[0094] (1.1) Ground pedestrian mortality risk model The risk of pedestrian death due to drone malfunction and subsequent collision should be considered from several angles. First, the probability of drone malfunction and collision should be taken into account; this is the likelihood of the drone crashing into the ground, and this value is related to the drone's inherent stability. Second, the probability of the drone colliding with a pedestrian should be considered, which is related to the population distribution and density on the ground. Finally, the probability of death resulting from a drone collision should also be considered, which is related to the impact kinetic energy of the drone at the time of the collision. In summary, the risk of pedestrian death is shown in the formula.

[0095] (12)

[0096] In the formula, C r_p The risk of death due to an unmanned vehicle colliding with a pedestrian on the ground; P crash Let be the probability of the drone failing in mid-air and crashing to the ground. Due to the high complexity and randomness of the drone crash process, this value is set to 10 to simplify the model and reduce the difficulty of solving the problem. -4 times / hour; P p_i This represents the total number of pedestrians in the urban area affected by the drone impact. R p_f This represents the probability of a pedestrian dying after being hit by a drone, also known as the pedestrian fatality rate.

[0097] For accurate calculation P p_iThe numerical value requires a comprehensive analysis of the drone's descent process. Typically, the drone's descent process can be categorized into ballistic descent, uncontrolled gliding, parachute descent, and active escape modes. This invention only considers the ballistic descent mode, where the drone falls along an approximate ballistic trajectory after losing most of its power. This descent process is entirely determined by the drone's aerodynamics. Based on the above analysis, the formula gives the total number of pedestrians on the ground affected by the drone impact. P p_i The calculation method.

[0098] (13)

[0099] In the above formula, a This is a weighting coefficient, which can be adjusted based on the terrain features at the drone crash site. N p The total number of pedestrians in the area hit by the drone can be calculated using the formula, where, S u_hit This refers to the area of ​​impact on the ground when the drone crashes. The population density within the ground grid affected by the drone impact can be calculated using a formula.

[0100] (14)

[0101] In the above formula, S u_hit The description covers the area of ​​impact on the ground when the drone crashed, such as... Figure 4 As shown in the figure, the blue area represents the area on the ground affected by the drone. This range is related to the drone's fall height, descent angle, three-dimensional dimensions of the drone, and the characteristics of pedestrians on the ground. The calculation method is shown in the formula.

[0102] In the formula, L u This indicates the maximum size of the drone. H p The average height of pedestrians on the ground is set to 1.75 meters in this model. Similarly, R p The average radius for pedestrians on the ground is set to 0.25 meters. This refers to the descent angle of the drone when it is unpowered.

[0103] (15)

[0104] The formula gives the descent angle The calculation method is described in this formula. v u The speed at which the drone descends. mu For the quality of drones, The air density at sea level is taken as 1.225. kg / m 3 , C D The drag coefficient, A u This indicates the windward area of ​​the drone. g This is the acceleration due to gravity, usually taken as 9.8. m / s 2 .

[0105] (16)

[0106] In addition, when quantifying the risk of pedestrian death on the ground, the impact kinetic energy when the drone hits the ground and the shading factor of different ground objects should also be calculated. The latter can be obtained from Table 3, while the former is related to the mass of the drone and the speed of the drone when it hits the ground. The calculation method is shown in the formula.

[0107] (17)

[0108] In the above formula, E k This refers to the impact kinetic energy when the drone collides with the ground. m u For the quality of drones, v hit It is the speed at which the drone hits the ground. The speed value is related to its acceleration during descent. As the drone is affected by both vertical drag and its own weight during descent, the former is related to the drag coefficient, the drone's frontal area, air density, and the drone's actual airspeed. The calculation formula is shown in the figure.

[0109] (18)

[0110] In the formula, F d For vertical resistance, C D The drag coefficient, A u This refers to the windward area of ​​the drone. The density of air is 1.225. kg / m 3 , v TAS Let be the true airspeed of the crashing drone; therefore, the drone's acceleration can be expressed by a formula.

[0111] (19)

[0112] In the formula, a For the acceleration of the drone, F g For the drone's own weight, F g = m u g Therefore, the crashed drone was t The velocity at the moment of impact with the ground can be calculated using a formula. v hit For drones t The speed at which it hits the ground at any moment. h This represents the altitude at which the drone crashed.

[0113] (20)

[0114] (1.2) Risk of death of occupants A drone malfunctioning in mid-air and colliding with a vehicle on the ground could pose a risk of death to the occupants. Figure 5 As shown, this risk is related to the probability of drone failure in the air and the probability of death of occupants in a vehicle after a drone collides with a ground vehicle. The formula provides a method for calculating the risk of death of occupants in a vehicle.

[0115] In the formula, C r_v For the risk of death of occupants in the vehicle, P crash The probability of a drone failing in mid-air and crashing into the ground. P v_i This represents the total number of ground vehicles in the urban area affected by drone impacts. R v_f The mortality rate of occupants in ground traffic accidents.

[0116] (twenty one)

[0117] P v_i The calculation method and P p_i Similarly, as shown in the formula, in this formula, S u_hit This refers to the area of ​​impact on the ground when the drone crashes. It refers to the density of vehicles on the ground.

[0118] (twenty two)

[0119] In the formula, R v_fThe mortality rate of vehicle occupants in ground traffic accidents can be calculated using the formula. E v_k This refers to the kinetic energy of a drone colliding with a vehicle on the ground. S f It is a masking factor. , yes S f = 0.5 is the impact energy required to achieve a 50% fatality rate. yes S f = 0 is the threshold of impact kinetic energy required to cause death to ground personnel.

[0120] (twenty three)

[0121] The formula gives the kinetic energy of the drone when it collides with a ground vehicle. E v_k The calculation method is as follows, in this formula, m u The mass of the drone when it collides with a ground vehicle. v u The speed at which the drone collides with the vehicle. v car By setting the road speed limit, the risk of death for occupants can ultimately be determined. C r_v .

[0122] (twenty four)

[0123] (2) Property loss risk cost model When drones operate in low-altitude urban airspace, they are highly susceptible to collisions with urban buildings, posing a risk of property damage. The formula provides a cost model for this property damage risk, in which... C r_b For the cost of property loss risk, This is the building type coefficient, with a value of [value missing]. It can be set according to different building types. The log-normal distribution function can be obtained from the formula. h b Indicates the height of the building. and These are the mean and standard deviation of the logarithmic variable, respectively.

[0124] (25)

[0125] (26)

[0126] The coefficients reflect the differences in economic losses caused by drone strikes to different types of buildings. Table 4 shows the coefficients for different building types in the city. The specific values, including those for the blue and red roofs. The values ​​are all 0.25, indicating that the economic losses caused by drone strikes to such buildings are relatively low; low-rise buildings (brown, gray, red) A value of 0.50 indicates that the economic losses following the damage are at a moderate level; industrial areas (silver-gray, white-blue) and high-rise buildings The values ​​were 0.80 and 1.00, respectively, indicating that the economic losses caused by collisions to industrial areas and high-rise buildings were significantly higher than those to other building types.

[0127] Table 4 Building Type Coefficient

[0128]

[0129] (3) Comprehensive risk cost model for UAV operation Based on the personnel death risk cost model and the property loss risk cost model derived from the formula, the formula presents a comprehensive risk cost model for UAV operation. Since the personnel death risk cost and the property loss risk cost have different orders of magnitude, this invention proposes a normalization coefficient. Its value range is (0, 1), which converts the two to the same order of magnitude.

[0130] (27)

[0131] In the above formula, C u The overall risk cost of drone operation; and These are the weighting coefficients for the risk of personal death and the risk of property loss, respectively. , + = 1; The normalization coefficient for risk cost is defined as follows: ,in C i_max It is the maximum value of all types of risk costs; C i These are the risks of personal death and property loss, respectively. C 1= C r_f , C 2= C r_b .

[0132] To construct a risk map of drone operations in urban areas, the urban airspace is divided into several voxel grids. Based on the comprehensive risk cost of drone operations obtained from the formula, a corresponding risk value is calculated for each voxel grid, such as... Figure 6 As shown, each voxel grid is defined as having a center point. C xyz The rectangular prism, and because the grid size in the multi-layer risk modeling of cities in this invention is 1 m ´ 1 m ´ 1 m Therefore, the resolution of the risk map needs to be converted, and the conversion method is shown in the formula. (28)

[0133] In this formula, C u,n_n_m ( i , j , k ) represents the average risk value of the voxel grid in the new risk map after resolution conversion, where the size of the voxel grid is . n ´ n ´ m ; i , j and k The new map voxel grid is in x , y and z Index of direction; a , b and c This is used to traverse all voxel locations of the original high-resolution risk map within the new grid to calculate the average risk value of the new grid.

Claims

1. A city unmanned aerial vehicle ground risk assessment method based on high-precision risk elements, characterized in that, Comprising the following steps: Step one, urban environment risk map element extraction and layer modeling First, the K-means clustering method is used to classify satellite images, and the classification results are fused with risk factor data to establish a refined risk factor layer by integrating multi-source data; Step two, unmanned aerial vehicle urban operation risk assessment and risk map construction First, two risk cost models are introduced, namely personnel death risk cost model and property loss risk cost model, which accurately quantify these two types of risks; Second, considering the personnel death risk and property loss risk, an unmanned aerial vehicle operation comprehensive risk cost model is established to realize unmanned aerial vehicle operation risk assessment; Finally, based on the unmanned aerial vehicle operation risk assessment results, a multi-layer superimposed unmanned aerial vehicle operation risk map construction method is adopted. 2.The high-precision risk element-based urban UAV ground risk assessment method according to claim 1, characterized in that, Step one specifically includes: (1) High-precision urban satellite image classification (1.1) K-Means clustering-based ground feature classification mechanism The Euclidean distance is used to calculate the distance of each pixel point to each cluster center u i The optimal clustering result is obtained by minimizing the objective function, and the mathematical expression of the objective function is as follows: ； J It is the objective function. K For the number of clusters, C i Indicates the first i The number of pixels in each cluster p ij Indicates the first i In the cluster, the first j Data points, u i Indicates the first i The centroid of a cluster; in the K-Means clustering method, the cluster center is... u i By minimizing the objective function J Perform iterative updates; (1.2) Similarity matching-based ground feature type identification mechanism Similarity matching is used to obtain the accurate matching of each cluster center and the corresponding ground feature type. Based on this, 15 types of ground features are defined from the selected satellite images; Using pixels u ’ i ( i = 1, 2, ¼, 15) represents land cover types, and the 15 defined land cover types constitute a set. U' ={ u' 1, u' 2, ¼, u' 15 The set of cluster centers obtained by the K-Means clustering algorithm is calculated. U With the defined set of land cover types U' The Euclidean distance between them yields the similarity matrix. D This is used to quantify the similarity between cluster centers and land cover types; In the similarity matrix D each element d ij represents the Euclidean distance between the first i cluster center and the first j feature type, and the calculation formula is d ij = || u i - u' j ||, the Euclidean distance obtained by calculation quantifies the spatial similarity between the cluster center and the feature type, thereby realizing accurate matching of the cluster center and the corresponding feature type, and the objective function of the similarity matching is: ； In the above formula, Z is the objective function, c ij = 1 indicates that the first i cluster center matches the first j feature type, c ij = 0 indicates that the first i cluster center does not match the first j feature type; (2) Urban risk map multi-layer modeling The urban satellite image ground feature classification results are combined with various risk factors to establish a multi-layer model of the urban risk map, including population density layer, ground vehicle layer, shelter factor layer and obstacle height layer, each layer corresponding to different urban geographical and social characteristics; (2.1) Population density layer N p represents the total population number within the urban ground grid, is the total area of the first land object in the grid; and i is the population density corresponding to the first land object type, A i is the total area of the first land object in the grid; and i is the population density corresponding to the first land object type, N p is the total area of the first land object in the grid; and A i is the population density corresponding to the first land object type, is the total area of the first land object in the grid; and ; A n'k representing the ground grid n the area of the k ground object type, the population density corresponding to the k ground object type, N n the total population of the ground grid n ; solved by the multiple linear regression method, in which the intercept term is zero, and finally the population density corresponding to each ground object type is obtained, that is, ; The WorldPop data is used to correct the population density of the same ground feature type in different urban areas, as shown in the formula: ； to correct the population density of the certain land feature type at different locations in the city, to correct the population density of the certain land feature type at different locations in the city, to correct the population density of the certain land feature type at different locations in the city, to correct the population density of the certain land feature type at different locations in the city, to correct the population density of the certain land feature type at different locations in the city, to correct the population density of the certain land feature type at different locations in the city, ; (2.2) Ground vehicle layer The normal distribution is used to simulate the ground vehicle density, as shown in the formula: ； the vehicle density for the urban ground traffic network, the average vehicle density for the urban ground traffic network, the standard deviation of the vehicle density within the urban ground traffic network, is a random value subject to a normal distribution ; (2.3) Shelter factor layer The personnel fatality rate based on the kinetic energy of unmanned aerial vehicle impact is introduced, as shown in the formula: ； R fatality is the lethality rate for a person, E k is the kinetic energy of the drone when it hits a person on the ground, S f is the obscuration factor, S f Î(0, 1], is S f is the impact energy required for a 50% lethality rate when = 0.5, is S f is the threshold impact kinetic energy required to cause a death of a person on the ground when = 0; The lethal rate under different shelter factors is calculated according to the formula S f wherein = 10 6 J , = 100 J ; (2.4) Obstacle height layer The 15 types of ground features are further divided into three types, the first type is open bare area with a shelter factor of zero, the second type is woody plants, and the third type is urban buildings; The first type of ground object is an open bare area, i.e. the height value is zero; for the second type of woody plants, the actual height of the woody plants in the urban environment is simulated by using a normal distribution to randomly generate the height value, and the height range of the woody plants is set to be between 3 meters and 8 meters, the mean value of the normal distribution is 5.5 meters, and the standard deviation is 1.25 meters, as shown in the following formula, h tree is the height of the woody plants in the urban area, is the average height of the woody plants, is the standard deviation of the height of the woody plants, h random is a random value subject to a normal distribution . ； The third type of urban buildings directly uses the building height data set and uses the random forest model to estimate the building height to generate the building height data. 3.The high-precision risk factor-based urban UAV ground risk assessment method according to claim 1, characterized in that, Step two specifically includes: (1) Personnel death risk cost model The following formula gives the personnel death risk cost model, in which the personnel death risk is further divided into ground pedestrian death risk and in-vehicle personnel death risk; ； In the above formula, C r_f is the cost of the risk of death to the person, C r_p is the cost of the risk of death to the person caused by the unmanned vehicle hitting the ground, C r_v is the cost of the risk of death to the person caused by the unmanned vehicle hitting the ground, (2) Property loss risk cost model The following formula gives the property loss risk cost model, where C r_b is the property loss risk cost, is the building type coefficient, which takes the value , is the lognormal distribution function, which is obtained from the following formula, h b represents the building height, and are the mean and standard deviation of the log variable, respectively; ； The coefficient reflects the differences in economic losses suffered by different types of buildings after being hit by drones. (3) Unmanned aerial vehicle operation comprehensive risk cost model Based on the personnel death risk cost model and the property loss risk cost model, the following formula shows the comprehensive risk cost model of the UAV operation; the normalization coefficient is proposed , whose value range is (0, 1), and both are converted to the same order of magnitude; ； In the above formula, C u is the comprehensive risk cost of the UAV operation; and are weight coefficients of the personnel death risk and the property loss risk, respectively, and , + = 1; is a normalization coefficient of the risk cost, defined as wherein C i_max is the maximum value of each type of risk cost; C i are the personnel death risk and the property loss risk, respectively, wherein C 1= C r_f , C 2= C r_b ; In order to construct the unmanned aerial vehicle city operation risk map, the city airspace is divided into a plurality of voxel grids, a corresponding risk value is calculated for each voxel grid according to the unmanned aerial vehicle operation comprehensive risk cost, each voxel grid is defined as a cuboid with a center point C xyz ; resolution conversion of the risk map is performed, and the conversion method is as shown in the formula: ； C u,n_n_m ( i , j , k ) are the average risk values of the voxel grid of the new risk map after resolution conversion, where the size of the voxel grid is n ´ n ´ m ; i , j and k are the indices of the new map voxel grid in x , y and z directions, respectively; a , b and c are used to iterate through all voxel positions of the original high-resolution risk map within the new grid to calculate the average risk values of the new grid.

Images