Flight path planning method for aircraft based on subdivision grid

Abstract

The present invention discloses a flight path planning method for an aircraft based on a dissected grid, which includes two major steps. First, based on the dissected grid, an environmental modeling of multiple radar detection areas is carried out to characterize the flight physical space of the aircraft. Then, in the physical space organized based on the dissected grid, an improved A* algorithm is used to plan the flight path of the aircraft. In the present invention, a new underlying representation framework of the three-dimensional dissected grid theory is adopted to effectively solve the problem of repeated modeling in the same space caused by the changing requirements of the spatial reference distance brought by different aircraft in the same space. And based on this theory, the traditional A* algorithm is improved to form a flight path planning method that is more adaptable to the background of aircraft penetration, more general, efficient and autonomous, and assists the pilot in making flight path decisions.

Description

Technical Field

[0001] The present invention relates to the technical field of trajectory planning, and particularly to an aircraft trajectory planning method based on a dissected grid. Background Art

[0002] With the continuous evolution of information-based warfare, various threats to aircraft in the aerospace battlefield have developed rapidly. This requires pilots to not only perform complex manual operations but also make judgment decisions such as trajectory selection when facing threats. However, due to physiological and psychological limitations, it has become increasingly difficult to rely solely on pilots to complete tasks, and there is an urgent need for a more general, efficient, and autonomous trajectory planning ability to assist in decision-making. The requirement for aircraft trajectory planning in the context of penetration in military applications is to quickly plan an effective trajectory that conforms to the combat intention for the aircraft before departure or during the mission after comprehensively considering factors such as the algorithm planning time, the final trajectory length, and the distribution of radar threat areas, and to achieve a penetration with controllable risks during this process.

[0003] Aircraft come in a wide variety of types with different characteristics. In traditional trajectory planning methods, the data organizational structure usually adopts a method of grid division based on a longitude-latitude grid or a local airspace rectangular coordinate system grid. In the target physical space, threat areas are discretely sampled and calculated according to a certain data organization method to support subsequent trajectory planning applications. Assuming that at a fixed time reference, when the flight speed of the aircraft changes and the spatial reference distance must be changed, the traditional method can only recalculate the threat area for the same environment at a new spatial scale and cannot effectively integrate aircraft with different speeds in the same method and framework for associated calculation, nor can it hierarchically associate different spatial organization granularities to meet the needs of multi-type aircraft trajectory planning. Summary of the Invention

[0004] Aiming at the above problems, the present invention aims to provide an aircraft trajectory planning method based on a dissected grid, which uses the dissected grid to solve the problem of airspace threat modeling, and on this basis, improves the A* algorithm for trajectory planning by using the characteristics of the grid, so as to obtain a more general, efficient, and autonomous trajectory planning method to assist pilots in making trajectory decisions.

[0005] In order to achieve the above object, the technical solution adopted by the present invention is as follows:

[0006] An aircraft trajectory planning method based on a dissected grid, characterized by including the following steps,

[0007] S1: Based on the dissected grid, perform multi-radar detection area environment modeling to represent the flight physical space of the aircraft;

[0008] S2: In the physical space organized by the dissected grid in step S1, use the improved A* algorithm to perform aircraft trajectory planning.

[0009] Further, the specific operations for multi-radar detection area environment modeling in step S1 include the following steps:

[0010] S101: Calculate the detection probability of a single radar;

[0011] S102: Calculate the combined detection probability of multiple radars in the multi-radar detection area based on the single-radar detection probability calculated in step S101;

[0012] S103: Based on the dissection grid organization, calculate the radar detection probability at any granularity level according to the combined detection probability of multiple radars;

[0013] S104: Use the correlation relationship between different granularity levels under the dissection grid organization to calculate the detection probability corresponding to the dissection volume element at the new granularity level;

[0014] S105: Perform visual modeling on the radar detection area based on the dissection grid according to the detection probability values within the corresponding dissection grid volume elements.

[0015] Further, the specific operations of step S101 include the following steps:

[0016] S1011: Express the classical single-radar detection probability as where represent the average detection probability and average signal power respectively, N represents the receiver noise power, b and n represent the detection threshold voltage value and the number of integrated pulses respectively;

[0017] S1012: Correlate and calculate the standard aircraft radar cross-section and the actual target radar cross-section to obtain where it is assumed that the cross-section of target i is the standard cross-section and the cross-section of target m is the actual cross-section, represent the single-radar average signal powers corresponding to i and m respectively, and represent the detection probabilities corresponding to i and m respectively;

[0018] S1013: According to the radar basic equation it can be obtained that where is the signal-to-noise ratio, is the maximum average radar cross-section of the target aircraft, R represents the Euclidean distance between the target and the radar, and K 0 is a constant;

[0019] S1014: When the detection probability of a radar with fixed parameters is set to 0.1, use σ mc and R to represent its corresponding target radar cross-section and maximum detection distance respectivelymc representation, then

[0020] S1015: Substitute the in step S1014 into the in step S1013, and the single radar detection probability can be obtained, and its calculation formula is

[0021] Furthermore, in step S102, assume that there are n radars with the same system in the physical space to be calculated, and they are independent of each other. Then their joint detection probability is In the rectangular coordinate system, assume that the coordinates of the radar are (x i , y i , 0), and the coordinates of the target are (x, y, z). Then the multi-radar joint detection probability can be expressed as

[0022] Furthermore, the specific operations of step S103 include the following steps

[0023] S1031: Use the dissection grid to divide the multi-radar detection area into regional grids, and regard the standard volume element as the size of the dissection grid equatorial standard volume element corresponding to the required grid granularity under the current parameters of the aircraft;

[0024] S1032: Use the volume element displacement operation method in the binary three-dimensional identification data of the dissection coding and the spatial volume element relationship calculation method to respectively obtain the three-dimensional dissection volume element quantity difference between the current position coordinate coding and the radar location coding, and then multiply by the length, width, and height of the basic volume element to obtain the coordinate difference in the multi-radar joint detection probability calculation formula in step S102;

[0025] S1033: Use the coordinate difference obtained in step S1032, combined with the multi-radar joint detection probability calculation formula in step S102, to calculate the radar detection probability at this granularity level.

[0026] Furthermore, in step S104, use the formula to calculate the detection probability size corresponding to the dissection volume element at the new granularity level. In the formula, code_L represents the coding identification corresponding to the dissection volume element block when the current granularity level is L; code_M represents the coding identification corresponding to each dissection volume element when dividing the airspace environment at the M level; P represents the joint radar detection probability value corresponding to this dissection coding; P code_L represents the probability that the aircraft is detected at least once when moving in the standard volume element at the L level, that is, the radar detection probability when the aircraft moves in the standard volume element at the L level; 8 L-MIt represents how many M-level dissection volume elements a basic dissection volume element at the L level under the three-dimensional octree organization consists of; j represents the corresponding volume element sequence; T represents the number of volume elements containing detection probability information among the M-level dissection volume elements included in the current basic dissection volume element structure at the L level.

[0027] Furthermore, the specific operations of step S2 include the following steps

[0028] S201: Set the constraint conditions of the aircraft and determine the corresponding grid level according to the required granularity;

[0029] S202: Determine the starting point and the ending point, and establish the open set and closed set data tables;

[0030] S203: Take the starting point of the aircraft as the parent node, and then determine the child nodes according to the change of the dissection coding bit based on the idea of variable step size for each segment, and calculate the f(n) value of each child node and put it into the open set;

[0031] S204: Find all child nodes corresponding to the minimum f(n) in the open set;

[0032] S205: If there are multiple child nodes corresponding to the minimum f(n) in the open set, continue to perform the secondary selection of the forward direction of the child nodes to find a better forward direction and its corresponding child nodes.

[0033] Furthermore, the specific operations of the method of searching all child nodes with variable step size for each segment in step S203 include the following steps

[0034] S2031: Calculate the heuristic cost factor h(n') of the parent node, and set the judgment value D in combination with the task background;

[0035] S2032: When the value of h(n') is less than the value of D, search all child nodes using a small step size and high-level grid;

[0036] S2033: When the value of h(n') is greater than the value of D, search all child nodes using a large step size and low-level grid.

[0037] Furthermore, the calculation method of f(n) of the child node in step S203 is

[0038]

[0039] In the formula, P code_L represents the radar detection probability when the aircraft is active in the standard volume element at the L level, and P(n) is the probability of safe passage of the corresponding child node.

[0040] Furthermore, the specific operations of the secondary selection of the forward direction of the child node in step S205 include the following steps

[0041] S2051: In the spatial division based on the spatial octree organizational structure, find the smallest octree structure that contains the target and the parent node, and use the formula to evaluate the risks of the eight sub-blocks. In the formula, q represents the eight sub-blocks of the smallest octree grid that contains the current position and the target position, and q = 1, 2, …, 8; N obstacle , N total represent the total number of threat voxels and the total number of all smallest granularity voxels that should be contained in the smallest granularity voxels respectively in the eight sectional voxel blocks at the current level under the smallest octree structure that contains the target and the parent node during the voxel modeling at the smallest granularity level;

[0042] S2052: In the smallest octree structure that contains the target and the parent node, select the sub-block with the largest P sparsity (q) value as the feasible direction, that is, the target point of the aircraft flight path.

[0043] The beneficial effects of the present invention are as follows:

[0044] 1. In the aircraft flight path planning method based on the sectional grid in the present invention, a new underlying representation framework of the three-dimensional sectional grid theory is adopted to effectively solve the problem of repeated modeling in the same space caused by the changing requirements of the spatial reference distance brought by different aircraft in the same space. Based on this theory, the traditional A* algorithm is improved to form a more adaptable, more general, efficient and autonomous flight path planning method to assist the pilot in making flight path decisions.

[0045] 2. In the present invention, based on the sectional grid, the multi-radar detection area environment is modeled to represent the flight physical space of the aircraft, and the radar threat area representation ability with multi-level association, adjustable levels and reusable repetition can be formed. Different-speed aircraft can be effectively integrated into the same method and framework for associated calculation, and the problem of repeated modeling of threat areas when different aircraft or different speeds of aircraft are involved can be effectively solved.

[0046] 3. In the present invention, the traditional A* algorithm is improved. First, in terms of improving the node search efficiency, a segmented variable step size search design is proposed to reduce the number of node calculations; second, the calculation method of the node f(n) is improved, and the detection probability of the child node is associated with the flight path search through risk controllable design, so that the entire algorithm search can consider passing through areas with less threat and balance between "risk" and "flight path length"; third, when the aircraft faces multiple feasible child node forward directions, the feasible directions are secondarily selected according to the sparsity to further reduce the number of invalid searches and node calculations. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 It is a schematic diagram of the radar threat area with a detection probability greater than 0.1 in the rectangular coordinate system for the present invention.

[0048] Figure 2 This is a schematic diagram of the grid quantity relationship based on the subdivision grid of the present invention.

[0049] Figure 3 This is the three-view diagram of the radar detection threat at the 13th level based on the subdivision grid of the present invention.

[0050] Figure 4 This is a schematic diagram of the level association of the radar detection threat at the 13th level and the 12th level based on the subdivision grid of the present invention.

[0051] Figure 5 These are the calculation elements of the Manhattan distance under the dissection organizational framework of the present invention.

[0052] Figure 6 This is a schematic diagram of the comparison of the dissection coding bits in the variable step size by segments of the present invention.

[0053] Figure 7 This is the applicable scenario of the horizontal plane sparsity of the present invention.

[0054] Figure 8 This is the corresponding relationship between the calculation of the horizontal plane sparsity and the moving direction of the present invention.

[0055] Figure 9 These are the feasible directions with the same cost in space of the present invention.

[0056] Figure 10 These are the sparsity calculation and corresponding directions in space of the present invention

[0057] Figure 11 This is a schematic diagram of the distribution of the radar joint detection threat area in the simulation experiment of the present invention.

[0058] Figure 12 These are the experimental results of the traditional method and the power exponents of 1 and 3 when n = 30 in Table 6 of the simulation experiment of the present invention.

[0059] Figure 13 These are the implementation effects of the algorithm under different radar threat area distribution scenarios when n = 30 in the simulation experiment of the present invention.

[0060] Figure 14 This is an application example of the improved algorithm when β = 1.5 in the simulation experiment of the present invention.

[0061] Figure 15 This is an application example of the improved algorithm when β = 1 in the simulation experiment of the present invention. Specific implementation manners

[0062] In order to enable those of ordinary skill in the art to better understand the technical solution of the present invention, the technical solution of the present invention will be further described below with reference to the accompanying drawings and embodiments.

[0063] An aircraft trajectory planning method based on a meshed grid, comprising the following steps,

[0064] S1: Based on the meshed grid, perform multi-radar detection area environmental modeling to characterize the flight physical space of the aircraft;

[0065] S2: In the physical space organized based on the meshed grid in step S1, use an improved A* algorithm to perform aircraft trajectory planning.

[0066] Specifically, the specific operations of the multi-radar detection area environmental modeling in step S1 include the following steps,

[0067] S101: Calculate the single-radar detection probability;

[0068] S102: According to the single-radar detection probability calculated in step S101, calculate the multi-radar joint detection probability in the multi-radar detection area;

[0069] S103: Based on the meshed grid organization, calculate the radar detection probability at any granularity level according to the multi-radar joint detection probability;

[0070] S104: Use the correlation relationship between different granularity levels under the meshed grid organization to calculate the detection probability corresponding to the meshed volume element at the new granularity level;

[0071] S105: Based on the detection probability values in the corresponding meshed grid volume element, perform visualization modeling on the radar detection area based on the meshed grid.

[0072] In the present invention, the method for calculating the single-radar detection probability is as follows:

[0073] S1011: Represent the classical single-radar detection probability as In the formula, respectively represent the average detection probability and the average signal power, N represents the receiver noise power, and b and n respectively represent the detection threshold voltage value and the number of integrated pulses;

[0074] S1012: To form a probability calculation method adapted to grid calculation, correlate and calculate the standard aircraft radar cross-section area and the actual target radar cross-section area. Assume that the cross-section area of target i is the standard cross-section area, and the cross-section area of target m is the actual cross-section area, respectively represent the single-radar average signal power corresponding to i and m. Then, in the same physical space, with the same background noise and propagation conditions, and the same radar performance, the detection probabilities corresponding to i and m can be respectively represented as and Taking the inverse of the natural logarithm of both equations and then dividing them, we can obtain In the formula, and respectively represent the detection probabilities corresponding to i and m;

[0075] S1013: According to the radar basic equation it can be obtained that in the same radar and the same physical space (with the noise level N unchanged), In the formula, is the signal-to-noise ratio, is the maximum average radar cross section of the target aircraft, R represents the Euclidean distance between the target and the radar, and K 0 is a constant related to factors such as the radar transmission power, antenna gain, and electromagnetic wave propagation environment; in the same physical space, these values can be considered unchanged.

[0076] S1014: When the detection probability of a fixed-parameter radar is set to 0.1, the corresponding target radar cross section and maximum detection distance are represented by σ mc and R mc respectively, then Ideally, in the selected physical space, for a detection radar located on the ground with known performance, its A value is unchanged and is a known number.

[0077] S1015: Substitute in step S1014 into in step S1013, and the calculation formula for the single-radar detection probability can be obtained as According to this calculation formula, it can be known that in the physical space, the average radar cross section of the aircraft the Euclidean distance R between the radar and the target i and the detection probability of the radar for the target at the corresponding position.

[0078] On this basis, the calculation method for the combined detection probability of multiple radars in the multi-radar detection area is as follows:

[0079] Assume that there are n radars with the same system in the physical space to be calculated, and they are independent of each other, then their combined detection probability is In the rectangular coordinate system, assume that the coordinates of the radar are (x i , y i , 0), and the coordinates of the target are (x, y, z), then the combined detection probability of multiple radars can be expressed as

[0080] Based on this calculation formula for the combined detection probability of multiple radars, assume that in 128×128×32 km 3In the physical space, there are 3 radars randomly distributed with the same parameters. When the detection probability of the radar is 0.1, let the value of A be 1 and the maximum detection distance be 30 km. The radar threat area with a detection probability greater than 0.1 in the rectangular coordinate system is shown in the appendix Figure 1 (Let L_grid represent the side length of the basic grid volume element, which is 1 km here), where the coordinates point to the lower left point of the coordinate, such as (1, 2, 1), and the color-coded diagram on the right represents the size of the detection probability value and the corresponding color.

[0081] On this basis, based on the hierarchical grid organization, the specific operations for calculating the radar detection probability at any granularity level include the following steps.

[0082] S1031: Use the hierarchical grid to divide the multi-radar detection area. According to the formula It can be known that assuming the values of A and are known, obtaining the three-dimensional coordinate difference in the rectangular coordinate system is similar to the solution of the Manhattan distance in the A* algorithm, and the target area divided by the hierarchical grid does not involve the influence of spatial curvature. Therefore, the standard volume element is regarded as the size of the equatorial standard volume element of the hierarchical grid corresponding to the grid granularity required under the current parameters of the aircraft.

[0083] S1032: Use the three-dimensional binary identification data of the hierarchical encoding and the volume element displacement operation method in the spatial volume element relationship calculation method to respectively obtain the three-dimensional hierarchical volume element quantity difference between the current position coordinate encoding and the radar position encoding, and then multiply by the length, width, and height of the basic volume element to obtain the coordinate difference in the multi-radar joint detection probability calculation formula in step S102, as shown in the appendix Figure 2 shown;

[0084] S1033: Use the coordinate difference obtained in step S1032 and combine it with the joint detection probability calculation formula to calculate the radar detection probability at this granularity level.

[0085] On this basis, when it is necessary to represent the radar threat area using the detection probability at different levels, the formula can be used to calculate the detection probability corresponding to the hierarchical volume element at the new granularity level. In the formula, code_L represents the encoding identifier corresponding to the volume element block at the current granularity level L; code_M represents the encoding identifier corresponding to each volume element at the M level when dividing the airspace environment; P represents the joint radar detection probability value corresponding to this hierarchical encoding; P code_L represents the probability that the aircraft is detected at least once when moving in the L-level standard volume element, that is, the radar detection probability when the aircraft moves in the L-level standard volume element; 8 L - MIt represents how many M-level dissection volume elements a basic dissection volume element at the L level under the three-dimensional octree organization consists of; j represents the corresponding volume element sequence; T represents the number of volume elements containing detection probability information among the M-level dissection volume elements included in the current basic dissection volume element structure at the L level.

[0086] There are two bases for this formula representation: one is based on the expression structure of "the actual space divided by the dissection grid - the unique identifier of the dissection code - the spatial inclusion information" formed by the dissection grid. On the basis of the corresponding coding relationship, the actual space and the corresponding information are in one-to-one correspondence; the other is to ensure that the exposure risk of the aircraft is controllable. When using the dissection grid to characterize the radar threat area, the probability of being detected by the radar in the entire dissection volume element should be indicated. This requires multiplying the inverted detection probabilities of each fine-grained level and then subtracting the result from 1 to represent the detection probability corresponding to the dissection volume element at the new granularity level.

[0087] On this basis, step S105 uses a dissection grid framework for visual modeling of the radar detection area based on the dissection grid. It is secondarily developed relying on the Cesium open-source platform in visualization to form a data interaction display platform, which can support the dissection calculation and representation of the physical space and the radar threat area in the flight path planning. The Cesium platform is a map engine written in JavaScript language using webGL rendering, used for map display in 3D, 2D, and 2.5D forms. After secondary development, it can interactively display the dissection grid division structure of the local area and the visual representation of the radar threat area. The software and hardware experimental environment is shown in Table 1 below.

[0088] Table 1 Operating Environment and Software and Hardware Configuration

[0089]

[0090] For the purpose of experimental demonstration, the parameters such as the airspace range and the number of radars are selected as shown in Table 2 below.

[0091] Table 2 Parameter Setting

[0092]

[0093] At the same time, due to computer performance limitations, as shown in the appendix Figure 3 where (a) is the top view, (b) is the front view, and (c) is the side view; the 13th-level dissection grid is selected as the bottom layer and the three views of the combined radar detection threat in the local airspace environment are shown based on the dissection grid framework and the above parameter information. The red depth represents the magnitude of the detection probability value, and the darker the color, the greater the detection probability value.

[0094] According to the different-level detection probability association method in step S104, in the appendix Figure 4In the middle, the 12th-level grid is used as a comparative example to show the grid-level association under the same environmental division. Among them, (a) is the schematic diagram of the 12th-level grid division, and (b) is the schematic diagram of the 13th-level grid division.

[0095] Furthermore, on the basis of the above-mentioned meshed grid, the specific operations for aircraft flight path planning using the improved A* algorithm include the following steps.

[0096] S201: Set the constraint conditions of the aircraft, and determine the corresponding grid level according to the required granularity.

[0097] To plan the flight path in the physical space organized by the meshed grid, the constraint conditions of the aircraft are set as follows:

[0098] (1) Aircraft speed constraint: Assume that the aircraft moves forward at a constant speed V during the mission time period. At the same time, in a certain mission under the meshed grid framework, the time step T for each type of aircraft to enter the next adjacent grid is the same. When using length to represent the minimum side length of the standard volume grid of a certain level, and the subscript such as L + 1 represents the (L + 1)th level, then the relationship among the three is length L ≥V×T≥length L+1 . In particular, when the aircraft type changes or its own flight speed changes, only the V value needs to be updated to find the next suitable level grid. For only one-time environmental modeling, there is also V≥V min , and find the minimum speed requirements for each type of aircraft or aircraft to be planned during the entire mission.

[0099] (2) Aircraft position expression constraint under appropriate grid granularity: When expressing the position and calculating the detection probability of the same aircraft at different levels, assume that the aircraft is a particle and is always located at the center position of the grid at that level.

[0100] (3) Aircraft climb angle / dive angle and minimum straight flight step convention: To adapt to the characteristics of grid movement, assume that the grid size in the simulation always satisfies the aircraft speed constraint in (1), and the aircraft can search and move forward in the airspace grid with a 26-neighborhood in three-dimensional space.

[0101] (4) Definition of Manhattan distance in the physical space organized by the meshed grid: As shown in the appendix Figure 5 shown, in the appendix Figure 5Among them, (a) are the three elements for calculating the Manhattan distance in the two-dimensional horizontal plane, and (b) are the three additional elements for calculating the Manhattan distance in the three-dimensional airspace. Based on the grid volume element near the equatorial surface at each level as the basic volume element, the calculation of the Manhattan distance is closely related to the size of the basic volume element at the corresponding level. In the horizontal two-dimensional plane, the Manhattan grid spacing is obtained by multiplying the ratio of its length, width, and diagonal length (with the length as the base) by 10 and then taking the floor value. In the three-dimensional airspace, the height dimension is added, and the ratio of the oblique diagonal length to the height side in different cases is multiplied by 10 and then the floor value is taken as the Manhattan grid spacing. Subsequently, according to the formula shown, when granularity conversion is required at different levels, the ratio of the lengths of the basic volume elements at the corresponding levels is converted into a cost ratio relationship for correlation relationship operations. Among them, radio L_M represents the ratio relationship of the M type edges of the dissection grid at the Lth level, and length L represents the length of the minimum side length of the regular volume element at the Lth level.

[0102] (5) Ratio relationship of grid volume elements in the simulation: As shown in Table 3 below, the range of the shortest side lengths at levels 7 - 18 is 217m - 44520m, which can meet the flight distance step constraints of different aircraft. At the same time, the side length relationships of the basic grid volume elements within these levels are approximately equal to 1:1:1. Therefore, in the simulation, the basic volume element is set as a cube volume element for algorithm simulation analysis.

[0103] Table 3 Comparison of the theoretical size and actual size of volume elements at levels 7 - 18 of the GeoSOT - 3D dissection framework

[0104]

[0105] (6) Convention on the maximum flight height of the aircraft: Taking a fighter jet as an example, the empirical value of its radar cross - section area is about 1 - 2m 2 , and the conventional flight height is 9 - 30km. Therefore, it is assumed that the maximum flight height of the aircraft in the present invention is 30km.

[0106] From the above - mentioned constraints, it also shows that the method in the present invention weakens the influence of the physical limitations of the aircraft on the algorithm search process. The core lies in improving the algorithm starting from the principle of the new underlying representation organization framework, so as to adapt to the requirements of mission - level trajectory planning and the characteristics of the aircraft penetration background.

[0107] S202: Determine the starting point and the ending point, and establish the open set and closed set data tables; The operation of this step is the same as that of the traditional A* algorithm, and will not be elaborated in the present invention.

[0108] S203: Take the starting point of the aircraft as the parent node, and then determine the child nodes according to the idea of segmented variable step size through the change of dissection coding bits, and calculate the f(n) value of each child node and put it into the open set;

[0109] In the context of mission-level trajectory planning for penetration, there are generally two characteristics: (1) The physical space involved is large, the distance from the takeoff point to the target point is far, and the span is wide; (2) To protect important targets, the radar threat areas are mainly distributed near the target point and are relatively concentrated. Based on the above two points, the present invention divides the global trajectory planning into two flight segments. The first flight segment is the cruise segment far from the target point, where the radar threat areas are few and sparsely distributed; the other flight segment is the penetration segment close to the target point, where the number of threats is large and the distribution is concentrated. This division is judged based on the magnitude of the heuristic cost factor h(n), and a judgment value D is set in combination with different mission backgrounds for variable step-size design in segments.

[0110] Calculate the heuristic cost factor h(n') of the parent node. When the value of h(n') is less than the value of D, it is considered that the aircraft enters the penetration segment where the radar threat areas are relatively concentrated, and a small step-size high-level grid is used to search all child nodes; when the value of h(n') is greater than the value of D, a large step-size low-level grid is used to search all child nodes. It should be noted here that in the traditional A* algorithm, the final movement cost value f(n) of the child node consists of two parts, f(n)=g(n)+h(n), where g(n) represents the actual movement path cost value spent from the starting point to the child node, and its calculation method is to add the actual movement path cost value accumulated from the starting point to the current point and the actual movement cost value spent from the current point to the child node; h(n) has nothing to do with the starting point and is the cost estimation value of the current position point to the target. Usually, there are Manhattan distance, Euclidean distance, and calculation methods combining the weights of the two. In this application, the calculation method of h(n') can adopt the existing Manhattan distance, Euclidean distance, or calculation methods combining the weights of the two, and will not be elaborated in the present invention.

[0111] Although the influence of curvature is not considered in the present invention, in combination with the actual situation, in order to avoid excessive changes in the inclusion relationship of the front and rear spaces caused by the influence of curvature, adjacent levels are uniformly used for variable step-size operation in segments, as shown in the appendix Figure 6 As shown, where (x, y, z) abstractly represent the coding arrays of a certain level in the latitude dimension, longitude dimension, and altitude dimension respectively.

[0112] In the appendix Figure 6 From left to right, the repeated (x, y, z) respectively represent the corresponding encodings in the longitude, latitude, and altitude dimensions of each level from the low level to the high level. In the bit comparison operation shown in the figure, the parent node encoding and the child node encoding respectively represent the current grid and the next neighborhood subdivision encoding under the 26-neighborhood search of the airspace. Assuming that the entire node encoding has n bits, when the aircraft uses a low-level large step-size grid encoding bit for comparison in the cruise segment, only the encoding comparison bit needs to be raised to the comparison of the (n - 3)-bit encoding, and its operation method using the change of the encoding bit is simple and easy.

[0113] Compared with the traditional search method, under the octree subdivision organizational structure, in the cruise stage, the improved algorithm can change from searching at least two small grids in the original eight small grids to searching a large octree grid, reducing the search resource consumption by more than half. Moreover, in the process of approaching the target point, due to the large-step search, the number of computing nodes in the radar threat area with a relatively simple distribution can be reduced, and it can quickly jump out of the meaningless local dilemmas. When the aircraft enters the penetration stage, small-step high-level grids are used for subdivision coding and comparison up to the nth bit to find a more reasonable feasible flight path in the scenario with dense threat distribution. The segmented variable-step design based on the characteristics of the subdivision grid can make the global flight path planning method of the aircraft improve efficiency while ensuring the rationality of the flight path.

[0114] On this basis, in order to ensure the rationality of the flight path, the detection probability values within the voxel of the subdivision grid corresponding to the child nodes are further integrated into a risk control coefficient, making the flight path selection of the aircraft closely related to the detection risks contained in the child nodes. On this basis, a risk-controllable design that can be self-defined is formed through the power exponent control factor acting on this coefficient. According to the detection probability of the current-level child nodes, the calculation method of f(n) of the child nodes in the present invention is

[0115]

[0116] In the formula, P code_L represents the radar detection probability when the aircraft is active in the standard voxel at the L level, and P(n) is the probability of safe passage of the corresponding child node.

[0117] Because P code_L , P(n) are radar detection probability values, both within the interval [0,1], and they cannot achieve good results when used as weights to affect the value of g(n). Therefore, a coefficient structure of 1 - ln(Pn) is designed. It takes the logarithm of the probability of safely passing through the child node and then subtracts it from 1 to obtain a positive feedback effect on the actual movement cost g(n), and at the same time expands the acting interval representation interval to [1, +∞). Finally, a β-th power is added outside 1 - ln(Pn) as a risk control factor to control the coefficient size, for the commander or pilot to further set the comprehensive influence degree of the detection threat on the final flight path of the aircraft in "exposure risk" and "flight path length" according to the mission requirements. Among them, the larger the β value, the more the path searched by the algorithm tends to be a safer (detoured) path.

[0118] S204: In the radar threat area of the subdivision grid organization, when the aircraft searches with a 26-neighborhood grid in the airspace, when a child node encounters a threat, there are usually multiple feasible child nodes with the same and minimum f(n) values in the open set. Therefore, it is necessary to first find all the child nodes corresponding to the minimum f(n) in the open set;

[0119] S205: If there are multiple child nodes corresponding to the minimum f(n) in the open set, continue to perform a secondary selection of the advancing direction of the child nodes to find a better advancing direction and its corresponding child nodes.

[0120] Specifically, S2051: In the space division based on the spatial octree organizational structure, find the minimum octree structure that contains the target and the parent node, and use the formula to evaluate the risks of the eight sub-blocks. In the formula, q represents the eight sub-blocks of the minimum octree grid that contains the current position and the target position, q = 1, 2,..., 8; N obstacle , N total represent the total number of threat voxels and the total number of all minimum-granularity voxels that should be contained in the minimum-granularity voxels respectively included in the eight sectional voxel blocks at the current level under the minimum octree structure that contains the target and the parent node during the minimum-granularity hierarchical grid voxel modeling.

[0121] Assume that the octree structure level that contains the current position and the target point is L, then N total = 2 L-1 ; and N obstacle can be quickly found through coding comparison, so as to relatively simply calculate the sparsity. By finding the direction with less potential risk among multiple feasible directions, the number of calculation nodes can be effectively reduced, and ineffective search can be avoided.

[0122] S2052: It can be seen from the formula that P sparsity (q) represents the sparsity of the eight sub-voxel blocks under the minimum octree structure. The larger its value, the fewer the number of grids with threats in the sub-blocks divided by the minimum-granularity grid size. Therefore, in the minimum octree structure that contains the target and the parent node, select the sub-block with the largest P sparsity (q) value as the feasible direction, that is, the target point of the aircraft flight path.

[0123] In the present invention, the scenarios that require sparsity judgment in the horizontal plane and the three-dimensional airspace are sorted out. The horizontal plane sparsity judgment scenarios are as shown in the appendix Figure 7 shown. In the appendix Figure 7 , the three scenarios (a), (b), and (c) respectively represent the feasible child nodes and the corresponding feasible directions when the current position and the target point are on the same vertical line, on the same horizontal line, and on the same oblique line. Among them, the red grid represents the current position point; the blue grid represents the target position point; the white grid represents the passable grid without detection threat; the gray grid represents the grid with detection threat. Taking the example that the target point and the current position are both in the Northern Hemisphere and the target point is to the right or upper right of the current position, as shown in the appendix Figure 8As shown (without indicating obstacle conditions), where the sparsity in the downward movement direction is determined by the sparsity of the voxel block where the current position is located; for the right, upward, and upper-right sparsity correspondence, it is consistent with the sparsity calculation position in (a) of Attachment Figure 8 .

[0124] In the legend shown in (b) of Attachment Figure 8 , the reason for only considering the neighborhood sparsity calculation in four directions is as follows: Combining the three scenarios shown in Attachment Figure 7 , in the Manhattan distance calculation method, when the target position is to the right or upper-right of the current position, there are only the four listed cases where the total movement cost f(n) of two feasible child nodes is equal. Therefore, the sparsity calculation only needs to judge based on the calculation results of these four directions. And when the feasible direction scenarios shown in Attachment Figure 7 appear in different vertical planes or upper horizontal planes, similar to mapping the respective P sparsity (n) values to the corresponding positions, which will not be elaborated here.

[0125] In the 26-neighborhood search in space, when the current position and the target position are not in the same plane, the direction situations that need to be considered are as shown in Attachment Figure 9 , where (a) shows the feasible directions with the same vertical position and f(n) cost, and (b) shows the feasible directions with the same oblique position and f(n) cost; in the figure, the gray grids represent the areas with radar detection threats, and the blue grid represents the target position. Except for the relationships shown in the figure, the rest can be simplified to relationships similar to the two-dimensional plane structure. Subsequently, as can be seen from the relationships shown in Attachment Figure 10 (in Attachment Figure 10 , (a) is the corresponding position of the sparsity, (b) is the corresponding direction of the vertical plane, and (c) is the corresponding direction of the inclined plane)), only need to map the position relationship of the sparsity in the octree structure to the corresponding feasible directions, and the logical structure remains unchanged.

[0126] Taking the encoding correspondence relationship in the Northern Hemisphere as an example, the corresponding position of the sparsity and the direction correspondence relationship under the octree are as shown, where the direction correspondence relationship of the horizontal plane is as shown in (b) of Attachment Figure 8 . And in Attachment Figure 10 , the directions corresponding to the sparsity shown in the spatial octree in the vertical plane and the inclined plane do not involve the direction judgment away from the target point. This is because theoretically, it can be inferred that among the 26 child nodes of the same parent node, moving to the child node away from the target point will only increase the corresponding final movement path cost f(n). And in the introduction of the basic principle of the A* algorithm, it can be found that in this case, the final costs of each child node cannot be the same, and according to the distribution characteristics of the radar-like spherical or hemispherical threat areas, it can be known that there are always feasible child nodes moving towards the target point.

[0127] Simulation experiment:

[0128] Generally, in the scenario where the distribution of the radar threat area is relatively complex, it is more meaningful for the aircraft to attempt to cross the risk area. To effectively compare the algorithm effects, in this simulation experiment, a cubic grid with a side length of 10 km is adopted, and in the scenario where the number of radars is between 10 and 30, the traditional A* algorithm (hereinafter referred to as the traditional algorithm) and the A* algorithm improved based on the dissected grid (hereinafter referred to as the improved algorithm) are compared and analyzed. The operating environment and software configuration of the experimental host are shown in Table 4 below.

[0129] Table 4 Simulation environment and software and hardware configuration

[0130]

[0131] Combined with the task background and the need for algorithm comparison and analysis, the regional conditions are set as shown in Table 5 below, where R represents the maximum detection radius of the radar that meets the conditions, and P represents the size of the combined detection probability of the radars.

[0132] Table 5 Regional condition settings

[0133]

[0134] Based on the above conditions, assuming that when the detection probability of the radar is 0.1, the maximum cross-sectional area of the aircraft is 1 m 2 , and the maximum detection range is 30 km. In this simulation experiment, it is assumed that aircraft with the same cross-sectional area cross the grid within the same time step T, and no further aircraft performance condition settings are made. At the same time, to compare the comprehensive performance of the traditional algorithm and the improved algorithm, the target point is set not to be covered by the radar threat in the simulation scenario, so that the traditional A* algorithm under the obstacle avoidance strategy can also find a feasible flight path. When the edge length of the volume element is 10 km and assuming n = 20, the distribution of the radar threat area is as shown in the appendix Figure 11 as follows.

[0135] Effect and comparative analysis of the improved A* algorithm

[0136] To reflect the degree of the final flight path exposure risk of the improved algorithm in the present invention, in addition to the three evaluation indicators of the planning time, the length of the final flight path, and the number of closed set nodes calculated, a risk coefficient is introduced to describe the final flight path risk degree of the improved algorithm. It is assumed that the risk coefficient is calculated by accumulating the combined detection threat probabilities of the radars for each grid node passed through in the final flight path, and its calculation method is In the formula, d(n) represents the magnitude of the accumulated detection probability risk value when the aircraft moves to the current position n, represents the probability that the aircraft is detected at least once when moving in the L-level standard volume element, and j represents its sequence position in the set of final flight path nodes.

[0137] 1. Analysis of the change in the magnitude of β power

[0138] To compare and analyze the influence of the magnitude of power β on the algorithm under different threat area distribution scenarios, in this simulation experiment, three values of β = 1, 2, and 3 are selected, and the number of randomly distributed radars n is 10, 20, and 30 for experiments, and then a comparative analysis is carried out with the traditional algorithm. The main results are shown in Table 6 below.

[0139] Table 6 Experimental results of the power change algorithm under different environmental complexities (L_grid = 10km)

[0140]

[0141] It should be noted that in the above data, the threat area distribution scenarios corresponding to different numbers of radars are the same, and there is no relationship between the threat area distributions of different numbers of radars. In different radar threat area distribution scenarios, the calculation time of each method is mainly positively correlated with the number scale of the closed set calculation nodes, and there is no inevitable correlation with the number of radars. Therefore, an increase in the number of radars does not necessarily lead to an increase in the method running time, because it is also affected by the radar threat area distribution.

[0142] According to the principle of the improved algorithm, the smaller the power value, the greater the possibility that the aircraft chooses to directly cross the risk area without detouring, and usually it can search for the end point faster, thus reducing the closed set calculation nodes, so the planning time is correspondingly reduced. However, from the performance of the improved algorithm with different powers when n = 10, it can be seen that there is a situation where the aircraft does not quickly find the target point when the power exponent is small, but instead calculates more nodes. Therefore, the inference about the influence of the power value on the aircraft planning efficiency is not absolute.

[0143] From the data comparison when n = 20 in the table, it can be seen that when β = 3, the final trajectory cost value f(n) of the improved algorithm is the same as that of the traditional algorithm, and the risk coefficient is 0. At the same time, although the number of closed set calculation nodes of the improved algorithm has decreased by nearly half, the algorithm planning time has not been reduced by half. This shows that although the improved algorithm reduces the calculation nodes, the increased algorithm logic structure makes the improvement effect of the algorithm planning efficiency not obvious. By comparing the results when n = 20, β = 1 and β = 3, it can be found that as the power increases, the trajectory tends to a farther and safer path, and the power adjustment effect takes effect. And the consistent experimental results between β = 1 and β = 2 can also lead to the conclusion that the power value should be reasonably set for different scenarios to play a role.

[0144] From the data comparison of n = 10 and n = 30 in the table, it can be seen that in the case of a more complex distribution of the radar threat area, the planning time of the improved algorithm is better than that of the traditional algorithm, and it is more adaptable to the scenario of a complex distribution of the radar threat area. To sum up, it can be seen that the improved algorithm in the present invention can find a shorter flight path under a certain exposure risk, and the power change adjustment function takes effect. At the same time, because the logical structure of the improved algorithm in the present invention is more complex, there is a phenomenon that the planning time is slightly higher than that of the traditional algorithm in the scenario of a relatively simple radar threat area distribution. However, in a more complex scenario, the improved algorithm in the present invention can achieve better results.

[0145] Appendix Figure 12 shows the results of the traditional method and the power exponents of 1 and 3 in Table 6 when n = 30. Among them, the final paths of the power exponents β = 1 and β = 2 are the same, only the number of closed-set nodes has a small difference, and they are not drawn again here. From the Figure 12 number of colored circles representing the closed-set nodes in (a) of the appendix, it can be seen that the number of closed-set nodes calculated by the traditional method is large and the calculation range involved is wide. (b) and (c) respectively show the final flight paths obtained by the improved algorithms with β = 1 and β = 3 and the closed-set nodes calculated by them. It can be found from Figures (b) and (c) that the final flight paths of the two are different, and most of the closed-set nodes are concentrated near the obstacles, and the numbers are similar, and the distribution of the closed-set nodes calculated by the improved algorithms is more reasonable than that of the traditional method.

[0146] 2. Applicability of the improved algorithm in the scenario of complex distribution of the radar threat area

[0147] To avoid special cases in the results shown in the appendix Figure 12 and Table 6, and further analyze and verify the applicability of the improved algorithm in the scenario of complex distribution of the radar threat area. In this simulation experiment, with the conditions of β = 1 and n = 30, the effect of the improved algorithm was randomly tested in four different scenarios of radar threat area distribution. The experimental results are shown in the appendix Figure 13 and Table 7 below. In the appendix Figure 13 , (a), (b), (c), and (d) respectively correspond to Scenario 1, Scenario 2, Scenario 3, and Scenario 4 in Table 7; among them, the yellow line segment represents the final flight path of the improved algorithm, the red line segment represents the final flight path of the traditional algorithm, and the colored circles represent the closed-set calculation nodes of the improved algorithm. The number of closed-set calculation nodes of the traditional algorithm is too large and is not shown here.

[0148] Table 7 Experimental results of the algorithm under different scenarios of radar threat area distribution when n = 30

[0149]

[0150] From Table 7 and the appendix Figure 13As can be seen from the results, in the case of a relatively complex distribution of the radar threat area, compared with the traditional algorithm, the improved A* algorithm based on the dissected grid can more effectively plan a feasible flight path through the threat area, with fewer closed-set calculation nodes and high operating efficiency. Moreover, in the environment shown in Scenario 2, the traditional algorithm cannot find a feasible flight path. It can be found that the improved A* algorithm based on the dissected grid can better adapt to the background of aircraft penetration.

[0151] At the same time, according to the comparison of the number of closed-set nodes calculated in Table 7, it can also be seen that compared with the traditional algorithm, in a scenario with a more complex environment, the improved algorithm greatly reduces the number of calculation nodes required to find the target point, thereby reducing the running time of the algorithm. The segmented variable-step design and the secondary selection of the forward direction of the child nodes effectively play their roles. And theoretically speaking, the segmented variable-step design only reduces the consumption of computing resources by half in the long-distance cruise section. In the attached Figure 13 In the part where the circles are shown to be relatively sparse in the distribution. And in this part of the interval, the traditional algorithm will not consume excessive computing resources either. Therefore, it can be concluded that the secondary selection design of the forward direction of the child nodes that plays a role when encountering obstacles plays an effective role.

[0152] To sum up, compared with the traditional A* algorithm, in the scenario where the distribution of the radar threat area is more complex, the comprehensive performance of the improved A* algorithm combined with the characteristics of the dissected grid framework in this chapter is better. Taking the data in Table 7 as an example, under the condition of n = 30, the number of closed-set calculation nodes of the improved algorithm is reduced by 50% - 90%, and the corresponding running time is reduced by 18% - 59%. And this method can further adjust between "distance" and "risk" through the selection of the power, can effectively plan a feasible flight path with a controllable exposure direction, and is more adaptable to the task-level planning aircraft penetration background under different requirements.

[0153] 3. Application Examples Based on the Cesium Platform

[0154] Set the application scenarios as shown in Table 8. Since the higher the level of the dissected grid division, the number of grids increases exponentially, and experiments cannot be carried out under the limitation of computer performance. Therefore, in this simulation experiment, the 12th and 13th levels are selected, and the number of randomly distributed radars n is 7 for the experiment. And because the 13th level is directly selected as the minimum granularity level, assuming that the aircraft still strictly meets various constraint conditions, the relevant parameters of the aircraft are not set here. The longitude and latitude distances of the standard volume elements of the 12th and 13th level grids are the same, and the height distances differ by a factor of two with an error of 1m. And the height distance (the shortest side) of the 12th level is about 14km, and the height distance (the shortest side) of the 13th level is about 7km. It is assumed that the influence of grid curvature on the flight path planning method is not considered in the dissected space of 0 - 100km.

[0155] Table 8 Parameter Settings of Application Scenarios

[0156]

[0157] Since the algorithm application example under the Cesium platform is a one-time simulation calculation, caching, and display, it cannot effectively demonstrate the secondary decision-making design of the advancing direction of the sub-nodes of the improved algorithm. Therefore, in this simulation experiment, only the final trajectories obtained with two power values of β = 1 and 1.5 in the same environment are selected for comparative demonstration. The application example of the improved algorithm when β = 1.5 is as shown in the appendix Figure 14 shown, where (a) is the front view, (b) is the top view, and (c) is the side view; the application example of the improved algorithm when β = 1 is as shown in the appendix Figure 15 shown, where (a) is the front view, (b) is the top view, (c) is the side view, and (d) is the illustration of the effective crossing of the risk area.

[0158] As shown in the appendix Figure 14 shown, when the improved algorithm sets β = 1.5, the final trajectory it obtains avoids the radar threat area and finds an absolutely safe trajectory. And in the final trajectory shown in the appendix Figure 15 shown, when β = 1 is set, after comprehensively evaluating the risk of the threat area and the trajectory distance, the improved algorithm selects a trajectory that is closer to the threat area, and it crosses part of the radar threat area with relatively low detection risk. From the results of the two illustrations, it can be found that the adjustment factor β of "exposure risk" and "trajectory length" is effectively applied, and the improved A* algorithm based on the triangulation grid can plan a feasible trajectory with controllable exposure risk for the aircraft.

[0159] The above shows and describes the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited by the above embodiments. What is described in the above embodiments and the specification only illustrates the principles of the present invention. Without departing from the spirit and scope of the present invention, the present invention will have various changes and improvements, and these changes and improvements all fall within the scope of the present invention claimed. The scope of protection claimed by the present invention is defined by the appended claims and their equivalents.

Claims

1. Aircraft trajectory planning method based on subdivision grid, Characterized in that, It includes the following steps, S1: Based on the subdivision grid, conduct multi-radar detection area environment modeling to characterize the flight physical space of the aircraft; S2: In the physical space organized by the subdivision grid in step S1, use the improved A* algorithm to plan the aircraft trajectory; The specific operations of step S2 include the following steps, S201: Set the constraint conditions of the aircraft and determine the corresponding grid level according to the required granularity; S202: Determine the starting point and the ending point, and establish open set and closed set data tables; S203: Take the starting point of the aircraft as the parent node, and then determine the child nodes according to the change of the subdivision coding bit based on the idea of variable step size for each segment, and calculate the f(n) value of each child node and put it into the open set; where f(n) is the final movement cost value of the child node; S204: Find all child nodes corresponding to the minimum f(n) in the open set; S205: If there are multiple child nodes corresponding to the minimum f(n) in the open set, continue to perform secondary selection of the forward direction of the child nodes to find a better forward direction and its corresponding child nodes; The specific operations of the method of searching all child nodes with variable step size for each segment in step S203 include the following steps, S2031: Calculate the heuristic cost factor h(n') of the parent node, and set the judgment value D in combination with the task background; S2032: When the value of h(n') is less than the value of D, use a small step size and high-level grid to search all child nodes; S2033: When the value of h(n') is greater than the value of D, use a large step size and low-level grid to search all child nodes; The calculation method of f(n) of the child node in step S203 is Wherein, P code_L represents the radar detection probability when the aircraft moves in the standard volume element at the L level, P(n) is the probability of safe passage of the corresponding child node; g(n) represents the actual moving path cost from the starting point to the child node, and g(n - 1) represents the actual moving path cost from the starting point to the (n - 1)-th child node; The specific operations of the secondary selection of the forward direction of the child node in step S205 include the following steps, S2051: In the spatial division based on the spatial octree organizational structure, find the smallest octree structure that contains the target and the parent node, and use the formula to evaluate the risks of the eight sub-blocks. In the formula, q represents the eight sub-blocks of the smallest octree grid that contains the current position and the target position, and q = 1, 2,..., 8; N obstacle , N total represent the total number of threat voxels and the total number of all the smallest granularity voxels that should be contained in the smallest granularity voxels respectively in the eight sectional voxel blocks at the current level under the smallest octree structure that contains the target and the parent node during the voxel modeling at the smallest granularity level; S2052: In the minimum octree structure containing the target and the parent node, select the sub-block with the largest P sparsity (q) value as the feasible direction, that is, the target point of the aircraft flight path.

2. The aircraft trajectory planning method based on subdivision grid according to claim 1, Characterized in that, The specific operations of the multi-radar detection area environment modeling in step S1 include the following steps, S101: Calculate the single radar detection probability; S102: According to the single radar detection probability calculated in step S101, calculate the multi-radar joint detection probability in the multi-radar detection area; S103: Based on the subdivision grid organization, calculate the radar detection probability at any granularity level according to the multi-radar joint detection probability; S104: Use the correlation relationship between different granularity levels under the subdivision grid organization to calculate the detection probability size corresponding to the subdivision volume element at the new granularity level; S105: Based on the detection probability value in the corresponding subdivision grid volume element, conduct visual modeling on the radar detection area based on the subdivision grid.

3. The aircraft trajectory planning method based on subdivision grid according to claim 2, Characterized in that, The specific operations of step S101 include the following steps, S1011: Express the classical single radar detection probability as where represent the average detection probability and average signal power respectively, N represents the receiver noise power, and b and n represent the detection threshold voltage value and the number of accumulated pulses respectively; S1012: Calculate the correlation between the standard aircraft radar cross-section and the actual target radar cross-section to obtain In the formula, assume that the cross-section of target i is the standard cross-section and the cross-section of target m is the actual cross-section, respectively represent the average single-radar signal power corresponding to i and m, and respectively represent the detection probabilities corresponding to i and m; S1013: According to the radar basic equation It can be obtained that In the formula, is the signal-to-noise ratio, is the maximum average radar cross section of the target aircraft, R represents the Euclidean distance between the target and the radar, and K 0 is a constant; S1014: When the detection probability of a fixed-parameter radar is set to 0.1, let the corresponding target radar cross section and maximum detection range be represented by σ mc and R mc respectively, then S1015: Substitute the in step S1014 into the in step S1013, and the calculation formula for the single radar detection probability is obtained as 4. The aircraft trajectory planning method based on subdivision grid according to claim 3, Characterized in that, In step S102, assume that there are n radars with the same system in the physical space to be calculated, and they are independent of each other. Then the combined detection probability is In a rectangular coordinate system, assume that the coordinates of the radar are (x i , y i , 0), and the coordinates of the target are (x, y, z). Then the combined detection probability of multiple radars can be expressed as 5. The aircraft trajectory planning method based on subdivision grid according to claim 4, Characterized in that, The specific operations of step S103 include the following steps, S1031: Divide the multi-radar detection area into regional grids using sub-divided grids, and regard the standard voxel as the size of the equatorial standard voxel of the sub-divided grid corresponding to the required grid granularity under the current parameters of the aircraft; S1032: Using the voxel displacement operation method in the binary three-dimensional identification data of sub-division coding and the spatial voxel relationship calculation method, calculate the three-dimensional sub-division voxel quantity difference between the current position coordinate coding and the radar location coding respectively, and then multiply by the length, width, and height of the basic voxel respectively to obtain the coordinate difference in the multi-radar joint detection probability calculation formula in step S102; S1033: Using the coordinate difference obtained in step S1032, combined with the multi-radar joint detection probability calculation formula in step S102, calculate the radar detection probability at this granularity level.

6. The aircraft trajectory planning method based on sub-divided grids according to claim 5, wherein: In step S104, the formula is used to calculate the detection probability corresponding to the dissected voxel at the new granularity level. In the formula, code_L represents the coding identifier corresponding to the dissected voxel block when the current granularity level is L; code_M represents the coding identifier corresponding to each sub-division voxel at the M level when dividing the airspace environment; P represents the joint radar detection probability value corresponding to this sub-division coding; P code_L represents the probability that the aircraft is detected at least once during its activities in the standard volume element at level L, that is, the radar detection probability of the aircraft during its activities in the standard volume element at level L; 8 L-M It represents how many M-level dissection volume elements a basic dissection volume element at the L level consists of under the three-dimensional octree organization; j represents the corresponding voxel sequence; T represents the number of voxels containing detection probability information among the M-level sub-division voxels included in the current L-level basic sub-division voxel structure.

Images