Urban low-altitude three-dimensional digital route network construction method and system

Abstract

The invention provides an urban low-altitude three-dimensional digital route network construction method and system, and relates to the technical field of low-altitude route construction, and the method comprises the steps: obtaining an airspace basic data set; performing airspace discretization on the airspace basic data set by adopting a discretization method; constructing a three-dimensional risk field; performing height layer division on a low-altitude airspace according to the type of an aircraft, and copying Delaunay network nodes in each height layer; generating an air route skeleton on the Delaunay network through potential energy gradient guidance; constructing a three-dimensional digital model based on the initial airway network, and establishing a digital twin system; adopting a multi-objective optimization algorithm to obtain an optimized air route network; performing path smoothing processing on the optimized air route network and verifying flight performance constraints; and updating the digital twin system to obtain the urban low-altitude three-dimensional digital airway network. According to the method, the urban low-altitude three-dimensional digital airway network which is reasonable in topological structure, high in generation efficiency, good in optimization performance and complete in digital expression can be constructed.

Description

Technical Field

[0001] This invention relates to the field of low-altitude airway construction technology, and in particular to a method and system for constructing a three-dimensional digital airway network for urban low-altitude airways. Background Technology

[0002] With the rapid development of UAV technology and urban air traffic, low-altitude aircraft are increasingly used in urban environments, and establishing a comprehensive urban low-altitude airway network has become a current research hotspot. Existing airway network construction methods mainly adopt technical approaches based on ground road network projection, obstacle Voronoi diagrams, and path search algorithms.

[0003] Chinese invention patent CN119516846A discloses a method for constructing a multi-level low-altitude airway network in complex urban environments. This method divides the 0-1000 meter low-altitude airspace into multiple altitude layers, assigns altitude layers to different types of aircraft, uses an octree structure for three-dimensional discretization of urban space, divides the airspace into a three-dimensional grid, designs a risk assessment function to evaluate the comprehensive risk of each grid unit, gradually constructs the airway network through a multi-level iterative approach, uses a hybrid iterative algorithm to plan and optimize airways at each altitude layer, and verifies path smoothing and flight performance using B-spline curves, ultimately generating a low-altitude airway network. While this patent achieves multi-level airspace division and multi-level airway network construction, its octree uniform grid discretization method lacks organic correlation with the spatial distribution of urban traffic demand. The multi-level iterative construction process requires generating multiple levels, such as skeleton airways, secondary airways, and feeder airways, resulting in low computational efficiency. Furthermore, the airway optimization stage mainly adopts single-objective or phased optimization strategies, making it difficult to achieve global collaborative balance among multiple optimization objectives. Summary of the Invention

[0004] In view of this, the present invention provides a method and system for constructing a three-dimensional digital airway network in urban low-altitude airspace. By improving the airspace discretization method, airway generation mechanism and optimization algorithm, it solves the technical problems in the prior art, such as insufficient correlation between airspace discretization results and traffic demand distribution, low efficiency of multi-level iterative construction, and lack of deep integration between multi-objective optimization algorithm and network topology. The invention constructs a three-dimensional digital airway network in urban low-altitude airspace with reasonable topology, high generation efficiency, good optimization performance and complete digital representation.

[0005] The technical solution of this invention is implemented as follows: On one hand, the present invention provides a method for constructing a three-dimensional digital airway network for urban low-altitude airspace, comprising: S1. Obtain traffic demand point distribution data, obstacle distribution data, and airspace control rules for urban low-altitude airspace to obtain the basic airspace dataset; S2. The weighted Voronoi-Delaunay discretization method is used to discretize the spatial basic dataset. The weight values ​​are determined by traffic demand and safety reservation coefficient to obtain the Voronoi cell set and the dual Delaunay network. S3. Based on the boundary curvature of the Voronoi unit, the flux density gradient of adjacent Voronoi units, and the shape characteristics of the Voronoi unit, a three-dimensional risk field is constructed to obtain the spatial three-dimensional risk field distribution. S4. Divide the low-altitude airspace into altitude layers according to the type of aircraft, and copy the Delaunay network nodes in each altitude layer to obtain a three-dimensional network node set. S5. Based on the three-dimensional risk field and the three-dimensional network node set, a potential field is generated by the potential field diffusion mechanism based on the graph Laplacian operator. The route skeleton is generated on the Delaunay network by the potential energy gradient guidance to obtain the initial route network. S6. Based on the initial route network, construct a three-dimensional digital model including a three-dimensional pipeline model, an interchange node model, and digital coding, and establish a digital twin system. The digital twin system records the operational status data of each route segment and node in real time. S7. Based on real-time data from the digital twin system, a multi-objective optimization algorithm is used to optimize the initial route network to obtain the optimized route network. S8. Perform path smoothing on the optimized route network and verify flight performance constraints to obtain route trajectories that meet dynamic requirements; S9. Based on the route trajectory update digital twin system, a three-dimensional digital route network for urban low-altitude airspace is obtained.

[0006] Preferably, step S2 includes: S21. Select a set of seed points in the airspace basic dataset. Each seed point corresponds to a traffic demand hotspot area center, important landmark building, or key point on the boundary of a controlled area. Set the seed point density according to the importance of the area. S22. Calculate the comprehensive weight value of each seed point based on traffic demand and safety reserve coefficient. The formula for calculating the weight value is as follows: Where v is the seed point index; This represents the overall weight value of the v-th seed point; Seed point location Normalized traffic demand at the location; The safety reserve coefficient for the v-th seed point; This is the adjustment coefficient; S23. Perform Voronoi subdivision of the spatial domain based on the weighted distance function to obtain a set of Voronoi cells. The weighted distance function is defined as: in, Let be any point in the airspace; Let v be the spatial coordinates of the v-th seed point; The distance is Euclidean. S24. Construct the dual network of the Voronoi graph to obtain the Delaunay network. The edges of the Delaunay network connect the seed points of adjacent Voronoi units. S25. Perform a safety screening on the Delaunay network, delete edges that cross no-fly zones, and retain edges that do not violate airspace control rules to obtain the Voronoi cell set and the dual Delaunay network.

[0007] Preferably, step S3 includes: S31. Integrate geographic information system data, real-time meteorological data, population density data and communication signal strength data to divide the urban low-altitude airspace into three-dimensional grid units. The grid unit size is 10 meters by 10 meters in the horizontal direction and 20 meters in the vertical direction. S32. Calculate the comprehensive risk assessment value for each grid cell. The comprehensive risk assessment function is defined as follows: Where q is the location point in the spatial domain; Let the obstacle risk function be used. For population density risk function; For meteorological conditions risk function; For communication quality risk function; The Voronoi topological risk function; These are the corresponding weighting coefficients; S33. Normalize the comprehensive risk assessment value to the range of [0,1] to obtain the three-dimensional risk field distribution of the airspace covering the entire city's low-altitude airspace.

[0008] Preferably, step S4 includes: S41. According to the aircraft classification standards and altitude management regulations, the low-altitude airspace from 0m to 1000m is divided into eleven altitude layers, with different altitude layers corresponding to different types of aircraft or flight missions. S42. For Delaunay network nodes, replicate and project them vertically according to eleven height layers to form a multi-layer three-dimensional route node network. The nodes in each height layer maintain their Voronoi topology on the horizontal plane. S43. For adjacent altitude layers, vertical connecting edges are set in traffic demand hotspot areas to allow aircraft to switch altitude layers under the premise of complying with the rules. The vertical connecting edges connect the node pairs corresponding to the horizontal positions in two adjacent altitude layers to obtain a three-dimensional network node set. S44. Map the calculated three-dimensional risk field onto the set of three-dimensional network nodes. For each node, its risk attribute is equal to the comprehensive risk assessment value of the grid cell corresponding to the three-dimensional spatial coordinates of that node. For each edge, its risk attribute is obtained by sampling along the edge path and calculating the weighted average value.

[0009] Preferably, step S5 includes: S51. Define the node potential function for each node in the set of nodes in the three-dimensional network: Where i is the node index; j is the neighbor node index; Let be the set of neighboring nodes of node i in the Delaunay network; For local potential energy, Let be the normalized traffic demand at node i. For the risk attribute value of node i, These are preset coefficients; Let be the local potential energy of neighbor node j; For the elements of the graph Laplace matrix; The Laplace diffusion coefficient is used. The integral of the risk gradient field within the Voronoi unit where node i is located; This represents the risk gradient impact coefficient. S52. Define the edge potential function for each edge in the set of nodes in the 3D network: in, The length of the edge between nodes i and j; The edge risk attribute between nodes i and j; The current flow of the edge between nodes i and j; Let i be the edge capacity between nodes i and j; These are the weighting coefficients; S53. The Dijkstra algorithm is used to calculate the path with the minimum potential energy from the starting point to the ending point on the three-dimensional network node set. The total potential energy of the path is equal to the sum of the potential energy of all nodes on the path and the potential energy of all edges. The path with the lowest potential energy is used as the main path and the path with the second lowest potential energy is used as the backup path. S54. Based on the traffic demand density of each region, the air route network is densified and refined. For regions where the normalized traffic demand is greater than a preset threshold, parallel feeder routes are inserted to form a three-level hierarchical structure: trunk routes connecting major urban areas, feeder routes covering the internal areas, and connecting routes connecting specific take-off and landing points, thus obtaining the initial air route network.

[0010] Preferably, the multi-objective optimization algorithm used in step S7 is an improved NSGA-II non-dominated sorting genetic algorithm, the specific steps of which include: S71. Generate an initial path population with a population size of N paths by utilizing the similarity of node embedding vectors in the low-dimensional feature space. S72. For each path in the current population, calculate three objective function values, including flight distance, energy consumption, and route congestion. At the same time, calculate the fitness function that includes Wasserstein distance. S73. Perform non-dominated sorting, divide the population into multiple front layers according to Pareto dominance, and calculate the crowding distance for the paths in each front layer. S74. Select the population based on the frontier level and crowding distance. Prioritize the retention of paths with lower frontier levels and paths with higher crowding within the same frontier level. Select the first M paths to enter the mating pool. S75. Perform crossover operation to generate offspring paths. Randomly select two parent paths from the mating pool, swap the node sequences of the two parent paths within the randomly selected interval to generate two offspring paths, and insert the shortest path segment to repair the break. S76. Perform mutation operation to mutate the child path with adaptive mutation probability. The mutation operation randomly selects a node in the child path and replaces it with a random neighbor node of that node in the Delaunay network. S77. Perform constraint checks on the generated child paths, checking for altitude layer constraints, no-fly zone constraints, and turning radius constraints. If the constraints are violated, apply a penalty to the fitness function. S78. Add the offspring paths to the population and enter the next iteration. When the change of the hypervolume index of the first front is less than the threshold or the maximum number of iterations is reached for ten consecutive generations, the iteration is terminated. Output the first front as the Pareto optimal path set. Select a suitable equilibrium solution from the Pareto front according to the aircraft mission preference to obtain the optimized route network.

[0011] Preferably, the fitness function including the Wasserstein distance in step S72 is defined as: in, A normalized multi-objective weighted sum; The topological influence weighting coefficient; The 2-Wasserstein distance between network traffic distributions before and after path selection; The Wasserstein distance is calculated in discrete form, and the formula is as follows: Where m and n are edge indices; For the transmission matrix; For the transmission matrix The element in the m-th row and n-th column represents the flow volume from edge m to edge n; Let m be the embedding vector of the node corresponding to edge m and edge n; The constraints are and as well as ,in Let $m$ be the flow rate before path selection and $n$ be the flow rate after path selection, respectively, and $k$ be the edge index used to sum the total flow rate over all edges.

[0012] Preferably, the node embedding vector in step S71 is calculated using a graph convolutional neural network, and the specific steps include: The node input feature vector is constructed, which includes spatial coordinates, potential energy value, risk value, demand, Minkowski functional vector of Voronoi element and topological feature vector. The Minkowski functional vector includes four components: volume of Voronoi element, half of surface area, average width and Euler characteristic number. Applying three layers of graph convolution operations to the Delaunay network, each graph convolution operation is defined as follows: Where b is the graph convolutional layer index, with values ​​from 1 to 3; It is a symmetric normalized adjacency matrix; This is the feature matrix of the (b-1)th layer; Let b be the learnable weight matrix of the b-th layer; It is the ReLU activation function; After three layers of graph convolution, a node embedding matrix is ​​obtained. The embedding vector of each node retains the topological location information of the node in the Delaunay network, the geometric information of the Voronoi unit, and the attribute information of the node.

[0013] Preferably, the height layer division includes: 0-20m, near-ground airspace; 20-50m, restricted airspace for consumer-grade aircraft; 50-120m, suitable airspace for consumer-grade aircraft; 120-150m, restricted airspace for light and small aircraft; 150-200m, suitable airspace for light and small aircraft; 200-230m, restricted airspace for medium-sized aircraft; 230-300m, suitable airspace for medium-sized aircraft; 300-400m, restricted airspace for large aircraft; 400-600m, suitable airspace for large aircraft; 600-700m, restricted airspace for manned aircraft at high speeds; 700-1000m, suitable airspace for manned aircraft at high speeds.

[0014] In addition, the present invention also provides a system for constructing a three-dimensional digital airway network for low-altitude urban areas, in order to implement the method described above.

[0015] The present invention has the following advantages over the prior art: This invention employs a weighted Voronoi-Delaunay discretization method, comprehensively determining weights based on traffic demand and safety reserve coefficients. This ensures a natural correlation between the generated route network topology and urban traffic demand distribution, improving the rationality of route layout compared to traditional uniform grid partitioning methods. A Voronoi topological risk function is introduced, comprehensively considering boundary curvature, cell volume, and flow gradient to achieve deep coupling between risk assessment and dynamic network topology characteristics, enabling routes to dynamically respond to the real-time operational status of the network. Based on the potential field diffusion mechanism of the graph Laplacian operator, global potential energy diffusion is achieved using the graph theory of Delaunay networks, avoiding the local optima problem of traditional local potential field methods and enabling route generation to perceive the global topology of the network. A graph convolutional neural network is introduced to extract node embedding vectors. The duality between the Voronoi graph and the Delaunay network is encoded into the feature space through the Minkowski functional of the Voronoi cell, allowing the optimization algorithm to simultaneously utilize network topological connectivity and spatial segmentation information. By using the Wasserstein distance metric to measure the change in network traffic distribution before and after path selection in the multi-objective optimization fitness function, the synergy between single aircraft path optimization and global network topology health is achieved, maintaining network load balance while meeting multiple objective requirements such as flight distance, energy consumption, and congestion. Attached Figure Description

[0016] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0017] Figure 1This is a flowchart of the method of the present invention; Figure 2 This is a schematic diagram of the tetrahedral network of the present invention; Figure 3 This is a schematic diagram of the route generation method of the present invention; Figure 4 This is a flowchart of the algorithm of the present invention; Figure 5 This is a diagram of the overall architecture of the three-dimensional digital airway network of the present invention. Detailed Implementation

[0018] The technical solutions of the present invention will be clearly and completely described below with reference to the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0019] like Figure 1 As shown, this invention provides a method for constructing a three-dimensional digital airway network in urban low-altitude airspace, including: S1, acquiring traffic demand point distribution data, obstacle distribution data, and airspace control rules for urban low-altitude airspace to obtain a basic airspace dataset; S2, using a weighted Voronoi-Delaunay discretization method to discretize the basic airspace dataset, where the weight values ​​are determined by a combination of traffic demand and safety reservation coefficients, to obtain a set of Voronoi cells and a dual Delaunay network; S3, constructing a three-dimensional risk field based on the boundary curvature of Voronoi cells, the flow density gradient of adjacent Voronoi cells, and the shape characteristics of Voronoi cells, to obtain the three-dimensional risk field distribution of the airspace; S4, dividing the low-altitude airspace into altitude layers according to aircraft type, and replicating Delaunay network nodes in each altitude layer to obtain... S5. Based on the 3D risk field and the 3D network node set, a potential field is generated using a potential field diffusion mechanism based on the graph Laplacian operator. The potential energy gradient guides the generation of the route skeleton on the Delaunay network, resulting in an initial route network. S6. Based on the initial route network, a 3D digital model is constructed, including a 3D pipeline model, an interchange node model, and digital coding. A digital twin system is established, which records the operational status data of each route segment and node in real time. S7. Based on the real-time data from the digital twin system, a multi-objective optimization algorithm is used to optimize the initial route network, resulting in an optimized route network. S8. The optimized route network undergoes path smoothing processing, and flight performance constraints are verified to obtain a route trajectory that meets dynamic requirements. S9. The digital twin system is updated based on the route trajectory to obtain a three-dimensional digital route network for urban low-altitude airspace.

[0020] In one embodiment of the present invention, step S1 includes: establishing a digital rule base for low-altitude flight that covers airspace control rules, aircraft classification standards, flight altitude restrictions, no-fly zone boundaries, flight speed restrictions, minimum safety intervals, and takeoff and landing operation specifications; collecting and integrating traffic demand data for urban low-altitude flights; identifying traffic demand hotspots in urban low-altitude airspace by analyzing historical flight records, logistics delivery demand, and emergency service demand; and assigning normalized demand values ​​to each key location, such as commercial center, hospital, fire station, and logistics hub. ,in A database of urban low-altitude traffic demand distribution is established using spatial location coordinates; urban take-off and landing point information is collected and organized, including the location coordinates, type classification, and capacity attributes of the take-off and landing points, to establish a take-off and landing point database; an aircraft performance parameter database is established, containing performance data such as power parameters, drag coefficients, mass parameters, turning radius, and climb and descent rates for different types of aircraft; the rule base and various databases are stored in a structured data format and support online updates to adapt to dynamic changes in control rules.

[0021] In one embodiment of the present invention, step S2 includes: S21. Select a set of seed points in the airspace basic dataset. Each seed point corresponds to a traffic demand hotspot area center, important landmark building, or key point on the boundary of a controlled area. Set the seed point density according to the importance of the area. Seed points are selected from the center of traffic demand hotspots, important landmarks, key points at the boundaries of controlled areas, and emergency service demand points. The density of seed points is determined according to the importance and complexity of the area. In densely built-up central business districts, the seed point spacing is set to 300m to 500m. In general urban areas, the seed point spacing is set to 800m to 1000m. In open urban fringe areas, the seed point spacing can be increased to 1500m to 2000m. S22. To address the technical problem of insufficient correlation between airspace discretization results and traffic demand distribution in the construction of urban low-altitude airway networks, a comprehensive weight value for each seed point is calculated based on traffic demand and safety reserve coefficient. The formula for calculating the weight value is as follows: Where v is the seed point index; This represents the overall weight value of the v-th seed point; Seed point location The normalized traffic demand at the location, with a value range of [0,1], is obtained from the demand distribution database in step S1; The safety reserve coefficient for the v-th seed point; This is the adjustment coefficient; Safety Reserve Coefficient The adjustment coefficient is calculated based on the density of surrounding obstacles and population distribution. Based on the dynamic setting of urban functional zones, commercial areas are weighted more heavily based on demand. A weight of 0.7 to 0.8 is used, with a focus on safety reserves for residential areas. The weight is set to 0.3 to 0.5. This weighting formula causes the Voronoi cells of high-weight seed points to automatically expand, forming a larger route influence range, while the Voronoi cells of low-weight seed points automatically shrink to restrict route passage. S23. Perform Voronoi subdivision of the spatial domain based on the weighted distance function to obtain a set of Voronoi cells. The weighted distance function is defined as: Where q is any point in the spatial domain; Let v be the spatial coordinates of the v-th seed point; The Euclidean distance is used; the introduction of weighted distance causes the Voronoi cells of high-weight seed points to automatically expand to form a larger route influence range, while the Voronoi cells of low-weight seed points automatically shrink to restrict route passage. Voronoi unit Defined as all points to the seed point The set of spatial points with the minimum weighted distance: The Voronoi unit boundary, as a natural boundary between adjacent influence domains, has natural obstacle avoidance characteristics. S24. Construct the dual network of the Voronoi graph to obtain the Delaunay network. The edges of the Delaunay network connect the seed points of adjacent Voronoi cells. The Delaunay network has the geometric property of maximizing the minimum angle, thus ensuring that the angles between the generated route segments are relatively gentle and avoiding sharp turns.

[0022] S25. Perform a safety screening on the Delaunay network, deleting edges that cross no-fly zones and retaining edges that do not violate airspace control rules, resulting in a set of Voronoi cells and the dual Delaunay network. The no-fly zone boundary is obtained from the rule base in step S1. Nodes in the skeleton network record their spatial three-dimensional coordinates, the volume attribute of their Voronoi cell, the shape complexity of their Voronoi cell, the list of connected edges, node type identifiers, and other geometric attributes. Edges in the skeleton network record the unique identifiers of the starting and ending nodes, the spatial length of the edge, the direction vector, the topological level of the edge in the Delaunay network, and other geometric attributes. In this embodiment, the schematic diagram of the constructed Delaunay tetrahedral network is shown below. Figure 2 As shown; This step uses the weighted Voronoi-Delaunay discretization method to naturally correlate the generated airway network topology with the distribution of urban traffic demand, improving the rationality of airway layout compared to the traditional uniform grid partitioning method.

[0023] In one embodiment of the present invention, step S3 includes: S31. Integrate geographic information system data, real-time meteorological data, population density data and communication signal strength data to divide the urban low-altitude airspace into three-dimensional grid units. The grid unit size is 10 meters by 10 meters in the horizontal direction and 20 meters in the vertical direction. S32. Calculate the comprehensive risk assessment value for each grid cell. The comprehensive risk assessment function is defined as follows: Where q is the location point in the spatial domain; Let the obstacle risk function be used. For population density risk function; For meteorological conditions risk function; For communication quality risk function; The Voronoi topological risk function; These are the corresponding weighting coefficients; Obstacle risk function Reflecting the threat level of static obstacles such as buildings, terrain, and trees, an inverse proportional function is used to map the distance between the aircraft and the nearest obstacle to a risk value; the closer the distance, the higher the risk. Population density risk function. To reflect the impact of ground population density on airway safety, dynamic population density data is processed using a spatiotemporal kriging interpolation algorithm and mapped to risk values; a meteorological condition risk function is also included. Taking into account dynamic meteorological factors such as wind speed, rainfall, visibility, and lightning, the comprehensive meteorological risk is calculated using the Mamdani fuzzy logic system; communication quality risk function. This model reflects the coverage strength and quality of communication signals. Based on the logarithmic distance path loss model and measured data, the signal strength at a given location is calculated and then mapped back to a risk value. This model is a classic model for radio wave propagation prediction. The back mapping means that the signal strength and risk value are negatively correlated, i.e., the stronger the signal strength, the lower the risk, and the weaker the signal strength, the higher the risk. The mapping function adopts a sigmoid function. Voronoi topological risk function This invention represents a core improvement to the Voronoi-Delaunay network structure, defined as follows: in Let q be the Euclidean distance from point q to the nearest Voronoi boundary; The safe threshold distance is set between 30m and 50m. The Frobenius norm of the Hessian matrix of the boundary distance field reflects the local curvature of the boundary. The greater the curvature, the more complex the boundary, where multiple Voronoi cells intersect, and the higher the risk of route conflict. Let q be the volume of the Voronoi cell to which point q belongs. If the volume is too small, it indicates that the airspace in this region is cramped and lacks maneuvering space. The shape complexity index is denoted by a, b, and c, which are the three side lengths of the smallest circumscribed cuboid of the Voronoi element. This index reflects the anisotropy of the element shape. An overly long or irregular shape indicates that the surrounding obstacles are complex. This represents the flux density gradient between adjacent Voronoi cells. Let q be the Voronoi elements adjacent to both sides of the boundary containing point q. Voronoi unit The aircraft flow density is such that when the flow differences between adjacent units are large, there is significant traffic flow convergence or divergence at the boundary, which significantly increases the risk. These are the weighting coefficients for the sub-items, set based on experimental statistics; Weighting coefficient satisfy Each weighting coefficient is dynamically adjusted using machine learning methods based on historical flight accident data and expert experience; S33. Normalize the comprehensive risk assessment value to the range of [0,1] to obtain the three-dimensional risk field distribution of the entire city's low-altitude airspace. Store the normalized risk field data in a spatial database for subsequent steps. The three-dimensional risk field established in this step deeply couples risk assessment with the dynamic characteristics of network topology through the Voronoi topological risk function. Compared with the traditional risk function that only considers the static distance to obstacles, it can dynamically respond to the real-time operating status of the airway network.

[0024] In one embodiment of the present invention, step S4 includes: S41. According to the aircraft classification standards and altitude management regulations, the low-altitude airspace from 0m to 1000m is divided into eleven altitude layers, with different altitude layers corresponding to different types of aircraft or flight missions. Specifically, it is divided as follows: 0-20m near-ground airspace, 20-50m isolated airspace for consumer-grade aircraft, 50-120m suitable airspace for consumer-grade aircraft, 120-150m isolated airspace for light and small aircraft, 150-200m suitable airspace for light and small aircraft, 200-230m isolated airspace for medium-sized aircraft, 230-300m suitable airspace for medium-sized aircraft, 300-400m isolated airspace for large aircraft, 400-600m suitable airspace for large aircraft, 600-700m isolated airspace for manned aircraft (hyperspeed), and 700-1000m suitable airspace for manned aircraft (hyperspeed). Physical isolation between different altitude levels is achieved through isolated airspace to avoid mutual interference between different types of aircraft; S42. For the Delaunay network nodes, replicate and project them vertically according to eleven altitude layers to form a multi-layered, three-dimensional route node network. , For altitude layer indexing, nodes within each altitude layer maintain their Voronoi topology on the horizontal plane. The vertical coordinates of nodes within each altitude layer are set to the midpoint altitude value of that altitude layer. The capacity and speed limit attributes of the connecting edges between nodes are adjusted according to the aircraft type and speed requirements of that altitude layer. The capacity attribute is calculated based on the airway width and the minimum aircraft spacing, representing the maximum number of aircraft that can pass through per unit time. S43. For adjacent altitude layers, vertical connecting edges are set in traffic demand hotspot areas, allowing aircraft to switch altitude layers under the premise of complying with rules. The vertical connecting edges connect the corresponding node pairs in the horizontal positions of two adjacent altitude layers, resulting in a three-dimensional network node set. The starting point of the vertical connecting edge is a node in the lower altitude layer, and the ending point is the corresponding node in the higher altitude layer directly above. The capacity of the vertical edge is set according to the safety interval requirements for altitude layer switching. The resulting multi-layer Delaunay network maintains the horizontal topology within each altitude layer while achieving three-dimensional interconnection between layers through vertical connecting edges. S44. Map the calculated three-dimensional risk field onto the set of three-dimensional network nodes. For each node, its risk attribute is equal to the comprehensive risk assessment value of the grid cell corresponding to the three-dimensional spatial coordinates of that node. For each edge, its risk attribute is obtained by sampling along the edge path and calculating the weighted average value. In this embodiment, for the height layer nodes Its risk attributes The three-dimensional spatial coordinates of this node The comprehensive risk assessment value corresponding to the grid cell, i.e. Retrieved from the risk field database; For the edge Including horizontal and vertical edges, their risk attributes Risk is calculated by sampling along the edge path, sampling multiple spatial points along the path with a sampling point spacing of 10m to 20m, querying the risk value of the corresponding raster cell for each sampling point, and calculating the weighted average as the risk attribute of the edge. in For the parameterized representation of edges, The weighting function is set so that the sampling point at the midpoint of the edge has the highest weight of 1.5, while the sampling points at both ends have a weight of 0. This weighting distribution reflects the greater impact of the midpoint of the edge on the overall route risk. After risk attribute mapping, each node and edge in the multi-layer Delaunay network has complete digital attributes such as spatial coordinates, Voronoi cell attributes, connectivity, capacity attributes, speed limit attributes, and risk attributes.

[0025] like Figure 3 As shown, in one embodiment of the present invention, based on the multi-layer three-dimensional Delaunay network generated in step S4, an improved artificial potential field theory is introduced to guide the automatic generation of routes, so as to solve the technical problem of the disconnect between the route generation mechanism and the network topology in the existing method; step S5 includes: S51. A distributed potential field model is established based on the node-edge topological characteristics of the Delaunay network. A node potential function is defined for each node in the Delaunay network, consisting of demand-attracting potential energy and risk-repelling potential energy. However, traditional potential functions are prone to getting trapped in local optima in complex topological networks. Therefore, a potential field diffusion mechanism based on the spectral characteristics of the Delaunay network is introduced. For nodes... Define the nodal potential function : Where i is the node index; j is the neighbor node index; Let be the set of neighboring nodes of node i in the Delaunay network; For local potential energy, Let be the normalized traffic demand at node i. For the risk attribute value of node i, These are preset coefficients; Let be the local potential energy of neighbor node j; For the elements of the graph Laplace matrix; The Laplace diffusion coefficient is used. The integral of the risk gradient field within the Voronoi unit where node i is located; This represents the risk gradient impact coefficient. In local potential energy, demand Generates attraction, risk value Repulsive effect occurs, coefficient Set according to task priority, generally take This strengthens the risk-repelling effect; Elements of the Graph Laplace matrix Defined as: in Let be the degree of node i. Laplace diffusion coefficient. The value ranges from 0.1 to 0.3, controlling the degree of diffusion of potential energy on the network; Risk gradient field integral Reflecting the spatial rate of change of risk in the region, a large risk gradient within a cell indicates the existence of a local risk abrupt change region. The potential energy increases accordingly to guide flight paths to avoid this region. The risk gradient influence coefficient... The value ranges from 0.5 to 1.0; S52. Define the edge potential function for each edge in the Delaunay network. The edge potential energy consists of three parts: spatial length potential energy, path risk potential energy, and congestion penalty potential energy. in, The length of the edge between nodes i and j; The edge risk attribute between nodes i and j; The current flow of the edge between nodes i and j; Let i be the edge capacity between nodes i and j; The weighting coefficient is typically taken as... The congestion penalty potential energy adopts an exponential form. When the flow at an edge approaches its capacity, the penalty term increases rapidly, effectively avoiding excessive concentration of air routes. In this embodiment, a potential field propagation constraint mechanism based on Voronoi units is introduced. This mechanism utilizes the spatial segmentation characteristics of the Voronoi diagram generated in step S2, ensuring that potential field energy propagates only between nodes belonging to the same or adjacent Voronoi units, and prohibiting long-range potential field interactions that cross multiple Voronoi unit boundaries. The potential field influence strength of node i on node j is... for: This constraint ensures that the potential field propagates only between adjacent Voronoi cells, avoiding unreasonable effects from a distance, and making the potential field distribution more consistent with the actual topology of urban airspace. S53. The Dijkstra algorithm is used to calculate the path with the minimum potential energy from the starting point to the ending point on the three-dimensional network node set. The total potential energy of the path is equal to the sum of the potential energy of all nodes on the path and the potential energy of all edges. The path with the lowest potential energy is used as the main path and the path with the second lowest potential energy is used as the backup path. The total potential energy of the path is defined as: The process of Dijkstra's potential minimum path algorithm is as follows: Initialize the distance from the starting point s. Distance to all other nodes Add all nodes to a priority queue Q; during the iteration, each time the unvisited node u with the lowest current potential energy is taken from queue Q. If u is the destination t, the process terminates and backtracks the path; otherwise, u is removed from Q. For each neighbor node v of u, calculate the total cost of reaching v through u. If the cost is less than the currently recorded Then update And record the predecessor node Iterate until the endpoint is retrieved or the queue is empty, then backtrack the predecessor pointer from the endpoint t. The path with minimum potential energy is obtained from the starting point s; S54. The airway network is refined and densified based on the traffic demand density of each region. For regions where the normalized traffic demand exceeds a preset threshold, parallel feeder routes are inserted, forming a three-tiered structure: trunk routes connecting major urban areas, feeder routes covering internal areas, and connecting routes connecting specific take-off and landing points, resulting in the initial airway network. For regions with high traffic demand, parallel feeder routes are inserted between trunk routes, with a spacing of 300m to 500m. Feeder routes are generated by creating new virtual nodes near trunk route nodes and connecting them to form paths parallel to the trunk routes. Each generated route records its path node sequence, total path potential value, risk value distribution along the route, expected capacity, and other digital attributes. This step introduces a potential field diffusion mechanism based on the graph Laplacian operator, and uses the graph theory of Delaunay networks to diffuse the global potential energy, avoiding the local optima problem of traditional local potential field methods, so that the route generation can sense the global topology of the network.

[0026] In one embodiment of the present invention, step S6 constructs a three-dimensional digital model based on the initial route network, including a three-dimensional pipeline model, an interchange node model, and digital coding, and establishes a digital twin system. The digital twin system records the operational status data such as flow rate, capacity, and risk value of each route segment and node in real time. A three-dimensional pipeline model of the airway network is constructed, expanding each airway path into an elliptical cylindrical pipeline structure. The centerline of the pipeline is the airway path curve, which is generated into a smooth three-dimensional spatial curve by cubic spline interpolation of the path node sequence. The cross-section of the pipeline is elliptical, with the major axis along the flight direction and the minor axis along the lateral deviation direction. The radius of the major axis is determined based on the aircraft's maximum wingspan plus a safety margin, and the radius of the minor axis is determined based on the allowable lateral deviation distance. The pipeline diameter for main airways is set to 80m to 100m, the pipeline diameter for branch airways is set to 50m to 60m, and the pipeline diameter for connecting airways is set to 30m to 40m.

[0027] Establish a digital coding system for air routes. Assign a unique digital identifier to each component in the air route network. The air route segment coding adopts a three-level structure: altitude layer code, area code, and air route number. The altitude layer code uses a two-digit alphanumeric combination to represent eleven altitude layers. The area code uses city grid numbers. The air route number is the sequential number of the air routes within that area. Node coding uses a combination format of node type identifier and spatial location hash value. Node types include categories such as ordinary nodes, traffic nodes, and take-off and landing points. The spatial location hash value is calculated from three-dimensional coordinates using the Geohash three-dimensional extended algorithm to ensure uniqueness. Take-off and landing point coding adopts a format of city code, functional area code, and point number. The coding system is compatible with the coding standards of the National Low Altitude Airspace Management Information System. A 3D digital model library of traffic nodes was constructed, and 3D models of grade-separated intersection nodes of different forms were designed based on parameters such as the number, direction, and altitude difference of merging routes at the nodes. For route intersections within the same altitude level, a planar roundabout node model was adopted, with the aircraft circling counterclockwise after entering the roundabout. For the convergence of routes at different altitude levels, a spiral grade-separated intersection node model was adopted, achieving a smooth transition of altitude levels through spiral ascending or descending ramps. A 3D model of the conical flight path at takeoff and landing points is constructed. The flight path adopts a frustum or inverted frustum structure, connecting the ground takeoff and landing points with the airway network. The conical flight path model supports real-time collision detection and conflict warning. When multiple aircraft use the conical flight path at the same takeoff and landing point simultaneously, the system calculates the minimum interval based on the aircraft's altitude and relative position and triggers a conflict warning. A route network graph database is established to store the topology and attribute data of the route network. Node objects in the database store attributes such as unique identifier, 3D coordinates, altitude layer, Voronoi cell number, potential energy value, risk value, demand, list of connecting edges, and node type. Edge objects store attributes such as unique identifier, starting node identifier, ending node identifier, edge length, direction vector, capacity, current flow, and risk value. A digital twin system for the airway network is established, linking a 3D digital airway model with real-time operational data to maintain a digital twin state vector for each airway segment and node. The state vector includes real-time status information such as the number of currently passing aircraft, the real-time position and speed of each aircraft, average flight speed, congestion index, most recent risk assessment value, communication quality indicators, and meteorological condition parameters. The digital twin system updates the state vector by collecting data from the aircraft's Automatic Dependent Surveillance-Broadcast (ADS-B) system and ground sensor networks every five to ten seconds, presenting the real-time operational status of the airway network through a 3D visualization interface.

[0028] like Figure 4 As shown, in one embodiment of the present invention, step S7, based on the three-dimensional digital route model and digital twin system established in step S6, uses the NSGA-II non-dominated sorting genetic algorithm based on graph embedding to perform refined route planning. The system receives flight path planning requests from aircraft, including information such as start and end points, aircraft type, expected takeoff time, and expected arrival time. It queries the aircraft performance parameter database established in step S1 for performance data such as power parameters, drag coefficient, mass, and turning radius for that aircraft type. It reads the topology of the flight path network from the graph database established in step S6, and retrieves attributes such as spatial coordinates, potential energy, risk value, and demand from each node, as well as attributes such as length and risk value from each edge. It also reads the current flow and capacity data of each node and edge in real time from the digital twin system, using these as inputs to the optimization algorithm. The nodes and edges of the Delaunay network are mapped to a low-dimensional feature space, and a graph neural network is used to learn the network's topology embedding. For each node, its input feature vector includes spatial coordinates, potential energy, risk value, demand, the Minkowski functional of the Voronoi unit, and topological features. The Minkowski functional vector contains four components: volume, half the surface area, average width, and Euler characteristic of the Voronoi unit, fully describing its geometry. The topological feature vector contains three components: node degree, clustering coefficient, and betweenness centrality, describing the topological position of a node in the network. Node degree is the number of edges connecting the node; the clustering coefficient reflects the connection density between the node's neighbors, calculated as the ratio of the actual number of edges between the node's neighbors to the possible number of edges; and betweenness centrality reflects the bridging effect of the node in the network, calculated as the proportion of the shortest paths passing through that node out of all shortest paths. The graph neural network contains three graph convolutional layers. Each graph convolutional operation uses a symmetrically normalized adjacency matrix to aggregate and transform node features. The mathematical expression of the graph convolutional operation is: Where b is the graph convolutional layer index, with values ​​from 1 to 3; It is a symmetric normalized adjacency matrix; This is the feature matrix of the (b-1)th layer; Let b be the learnable weight matrix of the b-th layer; It is the ReLU activation function; The hidden layer dimensions of the three-layer graph convolution are set to 128, 64, and 32, respectively. After the three-layer graph convolution, a 32-dimensional embedding vector is obtained for each node. The embedding vector retains the topological location information of the node in the Delaunay network, the geometric information of the Voronoi unit, and the attribute information of the node. The graph neural network is trained using supervised learning. The training dataset consists of historical flight paths and their quality evaluations. The loss function is the mean squared error of the path quality prediction, and the Adam optimizer is used for end-to-end training. A multi-objective optimization model is constructed, with optimization objectives including three dimensions: shortest flight distance, minimum energy consumption, and lowest route congestion; objective function 1 is the flight distance. ;in Let be the spatial length of the edge; objective function 2 is the energy consumption. Energy consumption at the edges Calculations based on the aircraft dynamics model: This formula comprehensively considers the aircraft's basic power consumption, air resistance power consumption, and potential energy changes caused by altitude variations, where... Based on power, air density, Where A is the drag coefficient, v is the frontal area, m is the flight speed, and g is the mass of the aircraft. For the change in height, For flight time; objective function 3 is the route congestion level. ;in The current traffic and capacity of the nodes are read in real time from the digital twin system. The current traffic and capacity are read in real time from the digital twin system, and the congestion target-guided algorithm selects the route segment with less traffic to achieve network load balancing; To avoid cascading congestion caused by drastic changes in network traffic distribution due to path selection during optimization, a network traffic distribution change metric based on Wasserstein distance is introduced. The fitness function for topology awareness is defined as follows: in For normalized multi-objective weighted sums, The average value of target k in the population is used for normalization. The weighting coefficients are set according to the aircraft mission type. The topological influence weighting coefficient ranges from 0.1 to 0.3. The 2-Wasserstein distance is the network traffic distribution before and after path selection. The Wasserstein distance measures the difference between two probability distributions; its physical meaning is the minimum transmission cost required to transform one probability distribution into another. The definition of the Wasserstein distance is: m and n are edge indices; For the transmission matrix; For the transmission matrix The element in the m-th row and n-th column represents the flow volume from edge m to edge n; Let m be the embedding vector of the node corresponding to edge m and edge n; The constraints are and as well as ,in Let $\mathbf{m}$ represent the flow of edge $m$ before path selection and the flow of edge $n$ after path selection, respectively, and $k$ be the edge index used to sum the total flow over all edges. This optimization problem can be solved using the network simplex method or the Sinkhorn iterative algorithm. An improved NSGA-II non-dominated sorting genetic algorithm is used to solve the multi-objective optimization problem. The algorithm first initializes the population and generates initial paths in the embedding space using the similarity of node embedding vectors. For the i-th initial path, starting from the starting point s, the next node is iteratively selected. Among all unvisited neighbors of the current node u, the selection probability of neighbor node v is calculated. , where h v ,h t Given the embedding vectors of node v and endpoint t, randomly select the next node to join the path according to this probability distribution, gradually expanding until the endpoint t is reached. Generate population size. The initial path is set, and the population size is determined based on the network complexity, typically taking [a certain value]. The range is 100 to 200; Entering the iterative optimization phase, the maximum number of iterations is set to 100 to 200 generations. Each iteration includes the following operations: First, calculate three objective function values ​​for each path in the current population. and fitness functions including Wasserstein distance Then, a non-dominated sort is performed, dividing the population into multiple front layers according to Pareto dominance. The paths in the first front layer represent the Pareto optimal set, where no two paths dominate each other. The crowding distance is calculated for each path within each front layer. The population is selected based on the frontier level and crowding distance. Paths with lower frontier levels are given priority for retention, and paths with higher crowding within the same frontier level are given priority for retention. The top 100 to 200 paths are selected to enter the mating pool. Perform a crossover operation to generate offspring paths. Randomly select two parent paths from the mating pool, randomly select two crossover points, and swap the node sequences between the two parent paths at the crossover points to generate two offspring paths. Since crossover may cause offspring path breaks, i.e., adjacent nodes are not connected in the Delaunay network, it is necessary to repair the breaks by inserting the shortest path segment connecting the two nodes at the breaks. The shortest path is calculated on the Delaunay network using Dijkstra's algorithm. Mutation operations increase population diversity. Offspring paths are mutated using an adaptive mutation probability, initially set to 0.3, which is gradually decreased with iterations. The mutation probability is calculated as follows: ,in Let g be the initial mutation probability, and g be the current iteration number. To achieve the maximum number of iterations, the mutation probability gradually decreases with each iteration, thus transitioning from a global search to a local search. The mutation operation randomly selects a node in the offspring path and replaces it with a random neighbor node of that node in the Delaunay network. The generated child paths are subjected to constraint checks. These checks include: altitude layer constraints (the path cannot cross the upper or lower boundaries of a given altitude layer); no-fly zone constraints (the path cannot cross no-fly zones defined in the rule base of step S1); and turning radius constraints (the radius of curvature at the path's turns cannot be less than the aircraft's minimum turning radius). If any constraints are violated, a penalty is applied to the fitness function. A large penalty coefficient ensures that paths violating the constraints are eliminated during the selection process. The offspring paths are added to the population for the next iteration. The algorithm is considered convergent when the change in the hypervolume index of the first front is less than the threshold of 0.01 for ten consecutive generations, or when the maximum number of iterations is reached, the iteration is terminated. The hypervolume index is defined as the size of the hypervolume dominated by the first front in the target space, and is a commonly used indicator to measure the quality of the Pareto front. The algorithm outputs the final first frontier as the Pareto optimal path set, which contains multiple non-dominant paths, each with an advantage on one or more objectives. Based on the aircraft's mission preferences, a suitable equilibrium solution is selected from the Pareto frontier. If the mission is time-sensitive, the shortest path is chosen; if energy conservation is required, the path with the lowest energy consumption is selected; if the network is congested, a path avoiding congested areas is selected. Alternatively, multi-criteria decision-making methods such as the ideal point method or the weighted Tchebycheff method can be used to select a comprehensive and balanced compromise solution from the Pareto frontier. The planned path results are updated in real time to the graph database and digital twin system established in step S6. The current traffic and remaining capacity of each edge traversed by the path are updated, as are the traffic and congestion index of each node traversed by the path. The state vector of the digital twin system is updated, including information such as the aircraft's mission number, position, speed, and estimated arrival time.

[0029] In one embodiment of the present invention, the path output by the optimization algorithm is a discrete sequence of nodes, which needs to be smoothed to generate a continuous smooth curve to meet the dynamic constraints of the aircraft. Cubic B-spline curves are used for path interpolation, with path nodes serving as control points to generate a smooth three-dimensional spatial curve. B-spline curves have the advantages of good local support and continuity, ensuring the continuity of the curve's position, first derivative, and second derivative. Node vectors are generated using uniform node vectors or a parameterization method based on chord length. Flight performance constraints are verified on the smoothed path curve. Multiple points are sampled along the path at intervals of ten to twenty meters, and four verifications are performed on each sampled point. Turning radius verification checks whether the minimum turning radius requirement of the aircraft is met by calculating the path curvature radius, which is read from the performance parameter database in step S1. Climb / descent rate verification checks whether the vertical velocity exceeds the maximum climb or descent rate of the aircraft. Velocity change verification checks whether the acceleration modulus exceeds the maximum acceleration of the aircraft. Safety margin verification checks whether the minimum distance from the calculated path to the airway boundary meets the safety margin requirements, generally requiring a minimum distance greater than or equal to 1.5 times the aircraft's wingspan. If any verification fails, an iterative adjustment strategy is employed. This involves adjusting the control point positions or increasing the number of control points using optimization algorithms. The objective function is to minimize the path adjustment while satisfying all constraints. This constraint optimization problem is solved using a sequential quadratic programming algorithm or the interior-point method, iteratively adjusting until all verification items pass. The smoothed path after successful verification is sent to the aircraft's flight control system as the final flight trajectory. The path data includes the parametric equations of the path curve, suggested flight speeds for each segment, timestamps for each key point, and other complete flight plan information.

[0030] like Figure 5 As shown, in one embodiment of the present invention, the output results of all the above processing steps are integrated to generate a complete urban low-altitude three-dimensional digital airway network system.

[0031] The airway network comprises four basic components in its spatial structure: horizontal segments, vertical segments, grade-separated intersections, and terminal routes. Horizontal segments consist of main and feeder routes at each altitude level. Main routes connect major urban transportation hubs, while feeder routes cover secondary demand points. Vertical segments connect different altitude levels and are located in traffic hotspots, allowing aircraft to change altitudes. Grade-separated intersections are located at the junctions of multiple airways, and traffic flow is managed using either planar roundabouts or spiral interchanges based on the intersection conditions. Terminal routes include transition routes from main routes to ground takeoff and landing points and conical takeoff and landing ramps. All components are uniformly identified and managed according to the coding system established in step S6.

[0032] A distributed route network database system is established. The graph database stores the route network topology and supports graph operations such as shortest path queries and subgraph matching. The spatial database stores 3D geometric model data, using a standard 3D modeling format. The time-series database records the real-time status of each component of the route network at a sampling frequency of five to ten seconds, including route segment traffic, node congestion index, aircraft position and speed, and weather conditions. These three databases are integrated through a unified data access interface.

[0033] A 3D visualization platform is established to render the airway network scene in real time within a browser. The platform displays the spatial layout of the airway network, distinguishing different types of airways with different colors and line widths, and marking intersections and takeoff / landing points with special symbols. The platform also displays the dynamic position of aircraft in real time, showing information such as their number, altitude, and speed, with their movement trajectories displayed as colored lines. Each airway segment is colored according to current traffic flow: green indicates unobstructed traffic, yellow indicates moderate traffic, and red indicates congestion.

[0034] Establish an intelligent management system for the airway network. A traffic flow monitoring module collects real-time traffic data and calculates a congestion index, triggering a congestion warning when the index exceeds a threshold. A dynamic capacity adjustment module dynamically adjusts airway capacity based on real-time traffic flow; when flow approaches its limit, it generates parallel temporary airways or expands the pipe diameter, automatically restoring capacity as flow decreases. A conflict detection and resolution module predicts future aircraft trajectories, detects potential conflicts, and resolves conflicts by adjusting speed, altitude, or path. An emergency response module handles emergencies, assesses affected airway segments and aircraft, and replans alternative routes or guides them to emergency landing points.

[0035] Through the above system integration, a complete urban low-altitude three-dimensional digital airway network system is formed. This system includes a three-dimensional spatial structure, a digital coding system, real-time status monitoring, and dynamic management functions, providing safe and efficient airway infrastructure for low-altitude aircraft.

[0036] In addition, the present invention also provides a system for constructing a three-dimensional digital airway network for low-altitude urban airspace, comprising: The data construction module is used to build a low-altitude airspace management rule base and an aircraft performance parameter base; The demand analysis module is used to identify demand points within the urban area and generate a continuous demand density field based on the location and intensity of the demand points. The candidate point generation module is used to generate a set of candidate waypoints for a particle system based on potential energy. By setting the potential energy field function, the particles are attracted in areas with high demand density and repelled near obstacles and no-fly zones. The particles move iteratively in the potential energy field until convergence to obtain the set of candidate waypoints. The altitude layer division module is used to divide low-altitude airspace into multiple altitude layers based on aircraft type, mission nature, and airspace management rules. The network generation module is used to construct the spatial network topology by performing Delaunay triangulation on candidate waypoints, calculate the Voronoi diagram through duality to obtain the spatial unit division, and map the potential energy field and risk field to the nodes and edges of the Delaunay network to obtain the initial waypoint network. The 3D modeling module is used to construct a 3D digital model based on the initial route network, including a 3D pipeline model, an interchange node model, and digital coding, and to establish a digital twin system. The digital twin system records the operational status data such as flow, capacity, and risk value of each route segment and node in real time. The route optimization module is used to optimize the initial route network using a multi-objective optimization algorithm. It reads the flow and capacity data of nodes and edges in real time from the digital twin system for the calculation of the route congestion objective function. The optimization objectives include the shortest flight distance, the minimum energy consumption, and the lowest route congestion, resulting in the optimized route network. The trajectory verification module is used to perform path smoothing on the optimized route network and verify flight performance constraints. Path smoothing adopts cubic B-spline curve interpolation, and flight performance constraints include turning radius constraints, climb and descent rate constraints, speed change constraints, and safety separation constraints, to obtain a route trajectory that meets the dynamic requirements. The system integration module is used to update the digital twin system based on airway trajectory, integrating three-dimensional spatial structure, digital coding system, real-time status monitoring and dynamic management functions to obtain a three-dimensional digital airway network for urban low-altitude airways.

[0037] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for constructing a three-dimensional digital airway network for urban low-altitude airspace, characterized in that, include: S1. Obtain traffic demand point distribution data, obstacle distribution data, and airspace control rules for urban low-altitude airspace to obtain the basic airspace dataset; S2. The weighted Voronoi-Delaunay discretization method is used to discretize the spatial basic dataset. The weight values ​​are determined by traffic demand and safety reservation coefficient to obtain the Voronoi cell set and the dual Delaunay network. S3. Based on the boundary curvature of the Voronoi unit, the flux density gradient of adjacent Voronoi units, and the shape characteristics of the Voronoi unit, a three-dimensional risk field is constructed to obtain the spatial three-dimensional risk field distribution. S4. Divide the low-altitude airspace into altitude layers according to the type of aircraft, and copy the Delaunay network nodes in each altitude layer to obtain a three-dimensional network node set. S5. Based on the three-dimensional risk field and the three-dimensional network node set, a potential field is generated by the potential field diffusion mechanism based on the graph Laplacian operator. The route skeleton is generated on the Delaunay network by the potential energy gradient guidance to obtain the initial route network. S6. Based on the initial route network, construct a three-dimensional digital model including a three-dimensional pipeline model, an interchange node model, and digital coding, and establish a digital twin system. The digital twin system records the operational status data of each route segment and node in real time. S7. Based on real-time data from the digital twin system, a multi-objective optimization algorithm is used to optimize the initial route network to obtain the optimized route network. S8. Perform path smoothing on the optimized route network and verify flight performance constraints to obtain route trajectories that meet dynamic requirements; S9. Based on the route trajectory update digital twin system, a three-dimensional digital route network for urban low-altitude airspace is obtained.

2. The method for constructing a three-dimensional digital airway network for urban low-altitude airspace according to claim 1, characterized in that, Step S2 includes: S21. Select a set of seed points in the airspace basic dataset. Each seed point corresponds to a traffic demand hotspot area center, important landmark building, or key point on the boundary of a controlled area. Set the seed point density according to the importance of the area. S22. Calculate the comprehensive weight value of each seed point based on traffic demand and safety reserve coefficient. The formula for calculating the weight value is as follows: Where v is the seed point index; This represents the overall weight value of the v-th seed point; Seed point location Normalized traffic demand at the location; The safety reserve coefficient for the v-th seed point; This is the adjustment coefficient; S23. Perform Voronoi subdivision of the spatial domain based on the weighted distance function to obtain a set of Voronoi cells. The weighted distance function is defined as: in, Let be any point in the airspace; Let v be the spatial coordinates of the v-th seed point; The distance is Euclidean. S24. Construct the dual network of the Voronoi graph to obtain the Delaunay network. The edges of the Delaunay network connect the seed points of adjacent Voronoi units. S25. Perform a safety screening on the Delaunay network, delete edges that cross no-fly zones, and retain edges that do not violate airspace control rules to obtain the Voronoi cell set and the dual Delaunay network.

3. The method for constructing a three-dimensional digital airway network for urban low-altitude airspace according to claim 1, characterized in that, Step S3 includes: S31. Integrate geographic information system data, real-time meteorological data, population density data and communication signal strength data to divide the urban low-altitude airspace into three-dimensional grid units. The grid unit size is 10 meters by 10 meters in the horizontal direction and 20 meters in the vertical direction. S32. Calculate the comprehensive risk assessment value for each grid cell. The comprehensive risk assessment function is defined as follows: Where q is the location point in the spatial domain; Let the obstacle risk function be used. For population density risk function; For meteorological conditions risk function; For communication quality risk function; The Voronoi topological risk function; These are the corresponding weighting coefficients; S33. Normalize the comprehensive risk assessment value to the range of [0,1] to obtain the three-dimensional risk field distribution of the airspace covering the entire city's low-altitude airspace.

4. The method for constructing a three-dimensional digital airway network for urban low-altitude airspace according to claim 1, characterized in that, Step S4 includes: S41. According to the aircraft classification standards and altitude management regulations, the low-altitude airspace from 0m to 1000m is divided into eleven altitude layers, with different altitude layers corresponding to different types of aircraft or flight missions. S42. For Delaunay network nodes, replicate and project them vertically according to eleven height layers to form a multi-layer three-dimensional route node network. The nodes in each height layer maintain their Voronoi topology on the horizontal plane. S43. For adjacent altitude layers, vertical connecting edges are set in traffic demand hotspot areas to allow aircraft to switch altitude layers under the premise of complying with the rules. The vertical connecting edges connect the node pairs corresponding to the horizontal positions in two adjacent altitude layers to obtain a three-dimensional network node set. S44. Map the calculated three-dimensional risk field onto the set of three-dimensional network nodes. For each node, its risk attribute is equal to the comprehensive risk assessment value of the grid cell corresponding to the three-dimensional spatial coordinates of that node. For each edge, its risk attribute is obtained by sampling along the edge path and calculating the weighted average value.

5. The method for constructing a three-dimensional digital airway network for urban low-altitude airspace according to claim 4, characterized in that, Step S5 includes: S51. Define the node potential function for each node in the set of nodes in the three-dimensional network: Where i is the node index; j is the neighbor node index; Let be the set of neighboring nodes of node i in the Delaunay network; For local potential energy, Let be the normalized traffic demand at node i. For the risk attribute value of node i, These are preset coefficients; Let be the local potential energy of neighbor node j; For the elements of the graph Laplace matrix; The Laplace diffusion coefficient is used. The integral of the risk gradient field within the Voronoi unit where node i is located; This represents the risk gradient impact coefficient. S52. Define the edge potential function for each edge in the set of nodes in the 3D network: in, The length of the edge between nodes i and j; The edge risk attribute between nodes i and j; The current flow of the edge between nodes i and j; Let i be the edge capacity between nodes i and j; These are the weighting coefficients; S53. The Dijkstra algorithm is used to calculate the path with the minimum potential energy from the starting point to the ending point on the three-dimensional network node set. The total potential energy of the path is equal to the sum of the potential energy of all nodes on the path and the potential energy of all edges. The path with the lowest potential energy is used as the main path and the path with the second lowest potential energy is used as the backup path. S54. Based on the traffic demand density of each region, the air route network is densified and refined. For regions where the normalized traffic demand is greater than a preset threshold, parallel feeder routes are inserted to form a three-level hierarchical structure: trunk routes connecting major urban areas, feeder routes covering the internal areas, and connecting routes connecting specific take-off and landing points, thus obtaining the initial air route network.

6. The method for constructing a three-dimensional digital airway network for urban low-altitude airspace according to claim 1, characterized in that, Step S7 uses an improved NSGA-II non-dominated sorting genetic algorithm, the specific steps of which include: S71. Generate an initial path population with a population size of N paths by utilizing the similarity of node embedding vectors in the low-dimensional feature space. S72. For each path in the current population, calculate three objective function values, including flight distance, energy consumption, and route congestion. At the same time, calculate the fitness function that includes Wasserstein distance. S73. Perform non-dominated sorting, divide the population into multiple front layers according to Pareto dominance, and calculate the crowding distance for the paths in each front layer. S74. Select the population based on the frontier level and crowding distance. Prioritize the retention of paths with lower frontier levels and paths with higher crowding within the same frontier level. Select the first M paths to enter the mating pool. S75. Perform crossover operation to generate offspring paths. Randomly select two parent paths from the mating pool, swap the node sequences of the two parent paths within the randomly selected interval to generate two offspring paths, and insert the shortest path segment to repair the break. S76. Perform mutation operation to mutate the child path with adaptive mutation probability. The mutation operation randomly selects a node in the child path and replaces it with a random neighbor node of that node in the Delaunay network. S77. Perform constraint checks on the generated child paths, checking for altitude layer constraints, no-fly zone constraints, and turning radius constraints. If the constraints are violated, apply a penalty to the fitness function. S78. Add the offspring paths to the population and enter the next iteration. When the change of the hypervolume index of the first front is less than the threshold or the maximum number of iterations is reached for ten consecutive generations, the iteration is terminated. Output the first front as the Pareto optimal path set. Select a suitable equilibrium solution from the Pareto front according to the aircraft mission preference to obtain the optimized route network.

7. The method for constructing a three-dimensional digital airway network for urban low-altitude airspace according to claim 6, characterized in that, The fitness function including the Wasserstein distance in step S72 is defined as follows: in, A normalized multi-objective weighted sum; The topological influence weighting coefficient; The 2-Wasserstein distance between network traffic distributions before and after path selection; The Wasserstein distance is calculated in discrete form, and the formula is as follows: Where m and n are edge indices; For the transmission matrix; For the transmission matrix The element in the m-th row and n-th column represents the flow volume from edge m to edge n; Let m be the embedding vector of the node corresponding to edge m and edge n; The constraints are and as well as ,in Let $m$ be the flow rate before path selection and $n$ be the flow rate after path selection, respectively, and $k$ be the edge index used to sum the total flow rate over all edges.

8. The method for constructing a three-dimensional digital airway network for urban low-altitude airspace according to claim 6, characterized in that, In step S71, the node embedding vector is calculated using a graph convolutional neural network. The specific steps include: The node input feature vector is constructed, which includes spatial coordinates, potential energy value, risk value, demand, Minkowski functional vector of Voronoi element and topological feature vector. The Minkowski functional vector includes four components: volume of Voronoi element, half of surface area, average width and Euler characteristic number. Applying three layers of graph convolution operations to the Delaunay network, each graph convolution operation is defined as follows: Where b is the graph convolutional layer index, with values ​​from 1 to 3; It is a symmetric normalized adjacency matrix; This is the feature matrix of the (b-1)th layer; Let b be the learnable weight matrix of the b-th layer; It is the ReLU activation function; After three layers of graph convolution, a node embedding matrix is ​​obtained. The embedding vector of each node retains the topological location information of the node in the Delaunay network, the geometric information of the Voronoi unit, and the attribute information of the node.

9. A method for constructing a three-dimensional digital airway network for urban low-altitude airspace according to claim 4, characterized in that, The height layer division includes: 0-20m, near-ground airspace; 20-50m, restricted airspace for consumer-grade aircraft; 50-120m, suitable airspace for consumer-grade aircraft; 120-150m, restricted airspace for light and small aircraft; 150-200m, suitable airspace for light and small aircraft; 200-230m, restricted airspace for medium-sized aircraft; 230-300m, suitable airspace for medium-sized aircraft; 300-400m, restricted airspace for large aircraft; 400-600m, suitable airspace for large aircraft; 600-700m, restricted airspace for manned aircraft at high speeds; 700-1000m, suitable airspace for manned aircraft at high speeds.

10. A system for constructing a three-dimensional digital airway network for low-altitude urban areas, characterized in that: The system is used to implement the method as described in any one of claims 1-9.

Images