Method and system for generating dynamic risk map of unmanned aerial vehicle operation on ground

Abstract

The invention provides a method and a system for generating a dynamic risk map of unmanned aerial vehicle operation on the ground, which provide detailed risk assessment for unmanned aerial vehicle flight, identify potential risks of a flight area, avoid a high-risk area, improve flight safety and provide a scientific basis for low-altitude flight path planning. And the damage to ground personnel and facilities caused by falling of the unmanned aerial vehicle or other accidents is reduced, so that the threat to the ground personnel and facilities caused by air crash accidents is reduced. According to the method and the system for generating the dynamic risk map of the unmanned aerial vehicle operation on the ground, the problems of insufficient multi-factor fusion, insufficient static data precision, poor dynamic adaptability and the like in an existing unmanned aerial vehicle risk assessment method can be efficiently solved.

Description

Technical Field

[0001] The present invention proposes a method and system for generating a dynamic risk map of UAV operation on the ground, relating to the field of UAV risk analysis. Background Art

[0002] Current drone flight risk assessment methods primarily rely on a single factor (such as population density or obstacle distribution) and lack the systematic integration of multi-source data, resulting in insufficient risk assessment accuracy. Traditional methods have the following problems:

[0003] (1) Qualitative methods lack quantitative depth: Existing frameworks (such as SORA) mainly provide qualitative analysis frameworks and lack quantitative research on key indicators such as risk probability and impact range, making it difficult to meet the needs of accurate safety assessment for large-scale applications of low-altitude economy.

[0004] (2) Insufficient accuracy of static data: Existing risk map generation methods are mostly based on simple matrices or low-resolution data (such as coarse-grained demographics), without combining high-precision geographic information layers, and are unable to dynamically reflect the spatiotemporal characteristics of the urban environment (such as changes in traffic conditions, wind speed and direction, and weather conditions).

[0005] (3) Insufficient integration of multiple factors: Most studies focus only on a single risk source (such as population density) and do not comprehensively consider the multidimensional characteristics of urban low-altitude scenes (such as no-fly zones and obstacle distribution). Summary of the Invention

[0006] In view of this, in order to fill the gaps and shortcomings of the existing technology, the present invention provides a method and system for generating a dynamic risk map of UAV operation on the ground, which is used to solve the problems of insufficient multi-factor fusion, insufficient static data accuracy, and poor dynamic adaptability in existing UAV risk assessment methods.

[0007] The present invention proposes a method and system for generating a dynamic risk map of UAV operations on the ground, which provides a detailed risk assessment for UAV flights, identifies potential risks in the flight area, avoids high-risk areas, improves flight safety, provides a scientific basis for low-altitude flight path planning, and reduces damage to ground personnel and facilities caused by UAV crashes or other accidents, thereby reducing the threat posed to ground personnel and facilities by crashes.

[0008] The present invention proposes a method and system for generating a dynamic risk map of the ground caused by UAV operation, including the following contents:

[0009] A method for generating a dynamic risk map of a UAV operating on the ground includes the following contents:

[0010] Step S1: Collect data and rasterize the collected data to form a data layer;

[0011] Step S2: constructing a ground risk assessment model, including calculating a population casualty risk value matrix based on the UAV failure and falling motion model;

[0012] Step S3: Fusion of multiple data layers to obtain a risk level matrix;

[0013] Step S4: Visualize the risk level matrix grid.

[0014] Furthermore, step S1 includes the following contents:

[0015] Step S11: collecting data including: collecting meteorological data, collecting population data, collecting building outline data and tree cover data, collecting traffic data, collecting no-fly zone data, collecting obstacle data, and collecting special area data;

[0016] Among them, the collection of meteorological data includes collecting precipitation, temperature, humidity, weather warnings, and real-time wind speed and direction data in the target area;

[0017] Wherein, collecting population data includes collecting population distribution data of the target area;

[0018] The collection of building outline data and tree cover data includes collecting the distribution of buildings and trees in the target area;

[0019] The collection of traffic condition data includes collecting the congestion level and road coordinates of the target area;

[0020] The collection of no-fly zone data includes collecting the no-fly zone range and the non-no-fly zone range of the target area;

[0021] Obstacle data collection includes collecting building heights and mountain terrain elevation data in the target area;

[0022] The collection of special area data includes collecting the distribution of special areas in the target area that are prone to causing secondary accidents, wherein the special areas that are prone to causing secondary accidents include but are not limited to gas stations, schools and hospitals.

[0023] Furthermore, step S1 also includes the following:

[0024] Step S12: Processing the data collected in step S11, including the following:

[0025] The collected meteorological data is rasterized, and corresponding risk levels are set according to the temperature, humidity, wind speed, precipitation and weather conditions of the grid cells to obtain a meteorological risk matrix;

[0026] The collected population distribution grid data is first corrected, where the population correction coefficient Ci The formula is as follows:

[0027]

[0028] Where i represents the target area, P i W represents the total population of the district or county where the target area is located, as obtained from the latest survey. i Indicates the total population of the district and county where the target area is located in the collected population distribution data;

[0029] The calculated population correction coefficient is multiplied by the population distribution grid data to obtain the corrected population distribution grid data, and then the population grid data is clipped to obtain the population distribution grid data of the target area, and the population distribution grid data of the target area is read to obtain the population data matrix;

[0030] Processing the collected building outline data and tree cover data includes rasterizing the building outline data and the tree cover data, clipping the building outline and tree cover raster data to obtain the building outline and tree cover raster data of the target area, and further calculating the percentage of the sum of the building cover area and the tree cover area in each grid cell to the total grid cell to obtain the shielding coefficient; and forming all the shielding coefficients into a shielding coefficient matrix;

[0031] The collected traffic condition data is rasterized, and the population correction coefficient of the grid unit is set according to the congestion level to generate a traffic congestion correction matrix;

[0032] The collected no-fly zone data is rasterized, and then the no-fly zone raster data is clipped to obtain the no-fly zone raster data of the target area, including setting the risk level of the rasterized no-fly zone cells in the target area to 4 and setting the risk level of the remaining rasterized non-no-fly zone cells to 0; further, a no-fly zone matrix is ​​generated by combining the no-fly zone cells and the non-no-fly zone cells;

[0033] The collected obstacle data is rasterized, including building heights and mountain obstacles. The rasterized building height and mountain obstacle data are then clipped to obtain the obstacle height data of the target area. If the UAV's flight altitude minus the obstacle height is less than 30 meters, the risk level of the grid unit is determined to be 4, and the obstacle matrix is ​​obtained.

[0034] The collected special area data is rasterized, including risk level marking for the special areas. The risk level of the special area is set according to the type of special area. If a grid cell contains a special area that is prone to secondary accidents, the risk level of the grid cell is set according to the highest risk level of the special area within the grid cell, resulting in a special area matrix.

[0035] Furthermore, step S2 includes the following contents:

[0036] Step S21: Based on the drone's crash motion model and aerodynamic model, the vertical and horizontal motion trajectories are simulated based on gravity, air resistance, and real-time wind speed and direction when the drone falls. The NED coordinate system is used to create the xyz axes, dividing the drone's falling motion into three directions. The x-axis points to the east, the y-axis points to the north, and the z-axis points to the ground. The expressions for the three directions of the drone's falling motion are as follows:

[0037]

[0038] Where m is the mass of the UAV; a x 、a y and a z The accelerations in three directions respectively; D x 、D y and D z They are the air resistance of the drone in three two different directions.

[0039] The formula for calculating air resistance is as follows:

[0040]

[0041] Among them C D is the air resistance coefficient; ρ a is the air density; v u is the flight speed of the UAV; A u is the windward area of ​​the UAV (i.e. the projected area perpendicular to the direction of movement).

[0042] Step S22: Without considering the wind force, the time t when the drone falls to the ground fall The calculation formula is as follows:

[0043]

[0044] Where m is the mass of the UAV; h is the initial flight altitude of the UAV; and g is the acceleration due to gravity.

[0045] x is the displacement of the UAV in the x-axis direction under the action of aerodynamic force. The expression of displacement x with respect to time t is as follows:

[0046]

[0047] where v x0 is the component of the initial velocity of the UAV in the x-axis direction; t fall Substituting the expression into the equation, we can calculate the displacement x;

[0048] Step S23: In the wind motion model, the motion of the UAV is only affected by wind force and can be modeled as:

[0049]

[0050] Where: x wind is the displacement of the UAV in the x-axis direction under the action of wind; v xwind is the component of wind speed in the x-axis direction; t fall is the time when the drone falls to the ground; therefore, the actual horizontal movement distance x of the drone during the failure descent phase d It is x and x wind The comprehensive manifestation of the two motion states, x d It can be expressed as:

[0051] x d =x+x wind

[0052] Similarly, the calculation method for the displacement of the UAV in the y-axis direction under the action of wind is the same as that of the UAV, and the y d =y+y wind .

[0053] Furthermore, step S2 includes the following contents:

[0054] Step S24: By simulating the falling process of the drone, a large number of random falling trajectories can be generated. Based on the distribution of the impact points of the drones in the area, a two-dimensional probability density function is generated using a two-dimensional Gaussian distribution;

[0055] The collision probability is calculated based on the two-dimensional probability density distribution of the impact point. The calculation formula is:

[0056]

[0057] Where: P collision is the probability that the drone crashes on the ground and collides with at least one person; PDF is a two-dimensional probability density function that describes the probability distribution of the location of the drone's impact point; ρ(x,y) is the population density at the ground geographic reference location (x,y); t(x,y) is the population density correction factor at the ground geographic reference location (x,y); A imp A is the impact area of ​​the drone hitting the ground, imp The formula includes the following:

[0058] A imp =2(r p +R UAV )h p cotθ+π(r p +R UAV ) 2

[0059] Where: r p and R UAV are the radius of the person being hit and the maximum radius of the falling drone respectively; h p is the average height of the people hit on the ground; θ is the falling angle of the drone, which is calculated as follows:

[0060]

[0061] Where: v u is the flight speed of the drone; C D is the air resistance coefficient; A u is the frontal area of ​​the UAV; ρ a is the air density.

[0062] Furthermore, step S2 also includes the following:

[0063] Step S25: The probability of ground personnel being hit and injured is related to the shielding effect of other obstacles such as trees and buildings on the ground and the kinetic energy of the drone when it hits the ground. The probability of ground personnel being hit and injured includes the following:

[0064]

[0065] Where: P fatality (x, y) is the probability of casualties at the geographical reference location (x, y), and (x, y) is within the influence range of the PDF of the drone hitting the ground; E imp is the kinetic energy of the drone when it hits the ground; P s (x, y) is the shielding coefficient corresponding to the area where (x, y) is located; α is the impact energy required to achieve a mortality rate of 50% when the shielding coefficient is 0.5; β is the impact energy required when P s The impact energy threshold required for death to occur when it approaches 0; E imp The calculation formula is as follows:

[0066]

[0067] Where: v imp Represents the speed of the drone when it hits the ground; v imp By analyzing the descent process of the drone, it was concluded that the specific calculation involves parameters such as air resistance and acceleration.

[0068] The formula for calculating air resistance is:

[0069]

[0070] The calculation formula of the UAV's acceleration a is:

[0071]

[0072] Where: G represents gravity, which can be expressed as G = mg. In the above formula, the speed v when the drone hits the ground imp The calculation formula is:

[0073]

[0074] Where h is the initial flight altitude of the UAV, the calculated v imp Substitute into the formula The impact kinetic energy of the drone when it hits the ground can be obtained;

[0075] Step S26: The calculation formula of the population casualty risk value matrix is ​​as follows:

[0076]

[0077] Where: the value of the element in R represents the probability of casualty, (x, y) is the ground geographic reference position; P event is the accident rate of drone crashing.

[0078] Furthermore, step S3 includes the following contents:

[0079] Step S31: Generate a comprehensive ground risk map by layer-by-layer superposition and the maximum risk value principle. The specific steps are as follows: First, read the data of the population casualty risk layer, special area layer, no-fly zone layer, obstacle layer, and meteorological risk layer, and unify the risk level matrix dimensions of all layers; then, using the maximum risk value principle, for the data stored in each grid cell in the map, select the maximum risk level among all layers to obtain the risk level matrix. The calculation formula is as follows:

[0080] RL(x,y)=max{F(x,y),O(x,y),S(x,y),R(x,y),W(x,y)}

[0081] Where RL is the risk level matrix; F, O, S, R, and W are the risk level matrices for the no-fly zone, obstacle layer, special area layer, population casualty risk layer, and meteorological risk layer, respectively. The specific meanings of the risk levels are: 0 indicates safety, 1 indicates low risk, 2 indicates medium risk, 3 indicates high risk, and 4 indicates that flying is prohibited in the area due to no-fly regulations, high obstacles, or weather conditions.

[0082] Furthermore, step S4 includes the following contents:

[0083] Step S41: Grid visualization of the risk level matrix, including the following contents:

[0084] First, verify the matrix dimensions, use different colors to map risk levels, adjust the coordinate axis offset to accurately align the grid lines with the matrix cells, and add black grid lines to strengthen the cell boundaries; finally, output a geographic heat map with a color bar.

[0085] According to a second aspect of the present invention, a system for generating a dynamic risk map of an unmanned aerial vehicle (UAV) operating on the ground comprises an electronic device, wherein the electronic device comprises a memory, a processor, and a computer program stored in the memory and executable on the processor, and is characterized in that when the processor executes the computer program, a method for generating a dynamic risk map of an unmanned aerial vehicle (UAV) operating on the ground is implemented as described in any one of the present inventions.

[0086] According to a third aspect of the present invention, a system for generating a dynamic risk map of an unmanned aerial vehicle operating on the ground includes a computer-readable storage medium storing a computer program, and is characterized in that when the computer program is executed by a processor, it implements a method for generating a dynamic risk map of an unmanned aerial vehicle operating on the ground as described in any one of the present inventions.

[0087] The present invention has the following advantages:

[0088] The method and system for generating a dynamic risk map of drone operations on the ground proposed in the present invention can use a multi-factor, multi-level quantitative model to consider both static factors of drone failure and crash events, such as obstacles and special areas, as well as dynamic factors such as real-time meteorological conditions, wind direction and speed, and traffic conditions. This can accurately quantify ground risks and reduce assessment errors. BRIEF DESCRIPTION OF THE DRAWINGS

[0089] Figure 1 It is a schematic diagram of the method and system for generating a dynamic risk map of the ground caused by the operation of a UAV in the present invention.

[0090] Figure 2 It is a two-dimensional probability density distribution diagram in the present invention.

[0091] Figure 3 Schematic diagram of the impact area when the UAV hits the ground in the present invention.

[0092] Figure 4 It is a schematic diagram of the research area in the present invention.

[0093] Figure 5 It is a schematic diagram of the obstacle layer in the present invention.

[0094] Figure 6 It is a schematic diagram of mapping different colors to different risk levels in the present invention. DETAILED DESCRIPTION

[0095] The technical solution of the present invention will be described in detail below with reference to the accompanying drawings.

[0096] It should be noted that the following detailed descriptions are illustrative and intended to provide further explanation of the present application. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by those skilled in the art to which the present application belongs.

[0097] It should be noted that the terms used herein are only for describing specific embodiments and are not intended to limit the exemplary embodiments according to the present application; as used herein, unless the context clearly indicates otherwise, the singular form is also intended to include the plural form, and it should be understood that when the terms "comprise" and / or "include" are used in this specification, they indicate the presence of features, steps, operations, devices, components and / or their combinations.

[0098] The present invention proposes a method and system for generating a dynamic risk map of the ground caused by UAV operation, including the following contents:

[0099] In one embodiment of the present invention, a method for generating a dynamic risk map of a drone operating on the ground includes the following:

[0100] Step S1: Collect data and rasterize the collected data to form a data layer;

[0101] Step S2: constructing a ground risk assessment model, including calculating a population casualty risk value matrix based on the UAV failure and falling motion model;

[0102] Step S3: Fusion of multiple data layers to obtain a risk level matrix;

[0103] Step S4: Visualize the risk level matrix grid.

[0104] In one embodiment of the present invention, step S1 includes the following contents:

[0105] Step S11: collecting data including: collecting wind data, collecting population data, collecting building vector data and tree cover data, collecting traffic data, collecting no-fly zone data, collecting obstacle data, and collecting special area data;

[0106] Among them, collecting wind data includes collecting real-time wind speed and wind direction data in the target area;

[0107] Among them, collecting population data includes collecting population distribution in the target area;

[0108] The collection of building quantity data and tree quantity data includes collecting the distribution of buildings and trees in the target area;

[0109] The collection of traffic condition data includes collecting the congestion level and road coordinates of the target area;

[0110] The collection of no-fly zone data includes collecting the no-fly zone range and the non-no-fly zone range of the target area;

[0111] Obstacle data collection includes collecting building heights and mountain terrain elevation data in the target area;

[0112] The collection of special area data includes collecting the distribution of special areas in the target area that are prone to causing secondary accidents, wherein the special areas that are prone to causing secondary accidents include but are not limited to gas stations, schools and hospitals.

[0113] In one embodiment of the present invention, step S1 further includes the following:

[0114] Step S12: Processing the data collected in step S11, including the following:

[0115] The collected meteorological data is rasterized, and the rasterized wind data is clipped according to the temperature, humidity, wind speed, precipitation and weather conditions of the grid cells to obtain a wind data matrix;

[0116] The collected population distribution grid data is first corrected, where the population correction coefficient C i The formula is as follows:

[0117]

[0118] Where i represents the target area, P i W represents the total population of the district or county where the target area is located, as obtained from the latest survey. i Indicates the total population of the district and county where the target area is located in the collected population distribution data;

[0119] The calculated population correction coefficient is multiplied by the population distribution raster data to obtain the corrected population distribution raster data, and then the population distribution raster data is clipped to obtain the population distribution raster data of the target area, and the population distribution raster data of the target area is read to obtain the population data matrix;

[0120] Processing the collected building outline data and tree cover data includes rasterizing the building outline data and the tree cover data, clipping the building outline and tree cover raster data to obtain the building outline and tree cover raster data of the target area, and further calculating the percentage of the sum of the building cover area and the tree cover area in each grid cell to the total grid cell to obtain the shielding coefficient; and forming all the shielding coefficients into a shielding coefficient matrix;

[0121] The collected traffic condition data is rasterized, and the population correction coefficient of the grid unit is set according to the congestion level to generate a traffic congestion correction matrix;

[0122] The collected no-fly zone data is rasterized, and then the no-fly zone raster data is clipped to obtain the no-fly zone raster data of the target area, including setting the risk level of the rasterized no-fly zone cells in the target area to 4 and setting the risk level of the remaining rasterized non-no-fly zone cells to 0; further, a no-fly zone matrix is ​​generated by combining the no-fly zone cells and the non-no-fly zone cells;

[0123] The collected obstacle data is rasterized, including building heights and mountain obstacles. The rasterized building height and mountain obstacle data are then clipped to obtain the obstacle height data of the target area. If the UAV's flight altitude minus the obstacle height is less than 30 meters, the risk level of the grid unit is determined to be 4, and the obstacle matrix is ​​obtained.

[0124] The collected special area data is rasterized, including risk level marking for the special areas. The risk level of the special area is set according to the type of special area. If a grid cell contains a special area that is prone to secondary accidents, the risk level of the grid cell is set according to the highest risk level of the special area within the grid cell, resulting in a special area matrix.

[0125] In one embodiment of the present invention, step S2 includes the following contents:

[0126] Step S21: Based on the drone's crash motion model and aerodynamic model, the vertical and horizontal motion trajectories are simulated based on gravity, air resistance, and real-time wind speed and direction when the drone falls. The NED coordinate system is used to create the xyz axes, dividing the drone's falling motion into three directions. The x-axis points to the east, the y-axis points to the north, and the z-axis points to the ground. The expressions for the three directions of the drone's falling motion are as follows:

[0127]

[0128] Where m is the mass of the UAV; a x 、a y and a z The accelerations in three directions respectively; D x 、D y and D z They are the air resistance of the drone in three two different directions.

[0129] The formula for calculating air resistance is as follows:

[0130]

[0131] Among them C D is the air resistance coefficient; ρ a is the air density; v u is the flight speed of the UAV; A u is the windward area of ​​the UAV (i.e. the projected area perpendicular to the direction of movement).

[0132] Step S22: Without considering the wind force, the time t when the drone falls to the ground fall The calculation formula is as follows:

[0133]

[0134] Where m is the mass of the UAV; h is the initial flight altitude of the UAV; and g is the acceleration due to gravity.

[0135] x is the displacement of the UAV in the x-axis direction under the action of aerodynamic force. The expression of displacement x with respect to time t is as follows:

[0136]

[0137] where v x0 is the component of the initial velocity of the UAV in the x-axis direction; t fall Substituting the expression into the equation, we can calculate the displacement x;

[0138] Step S23: In the wind motion model, the motion of the UAV is only affected by wind force and can be modeled as:

[0139]

[0140] Where: x wind is the displacement of the UAV in the x-axis direction under the action of wind; v xwind is the component of wind speed in the x-axis direction; t fall is the time when the drone falls to the ground; therefore, the actual horizontal movement distance x of the drone during the failure descent phase d It is x and x wind The comprehensive manifestation of the two motion states, x d It can be expressed as:

[0141] x d =x+x wind

[0142] Similarly, the calculation method for the displacement of the drone in the y-axis direction under the action of wind is the same.

[0143] In one embodiment of the present invention, step S2 includes the following contents:

[0144] Step S24: By simulating the falling process of the drone, a large number of random falling trajectories can be generated. Based on the distribution of the impact points of the drones in the area, a two-dimensional probability density function is generated using a two-dimensional Gaussian distribution;

[0145] The collision probability is calculated based on the two-dimensional probability density distribution of the impact point. The calculation formula is:

[0146]

[0147] Where: P collision is the probability that the drone crashes on the ground and collides with at least one person; PDF is a two-dimensional probability density function that describes the probability distribution of the location of the drone's impact point; ρ(x,y) is the population density at the ground geographic reference location (x,y); t(x,y) is the population density correction factor at the ground geographic reference location (x,y); A imp A is the impact area of ​​the drone hitting the ground, imp The formula includes the following:

[0148] A imp =2(r p +R UAV )h p cotθ+π(r p +R UAV ) 2

[0149] Where: r p and R UAV are the radius of the person being hit and the maximum radius of the falling drone respectively; h p is the average height of the people hit on the ground; θ is the falling angle of the drone, which is calculated as follows:

[0150]

[0151] Where: v u is the flight speed of the drone; C D is the air resistance coefficient; A u is the frontal area of ​​the UAV; ρ a is the air density.

[0152] In one embodiment of the present invention, step S2 further includes the following:

[0153] Step S25: The probability of ground personnel being hit and injured is related to the shielding effect of other obstacles such as trees and buildings on the ground and the kinetic energy of the drone when it hits the ground. The probability of ground personnel being hit and injured includes the following:

[0154]

[0155] Where: P fatality (x, y) is the probability of casualties at the geographical reference location (x, y), and (x, y) is within the influence range of the PDF of the drone hitting the ground; E imp is the kinetic energy of the drone when it hits the ground; P s (x, y) is the shielding coefficient corresponding to the area where (x, y) is located; α is the impact energy required to achieve a mortality rate of 50% when the shielding coefficient is 0.5; β is the impact energy required when P s The impact energy threshold required for death to occur when it approaches 0; E imp The calculation formula is as follows:

[0156]

[0157] Where: v imp Represents the speed of the drone when it hits the ground; v imp By analyzing the descent process of the drone, it was concluded that the specific calculation involves parameters such as air resistance and acceleration.

[0158] The formula for calculating air resistance is

[0159]

[0160] The calculation formula of the UAV's acceleration a is:

[0161]

[0162] Where: G represents gravity, which can be expressed as G = mg. In the above formula, the speed v when the drone hits the ground imp The calculation formula is

[0163]

[0164] Where h is the initial flight altitude of the UAV, the calculated v imp Substitute into the formula The impact kinetic energy of the drone when it hits the ground can be obtained;

[0165] Step S26: The calculation formula of the population casualty risk value matrix is ​​as follows:

[0166]

[0167] Where: the value of the element in R represents the probability of casualty, (x, y) is the ground geographic reference position; P event is the accident rate of drone crashing.

[0168] In one embodiment of the present invention, step S3 includes the following contents:

[0169] Step S31: Generate a comprehensive ground risk map by layer-by-layer superposition and the maximum risk value principle. The specific steps are as follows: First, read the data of the population casualty risk layer, special area layer, no-fly zone layer, obstacle layer, and meteorological risk layer, and unify the risk level matrix dimensions of all layers; then, using the maximum risk value principle, for the data stored in each grid cell in the map, select the maximum risk level among all layers to obtain the risk level matrix. The calculation formula is as follows:

[0170] RL(x,y)=max{F(x,y),O(x,y),S(x,y),R(x,y),W(x,y)}

[0171] Where RL is the risk level matrix; F, O, S, R, and W are the risk level matrices for the no-fly zone, obstacle layer, special area layer, population casualty risk layer, and meteorological risk layer, respectively. The specific meanings of the risk levels are: 0 indicates safety, 1 indicates low risk, 2 indicates medium risk, 3 indicates high risk, and 4 indicates that flying is prohibited in the area due to no-fly regulations, high obstacles, or weather conditions.

[0172] In one embodiment of the present invention, step S4 includes the following contents:

[0173] Step S41: Grid visualization of the risk level matrix, including the following contents:

[0174] First, verify the matrix dimensions, use different colors to map risk levels, adjust the coordinate axis offset to accurately align the grid lines with the matrix cells, and add black grid lines to strengthen the cell boundaries; finally, output a geographic heat map with a color bar.

[0175] In one embodiment of the present invention, according to the second aspect of the present invention, a system for generating a dynamic risk map of an unmanned aerial vehicle (UAV) operating on the ground comprises an electronic device, wherein the electronic device comprises a memory, a processor, and a computer program stored in the memory and executable on the processor, and is characterized in that when the processor executes the computer program, a method for generating a dynamic risk map of an unmanned aerial vehicle (UAV) operating on the ground as described in any one of the present inventions is implemented.

[0176] In one embodiment of the present invention, according to the third aspect of the present invention, a system for generating a dynamic risk map of an unmanned aerial vehicle operating on the ground includes a computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and is characterized in that when the computer program is executed by a processor, it implements a method for generating a dynamic risk map of an unmanned aerial vehicle operating on the ground as described in any one of the present inventions.

[0177] In addition to the above embodiments, the present invention also includes the following embodiments:

[0178] In one embodiment of the present invention, collecting wind data in step S1 includes the following:

[0179] Using the weather query API, obtain precipitation, temperature, humidity, weather warnings, and real-time wind speed and direction data for the target city. When precipitation occurs in the region, regional authorities issue a weather disaster warning, or a meteorological factor does not meet the set safety threshold, the risk level for all grids in the meteorological risk layer will be uniformly set to Level 4 (No Flight). If all weather data meets the set threshold range, the risk level for all grids in the meteorological risk layer will be set to 0. The set safety threshold range is shown in Table 1 (the safety threshold range can be adjusted accordingly based on the anti-interference capabilities of different drone models).

[0180] factor Safety threshold temperature 10℃～40℃ wind speed ≤ Level 5 wind (8-10.7m / s) humidity 30％～80％

[0181] Table 1 Safety threshold range of meteorological factors

[0182] In one embodiment of the present invention, the population data collected in step S1 includes the following:

[0183] The population density layer is based on WorldPop's population grid data for a certain year and the Nth census data of a certain country. The raster calculation tools of GIS software are used to count the population data of each district and county. The Nth national census data of a certain country are compared with the WorldPop data to calculate the population correction coefficient for each target area. i The calculation formula is as follows:

[0184]

[0185] Where: i represents a target area, P i represents the total population of target area i in the Nth national census data of a country, W i Represents the total population of target area i in the (N-1)th time period in the WorldPop population raster data for a particular year. Convert the calculated population correction coefficient to a raster format and multiply it with the WorldPop population raster data to obtain a more accurate corrected population distribution raster. Use the raster clipping tool in the GIS software to clip the population raster data for the target area. Then read the raster data and convert it into a matrix with the same dimensions as the risk map.

[0186] In one embodiment of the present invention, collecting building outline data and tree cover data in step S1 includes the following:

[0187] The shielding coefficient layer is based on the collected building outline data and tree cover data. The raster clipping tool of the GIS software is used to clip the building outline data and tree cover data of the target area. The two sets of data are then analyzed to calculate the proportion of building and tree cover area in each raster cell. The shielding coefficient of the raster cell is determined by referring to Table 2, and a shielding coefficient matrix with the same dimension as the risk map is generated.

[0188] Shielding factor Proportion of area covered by obstructions (buildings, trees) 0 <20％ 0.25 20% ≤ and < 40% 0.5 40% ≤ and < 60% 0.75 60% ≤ and < 80% 1 ≥80％

[0189] Table 2 Shielding coefficient reference table

[0190] In one embodiment of the present invention, collecting traffic condition data in step S1 includes the following:

[0191] The traffic condition layer retrieves real-time traffic condition data for a specified target area using a platform's rectangular area traffic situation query API. Each road condition data returned by the API includes the following: status: congestion level (0-3, corresponding to unknown, unblocked, slow, and congested); polyline: road coordinate point. The original rectangular area is split into grids with the same dimensions as the risk map. The grid in which the road coordinate point resides is determined based on the road coordinate point. The traffic condition level for each grid is recorded (if different roads have different levels, the highest value is taken). This is then mapped to a population density correction factor according to Table 3, generating a correction factor matrix with the same dimensions as the risk map.

[0192] Traffic condition level Population density correction factor 0: Unknown 1 1: Smooth 0.8 2: Slow down 1.2 3: Congestion 1.5

[0193] Table 3 Population density correction coefficient reference table

[0194] In one embodiment of the present invention, collecting no-fly zone data in step S1 includes the following:

[0195] The no-fly zone layer is based on the no-fly zone data for a specific drone system. GIS software is used to map the no-fly zones for drones and convert them into raster data. Using the GIS software's raster clipping tools, the no-fly zone data for the target area is clipped. The risk level for the no-fly zones is set to 4 (no flight allowed), while the remaining non-no-fly zones are set to 0 (safe), generating a matrix with the same dimensions as the risk map.

[0196] In one embodiment of the present invention, collecting obstacle data in step S1 includes the following:

[0197] The obstacle layer is based on building height data from the Geographic Remote Sensing Ecological Network and Copernicus terrain elevation data. Using the raster clipping tool in GIS software, we clipped the target area's data. If the height of an obstacle (such as a building or mountain) is greater than or equal to the drone's flight altitude, the risk level for that area is set to 4 (no flight allowed); otherwise, it is set to 0 (safe). This results in a matrix with the same dimensionality as the risk map.

[0198] In one embodiment of the present invention, collecting special area data in step S1 includes the following:

[0199] The special area layer uses the AutoNavi Map API to obtain POI data such as gas stations and hospitals. Table 4, defined in accordance with civil aviation safety regulations, is used to define risk levels for special areas. For example, the risk level for gas stations is 2 (medium risk). If multiple special areas exist within a grid unit (100m×100m), the highest risk level among these areas is used to create the special area layer.

[0200]

[0201] Table 4 Special area risk level reference table

[0202] In one embodiment of the present invention, the two-dimensional probability density distribution in step S2 is as follows: Figure 2 As shown, its meaning includes:

[0203] Figure 2 This model, based on a drone crash event, takes into account the uncertainties in aerodynamics (gravity, drag) and wind parameters (wind speed and direction). A two-dimensional probability density function is used to quantify the spatial probability distribution of a drone crash. For example, the red areas in the figure indicate the highest probability of a drone crash, while the dark blue areas indicate zero probability of a drone impact.

[0204] In one embodiment of the present invention, the impact area of ​​the drone hitting the ground in step S2 is as follows: Figure 3 shown.

[0205] In one embodiment of the present invention, the grid visualization of the risk level matrix in step S4 includes the following:

[0206] An area of ​​2km×2km in a certain region is selected as the analysis scene. The study area is divided into grids, such as Figure 4 As shown, a 20×20 grid is formed, and the side length of each grid unit is 100m.

[0207] Taking the obstacle layer as an example, when the drone's flight altitude is 60 meters, the risk level matrix of the obstacle layer is visualized as follows: Figure 5 shown.

[0208] In one embodiment of the present invention, the use of different colors to map risk levels in step S4 includes the following:

[0209] First, verify the matrix dimension and use different colors to map the risk levels, where black corresponds to risk level 4, purple corresponds to risk level 3, red corresponds to risk level 2, pink corresponds to risk level 1, and white corresponds to risk level 0. The final result is as follows Figure 6 shown.

[0210] The above are preferred embodiments of the present invention. Any changes made according to the technical solution of the present invention, as long as the resulting functions and effects do not exceed the scope of the technical solution of the present invention, shall fall within the scope of protection of the present invention.

Claims

1. A method for generating a dynamic risk map of UAV operation to the ground, characterized in that: The following steps are involved: Step S1: Collect data and rasterize the collected data to form a data layer; Step S2: constructing a ground risk assessment model, including calculating a population casualty risk value matrix based on the UAV failure and falling motion model; Step S3: Fusion of multiple data layers to obtain a risk level matrix; Step S4: Visualize the risk level matrix grid.

2. The method for generating a dynamic risk map of a UAV operating on the ground according to claim 1, characterized in that: Step S1 includes the following contents: Step S11: collecting data including: collecting meteorological data, collecting population data, collecting building vector data and tree cover data, collecting traffic data, collecting no-fly zone data, collecting obstacle vector data, and collecting special area data; Among them, the collection of meteorological data includes collecting precipitation, temperature, humidity, weather warnings, and real-time wind speed and direction data in the target area; Wherein, collecting population data includes collecting population distribution data of the target area; The collection of building outline data and tree cover data includes collecting the distribution of buildings and trees in the target area; The collection of traffic condition data includes collecting the congestion level and road coordinates of the target area; The collection of no-fly zone data includes collecting the no-fly zone range and the non-no-fly zone range of the target area; Obstacle data collection includes collecting building heights and mountain terrain elevation data in the target area; The collection of special area data includes collecting the distribution of special areas in the target area that are prone to causing secondary accidents, wherein the special areas that are prone to causing secondary accidents include gas stations, schools and hospitals.

3. The method for generating a dynamic risk map of a UAV operating on the ground according to claim 2, characterized in that: Step S1 also includes the following: Step S12: Processing the data collected in step S11, including the following: The collected meteorological data is rasterized, and corresponding risk levels are set according to the temperature, humidity, wind speed, precipitation and weather conditions of the grid cells to obtain a meteorological risk matrix; The collected population distribution grid data is first corrected, where the population correction coefficient C i The formula is as follows: Where i represents the target area, P i W represents the total population of the district or county where the target area is located according to the latest survey. i Indicates the total population of the district and county where the target area is located in the collected population distribution data; The calculated population correction coefficient is multiplied by the population distribution raster data to obtain the corrected population distribution raster data, and then the population distribution raster data is clipped to obtain the population distribution raster data of the target area, and the population distribution raster data of the target area is read to obtain the population data matrix; Processing the collected building outline data and tree cover data includes rasterizing the building outline data and the tree cover data, clipping the building outline and tree cover raster data to obtain the building outline and tree cover raster data of the target area, and further calculating the percentage of the sum of the building cover area and the tree cover area in each grid cell to the total grid cell to obtain the shielding coefficient; and forming all the shielding coefficients into a shielding coefficient matrix; The collected traffic condition data is rasterized, and the population correction coefficient of the grid unit is set according to the congestion level to generate a traffic congestion correction matrix; The collected no-fly zone data is rasterized, and then the no-fly zone raster data is clipped to obtain the no-fly zone raster data of the target area, including setting the risk level of the rasterized no-fly zone cells in the target area to 4 and setting the risk level of the remaining rasterized non-no-fly zone cells to 0; further, a no-fly zone matrix is ​​generated by combining the no-fly zone cells and the non-no-fly zone cells; The collected obstacle data is rasterized, including building heights and mountain obstacles. The rasterized building height and mountain obstacle data are then clipped to obtain the obstacle height data of the target area. If the UAV's flight altitude minus the obstacle height is less than 30 meters, the risk level of the grid unit is determined to be 4, and the obstacle matrix is ​​obtained. The collected special area data are rasterized, including risk level marking of special areas; the risk level of the special area is set according to the type of special area; when there is a special area in the grid unit that is prone to secondary accidents, the risk level of the grid unit is set according to the highest risk level of the special area in the grid unit to obtain the special area matrix.

4. The method for generating a dynamic risk map of a UAV operating on the ground according to claim 3, characterized in that: Step S2 includes the following contents: Step S21: Based on the drone's crash motion model and aerodynamic model, the vertical and horizontal motion trajectories are simulated based on gravity, air resistance, and real-time wind speed and direction when the drone falls. The NED coordinate system is used to create the xyz axes, dividing the drone's falling motion into three directions. The x-axis points to the east, the y-axis points to the north, and the z-axis points to the ground. The expressions for the three directions of the drone's falling motion are as follows: Where m is the mass of the UAV; a x 、a y and a z The accelerations in three directions respectively; D x 、D y and D z are the air resistance of the drone in three different directions; The formula for calculating air resistance is as follows: Among them C D is the air resistance coefficient; ρ a is the air density; v u is the flight speed of the UAV; A u is the windward area of ​​the UAV (i.e. the projected area perpendicular to the direction of movement); Step S22: Without considering the wind force, the time t when the drone falls to the ground fall The calculation formula is as follows: Where m is the mass of the drone; h is the initial flight altitude of the drone; g is the acceleration due to gravity; x is the displacement of the UAV in the x-axis direction under the action of aerodynamic force. The expression of displacement x with respect to time t is as follows: where v x0 is the component of the initial velocity of the UAV in the x-axis direction; t fall Substituting the expression x(t), the displacement x can be calculated; Step S23: In the wind motion model, the motion of the UAV is only affected by wind force and can be modeled as: Where: x wind is the displacement of the UAV in the x-axis direction under the action of wind; v xwind is the component of wind speed in the x-axis direction; t fall is the time when the drone falls to the ground; further, the actual horizontal movement distance x of the drone during the failure descent phase d It is x and x wind The comprehensive expression of the two motion states includes the following formula: x d =x+x wind Similarly, the calculation method for the displacement of the UAV in the y-axis direction under the action of wind is the same as that of the UAV, and the y d =y+y wind .

5. The method for generating a dynamic risk map of a UAV operating on the ground according to claim 4, characterized in that: Step S2 includes the following contents: Step S24: By simulating the falling process of the drone, a large number of random falling trajectories can be generated; according to the distribution of the impact points of the drones in the area, a two-dimensional probability density function is generated using a two-dimensional Gaussian distribution; The collision probability is calculated based on the two-dimensional probability density distribution of the impact point. The calculation formula is: Where: P collision is the probability that the drone crashes on the ground and collides with at least one person; PDF is a two-dimensional probability density function that describes the probability distribution of the location of the drone's impact point; ρ(x,y) is the population density at the ground geographic reference location (x,y); t(x,y) is the population density correction factor at the ground geographic reference location (x,y); A imp A is the impact area of ​​the drone hitting the ground, imp The formula includes the following: A imp =2(r p +R UAV )h p cotθ+π(r p +R UAV ) 2 Where: r p and R UAV are the radius of the person being hit and the maximum radius of the falling drone respectively; h p is the average height of the people hit on the ground; θ is the falling angle of the drone, which is calculated as follows: Where: v u is the flight speed of the drone; C D is the air resistance coefficient; A u is the frontal area of ​​the UAV; ρ a is the air density.

6. The method for generating a dynamic risk map of a UAV operating on the ground according to claim 5, characterized in that: Step S2 also includes the following: Step S25: The probability of ground personnel being hit and injured is related to the shielding effect of other obstacles such as trees and buildings on the ground and the kinetic energy of the drone when it hits the ground. The probability of ground personnel being hit and injured includes the following: Where: P fatality (x, y) is the probability of casualties at the geographical reference location (x, y), and (x, y) is within the influence range of the PDF of the drone hitting the ground; E imp is the kinetic energy of the drone when it hits the ground; P s (x, y) is the shielding coefficient corresponding to the area where (x, y) is located; α is the impact energy required to achieve a mortality rate of 50% when the shielding coefficient is 0.5; β is the impact energy required when P s The impact energy threshold required for death to occur when it approaches 0; E imp The calculation formula is as follows: Where: v imp Represents the speed of the drone when it hits the ground; v imp By analyzing the descent process of the drone, the specific calculation involves parameters such as air resistance and acceleration; The formula for calculating air resistance is: The calculation formula of the UAV's acceleration a is: Where: G represents gravity, which can be expressed as G = mg; in the above formula, the speed v when the drone hits the ground imp The calculation formula is: Where h is the initial flight altitude of the UAV, the calculated v imp Substitute into the formula The impact kinetic energy of the drone when it hits the ground can be obtained; Step S26: The calculation formula of the population casualty risk value matrix is ​​as follows: Where: the value of the element in R represents the probability of casualty, (x, y) is the ground geographic reference position; P event is the accident rate of drone crashing.

7. The method for generating a dynamic risk map of a UAV operating on the ground according to claim 6, characterized in that: Step S3 includes the following contents: Step S31: Generate a comprehensive ground risk map by layer-by-layer superposition and the principle of maximum risk value, including the following contents: First, we read the data of the population casualty risk layer, special area layer, no-fly zone layer, obstacle layer, and meteorological risk layer, and unify the risk level matrix dimensions of all layers. Using the maximum risk value principle, for the data stored in each grid cell in the map, we select the maximum risk level among all layers to obtain the risk level matrix. The calculation formula is as follows: RL(x,y)=max{F(x,y),O(x,y),S(x,y),R(x,y),W(x,y)} Where RL is the risk level matrix; F, O, S, R, and W are the risk level matrices for the no-fly zone, obstacle layer, special area layer, population casualty risk layer, and meteorological risk layer, respectively. The specific meanings of the risk levels are: 0 indicates safety, 1 indicates low risk, 2 indicates medium risk, 3 indicates high risk, and 4 indicates that flying is prohibited in the area due to no-fly regulations, high obstacles, or weather conditions.

8. The method for generating a dynamic risk map of a UAV operating on the ground according to claim 7, characterized in that: Step S4 includes the following contents: Step S41: Grid visualization of the risk level matrix, including the following contents: First, verify the matrix dimensions, use different colors to map different risk levels, adjust the coordinate axis offset to accurately align the grid lines with the matrix cells, and add black grid lines to strengthen the cell boundaries; finally, output a geographic heat map with a color bar.

9. A system for generating a dynamic risk map of a drone's operation to the ground, comprising an electronic device, wherein the electronic device comprises a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the computer program, it implements a method for generating a dynamic risk map of the ground for the operation of a drone as described in any one of claims 1 to 8.

10. A system for generating a dynamic risk map of a drone's operation on the ground, comprising a computer-readable storage medium storing a computer program, characterized in that: When the computer program is executed by a processor, it implements a method for generating a dynamic risk map of the ground for the operation of a drone as described in any one of claims 1 to 8.

Images