CN110160525B - Parallelization calculation method for changing flight paths of a large number of flights in dangerous weather based on discrete potential energy field - Google Patents

Parallelization calculation method for changing flight paths of a large number of flights in dangerous weather based on discrete potential energy field Download PDF

Info

Publication number
CN110160525B
CN110160525B CN201910283164.2A CN201910283164A CN110160525B CN 110160525 B CN110160525 B CN 110160525B CN 201910283164 A CN201910283164 A CN 201910283164A CN 110160525 B CN110160525 B CN 110160525B
Authority
CN
China
Prior art keywords
cell
flight
potential energy
value
density
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910283164.2A
Other languages
Chinese (zh)
Other versions
CN110160525A (en
Inventor
丁洋
郁嘉
严勇杰
李靖
付胜豪
肖雪飞
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
CETC 28 Research Institute
Original Assignee
CETC 28 Research Institute
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by CETC 28 Research Institute filed Critical CETC 28 Research Institute
Priority to CN201910283164.2A priority Critical patent/CN110160525B/en
Publication of CN110160525A publication Critical patent/CN110160525A/en
Application granted granted Critical
Publication of CN110160525B publication Critical patent/CN110160525B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01CMEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/20Instruments for performing navigational calculations

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Automation & Control Theory (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

A parallelization calculation method for a large number of flight diversion paths in dangerous weather based on a discrete potential energy field. The method comprises the steps of rasterizing an environment scene, extracting the range of a dangerous area, segmenting the environment scene, and constructing an environment inappropriateness field; constructing a speed field, establishing a flight density contribution function in the grid, and calculating the maximum speed which can be reached when the flight flies in the up, down, left and right directions; establishing a cost function, constructing functions related to three indexes of path length, time consumption and discomfort, and calculating the cost when a single flight selects a path to a target point; establishing a discrete potential energy field calculation function, carrying out discrete calculation on the potential energy field by using a unilateral differential equation, and taking the negative gradient direction of the potential energy field as a re-navigation optimal path; and (4) parallelizing potential energy field solution, adopting a partition operation rule, and solving the potential energy value of the cell in a parallel iteration manner until convergence. The invention effectively improves the calculation efficiency.

Description

Parallelization calculation method for changing flight paths of a large number of flights in dangerous weather based on discrete potential energy field
Technical Field
The invention belongs to flight diversion path planning technology, and particularly relates to a parallelization calculation method for a large number of flight diversion paths in dangerous weather based on a discrete potential energy field.
Background
Dangerous weather is an important reason threatening air transportation safety and causing flight delay, and when some airspaces or routes are influenced by the dangerous weather and cannot be used normally, the problem can be effectively solved by carrying out diversion on flights. The diversion is to arrange flights to temporarily select an unaffected route to bypass a limited flight area, and flight diversion path planning is a key problem of diversion.
The current research on flight re-voyage began in the last 90 s and after about 20 years of research, foreign and domestic scholars have proposed a variety of related algorithms, mainly including: the method comprises the steps of (1) a grid-based navigation path changing algorithm, namely a method for planning an air path by applying a path search algorithm in a network navigation environment model; the method comprises the steps of dividing areas influencing flight into different levels and giving response weight to the areas, and then applying a visual mode to determine a flight re-navigation path with the shortest weight distance; and based on the A-star search algorithm of the existing waypoints, performing path search from the waypoint set through the A-star algorithm. Although the above methods can effectively perform flight re-routing path planning, most of the methods are directed at a single flight, the calculation efficiency is low, the path planning problem when multiple flights or even a large number of flights are not considered, and meanwhile, the considered limiting factor is single and the expandability of the limiting factor is not realized.
Disclosure of Invention
Aiming at the problems, the invention provides a method for quickly generating a diversion path for a large number of flight flights under the influence of dangerous weather, and particularly provides a parallelization calculation method for the diversion path for the large number of flight flights under the influence of a discrete potential energy field.
The technical scheme of the invention is as follows: a parallelization calculation method for a large number of flight diversion paths in dangerous weather based on a discrete potential energy field comprises the following operation steps:
step 1, environment scene rasterization and discomfort field construction, namely extracting the range of a dangerous area for the dangerous area affected by weather, dividing the range of the dangerous area into polygonal areas, dividing the environment scene to form fine cells, calculating the number of cells covered by the dangerous area, and then calculating the discomfort value of each cell according to the flight density of the area;
step 2, constructing a speed field, establishing a flight density contribution function in the grid, determining the density contribution value of each flight to the cell and the density mapping value of each cell, then taking the density contribution value of the flight as a weight, calculating the weighted average speed in the cell, and decomposing the weighted average speed into the maximum speed which can be reached when the flight is flown towards the direction according to the upper direction, the lower direction, the left direction and the right direction;
step 3, establishing a cost function, namely establishing a function related to three indexes of path length, time consumption and discomfort to calculate the cost when a single flight selects a flight path to a target point, then discretizing the cost function, and establishing a unit cost function of the flight entering upper, lower, left and right adjacent cells from an initial cell according to rasterization processing in the step 1;
step 4, establishing a discrete potential energy field calculation function, namely establishing a search optimization problem of an optimal path from an initial position S to a target position G in a plane space according to the cost value of the path from the current position to the target position obtained in the step 3, calculating the potential energy value of any cell in the area, discretizing by using a one-side difference equation, and taking the direction of a negative gradient as a re-navigation optimal path;
and 5, parallelizing potential energy field solving, namely dividing all cells of the environment area into blocks tiles with predefined sizes by adopting a partition operation rule, and solving potential energy values of the cells in each tile in parallel so as to solve the potential energy field values in the cells and quickly generate the potential energy field by utilizing the calculation efficiency of the GPU.
Wherein, the environment scene rasterization and the discomfort field construction in the step 1 comprise the following steps:
step 1.1, dividing an environment scene into cells in the X direction and the Y direction respectively, converting a flight limited area into a polygon and mapping the polygon into the cells;
step 1.2, setting the more flights in a cell, the higher the discomfort value of the cell is, and when the density is higher than a certain threshold value, the discomfort value is positive and infinite; defining that if the discomfort value of the cell is positive infinity, the cell represents that the cell can not pass completely, and initializing the discomfort value of the cell belonging to the flight limited area to be positive infinity; the calculation of the discomfort value for a cell is given by the following equation:
Figure BDA0002022365820000021
wherein, the speed field in step 2 is constructed as follows:
step 2.1, establishing the density contribution of flights to the cells, and setting the distance between the flights and the center point of the cell C in the X-axis direction as delta X and the distance in the Y-axis direction as delta Y; the density contribution of the flight to cell A, B, C, D may be formulated as:
Figure BDA0002022365820000022
step 2.2, calculating the density contribution value of the cell, and adding the density values of flights having density contribution to the current cell according to the step 2.1 to obtain the density mapping value of the current cell, wherein the formula is as follows:
ρi,j=∑a∈Aρa
step 2.3, calculating the weighted average speed of the cells, taking the density contribution of the flight to the current cell as a weight, and calculating to obtain the weighted average speed in the (i, j) th cell, as described in the following formula:
Figure BDA0002022365820000031
step 2.4, calculating the maximum speed of the flight in the specific direction, wherein the calculation is related to the flight density, and the calculation formula is as follows:
Figure BDA0002022365820000032
in the step 3: the integral operation of the cost of the flight diversion path can be approximated by the sum of unit cost values of adjacent cells from the cell to the upper, lower, left and right, and the calculation mode of constructing the cost function is as follows:
Figure BDA0002022365820000033
Figure BDA0002022365820000034
the discrete potential energy field calculation function in the step 4 is constructed by the following steps:
step 4.1, constructing a discrete potential energy field calculation function, and performing discretization approximate solution according to a Upwind discrete method, wherein the potential energy field value calculation formula is as follows:
Figure BDA0002022365820000035
step 4.2, calculating the gradient of the potential energy field, according to the definition of the gradient, pointing the direction of the gradient to the direction in which the scalar field grows fastest, adopting a unilateral difference equation for discretization, and aiming at the gradient component in the X direction, the calculation method comprises the following steps:
Figure BDA0002022365820000036
in the step 5, the solution of the parallelization potential energy field is divided into three steps of algorithm definition, algorithm initialization and cyclic calculation, and the specific implementation process is as follows:
step 5.1, algorithm definition, namely setting t to represent a single tile for dividing the environment area into equal sizes, and setting an operation list L for storing the tiles needing to be updated;
step 5.2, algorithm initialization, namely uniformly dividing all cells into blocks tile with the same size, and moving all tiles containing target area cells into an operation list L;
and 5.3, circularly calculating, namely calculating the potential energy field value of each cell to be converged through multiple iterations to finish updating.
The specific implementation process of the loop calculation in the step 5.3 is as follows:
step 5.3.1, updating operation is carried out on each block t in the operation list L: calculating potential energy values of all cells in t circularly for n times; comparing the current result with the previous result after each calculation, and marking the cell as convergence if the difference value is smaller than a threshold value; traversing all the cells in the t, if all the cells are marked as convergence, marking the t as convergence, otherwise, continuing to perform the circulation in the previous step;
step 5.3.2, checking the four-way neighbor of the block t: calculating potential energy values of all the blocks t converged in the previous step once for all the adjacent blocks, and checking whether the potential energy values of the cells in the adjacent blocks of t are changed; adding the block t which is changed in the previous step into the operation list L again, and deleting the block t of which the potential energy value is not changed from the operation list L;
and 5.3.3, enabling the L to contain all unconverged tiles, and repeating the steps a and b until the L is empty.
The invention has the beneficial effects that: (1) a cost function suitable for a re-navigation scene is established by utilizing a potential energy field model concept, constraints such as fuel economy, a flight limited area and the like are considered, an extensible cost function is established, and the extensible cost function has the expansibility of a limiting condition; (2) the potential energy field model is a macroscopic calculation model, and the optimal path for diversion can be calculated for the diversion activities of a large number of flights by establishing the potential energy field; (3) a parallelized calculation algorithm is designed, the operation characteristics and efficiency of the GPU can be fully utilized, and the calculation performance is improved by about 5 times.
Drawings
FIG. 1 is a flow chart of the present invention;
FIG. 2 is a diagram of a discrete cell structure in accordance with the present invention;
FIG. 3 is a cell density mapping rule in accordance with the present invention;
FIG. 4 is a flow chart of a parallelized potential energy field solution algorithm in accordance with the present invention;
fig. 5 is a schematic diagram of a tile grid in the present invention.
Detailed Description
The following description of the embodiments of the present invention is provided in connection with the accompanying drawings and the examples, but the invention is not limited thereto.
The invention discloses a parallelization calculation method for a great number of flight diversion paths in dangerous weather based on a discrete potential energy field, which comprises the following operation steps:
step 1, environment scene rasterization and discomfort field construction, namely extracting the range of a dangerous area for the dangerous area affected by weather, dividing the range of the dangerous area into polygonal areas, dividing the environment scene to form fine cells, calculating the number of cells covered by the dangerous area, and then calculating the discomfort value of each cell according to the flight density of the area;
step 2, constructing a speed field, establishing a flight density contribution function in the grid, determining the density contribution value of each flight to the cell and the density mapping value of each cell, then taking the density contribution value of the flight as a weight, calculating the weighted average speed in the cell, and decomposing the weighted average speed into the maximum speed which can be reached when the flight is flown towards the direction according to the upper direction, the lower direction, the left direction and the right direction;
step 3, establishing a cost function, namely establishing a function related to three indexes of path length, time consumption and discomfort to calculate the cost when a single flight selects a flight path to a target point, then discretizing the cost function, and establishing a unit cost function of the flight entering upper, lower, left and right adjacent cells from an initial cell according to rasterization processing in the step 1;
step 4, establishing a discrete potential energy field calculation function, namely establishing a search optimization problem of an optimal path from an initial position S to a target position G in a plane space according to the cost value of the path from the current position to the target position obtained in the step 3, calculating the potential energy value of any cell in the area, discretizing by using a one-side difference equation, and taking the direction of a negative gradient as a re-navigation optimal path;
and 5, parallelizing potential energy field solving, namely dividing all cells of the environment area into blocks tiles with predefined sizes by adopting a partition operation rule, and solving potential energy values of the cells in each tile in parallel so as to solve the potential energy field values in the cells and quickly generate the potential energy field by utilizing the calculation efficiency of the GPU.
Wherein, the environment scene rasterization and the discomfort field construction in the step 1 comprise the following steps:
step 1.1, dividing an environment scene into cells in the X direction and the Y direction respectively, converting a flight limited area into a polygon and mapping the polygon into the cells;
step 1.2, setting the more flights in a cell, the higher the discomfort value of the cell is, and when the density is higher than a certain threshold value, the discomfort value is positive and infinite; defining that if the discomfort value of the cell is positive infinity, the cell represents that the cell can not pass completely, and initializing the discomfort value of the cell belonging to the flight limited area to be positive infinity; the calculation of the discomfort value for a cell is given by the following equation:
Figure BDA0002022365820000051
wherein, the speed field in step 2 is constructed as follows:
step 2.1, establishing the density contribution of flights to the cells, and setting the distance between the flights and the center point of the cell C in the X-axis direction as delta X and the distance in the Y-axis direction as delta Y; the density contribution of the flight to cell A, B, C, D may be formulated as:
Figure BDA0002022365820000061
step 2.2, calculating the density contribution value of the cell, and adding the density values of flights having density contribution to the current cell according to the step 2.1 to obtain the density mapping value of the current cell, wherein the formula is as follows:
ρi,j=∑a∈Aρa
step 2.3, calculating the weighted average speed of the cells, taking the density contribution of the flight to the current cell as a weight, and calculating to obtain the weighted average speed in the (i, j) th cell, as described in the following formula:
Figure BDA0002022365820000062
step 2.4, calculating the maximum speed of the flight in the specific direction, wherein the calculation is related to the flight density, and the calculation formula is as follows:
Figure BDA0002022365820000063
in the step 3: the integral operation of the cost of the flight diversion path can be approximated by the sum of unit cost values of adjacent cells from the cell to the upper, lower, left and right, and the calculation mode of constructing the cost function is as follows:
Figure BDA0002022365820000064
Figure BDA0002022365820000065
the discrete potential energy field calculation function in the step 4 is constructed by the following steps:
step 4.1, constructing a discrete potential energy field calculation function, and performing discretization approximate solution according to a Upwind discrete method, wherein the potential energy field value calculation formula is as follows:
Figure BDA0002022365820000071
step 4.2, calculating the gradient of the potential energy field, according to the definition of the gradient, pointing the direction of the gradient to the direction in which the scalar field grows fastest, adopting a unilateral difference equation for discretization, and aiming at the gradient component in the X direction, the calculation method comprises the following steps:
Figure BDA0002022365820000072
in the step 5, the solution of the parallelization potential energy field is divided into three steps of algorithm definition, algorithm initialization and cyclic calculation, and the specific implementation process is as follows:
step 5.1, algorithm definition, namely setting t to represent a single tile for dividing the environment area into equal sizes, and setting an operation list L for storing the tiles needing to be updated;
step 5.2, algorithm initialization, namely uniformly dividing all cells into blocks tile with the same size, and moving all tiles containing target area cells into an operation list L;
and 5.3, circularly calculating, namely calculating the potential energy field value of each cell to be converged through multiple iterations to finish updating.
The specific implementation process of the loop calculation in the step 5.3 is as follows:
step 5.3.1, updating operation is carried out on each block t in the operation list L: calculating potential energy values of all cells in t circularly for n times; comparing the current result with the previous result after each calculation, and marking the cell as convergence if the difference value is smaller than a threshold value; traversing all the cells in the t, if all the cells are marked as convergence, marking the t as convergence, otherwise, continuing to perform the circulation in the previous step;
step 5.3.2, checking the four-way neighbor of the block t: calculating potential energy values of all the blocks t converged in the previous step once for all the adjacent blocks, and checking whether the potential energy values of the cells in the adjacent blocks of t are changed; adding the block t which is changed in the previous step into the operation list L again, and deleting the block t of which the potential energy value is not changed from the operation list L;
and 5.3.3, enabling the L to contain all unconverged tiles, and repeating the steps a and b until the L is empty.
With reference to fig. 1, the method for calculating the parallelization diversion path based on the discrete potential energy field of the present invention includes the following steps:
firstly, performing scene rasterization on an environment scene, and establishing an uncomfortable field according to a flight limited area; dividing n cells in the X direction and the Y direction according to the range of the environment scene, and identifying each cell by using coordinates (i, j); each cell stores data such as density, speed, potential energy and the like, as shown in fig. 2, wherein the discomfort value represents the selection tendency of a certain group of flights with the same target for entering the cell, and the invention assumes that the discomfort degree of a cell is related to the flight density in the cell, namely, the more flights in the cell, the higher the discomfort value of the cell; secondly, grouping all flights to be calculated according to the diversion target cell; the specific implementation steps are as follows:
(1) and environment scene segmentation: dividing the environment scene into cells in the X direction and the Y direction respectively, converting the flight limited area into polygons and mapping the polygons into the cells, obviously, the smaller the area of the cells is, the finer the calculation result is; the larger the area of the cell is, the faster the calculation speed is, and in order to ensure the real-time performance of calculation, the cell with a proper size needs to be selected;
(2) constructing an inadaptation field: for the sake of simplifying the calculation, the invention assumes that the discomfort degree is only related to the flight density, i.e. the more flights in a cell, the higher the discomfort degree value of the cell, and when the density is higher than a certain threshold value, the discomfort degree value is positive and infinite; if the discomfort value of the cell is positive infinity, the cell represents that the cell cannot pass at all, the discomfort value of the cell belonging to the flight limited area is initialized to be positive infinity, and the calculation of the discomfort is given by the following formula:
Figure BDA0002022365820000081
in the formula: c. Ci,j-a discomfort value of cell (i, j); n isi,j-number of flights in cell (i, j); ε -flight number density threshold for a cell.
And secondly, as shown in fig. 3, establishing a correlation between flight density and speed, that is, establishing a density mapping rule, taking the density contribution of the flight to the current cell as a weight, calculating to obtain a weighted average speed of the flight in the (i, j) th cell, and decomposing the speed into a maximum speed which can be reached when all flights in the cell are influenced when the flights fly in four directions, namely, the maximum speed is called the flight limit speed of the flights in the four directions, namely, the flight limit speed of the cell (i, j). The method comprises the following specific steps:
(1) establishing the density contribution of the flights to the cells, wherein each flight has a separate density field, the density field takes the maximum value at the position of the flight and gradually decreases towards the periphery of the flight, and the distance between the flight and the center point of the cell C in the X-axis direction is set as deltax, and the distance in the Y-axis direction is set as deltay, as shown in the figure; the density contribution of the flight to cell A, B, C, D may be formulated as:
Figure BDA0002022365820000082
in the formula: rhoX-the density contribution of flights to cell X; Δ X is the distance between the flight and the center point of the unit where the flight is located in the X-axis direction; delta Y is the distance between the flight and the center point of the unit where the flight is located in the Y-axis direction; lambda-decay rate of density, lambda determines that the flight contributes no less to the density of the cell in which it is located
Figure BDA0002022365820000091
Density contribution to adjacent cells is not greater than
Figure BDA0002022365820000092
(2) Calculating the density contribution value of the cell, and calculating the density contribution of the flight to the cell which can be influenced by the flight according to a formula 3.3, so that the density values of the flights having the density contribution to the current cell are added to obtain the density mapping value of the current cell, wherein the formula is as follows:
ρi,j=∑a∈Aρa
in the formula: rhoi,j-a density map value of a cell (i, j); a-single flight; a-set of all flights that have a density contribution to cell (i, j); rhoa-the density contribution of flight a to cell (i, j);
(3) calculating the weighted average speed of the cells, taking the density contribution of the flight to the current cell as a weight, and calculating to obtain the weighted average speed in the (i, j) th cell, as described in the formula:
Figure BDA0002022365820000093
in the formula:
Figure BDA0002022365820000094
-weighted average velocity of cell (i, j); v. ofa-the speed of flight a;
(4) when the flight density in the adjacent cell corresponding to the flight direction is larger, the upper limit of the flight speed of the flight to the flight direction is the flight flow limit speed of the cell in the flight direction; when the density of the adjacent cells in the corresponding direction is lower, the flight in the direction can be free from interference, and the flight can keep the self ideal speed; if the density of the adjacent cells corresponding to the direction is at a middle value, the flight speed of the flight in the direction is obtained by linear interpolation. The formula of the influence speed of the density is obtained by the following steps:
Figure BDA0002022365820000095
in the formula:
Figure BDA0002022365820000096
-the achievable flight speed of the flight in cell (i, j) in direction d; v. ofmax-the desired flight speed of the flight; rhoi',j'-a density map value of a neighboring cell in the d-direction of cell (i, j); rhomin-a density map lower limit value for the cell; rhomax-a cell density map upper limit value.
Thirdly, establishing a cost function, wherein in a rasterized environment area, a path can be regarded as a series of ordered cell sequences, and flights can enter upper, lower, left and right adjacent cells from a certain cell along a certain direction, so that the integral operation of the cost of a single path can be approximated by the sum of unit cost values entering the upper, lower, left and right adjacent cells from the cell, namely, the calculation mode of constructing the cost function is as follows:
Figure BDA0002022365820000101
Figure BDA0002022365820000102
in the formula:
Figure BDA0002022365820000103
-a cost value from a cell (i, j) into a neighboring cell along the direction d;
Figure BDA0002022365820000104
-the achievable speed of flight in cell in direction d; c. Ci',j'-the discomfort index in the adjacent bin corresponding to the direction is positive infinity if the environment zone boundary is exceeded.
Fourthly, establishing a discrete potential energy field calculation function, constructing a scalar field p of the environmental scene, and taking a gradient descent curve of the scalar field p as an optimal path from the point A to the point B on the plane, wherein the specific steps are as follows:
(1) the construction function of the potential energy field can be expressed by the following formula:
Figure BDA0002022365820000105
in the formula: g-a set of cells representing a target area;
according to the Upwind discrete method, the above formula can be discretized to perform approximate solution, and the formula is as follows:
Figure BDA0002022365820000106
in the formula: pi,j-potential energy value of cell (i, j); px-the value of the neighbor potential in the X-axis direction; py-Y-axis direction neighbor potential value; e.g. of the typex-entry and PxThe same X-direction neighbor unit cost value; e.g. of the typey-entry and PyThe same Y-direction neighbor unit cost value. Wherein exAnd eyThe value of (A) is obtained by calculation in the step three;
(2) calculating the gradient of the potential energy field, wherein the gradient direction points to the direction in which the scalar field is fastest growing according to the definition of the gradient, so that the gradient of a unit cell (i, j) is calculated, potential energy values of the unit cell and four adjacent cells, namely an upper adjacent cell, a lower adjacent cell, a left adjacent cell and a right adjacent cell, need to be compared, a unilateral difference equation is adopted for discretization, and the calculation method is as follows aiming at the gradient component in the X direction:
Figure BDA0002022365820000111
the gradient component in the Y direction is calculated similarly.
Fifthly, numerical solution is carried out on the potential energy field by adopting a parallelization calculation algorithm, all cells of the environment area are divided into blocks (tiles) with predefined sizes by adopting a partition operation rule, wherein the solution of the potential energy value of the cell in each tile is carried out in parallel, and the specific steps are as follows:
(1) and algorithm definition:
a. setting t to represent a single tile which divides the environment area into equal size;
b. setting an operation list L for storing tiles needing to be updated;
(2) and algorithm initialization:
a. uniformly dividing all the cells into blocks tile with the same size;
b. moving all tiles containing target area cells into an operation list L;
(3) and circularly calculating:
a. updating operation is carried out on each block t in the operation list L: calculating potential energy values of all cells in t circularly for n times; and comparing the current result with the previous result after each calculation, and marking the cell as convergence if the difference value is less than the threshold value. And traversing all the cells in the t, if all the cells are marked as convergence, marking the t as convergence, and if not, continuing to perform the loop in the previous step.
b. Check block t four-way neighbor: calculating potential energy values of all the blocks t converged in the previous step once for all the adjacent blocks, and checking whether the potential energy values of the cells in the adjacent blocks of t are changed; adding the block t which is changed in the previous step into the operation list L again, and deleting the block t of which the potential energy value is not changed from the operation list L;
c. and c, repeating the steps a and b until the L is empty.
The invention has the following effects: (1) a cost function suitable for a re-navigation scene is established by utilizing a potential energy field model concept, constraints such as fuel economy, a flight limited area and the like are considered, an extensible cost function is established, and the extensible cost function has the expansibility of a limiting condition; (2) the potential energy field model is a macroscopic calculation model, and the optimal path for diversion can be calculated for the diversion activities of a large number of flights by establishing the potential energy field; (3) a parallelized calculation algorithm is designed, the operation characteristics and efficiency of the GPU can be fully utilized, and the calculation performance is improved by about 5 times. The method comprises the steps of rasterizing an environment scene, extracting the range of a dangerous area, segmenting the environment scene, and constructing an environment inappropriateness field; constructing a speed field, establishing a flight density contribution function in the grid, and calculating the maximum speed which can be reached when the flight flies in the up, down, left and right directions; establishing a cost function, constructing functions related to three indexes of path length, time consumption and discomfort, and calculating the cost when a single flight selects a path to a target point; establishing a discrete potential energy field calculation function, carrying out discrete calculation on the potential energy field by using a unilateral differential equation, and taking the negative gradient direction of the potential energy field as a re-navigation optimal path; and (4) parallelizing potential energy field solution, adopting a partition operation rule, and solving the potential energy value of the cell in a parallel iteration manner until convergence. The invention effectively improves the calculation efficiency.

Claims (1)

1. A parallelization calculation method for a large number of flight diversion paths in dangerous weather based on a discrete potential energy field is characterized by comprising the following operation steps:
step 1, environment scene rasterization and discomfort field construction, namely extracting the range of a dangerous area for the dangerous area affected by weather, dividing the range of the dangerous area into polygonal areas, dividing the environment scene to form fine cells, calculating the number of cells covered by the dangerous area, and then calculating the discomfort value of each cell according to the flight density of the area;
step 2, constructing a speed field, establishing a flight density contribution function in the grid, determining the density contribution value of each flight to the cell and the density mapping value of each cell, then taking the density contribution value of the flight as a weight, calculating the weighted average speed in the cell, and decomposing the weighted average speed into the maximum speed which can be reached when the flight is flown towards the direction according to the upper direction, the lower direction, the left direction and the right direction;
step 3, establishing a cost function, namely establishing a function related to three indexes of path length, time consumption and discomfort to calculate the cost when a single flight selects a flight path to a target point, then discretizing the cost function, and establishing a unit cost function of the flight entering upper, lower, left and right adjacent cells from an initial cell according to rasterization processing in the step 1;
step 4, establishing a discrete potential energy field calculation function, namely establishing a search optimization problem of an optimal path from an initial position S to a target position G in a plane space according to the cost value of the path from the current position to the target position obtained in the step 3, calculating the potential energy value of any cell in the area, discretizing by using a one-side difference equation, and taking the direction of a negative gradient as a re-navigation optimal path;
step 5, parallelizing potential energy field solving, namely adopting a partition operation rule to divide all cells of an environment area into blocks tiles with predefined sizes, and solving potential energy values of the cells in each tile in parallel so as to solve the potential energy field values in the cells and quickly generate a potential energy field by utilizing the calculation efficiency of a GPU;
specifically, the environment scene rasterization and discomfort field construction step in step 1 is as follows:
step 1.1, dividing an environment scene into cells in the X direction and the Y direction respectively, converting a flight limited area into a polygon and mapping the polygon into the cells;
step 1.2, setting the more flights in a cell, the higher the discomfort value of the cell is, and when the density is higher than a certain threshold value, the discomfort value is positive and infinite; defining that if the discomfort value of the cell is positive infinity, the cell represents that the cell can not pass completely, and initializing the discomfort value of the cell belonging to the flight limited area to be positive infinity; the calculation of the discomfort value for a cell is given by the following equation:
Figure FDA0003235698490000021
wherein, ci,jRepresenting a discomfort value for cell (i, j); n isi,jRepresenting the number of flights in cell (i, j); epsilon represents the flight number density threshold of the cell;
the velocity field in step 2 is constructed as follows:
step 2.1, establishing the density contribution of flights to the cells, and setting the distance between the flights and the center point of the cell C in the X-axis direction as delta X and the distance in the Y-axis direction as delta Y; the density contribution of the flight to cell A, B, C, D may be formulated as:
Figure FDA0003235698490000022
where ρ isA、ρB、ρC、ρDRepresents the density contribution of flights to cell A, B, C, D; Δ X represents the distance between the center point of the flight and the center point of the unit where the flight is located in the X-axis direction; Δ Y represents the distance between the center point of the flight and the center point of the unit where the flight is located in the Y-axis direction; lambda represents the attenuation speed of the density, and lambda determines that the flight contributes to the density of the cell in which the flight is positioned to be not less than
Figure FDA0003235698490000023
Density contribution to adjacent cells is not greater than
Figure FDA0003235698490000024
Step 2.2, calculating the density contribution value of the cell, and adding the density values of flights having density contribution to the current cell according to the step 2.1 to obtain the density mapping value of the current cell, wherein the formula is as follows:
ρi,j=∑a∈Aρa
where ρ isi,jA density map value representing a cell (i, j); a represents a single flight; a represents the set of all flights that have a density contribution to cell (i, j); rhoaRepresenting the density contribution value of flight a to cell (i, j);
step 2.3, calculating the weighted average speed of the cells, taking the density contribution of the flight to the current cell as a weight, and calculating to obtain the weighted average speed in the (i, j) th cell, as described in the following formula:
Figure FDA0003235698490000025
wherein the content of the first and second substances,
Figure FDA0003235698490000031
represents the weighted average velocity of cell (i, j); v. ofaRepresenting the speed of flight a;
step 2.4, calculating the maximum speed of the flight in the specific direction, wherein the calculation is related to the flight density, and the calculation formula is as follows:
Figure FDA0003235698490000032
wherein the content of the first and second substances,
Figure FDA0003235698490000033
representing the reachable flight speed of the flight in the cell (i, j) towards the direction d; v. ofmaxRepresenting the ideal flight speed of the flight; rhoi',j'Representing the density of cells adjacent in the d direction of cell (i, j)A degree mapping value; rhominA density map lower limit value representing a cell; rhomaxA density map upper limit value representing a cell;
in the step 3: the integral operation of the cost of the flight diversion path can be approximated by the sum of unit cost values of adjacent cells from the cell to the upper, lower, left and right, and the calculation mode of constructing the cost function is as follows:
Figure FDA0003235698490000034
Figure FDA0003235698490000035
wherein E isp' represents a cost value of the flight diversion path p;
Figure FDA0003235698490000036
represents the cost value from cell (i, j) into the neighboring cell along direction d;
Figure FDA0003235698490000037
representing the achievable speed of flight in the cell in the direction d; c. Ci',j'Representing the discomfort index in the adjacent grid in the corresponding direction, and if the discomfort index exceeds the boundary of the environment area, the discomfort index is positive infinity;
the discrete potential energy field calculation function in the step 4 is constructed by the following steps:
step 4.1, constructing a discrete potential energy field calculation function, and performing discretization approximate solution according to a Upwind discrete method, wherein the potential energy field value calculation formula is as follows:
Figure FDA0003235698490000038
wherein, Pi,jRepresents the potential energy value of the cell (i, j); pxRepresenting the potential energy value of the adjacent lattice in the X-axis direction; pyTo representThe adjacent lattice potential energy value in the Y-axis direction; e.g. of the typexIndicating entry with PxThe same X-direction neighbor unit cost value; e.g. of the typeyIndicating entry with PyThe same adjacent lattice unit cost value in the Y direction; wherein e isxAnd eyThe value of (A) is obtained by calculation in the step three;
step 4.2, calculating the gradient of the potential energy field, according to the definition of the gradient, pointing the direction of the gradient to the direction in which the scalar field grows fastest, adopting a unilateral difference equation for discretization, and aiming at the gradient component in the X direction, the calculation method comprises the following steps:
Figure FDA0003235698490000041
in the step 5, the solution of the parallelization potential energy field is divided into three steps of algorithm definition, algorithm initialization and cyclic calculation, and the specific implementation process is as follows:
step 5.1, algorithm definition, namely setting t to represent a single tile for carrying out equal-size division on the environment area, and setting an operation list L for storing the tiles needing to be updated;
step 5.2, algorithm initialization, namely uniformly dividing all cells into blocks tile with the same size, and moving all tiles containing target area cells into an operation list L;
step 5.3, circularly calculating, namely calculating the potential energy field value of each cell to be converged through multiple iterations to complete updating;
the specific implementation process of the loop calculation in the step 5.3 is as follows:
step 5.3.1, updating operation is carried out on each block t in the operation list L: calculating potential energy values of all cells in t circularly for n times; comparing the current result with the previous result after each calculation, and marking the cell as convergence if the difference value is smaller than a threshold value; traversing all the cells in the t, if all the cells are marked as convergence, marking the t as convergence, otherwise, continuing to perform the circulation in the previous step;
step 5.3.2, checking the four-way neighbor of the block t: calculating potential energy values of all the blocks t converged in the previous step once for all the adjacent blocks, and checking whether the potential energy values of the cells in the adjacent blocks of t are changed; adding the block t which is changed in the previous step into the operation list L again, and deleting the block t of which the potential energy value is not changed from the operation list L;
and 5.3.3, enabling the L to contain all unconverged tiles, and repeating the steps a and b until the L is empty.
CN201910283164.2A 2019-04-10 2019-04-10 Parallelization calculation method for changing flight paths of a large number of flights in dangerous weather based on discrete potential energy field Active CN110160525B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910283164.2A CN110160525B (en) 2019-04-10 2019-04-10 Parallelization calculation method for changing flight paths of a large number of flights in dangerous weather based on discrete potential energy field

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910283164.2A CN110160525B (en) 2019-04-10 2019-04-10 Parallelization calculation method for changing flight paths of a large number of flights in dangerous weather based on discrete potential energy field

Publications (2)

Publication Number Publication Date
CN110160525A CN110160525A (en) 2019-08-23
CN110160525B true CN110160525B (en) 2022-02-11

Family

ID=67638547

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910283164.2A Active CN110160525B (en) 2019-04-10 2019-04-10 Parallelization calculation method for changing flight paths of a large number of flights in dangerous weather based on discrete potential energy field

Country Status (1)

Country Link
CN (1) CN110160525B (en)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112396870A (en) * 2020-10-14 2021-02-23 华北理工大学 Flight forbidden zone setting method based on coordinate sorting Graham-scan
CN115686072B (en) * 2022-12-29 2023-03-14 中国电子科技集团公司第二十八研究所 Automatic generation method of unmanned aerial vehicle group safe operation air line based on spatial grid

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2690702B2 (en) * 1994-09-13 1997-12-17 日本電気株式会社 Self-organizing device
CN107480096B (en) * 2017-08-21 2020-03-31 西安交通大学 High-speed parallel computing method in large-scale group simulation
CN109459026B (en) * 2018-11-08 2021-01-01 北京理工大学 Multi-moving-body collaborative full-coverage path planning method

Also Published As

Publication number Publication date
CN110160525A (en) 2019-08-23

Similar Documents

Publication Publication Date Title
CN107145161B (en) Flight path planning method and device for unmanned aerial vehicle to access multiple target points
CN107103164B (en) Distribution method and device for unmanned aerial vehicle to execute multiple tasks
CN109144102A (en) A kind of Path Planning for UAV based on improvement bat algorithm
Chakrabarty et al. Flight path planning for uav atmospheric energy harvesting using heuristic search
CN106225800B (en) Environmentally friendly vehicle navigation path construction method based on real-time road condition information
CN110262563A (en) Multiple no-manned plane collaboratively searching mesh calibration method waterborne
CN106444835A (en) Underwater vehicle three-dimensional path planning method based on Lazy Theta satellite and particle swarm hybrid algorithm
CN102880186A (en) Flight path planning method based on sparse A* algorithm and genetic algorithm
CN110160525B (en) Parallelization calculation method for changing flight paths of a large number of flights in dangerous weather based on discrete potential energy field
CN108388250A (en) A kind of unmanned surface vehicle paths planning method based on adaptive cuckoo searching algorithm
CN109032167A (en) Unmanned plane paths planning method based on Parallel Heuristic Algorithm
CN111121784B (en) Unmanned reconnaissance aircraft route planning method
CN114812564B (en) Multi-target unmanned aerial vehicle path planning method based on urban dynamic space-time risk analysis
CN110986954B (en) Military transport plane route planning method based on gray wolf optimization algorithm
CN104504198A (en) Airway network topology designing method based on double-layered co-evolution
CN113610312A (en) Ship navigation real-time optimal route planning method based on improved genetic algorithm
CN115713856A (en) Vehicle path planning method based on traffic flow prediction and actual road conditions
CN108445894A (en) A kind of secondary paths planning method considering unmanned boat movenent performance
CN115454132A (en) Unmanned aerial vehicle inspection three-dimensional path planning method and system for mountain photovoltaic power station
CN109445462B (en) Unmanned aerial vehicle robust path planning method in uncertain environment
CN113256013B (en) Intelligent vehicle path searching method under environmental constraint
CN107229998A (en) A kind of autonomous pathfinding strategy process of unmanned plane
CN107480096B (en) High-speed parallel computing method in large-scale group simulation
Ghorpade Airspace configuration model using swarm intelligence based graph partitioning
CN107590297A (en) The online planing method of aircraft reentry trajectory based on particle swarm optimization algorithm

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
CB02 Change of applicant information
CB02 Change of applicant information

Address after: 210046 No.1, Lingshan South Road, Qixia District, Nanjing City, Jiangsu Province

Applicant after: THE 28TH RESEARCH INSTITUTE OF CHINA ELECTRONICS TECHNOLOGY Group Corp.

Address before: 210001 No.1, alfalfa Garden East Street, Qinhuai District, Nanjing City, Jiangsu Province

Applicant before: THE 28TH RESEARCH INSTITUTE OF CHINA ELECTRONICS TECHNOLOGY Group Corp.

GR01 Patent grant
GR01 Patent grant