Urban low-altitude corridor-oriented flying car fine route generation method and system

Abstract

The application relates to the technical field of low-altitude aircraft route planning and urban air traffic, and provides a fine route generation method and system for a flying car in a city low-altitude corridor. The method establishes a three-dimensional dynamic map by loading static no-fly zones, three-dimensional building models, real-time meteorological wind fields and multi-source data of take-off and landing point capacity, calculates the no-fly zone index potential field obstacle avoidance cost, building position classification wind field cost and queuing theory take-off and landing point capacity cost of each grid unit, and fuses to generate a multi-constraint cost map, adopts a hybrid A-star search algorithm to search an optimal path, adopts a B-spline curve smoothing process and combines a dynamic window method to execute local obstacle avoidance and generate a space-time trajectory point sequence, and triggers dynamic re-planning when a sudden no-fly zone or insufficient take-off and landing point capacity is detected. The method has the beneficial effects of high route safety, low energy consumption, good flight stability and strong execution robustness.

Description

Technical Field

[0001] This application relates to the fields of low-altitude aircraft route planning and urban air traffic technology, and in particular to a method and system for generating refined flight routes for flying cars oriented towards urban low-altitude corridors. Background Technology

[0002] Existing methods for generating flight routes for flying cars in urban low-altitude corridors have several shortcomings, including: 1. Existing methods mostly use two-dimensional planar path planning, which do not fully consider the impact of the three-dimensional building form of the city and the differences in wind field at different altitudes on flight energy consumption and stability. When the flight path is flying between urban buildings, it is easy to encounter strong winds in the funnel effect and turbulence in the leeward area, resulting in high energy consumption and poor flight stability. 2. Existing methods for obstacle avoidance in no-fly zones mostly employ simple hard constraints, only setting the interior of the no-fly zone as impassable, without setting safe buffer distances and gradient costs at the boundaries of the no-fly zone. This results in flight paths being too close to the boundaries of the no-fly zone, leading to insufficient safety margins. 3. Existing methods do not take into account the real-time remaining capacity and queuing time of take-off and landing points. They assume that take-off and landing points are available when planning routes. In actual execution, the saturation of take-off and landing point capacity often leads to circling and waiting in the air or forced diversion, resulting in poor route execution feasibility. 4. Existing methods for path smoothing often use simple polylines or circular arcs for transition, without combining the maximum curvature constraint of the flying car for smoothing. As a result, the generated flight path has abrupt curvature changes, and the flying car cannot track it accurately. 5. Existing methods lack an online dynamic replanning mechanism, and cannot adjust the flight path in a timely manner when sudden no-fly zones or changes in the capacity of take-off and landing points occur during flight, which may easily lead to mission failure or safety risks. 6. Existing methods do not consider the switching characteristics between the vertical take-off and landing mode and the cruise mode of the flying car. The trajectory is generated by a single dynamic model throughout the process. The trajectory is discontinuous at the mode switching point, which poses a risk of attitude change in actual execution.

[0003] Therefore, a method is urgently needed to solve at least one of the above problems. Summary of the Invention

[0004] This application provides a method and system for generating refined flight routes for flying cars in urban low-altitude corridors, aiming to solve the various shortcomings of existing methods for generating flight routes for flying cars in urban low-altitude corridors.

[0005] In a first aspect, embodiments of this application provide a method for generating refined flight routes for flying cars in urban low-altitude corridors, the method comprising: Load static no-fly zones, 3D building models, real-time meteorological wind fields, and multi-source data on take-off and landing point capacity of urban low-altitude corridors, establish a 3D dynamic map, and divide it into basic grid units; The no-fly zone obstacle avoidance cost, building wind field cost, and take-off and landing point temporal capacity cost are calculated for each grid cell. These costs are then fused to generate a multi-constraint cost map. The no-fly zone obstacle avoidance cost is calculated by setting the grid inside the no-fly zone as impassable and calculating the cost of the grid within the safe buffer distance according to the exponential potential field. The building wind field cost is calculated by classifying the grid relative to the building position and calculating the energy cost and stability cost. The take-off and landing point temporal capacity cost is calculated by queuing theory model to calculate the waiting time cost. The optimal path is searched on the multi-constraint cost map using a hybrid A* search algorithm, and the optimal path is smoothed using B-spline curves. Local obstacle avoidance is performed in combination with the dynamic window method to generate a spatiotemporal trajectory point sequence. When a sudden no-fly zone or take-off and landing point capacity is detected to be lower than the preset threshold, the data is reloaded and the route replanning is performed.

[0006] Secondly, this application provides a refined flight path generation system for flying cars oriented towards urban low-altitude corridors, the system comprising: The data loading unit is used to load static no-fly zones, 3D building models, real-time meteorological wind fields, and multi-source data on take-off and landing point capacity for urban low-altitude corridors, to create a 3D dynamic map and divide it into basic grid units. The cost calculation unit is used to calculate the no-fly zone obstacle avoidance cost, building wind field cost, and take-off and landing point temporal capacity cost for each grid cell, and merge them to generate a multi-constraint cost map. The no-fly zone obstacle avoidance cost sets the grid inside the no-fly zone as impassable and calculates the cost of the grid within the safe buffer distance according to the exponential potential field. The building wind field cost calculates the energy cost and stability cost according to the grid relative to the building position. The take-off and landing point temporal capacity cost calculates the waiting time cost according to the queuing theory model. The route planning unit is used to search for the optimal path on a multi-constraint cost map using a hybrid A* search algorithm, smooth the optimal path using B-spline curves, perform local obstacle avoidance using a dynamic window method, and generate a spatiotemporal trajectory point sequence. When a sudden no-fly zone or take-off and landing point capacity is detected to be lower than a preset threshold, the data is reloaded and route replanning is performed.

[0007] The method provided in this application loads multi-source data, including static no-fly zones, 3D building models, real-time meteorological wind fields, and take-off and landing point capacity, into a unified coordinate system to establish a 3D dynamic map and divide it into basic grid units. This method breaks through the limitations of traditional 2D planar planning and fully considers the constraints of urban 3D space and dynamic environmental changes.

[0008] It should be understood that the above general description and the following detailed description are exemplary and explanatory only, and do not limit this application. Attached Figure Description

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

[0010] Figure 1 This is a schematic flowchart illustrating the steps of a method for generating refined flight routes for flying cars in urban low-altitude corridors, as provided in an embodiment of this application. Figure 2 This is a complete flowchart of an overall method provided in an embodiment of this application; Figure 3 This is a flowchart of a building wind field cost calculation method provided in one embodiment of this application; Figure 4 This is a flowchart illustrating the calculation of obstacle avoidance costs in a no-fly zone, provided in one embodiment of this application. Figure 5 This is a flowchart of a takeoff and landing point timing capacity cost calculation provided in one embodiment of this application; Figure 6 This is a schematic block diagram of the structure of a refined flight route generation system for flying cars in urban low-altitude corridors provided in one embodiment of this application; Figure 7 This is a schematic block diagram of the structure of a computer device provided in an embodiment of this application.

[0011] It should be understood that the above general description and the following detailed description are exemplary and explanatory only, and do not limit this application. Detailed Implementation

[0012] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0013] The flowchart shown in the attached diagram is for illustrative purposes only and does not necessarily include all content and operations / steps, nor does it necessarily have to be performed in the order described. For example, some operations / steps can be broken down, combined, or partially merged, so the actual execution order may change depending on the actual situation.

[0014] It should be understood that, in order to clearly describe the technical solutions of the embodiments of the present invention, the terms "first" and "second" are used in the embodiments of the present invention to distinguish identical or similar items with essentially the same function and effect. Those skilled in the art will understand that the terms "first" and "second" do not limit the quantity or execution order, and the terms "first" and "second" are not necessarily different.

[0015] It should be understood that the terminology used in this specification is for the purpose of describing particular embodiments only and is not intended to limit the scope of the application. As used in this specification and the appended claims, the singular forms “a,” “an,” and “the” are intended to include the plural forms unless the context clearly indicates otherwise.

[0016] It should also be understood that the term “and / or” as used in this application specification and the appended claims means any combination of one or more of the associated listed items and all possible combinations, and includes such combinations.

[0017] Please refer to Figure 1 This application provides a method for generating refined flight routes for flying cars in urban low-altitude corridors. The provided method includes steps S101 to S103, detailed below: Step S101. Load the static no-fly zone, 3D building model, real-time meteorological wind field and take-off and landing point capacity multi-source data of the urban low-altitude corridor, establish a 3D dynamic map and divide it into basic grid units.

[0018] Specifically, this involves loading multi-source data on static no-fly zones in urban low-altitude corridors, 3D building models, real-time meteorological wind fields, and take-off and landing point capacity. Static no-fly zone data is loaded from the no-fly zone database of the airspace control department, including permanent no-fly zones (such as airport airspace, military control zones, airspace above government offices, and core airspace above densely populated areas) and temporary no-fly zones (such as activity control zones, construction airspace, and temporary flight restriction zones). Each no-fly zone record includes spatial boundaries (latitude and longitude coordinates or planar coordinates of the polygon vertices), altitude range (lower and upper limits, in meters), effective time (start and end times; permanent no-fly zones are effective 24 hours a day), no-fly type (no overflying, altitude restriction, speed restriction, etc.), and control level. The 3D building model data is loaded from the city's geographic information system or building information model database. It includes the geographical location of each building (coordinates of the vertices of the base polygon), height (elevation of the building's top), shape parameters (building shape type such as rectangle, L-shape, arc, facade angle, roof type such as flat roof, pitched roof), and building material properties (surface roughness coefficient used for wind flow calculations). The spatial accuracy of the building model typically reaches the meter level, covering all buildings within a certain range on both sides of the city's low-altitude corridor. Real-time meteorological wind field data is loaded from meteorological service interfaces or airborne meteorological sensors. It includes wind direction (angle, 0 degrees clockwise from true north), wind speed (meters per second), and turbulence intensity (a percentage of turbulence intensity, representing the amplitude of wind speed fluctuations) at different altitude levels (such as 100 meters, 200 meters, 300 meters, 500 meters, etc., the low-altitude corridor coverage heights). The temporal resolution of the wind field data is typically at the minute or hour level, and the spatial resolution is gridded (e.g., a 1 km × 1 km grid). Wind field parameters at any location and altitude can be obtained through spatiotemporal interpolation. Takeoff and landing point capacity data is loaded in real time from the takeoff and landing point management system. This includes the geographical location (latitude, longitude, and altitude) of each takeoff and landing point (vertical takeoff and landing field, rooftop apron, and ground takeoff and landing area), remaining capacity (number of currently available takeoff and landing slots), queue length (number of flying cars waiting to take off or land), average service time (average time per takeoff and landing, in minutes), and takeoff and landing point type (vertical takeoff and landing field, multi-rotor apron, and fixed-wing runway mixed field). The takeoff and landing point data is typically updated at the second or minute level to ensure that the latest capacity status is used when planning flight routes. After loading multi-source data, data cleaning and format unification are performed, including coordinate system unification (all data is converted to the same plane coordinate system or geographic coordinate system), timestamp alignment, missing data interpolation, and outlier filtering, providing clean and consistent input data for subsequent 3D dynamic map construction.

[0019] A 3D dynamic map is created and divided into basic grid units. The 3D dynamic map is bounded by the spatial range of the urban low-altitude corridor, covering a certain width on both sides of the corridor horizontally (e.g., 500 meters to 1 kilometer on each side of the corridor's centerline), and vertically covering the altitude from the ground to the upper limit of the low-altitude corridor (e.g., 120 meters to 500 meters, set according to urban low-altitude control regulations). The 3D dynamic map is divided into basic grid units according to the set spatial resolution. The horizontal resolution (length and width of the grid unit) is set according to the flight path accuracy requirements (e.g., 5 meters, 10 meters, or 20 meters; higher accuracy requires higher resolution and more computation). The vertical resolution (height of the grid unit) is set according to the differences in wind fields at different altitudes and the accuracy of flight altitude (e.g., 10 meters, 20 meters, or 50 meters). Each basic grid unit is a 3D cube (or cuboid) with a unique spatial index number (e.g., the 3D coordinate index (i,j,k) corresponds to the horizontal grid point (i,j) and the altitude layer k). Each basic grid cell stores multi-source attribute data for that location, including: no-fly zone attributes (whether it is located inside a no-fly zone, its no-fly zone number, and distance to the nearest no-fly zone boundary), building attributes (whether it is located inside a building, distance to the nearest building, and the location type relative to the building, such as open area / canyon area / roof area / leeward area), wind field attributes (wind direction, wind speed, and turbulence at the center of the grid, obtained through 3D spatiotemporal interpolation of the wind field grid data), and take-off and landing point attributes (distance to the nearest take-off and landing point, the remaining capacity of the take-off and landing point, and the queuing time). The 3D dynamic map supports dynamic updates: when real-time meteorological wind field data is updated, the wind field attributes of all grid cells are refreshed; when take-off and landing point capacity data is updated, the take-off and landing point attributes of the relevant grid cells are refreshed; when a temporary no-fly zone notice is received, the no-fly zone attributes of the relevant grid cells are updated. This dynamic update mechanism ensures that the 3D dynamic map always reflects the latest environmental status, providing a real-time and accurate environmental model for subsequent cost calculations and route searches. Through a complete process of multi-source data loading and cleaning, three-dimensional spatial grid division, multi-attribute storage and dynamic updating of grid units, an accurate, real-time, and multi-dimensional three-dimensional dynamic map of urban low-altitude corridors was established, laying a reliable environmental modeling foundation for subsequent multi-constraint cost calculation and route search.

[0020] Step S102. Calculate the no-fly zone obstacle avoidance cost, building wind field cost, and take-off and landing point temporal capacity cost for each grid cell, and fuse them to generate a multi-constraint cost map. The no-fly zone obstacle avoidance cost sets the grid inside the no-fly zone as impassable, and calculates the cost of the grid within the safe buffer distance according to the exponential potential field. The building wind field cost calculates the energy cost and stability cost according to the grid relative to the building position. The take-off and landing point temporal capacity cost calculates the waiting time cost according to the queuing theory model.

[0021] Specifically, the obstacle avoidance cost in no-fly zones is calculated for each grid cell. Each basic grid cell in the 3D dynamic map is traversed, and the obstacle avoidance cost calculation for each grid cell is performed. First, it is determined whether the center point of the grid cell is located inside any no-fly zone: the 3D coordinates (horizontal and altitude) of the grid cell center point are compared with the spatial boundaries (horizontal polygons and altitude ranges) of each no-fly zone using an inclusion test (point within the polygon test + altitude range judgment). If the center point is located within the horizontal boundary of any no-fly zone, the altitude is within the altitude range of that no-fly zone, and the current time is within the effective time of that no-fly zone, then the grid cell is determined to be inside the no-fly zone. The obstacle avoidance cost is set to positive infinity (or a maximum value such as 10 to the power of 9), indicating that the grid is absolutely impassable, and the search algorithm will not choose a path through that grid. If a grid cell is not located within any no-fly zone, calculate the shortest Euclidean distance d from the center point of that grid cell to the nearest no-fly zone boundary: Traverse all no-fly zones, calculate the shortest distance (point-to-line segment distance) from the grid center point to each side of the horizontal boundary polygon of each no-fly zone and the distance to the altitude range boundary, and take the shortest distance in three-dimensional space as d, in meters. Set the safety buffer distance D_safe as a configurable parameter (set according to the positioning accuracy, control accuracy, and safety margin requirements of the flying car, such as 20 meters, 50 meters, or 100 meters; higher safety requirements require a larger D_safe). Determine the relationship between d and D_safe: If d is less than D_safe (the grid is within the safe buffer zone of the no-fly zone boundary), then calculate the obstacle avoidance cost of the no-fly zone according to the exponential potential field. The cost formula is C_no_fly=A×exp(-σ×d), where A is the maximum cost coefficient (when d approaches 0, the cost approaches A. A is set to a large value, such as 1000, to ensure that the cost of the grid adjacent to the no-fly zone boundary is very high), and σ is the attenuation factor (controlling the exponential attenuation rate. The larger σ is, the faster the cost attenuates with increasing distance. For example, σ=0.05 to 0.1 per meter corresponds to the cost at D_safe attenuating to about 1% to 5% of A). The characteristic of the exponential potential field is that the closer to the no-fly zone boundary, the higher the cost, and the cost decreases exponentially with increasing distance, forming a smooth gradient cost field, guiding the flight path to naturally move away from the no-fly zone boundary without causing a sudden change in cost. If d is greater than or equal to D_safe (the grid is outside the safe buffer zone), then the obstacle avoidance cost of the no-fly zone is set to zero, indicating that the grid is not affected by the no-fly zone constraint. The combination of hard constraints (no passage within the no-fly zone) and soft constraints (exponential potential field within the buffer zone) in the obstacle avoidance cost of the no-fly zone not only strictly ensures that the flight path does not enter the no-fly zone, but also guides the flight path to maintain a safe distance from the no-fly zone through gradient cost, thus avoiding the problem of insufficient safety margin when the flight path is close to the boundary of the no-fly zone in the traditional hard constraint method.

[0022] Calculate the building wind field cost for each grid cell. Traverse each basic grid cell in the 3D dynamic map and perform building wind field cost calculation for each grid cell. First, determine the location type of the grid cell relative to the building: Based on the spatial relationship between the center point of the grid cell and the surrounding buildings (distance to the nearest building, orientation relative to the building facade, and relationship between grid height and building height), the location type is divided into four categories: Open area (grid is far from all buildings, the distance to the nearest building is greater than 2 to 3 times the building height, and it is not affected by building bypass flow), Building canyon area (grid is located above a street or passage between two or more high-rise buildings, and the height of the buildings on both sides is significantly higher than the grid height, forming a canyon terrain), Building roof area (grid is located near the top of the building, and its height is within 0 to 1 times the building height above the building top, and it is affected by airflow separation and downwash from the building top), Building leeward area (grid is located on the leeward side of the building (relative to the prevailing wind direction), and its distance from the building facade is within 1 to 3 times the building height, and it is affected by building wake and turbulence). For different location types, corresponding physical models are used to calculate wind field costs: In open areas, the building wind field cost is calculated based on the basic atmospheric wind field, that is, the real-time meteorological wind field parameters (wind direction, wind speed) of the grid are used directly without building disturbance correction. The wind field cost is the sum of the basic energy consumption cost and the basic stability cost (the value is small and serves as the cost benchmark). In building canyon areas, the wind speed amplification factor and energy cost are calculated based on the funnel effect model: the funnel effect refers to the phenomenon that airflow is accelerated in a narrow channel due to the reduction of cross-sectional area. The wind speed amplification factor is calculated based on the ratio of canyon width (distance between buildings on both sides) to building height. The formula can be expressed as k_venturi=1+α×(H / W), where H is the building height, W is the canyon width, and α is the funnel effect coefficient (such as 0.3 to 0.5). The amplified actual wind speed v_real=k_venturi×v_base (basic wind speed); the angle between the actual wind speed and the heading of the flying car determines whether it is a tailwind, headwind, or crosswind. Tailwind (angle between wind direction and heading less than 45 degrees) reduces the ground speed requirement of the flying car and reduces energy consumption, with a negative energy cost (benefit); headwind (angle greater than 135 degrees) increases energy consumption, with a positive energy cost that increases with increasing wind speed; crosswind (angle between 45 and 135 degrees) increases the energy consumption for attitude correction and the difficulty of lateral control of the flying car, with a moderately positive energy cost; energy cost C_energy=k_e×v_real×cos(θ) (θ is the angle between wind direction and heading, cos is positive and the cost is negative (benefit) when tailwind, and negative and the cost is positive when headwind), k_e is the energy cost coefficient.The vertical wind shear intensity and stability cost for the building roof area are calculated using a downwash model: airflow separation at the building roof generates a downward downwash and a vertical wind speed gradient (vertical wind shear). The vertical wind shear intensity is calculated based on the distance difference between the grid height and the building roof; the closer the distance, the stronger the wind shear. The formula can be expressed as S_windshear=β×(H_building / (h_grid-H_building)), where h_grid is the grid height, H_building is the building height, and β is the wind shear coefficient. The stability cost C_stability=k_s×S_windshear, where k_s is the stability cost coefficient. Vertical wind shear can cause altitude control disturbances and attitude instability in flying cars, and the cost increases with the increase of wind shear intensity. The turbulent kinetic energy and stability cost in the leeward area of ​​a building are calculated using a high-turbulence model. A Karman vortex street and turbulent wake are formed on the leeward side of the building. The turbulent kinetic energy is calculated based on the distance from the grid to the leeward side of the building and the building height; the closer the distance and the taller the building, the greater the turbulent kinetic energy. The formula can be expressed as TKE=γ×H_building×exp(-δ×d_leeward), where d_leeward is the distance to the leeward side of the building, γ is the turbulence intensity coefficient, and δ is the turbulence attenuation coefficient. The stability cost C_stability=k_t×TKE, where k_t is the turbulence cost coefficient. High turbulence can cause attitude jitter and trajectory deviation in flying cars, and the cost increases with the increase of turbulent kinetic energy. The comprehensive cost of the building wind field C_wind=C_energy+C_stability (a weighted sum of energy cost and stability cost, with configurable weights for each). The building wind field cost for each grid cell includes both an energy gain coefficient (negative when downwind, representing energy savings) and a turbulence risk coefficient (stability cost, representing a decrease in flight stability). By classifying the relative positions of buildings and quantitatively calculating physical models such as the funnel effect, downwash, and high turbulence, we have achieved accurate modeling of the impact of urban three-dimensional building morphology on the energy consumption and stability of flying cars. This enables flight path search to actively select flight paths with low energy consumption and good stability (such as downwind canyon areas, avoiding leeward turbulence areas and roof wind shear areas).

[0023] Calculate the takeoff and landing time-series capacity cost for each grid cell. Traverse each basic grid cell in the 3D dynamic map and perform takeoff and landing time-series capacity cost calculation for each grid cell. First, determine the destination takeoff and landing point or the nearest takeoff and landing point corresponding to the grid cell: for grids near the end of the route, associate them with the destination takeoff and landing point specified by the task; for grids in the middle of the route, associate them with the nearest available takeoff and landing point (as an emergency backup point reference). Obtain the remaining capacity C_rem (number of currently available takeoff and landing slots, integer), queue length L_q (number of flying cars currently waiting for takeoff and landing, integer), and average service time T_service (average time for a single takeoff and landing, in minutes or seconds) of the takeoff and landing point in real time through the application programming interface of the takeoff and landing point management system. Based on the M / M / 1 queuing model (a single-server Poisson arrival and exponential service time queuing model, applicable to single-runway or single-landing scenarios; for multi-runway scenarios, it can be extended to the M / M / c model), the estimated waiting time T_wait is calculated as follows: In the M / M / 1 model, the average waiting time (queue waiting time, excluding service time) is calculated as T_wait = L_q × T_service (Little's Law, the average number of customers in the system equals the average arrival rate multiplied by the average dwell time, simplified here to queue length multiplied by the average service time). A more accurate M / M / 1 waiting time formula is T_wait = (λ × T_service²) / (1 - λ × T_service), where λ is the arrival rate (flights / unit time), which can be inferred from the queue length and service time. For the multi-server M / M / c model, the waiting time formula is more complex, requiring the calculation of the probability that all c servers are busy (Erlang C formula) and then multiplying it by the average service time and dividing by (c - utilization rate). This method uses a simplified but practical T_wait=L_q×T_service as the expected waiting time, with the unit being consistent with T_service.The waiting time T_wait is compared with a preset threshold T_max (a configurable parameter, such as 5 minutes or 10 minutes, representing the maximum acceptable waiting time): If T_wait is less than or equal to T_max, the takeoff and landing point time-series capacity cost C_capacity = w4 × T_wait, where w4 is a time cost weighting coefficient (linearly mapping the waiting time to a cost value, e.g., w4 = 10 cost units / minute, the longer the waiting time, the higher the cost). This linear mapping ensures that the cost is proportional to the time when the waiting time is within an acceptable range; if T_wait is greater than T_max (the waiting time exceeds the acceptable threshold, and the takeoff and landing point capacity is saturated), then nearby alternative takeoff and landing points are searched: [The text abruptly ends here, likely due to an incomplete sentence or missing information.] The system searches for all available takeoff and landing points within a center and a radius of d_max (configurable, e.g., 3 km, 5 km), filtering for backup takeoff and landing points that meet the criteria of sufficient capacity (remaining capacity greater than 0 and queuing time less than T_max). If a backup takeoff and landing point exists, the destination is changed to that backup takeoff and landing point, and the capacity cost is recalculated (the cost may include the increased flight distance cost due to changing the takeoff and landing point). If no available backup takeoff and landing point exists (all nearby takeoff and landing points are saturated), a forced waiting strategy is implemented, setting the takeoff and landing point temporal capacity cost to positive infinity (indicating that the takeoff and landing point is currently unavailable, and the route search needs to consider circling and waiting or selecting other options), while outputting a departure circling and waiting strategy suggestion. The takeoff and landing point temporal capacity cost transforms the real-time capacity status and queuing time of the takeoff and landing point into a cost item for route search, enabling route planning to proactively avoid saturated takeoff and landing points or select takeoff and landing points with shorter waiting times, ensuring the temporal feasibility of the route and avoiding circling and waiting or forced diversion due to takeoff and landing point capacity saturation during actual execution.

[0024] A multi-constraint cost map is generated through fusion. Corresponding fusion weight coefficients are configured for the obstacle avoidance cost in the no-fly zone, the building wind farm cost, and the take-off and landing point time-series capacity cost: No-fly zone obstacle avoidance cost weight w1 (usually large, such as 1.0 to 2.0, because safety constraints have the highest priority), building wind farm cost weight w2 (such as 0.5 to 1.0, balancing energy consumption and stability optimization), and take-off and landing point time-series capacity cost weight w3 (such as 0.3 to 0.8, time-series feasibility constraints). These weight coefficients can be dynamically adjusted according to task priorities (such as safety priority, economic priority, and timeliness priority). The three types of costs for each basic grid cell are weighted and summed according to the fusion weight coefficients to generate the comprehensive cost of each grid cell: C_total = w1 × C_no_fly + w2 × C_wind + w3 × C_capacity. If any cost is positive infinity (e.g., the grid inside the no-fly zone or the take-off and landing point is completely unusable), then the comprehensive cost is also positive infinity (impassable). The comprehensive cost of all basic grid units is combined according to spatial location (3D index (i,j,k)) to generate a multi-constraint cost map. The cost map has a 3D grid structure, with each grid storing the comprehensive cost and the cost of each component (for easy debugging and analysis). The multi-constraint cost map is the search space for the subsequent hybrid A* search algorithm. The search algorithm finds the path with the minimum cumulative comprehensive cost from the starting point to the ending point on the cost map, thereby simultaneously satisfying the no-fly zone safety constraints, building wind field energy consumption and stability optimization, and take-off and landing point timing feasibility constraints. Through independent and accurate calculation of the three types of costs (no-fly zone exponential potential field, building wind field physical model, and take-off and landing point queuing theory) and weighted fusion, a multi-constraint cost map that comprehensively reflects the multiple constraints of urban low-altitude corridors is generated, providing an accurate and multi-dimensional cost space for route search.

[0025] Step S103. Use the hybrid A* search algorithm to search for the optimal path on the multi-constraint cost map, use B-spline curves to smooth the optimal path, combine the dynamic window method to perform local obstacle avoidance, and generate a spatiotemporal trajectory point sequence; when a sudden no-fly zone or take-off and landing point capacity is detected to be lower than the preset threshold, reload the data and perform route replanning.

[0026] Specifically, a hybrid A* search algorithm is used to search for the optimal path on a multi-constraint cost graph. The hybrid A* algorithm is a hybrid search algorithm combining the A* algorithm with kinematic constraints. It introduces the kinematic states (position, velocity, heading angle) of the flying car into the grid state space of traditional A*, ensuring that the generated path naturally satisfies the flying car's kinematic constraints (minimum turning radius, maximum climb / descent rate, speed range), avoiding unexecutable sharp turns or excessive gradients in paths generated by traditional A*. Search initialization: The takeoff point in the flight mission parameters is used as the search starting point (initial state includes the starting point's three-dimensional coordinates, initial velocity, and initial heading angle), and the landing point is used as the search endpoint. The open node set (priority queue, sorted by comprehensive cost value) and the closed node set (expanded node set) are initialized. The starting point is added to the open node set. The actual cost of the starting point g=0, the heuristically estimated cost h=heuristic function(start point, endpoint), and the comprehensive cost f=g+h.Search loop: Select the node with the smallest comprehensive cost f from the open node set as the current expansion node, and move the current expansion node from the open node set to the closed node set; determine whether the current expansion node is within the endpoint neighborhood (3D distance less than a preset threshold, such as 10 meters), if so, the search is successful and the loop is exited; generate the neighboring nodes of the current expansion node according to the kinematic constraints of the flying car: the neighboring nodes are not simple grid adjacencies (such as 26 neighbors), but are generated by forward simulation based on the kinematic model of the flying car (such as a bicycle model or a six-degree-of-freedom model) from the current state. The simulation inputs include changes in heading angle (left turn, straight, right turn, constrained by the maximum turning angle), speed changes (acceleration, deceleration, constant speed, constrained by the maximum acceleration), and altitude changes (climb, level flight, descent, constrained by the maximum climb / descent rate). Each input combination is simulated forward for one time step (e.g., 1 second) to generate a new state (new position, new speed, new heading angle) for a neighboring node. The number of neighboring nodes is determined by the number of input combinations (e.g., 3×3×3=27 neighboring nodes). For each neighboring node, the multi-constraint cost graph network corresponding to its position is checked. The comprehensive cost of a grid is calculated. If the cost is positive infinity (inaccessible, such as within a no-fly zone), the neighboring node is discarded. Otherwise, the actual cost of the neighboring node is calculated as g_new = g_current + movement cost (the movement cost includes the base cost corresponding to the path length plus the comprehensive cost integral of the grids traversed, which can be simplified to the weighted average of the comprehensive costs of each grid on the line segment from the current node to the neighboring node multiplied by the distance). The heuristic estimated cost h_new = heuristic function(neighboring node, destination) is calculated. The heuristic function uses three-dimensional Euclidean distance multiplied by the base cost per unit distance (along with the addition of altitude difference penalty and turning prediction cost) to ensure that the heuristic value is the lower bound of the actual cost (ensuring A* optimality). The comprehensive cost f_new = g_new + h_new. If the neighboring node is not in the closed node set and (not in the open node set or f_new is less than the existing f value of the node in the open node set), the neighboring node is added to or updated to the open node set, and its parent node is recorded as the current expanded node. The above loop of selecting the node with the smallest f, expanding the neighborhood, and updating the open set is repeated until the node is expanded to the destination or the open node set is empty (search failure). After a successful search, the system backtracks from the search endpoint back to the search starting point along the parent node pointers, generating an optimal path node sequence. Each node contains three-dimensional coordinates, velocity, heading angle, and a timestamp (calculated from the cumulative time step). The path node sequence is a discrete sequence of state points that satisfies the kinematic constraints of the flying car and minimizes the cumulative comprehensive cost. The hybrid A* search ensures the executability of the path through neighborhood generation of kinematic constraints, ensures the comprehensive optimization of the path in terms of safety, energy consumption, stability, and temporal feasibility through cost integration of multi-constraint cost graphs, and ensures search efficiency through heuristic function guidance.

[0027] B-spline curves are used to smooth the optimal path. Nodes in the optimal path node sequence are used as control points for the B-spline curve. A B-spline curve is a piecewise polynomial curve with local controllability (moving a control point only affects a few local curve segments) and high-order continuity (C2 continuity, i.e., position, first derivative, and second derivative are all continuous), making it suitable for smoothing discrete path points. The order of the B-spline curve is determined based on the number of control points: cubic B-splines (order 3, i.e., a 4th-order polynomial, C2 continuity) are typically used. When the number of control points is small (e.g., less than 4), the order can be reduced or the number of control points increased. Weighted interpolation is performed on the control points according to the B-spline basis function: the parametric equation of the B-spline curve is P(t)=Σ(N_{i,p}(t)×P_i), where P_i is the coordinate of the i-th control point, N_{i,p}(t) is the i-th p-th B-spline basis function (calculated from the node vector using the De Boer recursive formula), and t is the curve parameter (from 0 to 1); the parameter t is sampled at a set interval (e.g., 0.01) within the range [0,1], and the coordinates of the points on the curve corresponding to each t are calculated to generate a continuous and smooth B-spline curve. The curve passes through or approximates the control points (if the open curve does not pass through the first and last control points, it can be adjusted to a clamped B-spline so that the curve passes through the first and last control points). Curvature continuity is checked on the B-spline curve: the curvature of each sampling point is calculated along the curve (curvature=|r'(t)×r''(t)| / |r'(t)|). 3 (where r' is the first derivative and r'' is the second derivative). The curvature is compared with the maximum turning curvature threshold of the flying car (maximum turning curvature = 1 / minimum turning radius, the minimum turning radius is determined by the flying car's performance parameters, such as 50 meters corresponding to a maximum curvature of 0.02 per meter). When the curvature of a segment of the B-spline curve exceeds the maximum turning curvature threshold of the flying car, it indicates that the curve is turning too sharply and the flying car cannot track it accurately. The position of the corresponding control point needs to be adjusted (moving the control point to the outside of the turn or increasing the control point density), the B-spline curve is regenerated, and the curvature is checked again. Iterative adjustments are made until the curvature of the entire B-spline curve satisfies the kinematic constraints of the flying car (curvature is everywhere less than or equal to the maximum turning curvature). Curvature checking and adjustment ensure that the smoothed curve is within the trackable range of the flying car, avoiding drastic attitude changes caused by sudden curvature changes. The smoothed B-spline curve is a continuous, smooth, C2 continuous three-dimensional path with curvature satisfying constraints, which is more suitable for accurate tracking of the flying car than the original discrete node sequence.

[0028] Local obstacle avoidance is performed using the dynamic window method, generating a spatiotemporal trajectory point sequence. The dynamic window method (DWA) is a velocity-space-based local obstacle avoidance algorithm that considers current velocity and acceleration constraints in real time during the flying car's trajectory execution, searching for the optimal velocity command within a dynamic velocity window to achieve local dynamic obstacle avoidance and trajectory tracking. Samples are taken along a smoothed B-spline curve at set time intervals (e.g., 0.5 seconds or 1 second) to generate an initial spatiotemporal trajectory point sequence. Each spatiotemporal trajectory point contains spatial coordinates (x, y, z), a timestamp t, and velocity information (velocity magnitude and direction). This initial trajectory point sequence serves as the flying car's reference trajectory. During the flight of the flying car, a dynamic speed window is determined based on the current speed v_current and acceleration constraints (maximum acceleration a_max, maximum deceleration a_min). The achievable speed range is [v_current - a_min × Δt, v_current + a_max × Δt] (limited by both maximum and minimum speed v_max), and the achievable heading angle change rate range is [-ω_max, ω_max] (ω_max is the maximum turning angular velocity, determined by the maximum turning radius and speed). The speed window is a two-dimensional combination space of speed magnitude and heading angle change rate. Within the dynamic speed window, all candidate speed combinations are traversed (discrete at a set resolution, such as a speed step of 0.5 m / s and an angular velocity step of 1 degree / s). For each candidate speed combination, a forward simulation is performed: using the current state as the initial condition and the candidate speed magnitude and heading angle change rate as input, the motion trajectory within the set prediction time domain (e.g., 3 to 5 seconds) is simulated forward to obtain a series of position points on the predicted trajectory. Calculate an evaluation score for each predicted trajectory. Evaluation dimensions include: minimum distance to obstacles (the minimum distance from each point on the predicted trajectory to dynamic obstacles (such as other flying cars, temporary obstacles, or sudden no-fly zone boundaries; a larger distance is better, and if the distance is less than a safety threshold, the candidate speed is discarded), deviation from the B-spline reference curve (the vertical distance from the predicted trajectory endpoint to the reference curve; a smaller deviation is better), speed smoothness (the change between the candidate speed and the current speed; a smaller change is better), and progress towards the target (the decrease in distance from the predicted trajectory endpoint to the endpoint; a larger decrease is better). The evaluation values ​​of each dimension are weighted and summed to obtain a comprehensive score for the candidate speeds. The candidate speed combination with the highest comprehensive score (maximum distance to obstacles, minimum deviation from the reference curve, good smoothness, and fast progress) is selected as the current execution speed. The actual motion state of the flying car is updated based on the current execution speed, and the updated state points are added to the spatiotemporal trajectory point sequence (replacing or supplementing the initial reference trajectory points) to form the actual execution spatiotemporal trajectory point sequence.The dynamic window method operates at high frequencies (e.g., 10 Hz to 20 Hz), continuously sensing local dynamic obstacles and adjusting speed commands. This enables the avoidance of local dynamic obstacles while tracking a reference trajectory, generating a spatiotemporal trajectory point sequence that is both globally optimal and locally safe and feasible. The final spatiotemporal trajectory point sequence (containing the spatial coordinates, timestamp, velocity, and heading angle of each point) is output for the flying car's flight control system to execute.

[0029] When a sudden no-fly zone or landing point capacity falls below a preset threshold is detected, the data is reloaded and route replanning is executed. During route execution, the system receives airspace control notices (via the aviation information service interface or pushed by the control department) and landing point status updates (via the landing point management system) in real time. It continuously monitors for sudden no-fly zones (newly issued temporary no-fly zones, active control zones, and emergency control airspace, the spatial scope and effective time of which may affect the currently executed route) or landing point remaining capacity falling below a preset capacity threshold (e.g., the remaining capacity of the destination landing point drops to 0 or the queuing time exceeds T_max). When a sudden no-fly zone or landing point remaining capacity falls below the preset capacity threshold is detected, the route replanning process is triggered: first, the current route execution is frozen (the flying car remains in its current state, hovering or leveling while waiting), and the flying car's current position (3D coordinates), current speed (speed magnitude and heading angle), and remaining battery status (remaining battery percentage or remaining flight time) are recorded. These statuses serve as the initial conditions for replanning. The updated multi-source data is reloaded, including updated static and temporary no-fly zone data (newly added sudden no-fly zones), 3D building model data (usually unchanged, but updated if there are changes such as construction), real-time meteorological wind field data (latest wind field), and real-time take-off and landing point status data (latest capacity and queuing information). The updated data is then integrated into the 3D dynamic map, and the multi-constraint cost map is refreshed (grid costs affected by sudden no-fly zones are updated to be impassable, and grid costs affected by changes in take-off and landing point capacity are updated). Using the flying car's current location as the new search starting point, and the original landing point (if still available) or a backup take-off and landing point (if the original landing point's capacity is unavailable, a nearby available backup take-off and landing point is selected) as the new search endpoint, the flight path planning is re-executed (including a complete process of hybrid A* search, B-spline smoothing, and dynamic window local obstacle avoidance), generating a new spatiotemporal trajectory point sequence from the current location to the new endpoint. The re-planned spatiotemporal trajectory point sequence is then sent to the flying car for execution, and the flying car smoothly transitions from its current state to the new trajectory (the trajectory connection algorithm ensures continuity of position, speed, and acceleration), restoring the flight path execution state. The replanning process should be as short as possible (typically requiring completion within a few seconds to tens of seconds) to reduce flying car waiting time. Incremental replanning strategies (only re-searching affected segments, retaining unaffected flown segments and unaffected future segments) can be used to accelerate the replanning process. Through a complete closed loop of real-time environmental monitoring, trigger condition detection, state freezing, data updating, incremental replanning, trajectory connection, and execution recovery, online adaptive adjustment of flight routes under dynamic environmental changes (sudden no-fly zones, changes in takeoff and landing point capacity) is achieved, ensuring mission continuity and safety.By employing a complete process that combines A* kinematic constraint search, B-spline curvature constraint smoothing, dynamic window local obstacle avoidance, and online dynamic replanning, a spatiotemporal trajectory point sequence that is globally optimal, locally safe, accurately trackable, and dynamically adaptive is generated. This fully meets the combined requirements of safety, economy, feasibility, and accuracy for generating flight routes for flying cars in urban low-altitude corridors.

[0030] In some embodiments, the process of loading multi-source data on static no-fly zones, 3D building models, real-time meteorological wind fields, and take-off and landing point capacity for urban low-altitude corridors, establishing a 3D dynamic map, and dividing it into basic grid units includes: loading a static no-fly zone database and temporary no-fly zone notices, extracting the spatial boundaries, altitude range, and effective time of the no-fly zones; loading urban 3D building outline models, extracting the geographical location, height, and shape parameters of the buildings; loading real-time meteorological wind field data, extracting wind direction, wind speed, and turbulence parameters at different altitude levels; loading real-time status data of take-off and landing points, extracting the geographical location, remaining capacity, queue length, and average service time parameters of the take-off and landing points; and integrating the static no-fly zone data, 3D building model data, real-time meteorological wind field data, and real-time status data of take-off and landing points into a unified coordinate system to establish a 3D dynamic map, and dividing the 3D dynamic map into basic grid units according to a set spatial resolution.

[0031] Load the static no-fly zone database and temporary no-fly zone notices, and extract the spatial boundaries, altitude range, and effective time of the no-fly zones. The static no-fly zone database is loaded from the geographic information service interface of the airspace control department. Each no-fly zone record in the database includes: a unique identifier number for the no-fly zone, the name of the no-fly zone (e.g., the airspace above an airport, a military management zone, the airspace above a government agency, the core airspace above a densely populated area), the type of no-fly zone (permanent no-fly zone, restricted zone, danger zone, controlled zone), the spatial boundary (the latitude and longitude coordinate sequence of the polygon vertices, the number of vertices varies from 3 to dozens depending on the complexity of the boundary, the coordinate system is WGS84 or CGCS2000 geographic coordinate system), the altitude range (lower and upper limits, in meters, relative to mean sea level or relative to ground level, the lower limit is usually 0 meters above ground level, the upper limit is set according to control requirements such as 120 meters, 500 meters, 1000 meters, etc.), the effective time (permanent no-fly zones are effective all day and all year round, temporary no-fly zones have a clear start and end time, accurate to the minute), the control level (no overflying, altitude restriction, speed restriction, entry requires application), and remarks. Temporary no-fly zone notices are loaded in real-time from aviation information service interfaces or air traffic control departments. The notices include the spatial boundaries, altitude range, effective and expiration times, and reasons for establishment (e.g., large-scale events, air shows, construction hoisting, emergency rescue, VIP security). Temporary no-fly zones are time-sensitive; they must be included within the no-fly zone restrictions during their effective period and will automatically be lifted upon expiration. During loading, the no-fly zone data is parsed and format-converted: the polygon vertices of the geographic coordinate system are converted to a planar coordinate system (e.g., UTM projection or a local urban planar coordinate system, or a northeast-sky coordinate system with the flight path origin as the origin), facilitating subsequent distance calculations and grid inclusion testing; the altitude is standardized to an elevation relative to mean sea level; and the effective time is converted to a unified timestamp format (Unix timestamp or ISO8601 format). The parsed no-fly zone data is stored in a no-fly zone spatial index structure (e.g., an R-tree or quadtree spatial index), supporting quick queries of no-fly zones for a specific point or area, improving the efficiency of subsequent inclusion testing and distance calculations.

[0032] Load a 3D urban building outline model and extract the geographical location, height, and shape parameters of the buildings. The urban 3D building outline model is loaded from a city geographic information system database, a building information model database, or an oblique photogrammetry 3D reconstruction model. The model covers all buildings within a certain range on both sides of the city's low-altitude corridor (e.g., 1 km on each side of the corridor's centerline), with data accuracy reaching meter-level or sub-meter-level. Each building's record includes: a unique building identifier, building name or address, geographical location (planar coordinates of the vertices of the base polygon, describing the building's projected outline on the horizontal plane; the number of polygon vertices varies from 4 (rectangular buildings) to dozens (complex-shaped buildings) depending on the building shape), height parameters (elevation of the building's top, in meters, with an accuracy of 0.1 meters; the elevation of the building's base can also be included for calculating relative height), and shape parameters (building shape type classification: rectangular, L-shaped, U-shaped, arc-shaped, circular, irregular; facade tilt). Angles (angle of inclination of vertical or sloping facades); roof type (flat roof, pitched roof, dome, irregular roof); roof accessories (location and height of antennas, cooling towers, elevator machine rooms, etc.); building material properties (surface roughness coefficient, used for wind flow calculations; roughness of concrete buildings is 0.01 to 0.05, glass curtain walls are 0.005 to 0.01, and metal roofs are 0.001 to 0.005); and building function properties (residential, commercial, office, industrial, public facilities, used to determine the pedestrian density and safety requirements around the building). During loading, the building model is simplified and preprocessed: for overly complex building outlines (too many vertices), the Douglas-Puk algorithm is used to simplify the polygons, reducing the number of vertices while maintaining the outline shape, thus reducing subsequent computational load; the building's base polygons are combined with its height to construct a 3D building bounding box or 3D polyhedron model for spatial inclusion testing and distance calculation; a building spatial index (R-tree or 3D mesh index) is established to support quick queries of buildings near a given location. Three-dimensional building models are the foundation for calculating building wind field costs, and their accuracy directly affects the accuracy of wind field costs.

[0033] Real-time meteorological wind field data is loaded, and wind direction, wind speed, and turbulence parameters at different altitude levels are extracted. The real-time meteorological wind field data is loaded from meteorological service interfaces (such as numerical weather prediction services and urban microscale meteorological simulation services) or airborne meteorological sensors (wind speed and direction sensors mounted on flying cars and atmospheric data computers). The wind field data is three-dimensional gridded data: the horizontal direction is a regular grid (such as 1 km × 1 km or 500 m × 500 m resolution), the vertical direction is multiple altitude levels (such as 50 m, 100 m, 150 m, 200 m, 300 m, 500 m, etc., covering the entire altitude range of the low-altitude corridor), and the time dimension is multiple forecast times (such as the current time and the next 1 hour, 2 hours, 3 hours, etc., with a time resolution of 15 minutes or 1 hour). The wind field parameters for each grid point, each altitude layer, and each moment include: wind direction (angle, 0 to 360 degrees, with due north being 0 degrees and rotating clockwise; wind direction refers to the direction from which the wind is coming), wind speed (meters per second, positive value), and turbulence intensity (a percentage of turbulence intensity, representing the proportion of wind speed fluctuations to the average wind speed, such as 5%, 10%, 20%; higher turbulence intensity indicates more unstable winds). It may also include auxiliary parameters such as air temperature, air pressure, and air density (used for aerodynamic calculations). During loading, spatiotemporal interpolation is performed on the wind field data: For any grid cell center location (arbitrary horizontal coordinates and height) and any time in the 3D dynamic map, the wind direction, wind speed, and turbulence parameters at that location, height, and time are obtained from the grid wind field data through trilinear interpolation (bilinear interpolation in the horizontal direction + linear interpolation in the vertical direction) or a more precise inverse distance weighted interpolation; for height layers not covered by the wind field data, estimation is performed through vertical extrapolation (e.g., power-law wind profile v(h)=v_ref×(h / h_ref)^α, where α is the wind shear index, and α is approximately 0.2 to 0.4 in urban areas); interpolation ensures that each grid cell in the 3D dynamic map has corresponding wind field parameters. Real-time wind field data is updated at a set frequency (e.g., every 15 minutes or every hour), and after the update, the wind field attributes and building wind field costs of all grid cells in the 3D dynamic map are refreshed.

[0034] Load real-time status data of takeoff and landing points, and extract parameters such as geographical location, remaining capacity, queue length, and average service time. Real-time status data is loaded in real-time from the takeoff and landing point management system (vertical takeoff and landing field operation management platform, urban low-altitude traffic management system) via an application programming interface. The data update frequency is typically on the order of seconds or minutes (e.g., every 10 seconds or every minute) to ensure that the latest takeoff and landing point capacity status is used during route planning. Each takeoff and landing point record includes: a unique identifier, name or address, geographical location (latitude and longitude coordinates and altitude, converted from planar coordinates for distance calculation), takeoff and landing point type (vertical takeoff and landing field (dedicated to multi-rotor / eVTOL), rooftop helipad, ground takeoff and landing area, fixed-wing runway mixed field, helicopter helipad), takeoff and landing point capacity parameters (total number of takeoff and landing slots (number of flying cars that can be served simultaneously), number of currently available takeoff and landing slots (remaining capacity C_rem), number of currently occupied slots), and queue status parameters (number of flying cars currently waiting for takeoff and landing (queue length)). The system includes parameters for takeoff and landing point data, such as length L_q, estimated waiting time for each flying car in the queue, service efficiency parameters (average service time per takeoff and landing T_service, in seconds or minutes, including the entire process of landing, taxiing, passenger boarding and alighting, refueling, and takeoff; distinguishing between landing and takeoff service times), operational status parameters (open / closed, weather suitability (e.g., closed during strong winds or thunderstorms), equipment status (normal / maintenance / fault)), and refueling facility parameters (whether charging piles or battery swapping stations are available, refueling power, and current refueling queue length, used for mission planning requiring refueling). During loading, the system performs validity checks on takeoff and landing point data: filtering out closed or faulty takeoff and landing points (unavailable), and performing reasonableness checks on remaining capacity and queue length (remaining capacity cannot be negative, and queue length cannot exceed the system's maximum capacity). Abnormal data is marked and replaced with data from the previous moment. Takeoff and landing point data is stored in a spatial index structure, supporting quick location-based queries for nearby available takeoff and landing points. Real-time takeoff and landing point status data is the foundation for calculating takeoff and landing point time-series capacity costs; its real-time performance directly affects the accuracy of route time-series feasibility.

[0035] Static no-fly zone data, 3D building model data, real-time meteorological wind field data, and real-time take-off and landing point status data are integrated into a unified coordinate system to create a 3D dynamic map. The 3D dynamic map is then divided into basic grid units according to a set spatial resolution. Multi-source data are integrated into a unified coordinate system: a unified planar coordinate system (such as a local planar coordinate system with the origin of the flight path or city center as the origin, or a UTM projection coordinate system) and a unified elevation datum (mean sea level elevation) are selected. The polygonal boundaries of the static no-fly zone, the bottom outline and height of the 3D buildings, and the geographical locations of the take-off and landing points are transformed from the original geographic coordinate system (latitude and longitude) to the unified planar coordinate system through map projection transformation (such as Gauss-Kruger projection or Universal Transverse Mercator projection), with the elevation unified as altitude. The grid coordinates of the meteorological wind field are also transformed to the unified planar coordinate system, with the elevation layer unified as altitude. After transformation, all data can be directly used for spatial calculations (distance calculation, inclusion testing, interpolation) under the same spatial reference system. Establish a three-dimensional dynamic map: The spatial range of the three-dimensional dynamic map is set according to the mission area. In the horizontal direction, it covers the flight corridor from the take-off point to the landing point and a certain width on both sides (such as 500 meters to 2 kilometers on both sides of the center line of the corridor, to ensure that all possible detour paths and nearby take-off and landing points are included). In the vertical direction, it covers the altitude from the ground (0 meters above sea level or the lowest altitude of the mission area) to the upper limit of the low-altitude corridor (such as 120 meters, 300 meters or 500 meters, set according to the city's low-altitude control regulations and the performance of the flying car). Dividing into basic grid units: The 3D dynamic map is divided into regular 3D grids according to the set spatial resolution. The horizontal resolution is Δx = Δy (the horizontal side length of the grid unit, set according to the balance between flight path accuracy requirements and computational load, such as 5 meters, 10 meters, and 20 meters; 5 meters is used in urban core areas with high accuracy requirements, and 20 meters is used in open suburbs), and the vertical resolution is Δz (the height of the grid unit, set according to the differences in wind field at different altitudes and flight altitude accuracy, such as 10 meters, 20 meters, and 50 meters; 10 meters is used in urban canyon areas where the wind field changes greatly with altitude, and 50 meters is used in open areas). Each basic grid unit is a 3D cube with a unique 3D index number (i,j,k), where i is the grid point index in the x-direction, j is the grid point index in the y-direction, and k is the grid point index in the z-direction (altitude). The mapping relationship between the index and spatial coordinates is x = x0 + i × Δx, y = y0 + j × Δy, z = z0 + k × Δz (x0, y0, z0 are the coordinates of the map origin).Each grid cell stores the multi-source fusion attributes for that location: no-fly zone attributes (whether it is located inside a no-fly zone, the no-fly zone number, and the distance d to the nearest no-fly zone boundary), building attributes (whether it is located inside a building, the distance to the nearest building, the location type relative to the building (open area / canyon area / roof area / leeward area), and the height and spacing of surrounding buildings), wind field attributes (wind direction, wind speed, and turbulence at the center of the grid, obtained through three-dimensional spatiotemporal interpolation of wind field grid data), and take-off and landing point attributes (distance to the nearest take-off and landing point, the remaining capacity of the take-off and landing point, queue length, average service time, and estimated waiting time). The 3D dynamic map supports a dynamic update mechanism: when meteorological wind field data is updated, the wind field attributes and building wind field costs of all grids are refreshed; when the take-off and landing point status is updated, the take-off and landing point attributes and capacity costs of the relevant grids are refreshed; when a temporary no-fly zone notice is received, the no-fly zone attributes and no-fly zone costs of the relevant grids are updated; when the building model changes (such as adding a new construction building), the building attributes and wind field costs of the relevant grids are updated. Dynamic updates are implemented incrementally (updating only the affected grid areas), avoiding the performance overhead of recalculating the entire map. Through a complete process of unified multi-source data coordinates, 3D spatial range setting, regular grid division, multi-attribute storage of grids, and dynamic incremental updates, an accurate, real-time, multi-dimensional, and dynamically evolving 3D dynamic map of urban low-altitude corridors has been established, providing a reliable environmental modeling foundation for subsequent multi-constraint cost calculations and flight path searches.

[0036] In some embodiments, the calculation of the no-fly zone obstacle avoidance cost, building wind field cost, and take-off and landing point temporal capacity cost for each grid cell, and the fusion to generate a multi-constraint cost map, includes: traversing each basic grid cell in the three-dimensional dynamic map, calculating the no-fly zone obstacle avoidance cost, building wind field cost, and take-off and landing point temporal capacity cost for each basic grid cell respectively; configuring corresponding fusion weight coefficients for the no-fly zone obstacle avoidance cost, building wind field cost, and take-off and landing point temporal capacity cost respectively; performing a weighted summation of the no-fly zone obstacle avoidance cost, building wind field cost, and take-off and landing point temporal capacity cost for each basic grid cell according to the fusion weight coefficients to generate the comprehensive cost value of each basic grid cell; and combining the comprehensive cost values ​​of all basic grid cells according to their spatial location to generate a multi-constraint cost map.

[0037] The algorithm iterates through each basic grid cell in the 3D dynamic map, calculating the no-fly zone obstacle avoidance cost, building wind field cost, and takeoff and landing point temporal capacity cost for each basic grid cell. The traversal uses a 3D nested loop (i from 0 to Nx-1, j from 0 to Ny-1, k from 0 to Nz-1, where Nx, Ny, and Nz are the number of grid cells in the x, y, and z directions, respectively), calculating the three types of costs for each grid cell sequentially. No-fly zone obstacle avoidance cost calculation: For each grid cell center point (x_c, y_c, z_c), firstly, quickly query the list of no-fly zones that may contain the point using the no-fly zone spatial index. For each no-fly zone in the list, perform point testing within the polygon (horizontal direction: ray method or orbital number method to determine if the point is inside the horizontal polygon of the no-fly zone) + altitude range judgment (whether z_c is between the lower and upper limits of the no-fly zone) + effective time judgment (whether the current time is within the effective time of the no-fly zone). If any no-fly zone meets all the conditions, the grid is determined to be inside the no-fly zone, and the cost is set to positive infinity; if it is not inside any no-fly zone, calculate the shortest distance from the point to the boundaries of all no-fly zones (the distance from the point to the line segment on each side of the horizontal polygon of each no-fly zone + the distance to the altitude boundary, taking the three-dimensional shortest distance d). Calculate the cost based on the relationship between d and the safety buffer distance D_safe using an exponential potential field or a zero value. Building wind field cost calculation: For each grid cell, the building list within a certain range (e.g., 3 times the maximum building height) is queried using the building spatial index. The location type (open area / canyon area / roof area / leeward area) is determined based on the spatial relationship between the grid point and surrounding buildings (distance to the nearest building, relationship between grid height and building height, and relationship with the prevailing wind direction relative to the building facade). Energy and stability costs are calculated based on the corresponding physical model (basic atmospheric wind field / pipeline effect / downwash / high turbulence), and summed to obtain the comprehensive building wind field cost. Take-off and landing point time-series capacity cost calculation: For each grid cell, associated take-off and landing points are determined (grids near the route endpoint are associated with destination take-off and landing points; grids along the route are associated with the nearest available take-off and landing points as alternative landing references). The remaining capacity, queue length, and average service time of the take-off and landing point are obtained. The expected waiting time T_wait is calculated using the M / M / 1 queuing theory model. The capacity cost (linear mapping, alternative take-off and landing point search, or positive infinity forced waiting) is calculated based on the relationship between T_wait and the threshold T_max. The three types of cost calculations are independent and can be executed in parallel (e.g., multi-threaded or GPU parallel computation of different grid regions), improving computational efficiency. The calculation results are stored in the attribute structure of each grid cell, containing the three individual cost values ​​and the combined cost value of subsequent calculations.

[0038] Assign corresponding fusion weight coefficients to the no-fly zone obstacle avoidance cost, building wind field cost, and takeoff and landing point time-series capacity cost. The fusion weight coefficient determines the relative importance of the three types of costs in the overall cost; the larger the weight, the higher the priority of that type of constraint in route optimization. The weight coefficients are configured according to the task type and optimization objective: No-fly zone obstacle avoidance cost weight w1: safety constraint, usually set to the highest weight (e.g., 1.0 to 3.0), because flight safety is the most basic requirement, and violations in no-fly zones will lead to serious safety consequences and legal penalties; in safety-priority tasks (such as manned flight, medical emergency), w1 takes a high value (2.0 to 3.0), and in scenarios such as cargo transportation where the safety margin can be appropriately reduced, w1 takes a standard value (1.0 to 1.5). Building wind farm cost weight w2: Energy consumption and stability optimization constraint, usually set to a medium weight (e.g., 0.5 to 1.5), because energy consumption affects range and operating costs, and stability affects passenger comfort and flight safety. In economically prioritized missions (e.g., logistics delivery, long-range flights), w2 takes a high value (1.0 to 1.5) to prioritize reducing energy consumption. In comfort-priority manned missions, the stability cost portion of w2 takes a high value to prioritize avoiding turbulence. Take-off and landing point timing capacity cost weight w3: Timing feasibility constraint, usually set to a medium weight (e.g., 0.3 to 1.0), because take-off and landing point waiting time affects mission timeliness and operational efficiency. In time-priority missions (e.g., emergency rescue, time-limited delivery), w3 takes a high value (0.8 to 1.0) to prioritize take-off and landing points with shorter waiting times. In non-emergency missions, w3 takes a standard value (0.3 to 0.5). The weighting coefficients can be configured as a globally uniform value or as dynamic weights that vary with spatial regions (e.g., w1 automatically increases in densely populated urban core areas with no-fly zones, w2 automatically increases in areas with densely built-up high-rise buildings, and w3 automatically increases in busy take-off and landing point areas). Dynamic weights are automatically adjusted through regional attribute detection. The sum of the weighting coefficients does not require normalization (because the dimensions and numerical ranges of various costs differ; the role of weighting is to adjust relative importance rather than normalize), but testing is necessary to ensure that the generated flight paths under the weighting configuration meet expectations (e.g., the flight path will not be overly circuitous due to excessively high w2, nor will take-off and landing point saturation be ignored due to excessively low w3). The weighting configuration is stored in the system parameter file, supporting quick selection of preset weighting schemes (safety priority scheme, economy priority scheme, timeliness priority scheme, balanced scheme) based on the task type before task startup.

[0039] The no-fly zone obstacle avoidance cost, building wind field cost, and take-off and landing point temporal capacity cost of each basic grid cell are weighted and summed according to the fusion weight coefficients to generate the comprehensive cost of each basic grid cell. The comprehensive cost calculation formula is: C_total(i,j,k)=w1×C_no_fly(i,j,k)+w2×C_wind(i,j,k)+w3×C_capacity(i,j,k), where C_no_fly is the no-fly zone obstacle avoidance cost, C_wind is the building wind field cost, C_capacity is the take-off and landing point temporal capacity cost, and w1, w2, and w3 are the corresponding fusion weight coefficients. When calculating, positive infinity values ​​need to be handled: If any component cost is positive infinity (e.g., C_no_fly=+∞ for a grid inside a no-fly zone, or C_capacity=+∞ for a take-off and landing point that is completely unusable), then the total cost C_total is also set to positive infinity (indicating that the grid is absolutely impassable, and the search algorithm will not choose a path through the grid), to avoid positive infinity participating in arithmetic operations and causing numerical overflow; if all component costs are finite values, then they are summed normally according to the formula. The total cost value is a non-negative real number (the cost of a no-fly zone is non-negative, the energy gain in the wind field cost can be negative but is usually non-negative after summing, and the capacity cost is non-negative, and is non-negative overall after weighting), and the value ranges from zero (a completely safe, low-energy, zero-wait ideal grid) to positive infinity (impassable grid). After calculating the overall cost for each grid cell, the result is stored in the attribute structure of that grid, while retaining the costs of three sub-items (for subsequent debugging, analysis, and visualization; for example, heatmaps of no-fly zone costs, wind farm costs, and capacity costs can be rendered separately to check the rationality of each cost). The calculation of the overall cost can be completed simultaneously with the calculation of the sub-items during grid traversal (the overall cost is calculated immediately after the sub-item calculation is completed for each grid), or the overall cost can be calculated uniformly after all sub-item calculations are completed (for parallelization). The results are consistent in both methods.

[0040] The comprehensive cost value of all basic grid cells is combined according to spatial location to generate a multi-constraint cost map. The multi-constraint cost map is a three-dimensional grid data structure, completely consistent with the grid division of the three-dimensional dynamic map (the same Nx×Ny×Nz grid, the same spatial resolution and coordinate mapping). Each grid cell stores the comprehensive cost value C_total(i,j,k) and three component cost values ​​for that location. The cost map is stored in memory as a three-dimensional array (data type is floating-point, such as float32, positive infinity is represented by the maximum value of float32 or a special marker value), supporting fast random access to the cost value of any grid by three-dimensional index (the cost of frequently querying neighboring grids in the search algorithm). The cost map also provides spatial indexing and region query functions: it can query all grid cost values ​​within a certain spatial region (such as within a cube or sphere) for path cost integration calculation; it can query the set of traversable grids with cost values ​​below a certain threshold for preprocessing the search space. The multi-constraint cost map serves as the search space and cost basis for the hybrid A* search algorithm. The algorithm searches the cost map from the starting point to the endpoint for the path with the minimum cumulative comprehensive cost. The total cost of the path is the sum of the comprehensive costs of each grid it passes through (or the integral along the path segment), ensuring that the generated flight path simultaneously satisfies no-fly zone safety constraints (avoiding high-cost and impassable grids), building wind field energy consumption and stability optimization (selecting low-wind-field-cost grids, i.e., downwind and low-turbulence areas), and take-off and landing point temporal feasibility constraints (selecting low-capacity-cost grids, i.e., take-off and landing points with short waiting times). The cost map supports dynamic updates: when any attribute of the 3D dynamic map is updated (no-fly zone change, wind field update, take-off and landing point capacity change), the component costs and comprehensive costs of the corresponding grids are recalculated, the cost map is partially refreshed, and the search algorithm can perform replanning based on the updated cost map. Through a complete process of independent calculation of three types of costs, weight coefficient configuration, weighted summation, and 3D grid combination storage, a multi-constraint cost map that comprehensively reflects the multiple constraints of urban low-altitude corridors is generated, providing an accurate, multi-dimensional, and dynamically updatable cost space for flight path search.

[0041] In some embodiments, the obstacle avoidance cost in the no-fly zone includes setting the grid inside the no-fly zone as impassable, calculating the cost of the grid within the safe buffer distance using an exponential potential field, calculating the energy cost and stability cost of the building wind field according to the grid's relative building position, and calculating the waiting time cost of the take-off and landing point time-series capacity using a queuing theory model. This includes: determining whether the center point of the basic grid cell is located inside the no-fly zone; if it is, setting the obstacle avoidance cost to impassable; if it is not located inside the no-fly zone, calculating the distance from the basic grid cell to the nearest no-fly zone boundary; if the distance is less than the safe buffer distance, calculating the obstacle avoidance cost according to an exponential potential field; and if the distance is greater than or equal to the safe buffer distance, setting the obstacle avoidance cost to impassable. Zero; Determine the location type of the basic grid cell relative to the building. When the location type is an open area, calculate the building wind field cost according to the basic atmospheric wind field. When the location type is a building canyon area, calculate the wind speed amplification factor and energy cost according to the septum effect model. When the location type is a building roof area, calculate the vertical shear intensity and stability cost according to the downwash flow model. When the location type is a building leeward area, calculate the turbulent kinetic energy and stability cost according to the high turbulence model. Based on the queue length and average service time of the take-off and landing points, calculate the waiting time of the take-off and landing points according to the queuing theory model. Linearly map the waiting time to the take-off and landing point temporal capacity cost. When the waiting time exceeds the preset threshold, search for nearby alternative take-off and landing points and recalculate the take-off and landing point temporal capacity cost.

[0042] The algorithm determines whether the center point of the basic grid cell is located inside the no-fly zone. If it is, the obstacle avoidance cost is set to impassable. If it is not located inside the no-fly zone, the algorithm calculates the distance from the basic grid cell to the nearest no-fly zone boundary. If the distance is less than the safe buffer distance, the obstacle avoidance cost is calculated using an exponential potential field. If the distance is greater than or equal to the safe buffer distance, the obstacle avoidance cost is set to zero. The obstacle avoidance cost in the no-fly zone employs a two-level calculation mechanism combining hard and soft constraints. Hard Constraints (No Passage Within the No-Fly Zone): For each grid cell center point (x_c, y_c, z_c), a list of candidate no-fly zones containing that point is quickly retrieved using the no-fly zone spatial index (R-tree). For each no-fly zone in the candidate list, a triple inclusion test is performed sequentially: Horizontal Inclusion Test: The ray casting method (radiating a horizontal ray from the point to the right and counting the number of intersections with the polygon edges; an odd number indicates the point is inside the polygon, and an even number indicates it is outside) or the orbital number method (calculating the total angle around the point on the polygon boundary; 2π indicates inside, and 0 indicates outside) is used to determine whether (x_c, y_c) is inside the horizontal polygon boundary of the no-fly zone. For no-fly zones with holes (polygons with inner rings), it is also necessary to determine whether the point is inside any hole; Altitude Inclusion Test: It is determined whether z_c is greater than or equal to the lower limit altitude of the no-fly zone and less than or equal to the upper limit altitude of the no-fly zone (the no-fly zone is a cylinder within the altitude range); Time Validity Test: It is determined whether the current system time is greater than or equal to the start time of the no-fly zone's effect and less than or equal to the end time of the effect (this test is always true for permanent no-fly zones). If any no-fly zone passes all triple tests, the grid cell is determined to be inside the no-fly zone. The obstacle avoidance cost C_no_fly for the no-fly zone is set to positive infinity (represented by the maximum floating-point number 1e30 or a special notation), meaning it is impassable. This grid cell will not be selected as a path node in the search algorithm, fundamentally preventing flight paths from entering the no-fly zone. Soft constraint (exponential potential field within a safe buffer distance): If a grid cell is not located inside any no-fly zone, the shortest three-dimensional Euclidean distance d from the center point of that grid cell to the nearest no-fly zone boundary is calculated. Distance calculation method: Traverse all valid no-fly zones (no-fly zones in effect at the current time), and calculate the shortest distance from a point to its three-dimensional boundary for each no-fly zone: Calculate the shortest distance from the point to each side of the polygon in the horizontal direction (the formula for distance from a point to a line segment: for line segment AB and point P, distance = |(PA)×(BA)| / |BA| if the projection is on the line segment, otherwise min(|PA|,|PB|)), and take the minimum horizontal distance among all sides of the polygon; Calculate the distance from the point to the upper and lower altitude limits in the vertical direction (min(|z_c-z_low|,|z_c-z_high|) if z_c is outside the altitude range, otherwise the vertical distance is 0); The three-dimensional shortest distance is the square root of the sum of the squares of the horizontal and vertical distances (√(d_horiz²+d_vert²)); Take the minimum three-dimensional distance among all no-fly zones as d.Set the safety buffer distance D_safe as a configurable parameter, which is comprehensively determined according to the navigation and positioning accuracy of the flying car (such as a GNSS positioning error of 1 to 3 meters), the flight control accuracy (track tracking error of 1 to 2 meters), and safety margin requirements (high margin is required for manned flight), with a typical value of 20 meters to 100 meters. A larger value such as 50 to 100 meters is taken when no-fly zones are dense in the core urban area, and a smaller value such as 20 to 30 meters is taken in open suburbs. Judge the relationship between d and D_safe: if d < D_safe (the grid is within the safety buffer range of the no-fly zone boundary, which is not in the no-fly zone but too close with insufficient safety margin), the no-fly zone obstacle avoidance cost is calculated according to the exponential potential field as C_no_fly = A×exp(-σ×d), where A is the maximum cost coefficient (when d→0, the cost→A, and A is set to a relatively large value such as 500 to 2000, to ensure that the cost of the grid adjacent to the no-fly zone boundary is very high, so that the search algorithm naturally avoids it), and σ is the attenuation factor (controlling the exponential attenuation speed; the larger σ is, the faster the cost attenuates with the increase of distance, with a typical value of σ = 0.02 to 0.1 per meter, corresponding to the cost at d=D_safe attenuating to A×exp(-σ×D_safe); for example, when σ=0.05 and D_safe=50 meters, the cost at the boundary is A×exp(-2.5) ≈ 0.082A, that is, about 8% of the maximum cost, forming a smooth gradient transition). The mathematical characteristic of the exponential potential field is that the cost is the highest (=A) at d=0, decreases monotonically and exponentially with the increase of d, and has a continuous derivative (the cost field is smooth without mutation), so that the route generated by the search algorithm naturally keeps away from the no-fly zone boundary and will not cause path oscillation due to cost mutation. If d ≥ D_safe (the grid is outside the safety buffer range, with a sufficient distance from the no-fly zone and sufficient safety margin), the no-fly zone obstacle avoidance cost is set to zero, indicating that the grid is not affected by the no-fly zone constraint. Through the combination of hard constraints (absolutely impassable inside) and soft constraints (exponential gradient cost within the buffer range), it not only strictly ensures that the route does not enter the no-fly zone, but also guides the route to maintain a safe distance from the no-fly zone through the gradient cost, solving the problem of insufficient safety margin in the traditional hard constraint method where the route flies close to the no-fly zone boundary.

[0043] The system determines the location type of the basic grid cells relative to the building. For open areas, the building wind field cost is calculated using the basic atmospheric wind field. For canyons, the wind speed amplification factor and energy cost are calculated using the funnel effect model. For rooftop areas, the vertical shear intensity and stability cost are calculated using the downwash model. For leeward areas, the turbulent kinetic energy and stability cost are calculated using the high-turbulence model. The building wind field cost is quantitatively calculated based on the spatial location of the grid cells relative to the building, using the corresponding aerodynamic physical model. This enables accurate modeling of the impact of urban 3D building morphology on the energy consumption and stability of flying cars. Location type determination: For the center point of each grid cell, query all buildings within a certain range (e.g., 3 times the height of the tallest surrounding building) using the building space index. Calculate the horizontal distance d_building from the grid point to the nearest building, the relationship between the grid height z_grid and the nearest building height H_building, and the relationship between the grid point's orientation relative to the building facade (which side of the building it's on) and the prevailing wind direction (the angle between the incoming flow direction and the building facade). Based on this, determine the location type: Open area: d_building is greater than 2 to 3 times H_building (the grid is far from all buildings, and the influence of building bypass flow is negligible); Building canyon area: the grid is located between two or more high-rise buildings. Above the passage between buildings, the height of both buildings is greater than z_grid and the building spacing W is less than 2 times H_building (forming a narrow canyon); in the building roof area: z_grid is 0 to 1 times H_building above H_building and its horizontal position is within the projection range of the building's top surface or adjacent to the edge (affected by airflow separation at the building's top); in the building leeward area: the grid is located on the leeward side of the building (relative to the prevailing wind direction, i.e., the building is upstream of the incoming flow direction and the grid is downstream), with a horizontal distance from the building facade between 0.5 and 3 times H_building, and a height between 0 and 1.5 times H_building (affected by building wake turbulence). If the grid simultaneously satisfies multiple types of characteristics (such as being both in a canyon and on the leeward side of a building), the type with the most significant impact is taken (or multiple costs are calculated superimposed and the maximum value is taken).

[0044] In open areas, the building wind field cost is calculated based on the basic atmospheric wind field: real-time meteorological wind field parameters (wind direction θ_wind, wind speed v_wind, turbulence intensity I_turb) obtained by spatiotemporal interpolation at the grid location are directly used, without building disturbance correction; the energy cost C_energy_open=k_e×v_wind×cos(θ_wind-θ_heading), where θ_heading is the heading angle of the flying car (the tangent direction of the flight path at this grid), and cos(θ_wind-θ_heading) is the cosine of the angle between the wind direction and the heading (when the wind is downwind, the angle is small, cos is positive, and the energy cost is negative, i.e., energy consumption benefit). When facing a headwind, the angle is large, cosine is negative, and energy cost is positive, meaning energy consumption increases. When facing a crosswind, cosine is close to zero, and energy cost is small. k_e is the energy cost coefficient (converting wind speed × angle projection into cost value, e.g., k_e = 2 cost units / (m / s)); stability cost C_stability_open = k_s × I_turb × v_wind, where turbulence intensity is multiplied by wind speed to represent the absolute value of turbulence intensity, and k_s is the stability cost coefficient; open area wind field cost C_wind_open = C_energy_open + C_stability_open, which is usually small and serves as a benchmark reference for building wind field cost.

[0045] The wind speed amplification factor and energy cost in urban canyon areas are calculated according to the venturi effect model: The venturi effect refers to the phenomenon that airflow velocity increases and static pressure decreases in narrow channels due to the reduced cross-sectional area. It is the main reason why wind speeds in urban canyons are significantly higher than in open areas. The wind speed amplification factor k_venturi is calculated based on the canyon's geometric parameters: k_venturi = 1 + α × (H_avg / W_canyon) × f(θ_flow), where H_avg is the average height of buildings on both sides of the canyon, W_canyon is the canyon width (distance between building facades on both sides), H_avg / W_canyon is the aspect ratio (the larger the aspect ratio, the stronger the funneling effect; typical urban canyons have an aspect ratio of 1 to 3), α is the funneling effect coefficient (an empirical constant, typically 0.3 to 0.6, which can be calibrated through wind tunnel experiments or numerical simulations), and f(θ_flow) is the wind direction correction function (f = 1 when the wind direction is parallel to the canyon's orientation, the funneling effect is strongest when it is perpendicular to the canyon's orientation, f = 0 when it is perpendicular, and interpolation is performed using sin or cos, such as f = |sin(θ_wind - θ_canyon)|, where θ_canyon is the canyon's orientation angle); the amplified actual wind speed within the canyon is v_canyon = k_venturi × v_wind (base wind speed). The energy cost C_energy_canyon = k_e × v_canyon × cos(θ_wind - θ_heading) (same as the formula for open areas but using amplified wind speeds; energy gains are greater with the wind and energy losses are greater with the headwind, guiding flight paths to choose the windward direction in canyons); the stability cost in canyon areas is relatively small (although airflow accelerates in canyons, its directionality is strong and turbulence is relatively low, unless there are bends), C_stability_canyon = k_s × I_turb × v_canyon × 0.5 (the coefficient 0.5 indicates that the turbulence impact in canyon areas is lower than in leeward areas); the wind field cost in canyon areas C_wind_canyon = C_energy_canyon + C_stability_canyon. The canyon effect model quantitatively characterizes the impact of amplified wind speeds in urban canyons on the energy consumption of flying cars, enabling flight path searches to actively choose windward canyons (higher energy gains) and avoid headwind canyons (higher energy costs).

[0046] The vertical wind shear intensity and stability cost in the building roof area are calculated according to the downwash model: the airflow at the top of the building is separated due to the building's obstruction, forming a downward downwash (airflow deflected downward) and strong vertical wind shear (rapid change in wind speed in the vertical direction) above the building roof. The downwash and vertical wind shear will cause altitude control disturbances, lift abrupt changes and attitude instability of flying cars, which are the main stability risks of flying in the roof area. The vertical wind shear intensity S_windshear is calculated based on the distance difference between the grid height and the top of the building: S_windshear=β×(H_building / Δh)×g(θ_flow), where Δh=z_grid-H_building is the height difference between the grid height and the top of the building (the smaller Δh is, the closer to the roof and the stronger the wind shear; Δh is significant in the range of 0 to 1 times H_building), H_building / Δh is the ratio of the building height to the distance difference (characterizing the wind shear intensity; the closer to the roof, the larger the ratio and the stronger the wind shear), β is the wind shear coefficient (an empirical constant, typically 0.1 to 0.3, which can be calibrated through wind tunnel experiments), and g(θ_flow) is the wind direction correction (the roof separation is strongest when the wind direction is perpendicular to the building facade, g=1; when it is parallel, g=0.3); when Δh≤0 (the grid is inside the building or on the top surface), S_windshear is set to a maximum value (impassable or high cost). The stability cost C_stability_roof = k_s × S_windshear (where k_s is the stability cost coefficient, and the stronger the wind shear, the higher the stability cost). The energy cost in the roof area is relatively small (the main risk is stability rather than energy consumption), C_energy_roof = k_e × v_wind × cos(θ_wind - θ_heading) × 0.3 (the coefficient 0.3 indicates that the energy consumption impact of the roof area is lower than the stability impact); the wind field cost in the roof area C_wind_roof = C_energy_roof + C_stability_roof. The downwash model quantitatively characterizes the impact of vertical wind shear above building roofs on flight stability, enabling flight path searches to actively avoid areas near building roofs (especially the altitude layer immediately adjacent to the roof) and select altitude layers with weaker wind shear for flight.

[0047] The turbulent kinetic energy and stability cost in the leeward area of ​​a building are calculated using a high-turbulence model: On the leeward side (wake region), airflow separation creates a Karman vortex street and a strong turbulent wake. High turbulent kinetic energy, large wind speed fluctuations, and periodic vortex shedding can lead to attitude wobbling, trajectory deviation, and structural fatigue in flying vehicles, representing the main stability risk for flight in the leeward area. Turbulent kinetic energy (TKE) is calculated based on the distance from the grid to the leeward side of the building and the building height: TKE = γ × H_building × exp(-δ × d_leeward) × h(θ_flow), where d_leeward is the horizontal distance from the grid to the leeward side of the building (the downstream facade of the building, relative to the prevailing wind direction), H_building is the building height, γ is the turbulence intensity coefficient (an empirical constant, typically 0.5 to 2.0 m² / s² / m² building height, calibrated through wind tunnel experiments or large eddy simulation), exp... (-δ×d_leeward) represents the exponential decay of turbulence along distance (δ is the decay coefficient, typically 0.3 to 1.0 per meter; turbulence gradually decays with increasing distance in the leeward area, usually returning to the background turbulence level at 3 to 5 times the building height), h(θ_flow) is the wind direction correction (the wake is strongest when the wind is perpendicular to the building facade, h=1; it weakens when the wind is oblique; and there is no obvious leeward area when the wind is parallel, h=0.1); TKE is in square meters per square second (the dimension of turbulent kinetic energy), typically 1 to 5 square meters per square second on the immediate leeward side, decaying with distance. The stability cost C_stability_leeward = k_t × TKE (where k_t is the turbulence cost coefficient, converting turbulent kinetic energy into cost value, e.g., k_t = 10 cost units / (m² / s²), the higher the turbulent kinetic energy, the higher the stability cost). The energy cost in the leeward region is moderate (the average wind speed in the wake region decreases but fluctuates greatly, resulting in high energy consumption uncertainty), C_energy_leeward = k_e × v_wind × cos(θ_wind - θ_heading) × 0.5 + k_e × sqrt(TKE) (adding an energy consumption uncertainty cost related to turbulence intensity); the wind field cost in the leeward region C_wind_leeward = C_energy_leeward + C_stability_leeward. The high-turbulence model quantitatively characterizes the impact of turbulent wake on the leeward side of buildings on flight stability, enabling flight path searches to actively avoid the leeward region of buildings (especially the high-turbulence region immediately adjacent to the leeward side) and select areas with weaker turbulence for flight.

[0048] The comprehensive cost of building wind fields is the sum of energy cost and stability cost for the corresponding location type (C_wind = C_energy + C_stability). The building wind field cost of each grid cell simultaneously reflects the energy consumption level (energy cost, with negative benefits for tailwinds and positive costs for headwinds) and flight stability risk (stability cost, with higher costs for stronger turbulence and wind shear). Through quantitative calculations using four location classifications and corresponding physical models (basic atmospheric wind field, funnel effect, downwash vertical wind shear, and high-turbulence wake), accurate and regional modeling of the impact of urban 3D building morphology on the energy consumption and stability of flying cars is achieved. This enables flight path search to proactively select flight paths with low energy consumption and good stability (such as tailwind valleys, open areas, avoiding roof wind shear areas and leeward turbulence areas), significantly reducing the energy consumption and stability risks of inter-building flights in the city.

[0049] Based on the queue length and average service time of takeoff and landing points, the waiting time at takeoff and landing points is calculated according to the queuing theory model. The waiting time is linearly mapped to the time-series capacity cost of takeoff and landing points. When the waiting time exceeds a preset threshold, nearby alternative takeoff and landing points are searched and the time-series capacity cost of takeoff and landing points is recalculated. The time-series capacity cost of takeoff and landing points is quantitatively calculated based on the queuing theory model, enabling the real-time capacity status and queue delay of takeoff and landing points to be proactively considered during route planning, ensuring the time-series feasibility of routes. The queuing theory model calculates the waiting time: For associated takeoff and landing points (destination takeoff and landing points or nearest alternative takeoff and landing points), its real-time status parameters are obtained: queue length L_q (the number of flying cars currently waiting to take off or land at this takeoff and landing point, including air waiting and ground waiting), average service time T_service (the average full-process service time of a single flying car from landing to takeoff, including approach, landing, taxiing, passenger loading and unloading, refueling, takeoff and departure, etc., in seconds or minutes, updated in real time according to the takeoff and landing point operation statistics), and remaining capacity C_rem (the number of currently available takeoff and landing slots). The estimated waiting time is calculated using the M / M / 1 queuing model (a single-server queuing model that assumes the arrival process is a Poisson process and the service time is an exponential distribution, applicable to single-runway or single-group take-off and landing points; for multi-runway take-off and landing points, it is extended to the M / M / c model, where c is the number of parallel servers, i.e., the number of runways or the number of take-off and landing positions) of queuing theory. In the M / M / 1 model, according to Little's Law (the average number of customers in the system equals the average arrival rate multiplied by the average dwell time), the simplified formula for calculating the average queuing wait time (excluding service time) is T_wait=L_q×T_service (queue length multiplied by average service time, physically meaning the average time each of the previous L_q aircraft occupied T_service, and the total time the current aircraft needs to wait). The more accurate M / M / 1 waiting time formula is T_wait=(λ×T_service²) / (1-λ×T_service), where λ is the arrival rate (flights / unit time), which can be inferred from the queue length and service time (λ=L_q / (T_wait+T_service), which needs to be solved iteratively). When the system utilization ρ=λ×T_service approaches 1, the waiting time increases sharply. For the multi-server (M / M / c) model, the average waiting time is calculated as T_wait = (P_q × T_service) / (c × (1-ρ)), where P_q is the probability that all c servers are busy (calculated using the Erlang C formula), and ρ = λ × T_service / c is the utilization rate of a single server. This method uses the simplified but practical T_wait = L_q × T_service as the estimated waiting time (simple calculation, good real-time performance, and accuracy meeting route planning requirements). For scenarios requiring high precision, the more precise M / M / c formula can be used.

[0050] Linearly map waiting time to takeoff and landing time-series capacity cost: Set a preset waiting time threshold T_max (configurable parameter, representing the maximum acceptable takeoff and landing waiting time for the task, such as 5 minutes, 10 minutes, 15 minutes; for emergency tasks, T_max is small, such as 3 to 5 minutes; for non-emergency tasks, T_max is large, such as 10 to 15 minutes). If T_wait ≤ T_max (waiting time is within the acceptable range), then the takeoff and landing time-series capacity cost C_capacity = w4 × T_wait, where w4 is the time cost weighting coefficient (linearly mapping waiting time to cost value, such as w4 = 5 to 20 cost units / minute; the longer the waiting time, the higher the cost; linear mapping ensures that within the acceptable range, cost is proportional to time, and the search algorithm prioritizes takeoff and landing points with shorter waiting times); the physical meaning of linear mapping is that each minute of waiting time corresponds to a fixed operating cost (time cost, fuel cost, passenger time value); the longer the waiting time, the higher the total cost.

[0051] When the waiting time exceeds a preset threshold, nearby alternative landing points are searched and the timing capacity cost of the landing points is recalculated: If T_wait > T_max (the waiting time at the destination landing point exceeds the acceptable threshold, the landing point capacity is saturated, or the queue is too long), the alternative landing point search mechanism is triggered. Search range: All available landing points (status: open, equipment is normal, weather is suitable) are searched within a circular area centered on the destination with a radius d_max (configurable, such as 3 km, 5 km, 10 km, representing the maximum acceptable distance of the alternative landing point from the destination; if the distance is too far, the ground connection time is too long, and the alternative landing is meaningless). Filtering conditions: For each alternative landing point within the search range, its real-time status is obtained, and alternative landing points that meet the sufficient capacity condition are filtered: remaining capacity C_rem > 0 (there are available slots) and the estimated waiting time T_wait_backup calculated according to queuing theory ≤ T_max (the waiting time is within the acceptable range), while also satisfying the distance increase constraint d_backup ≤ d_max (the ground distance from the alternative landing point to the destination is within the acceptable range, which can be evaluated in conjunction with ground traffic time). If a suitable alternative take-off and landing point exists, the optimal alternative take-off and landing point is selected (based on a weighted score, considering factors such as shortest waiting time, closest distance to the destination, and matching refueling facilities). The destination is changed to this alternative take-off and landing point, and the capacity cost C_capacity=w4×T_wait_backup+w_dist×d_backup (adding the increased ground connection distance cost due to the change of take-off and landing point, where w_dist is the distance cost coefficient) is recalculated. The alternative landing point information is then marked in the output (for task scheduling and passenger notification). If no alternative take-off and landing points are available (all nearby take-off and landing points are saturated, waiting time exceeds the threshold, or are unavailable), a forced waiting strategy is executed: the temporal capacity cost of the take-off and landing point is set to positive infinity (C_capacity=+∞, indicating that the destination is currently unavailable and other options need to be considered in the route search), and a departure air circling waiting strategy suggestion is output (it is suggested that the flying car circle in the designated waiting airspace and wait for the take-off and landing point capacity to be released before landing; circling waiting needs to consider the remaining range and airspace constraints); the positive infinity cost makes the search algorithm either choose another time window (adjust the estimated arrival time to avoid peak hours) or reflect the unavailability of the take-off and landing point in the cost graph during planning. Through a complete mechanism of queuing theory waiting time calculation, linear time cost mapping, alternative take-off and landing point search, and forced waiting strategy, the real-time capacity status and queuing delay of take-off and landing points are transformed into cost terms for route search. This enables route planning to proactively avoid take-off and landing points with saturated capacity, select take-off and landing points with short waiting times, or reasonable alternative landing points. This ensures the temporal feasibility of routes under take-off and landing point capacity constraints and avoids in-flight circling or forced diversions caused by take-off and landing point capacity saturation in actual execution.

[0052] In some embodiments, the method of using a hybrid A* search algorithm to search for the optimal path on a multi-constraint cost graph includes: taking the takeoff point in the flight mission parameters as the search start point and the landing point in the flight mission parameters as the search end point, and initializing an open node set and a closed node set; selecting the node with the minimum comprehensive cost value from the open node set as the current expansion node, and moving the current expansion node from the open node set to the closed node set; generating neighboring nodes of the current expansion node according to the kinematic constraints of the flying car, calculating the actual cost and heuristic estimated cost of each neighboring node, and adding the actual cost and heuristic estimated cost to obtain the comprehensive cost value of the neighboring node; adding neighboring nodes not in the closed node set to the open node set, and repeatedly selecting the node with the minimum comprehensive cost value for expansion until it is expanded to the search end point; and backtracking from the search end point along the parent node pointer back to the search start point to generate an optimal path node sequence.

[0053] The search starts with the takeoff point from the flight mission parameters and ends with the landing point, initializing the open and closed node sets. During the initialization phase, key information is extracted from the flight mission parameters: the three-dimensional coordinates of the takeoff point (x_start, y_start, z_start, typically the takeoff position of a certain takeoff / landing point, with altitude equal to the takeoff / landing point elevation plus the safe takeoff altitude), the three-dimensional coordinates of the landing point (x_goal, y_goal, z_goal, the landing position of the destination takeoff / landing point), the flying car's performance parameters (maximum speed v_max, minimum speed v_min, maximum acceleration a_max, maximum deceleration a_min, maximum turning angular velocity ω_max, maximum climb rate climb_max, maximum descent rate descend_max, minimum turning radius R_min, these parameters determine the boundaries of kinematic constraints), and mission constraints (expected arrival time window, required transit points, no-fly zone avoidance priority, etc.). The initial state, starting from the takeoff point, is defined as follows: the state vector includes a 3D position (x, y, z), a velocity vector (vx, vy, vz) or velocity magnitude and heading angle (v, θ_heading, φ_climb), and a timestamp t. The position of the starting state is the takeoff point coordinates, the velocity is the initial climb speed or hovering speed (transitioning to cruise speed after takeoff in vertical takeoff and landing mode), the heading angle is the initial flight path direction, and the timestamp t=0. The landing point is defined as the search endpoint. The endpoint determination criteria are that the 3D distance between the search node's position and the landing point is less than a preset endpoint neighborhood threshold (e.g., 10 meters, 20 meters, to ensure endpoint accuracy) and the velocity state meets landing requirements (speed reduced to landing speed, altitude matching the landing point altitude). Initialize the open and closed node sets: The open node set uses a priority queue (min-heap) data structure, sorted by the comprehensive cost f of each node in ascending order, supporting fast retrieval of the node with the smallest f (O(logn)) and insertion / update of nodes (O(logn)). Initially, the starting state is added to the open node set, the actual cost g of the starting point is 0 (the cumulative cost from the starting point to the starting point is zero), the heuristically estimated cost h = heuristic function (starting state, ending point), and the comprehensive cost f = g + h = h. The closed node set uses a hash table or a three-dimensional grid marking array to store the states of nodes that have been expanded (discreteized states of position + velocity + heading) to avoid redundant expansion; it is initially empty. Simultaneously, the parent node pointer table (storing the parent node reference of each node for final path backtracking) and the node state table (storing the optimal g value of each state to determine if a newly arrived state is better) are initialized. The open node set, closed node set, parent pointer table, and state table together constitute the core data structure of the hybrid A* search.

[0054] The node with the lowest overall cost is selected from the open node set as the current expansion node, and then moved from the open node set to the closed node set. Each iteration of the main search loop first performs node selection: the node with the lowest overall cost *f* from the top of the open node set (priority queue) is selected as the current expansion node (this is the core of the A* algorithm: always expand the node that seems most likely to reach the destination; in *f* = *g* + *h*, *g* is the actual cost already incurred, and *h* is the estimated cost to reach the destination; the smallest *f* means the smallest estimated total cost). The current node is removed from the open node set (popped from the top of the heap), and its state (discrete identifier of position, velocity, and heading) is added to the closed node set (marked as expanded; if a neighboring state is already in the closed set during subsequent neighborhood generation, it is skipped to avoid duplicate expansion). After selection, a termination check is performed: it is checked whether the current node's position is within the destination's neighborhood (using the three-dimensional Euclidean distance *(x-x_goal)*). 2 +(y-y_goal) 2 +(z-z_goal)²) 0.5 If the destination threshold is met and the speed status satisfies the landing requirements, the search is successful, the main loop is exited, and the path backtracking phase begins. If the open node set is empty (all reachable nodes have been expanded but the destination has not been reached), the search fails, and a "no feasible path" message is output (the task parameters or constraints need to be adjusted and the search retried). Otherwise, the current node set is expanded. Moving the current node set into the closed set ensures that each state is expanded only once (when the same state is reached via other paths, if the g value of the new path is not better than the recorded g value, it is not expanded again; the optimality of A* guarantees that the g value when a state is first expanded is the optimal actual cost for that state), thus avoiding search space explosion. Efficient implementation of priority queues (such as binary heaps or Fibonacci heaps) ensures the real-time performance of selection operations for a large number of nodes (tens to hundreds of thousands of nodes in urban low-altitude scenarios).

[0055] The neighboring nodes of the current expanding node are generated according to the kinematic constraints of the flying car. The actual cost and heuristic estimated cost of each neighboring node are calculated, and the actual cost and heuristic estimated cost are added together to obtain the comprehensive cost of the neighboring node. The generation of neighboring nodes is the core difference between hybrid A* and traditional A*: traditional A* generates fixed neighbors (such as 26 neighbors) on the grid, without considering kinematic constraints, and the generated paths may be unexecutable; hybrid A* generates neighboring nodes based on the kinematic model of the flying car by simulating forward from the current state, so that the generated paths naturally satisfy the kinematic constraints. Kinematic model: A simplified dynamic model of the flying car (such as a three-degree-of-freedom point mass model or a six-degree-of-freedom model) is adopted. The state variables are position (x, y, z), velocity magnitude v, heading angle θ (horizontal direction angle), and climb angle φ (vertical direction angle). The control inputs are acceleration a (longitudinal direction, range [a_min, a_max]), heading angular velocity ω (horizontal turning, range [-ω_max, ω_max], ω_max=v / R_min is determined by the minimum turning radius and velocity), and climb angular velocity γ (vertical direction, range [-descend_max, climb_max]). Forward simulation: Starting from the current state, the control input space is discretized and sampled (e.g., acceleration is sampled in 3 levels: {deceleration, constant speed, acceleration}; yaw rate is sampled in 3 levels: {left turn, straight, right turn}; climb rate is sampled in 3 levels: {descent, level flight, climb}, for a total of 3×3×3=27 control input combinations; for higher accuracy requirements, the number of levels can be increased, such as 5 levels, but the computational load will increase). For each control input combination, kinematic integration (Euler method or Runge-Kutta method) is performed with a fixed time step Δt (e.g., 1 second, 0.5 seconds; the smaller the step, the smoother the path but the more nodes): position update x _new=x+v×cos(θ)×cos(φ)×Δt, y_new=y+v×sin(θ)×cos(φ)×Δt, z_new=z+v×sin(φ)×Δt; velocity update v_new=v+a×Δt (clamped to [v_min,v_max]); heading update θ_new=θ+ω×Δt; climb angle update φ_new=φ+γ×Δt (clamped to [-max_descend_angle,max_climb_angle]); time update t_new=t+Δt. Each combination of control inputs generates a new state (position, velocity, heading, climb angle, time) for one neighboring node, for a total of 27 neighboring nodes (or fewer than 27 after deduplication if some input combinations are clamped when the current velocity / altitude is close to the boundary).Perform a feasibility check on each neighboring node: Check if the comprehensive cost of the multi-constraint cost graph mesh corresponding to the location is positive infinity (inaccessible, such as inside a no-fly zone, inside a building, or at a completely unusable takeoff and landing point); if so, discard the neighboring node. Check if the speed, heading, and climb angle are within the flying car's performance limits (already clamped during forward simulation but need confirmation); check if the location is within the 3D dynamic map area (discard if outside the map boundary). Neighboring nodes that pass the checks proceed to cost calculation.

[0056] Actual cost calculation: The actual cost of a neighboring node, g_new = g_current + movement cost Δg. The movement cost Δg is the comprehensive cost incurred by the path segment from current to the neighboring node. It is calculated by performing a line integral along the line segment from current to the neighboring node (or multiple sampling points of the forward simulation trajectory) on the comprehensive cost map of the multi-constraint cost map. This is simplified to the average comprehensive cost of sampling several points (e.g., 5 to 10 points) at equal intervals along the line segment multiplied by the line segment length (distance), i.e., Δg = (ΣC_total(p_i) / n) × d_segment, where p_i is the sampling point on the line segment, n is the number of sampling points, and d_segment is the three-dimensional distance from current to the neighboring node. The movement cost includes both time cost (Δt multiplied by a time weight, reflecting flight time cost) and energy cost (energy consumption can be estimated based on speed, acceleration, and climb angle, and added to Δg). g_new is the cumulative actual cost from the starting point to the neighboring node. Heuristic cost estimation calculation: h_new = heuristic function(neighborhood node state, destination state). The heuristic function uses the 3D Euclidean distance multiplied by the base cost per unit distance as a lower bound estimate, i.e., h = d_3D × c_base, where d_3D = √((x_new - x_goal)² + (y_new - y_goal)² + (z_new - z_goal)²) is the 3D straight-line distance from the neighboring node to the destination, and c_base is the base cost per unit distance (e.g., 1 cost unit / meter, which can be determined based on the average grid cost). Price adjustment); to improve heuristic accuracy (reduce the number of search nodes), a height difference penalty (the distance weight in the z-direction is higher than that in the horizontal direction because vertical maneuvering consumes more energy), turn prediction cost (predict the turn cost based on the angle between the current heading and the destination direction), and no-fly zone detour prediction (predict the detour cost by quickly querying whether there is a no-fly zone blocking the destination direction); the heuristic function must satisfy admissibility (h does not exceed the minimum actual cost to the destination, i.e., h is a lower bound) to ensure the optimality of the A* algorithm (the path when first reaching the destination is the optimal path). The comprehensive cost f_new = g_new + h_new (actual cost paid + estimated future cost, total cost estimate).

[0057] Add neighboring nodes not in the closed node set to the open node set, and repeatedly select the node with the lowest comprehensive cost for expansion until the search endpoint is reached. For each neighboring node that passes the feasibility check, check whether its state (discrete identifier of position + velocity + heading, with position discretized by grid resolution, velocity by gear, and heading by angle step size to form a state key) is already in the closed node set: if it is already in the closed node set, it means that the state has been expanded by a better path (A* optimality guarantees that g is optimal during the first expansion), and skip the neighboring node; if it is not in the closed node set, further check whether it is already in the open node set: if it is not in the open node set, add the neighboring node to the open node set (insertion priority queue), record its parent node as current, and record its g, h, and f values; if it is already in the open node set, compare the newly calculated g_new with the g value already recorded for the node in the open node set. If g_new is smaller (the new path is better), update the g, f values ​​and parent node pointer of the node (perform a decrement key operation in the priority queue), and if g_new is not smaller, keep the original record unchanged. This update mechanism ensures that each state in the open set always retains the optimal (minimum g) path to that state. After processing all neighboring nodes of the current state, the process returns to the beginning of the main loop, selecting a new minimum f node from the open node set for expansion. This process iterates until the search reaches the endpoint (successful search) or the open set becomes empty (f-f). The number of nodes during the search process can be very large (hundreds of thousands of nodes in complex urban low-altitude scenarios). Efficient data structures (priority queues using binary heaps, state tables using hash tables or three-dimensional arrays) and state discretization (position by grid, speed and heading by gear) are used to control memory and computation. Typical urban route searches are completed within seconds to tens of seconds (meeting the real-time requirements of route planning).

[0058] The optimal path node sequence is generated by backtracking from the search endpoint to the search starting point along the parent node pointer. After a successful search (current is in the neighborhood of the endpoint and the landing conditions are met), starting from the endpoint node current, backtracking level by level along the parent node pointer: current→parent(current)→parent(parent(current))→...→starting point. At each node, the state (position, velocity, heading, climb angle, timestamp, cumulative cost g) of that node is added to the path list until the starting point (the node with an empty parent pointer) is reached. The node list obtained from the backtracking is then reversed in chronological order (from starting point to endpoint) to generate the optimal path node sequence. Each node in the sequence contains: three-dimensional spatial coordinates (x, y, z), velocity magnitude v and heading angle θ, climb angle φ, timestamp t, and cumulative actual cost g. The path node sequence is a discrete sequence of state points, and the adjacent nodes are kinematically feasible trajectory segments simulated in the forward motion (satisfying velocity, acceleration, turning, and climb constraints). Therefore, the entire path naturally satisfies the kinematic constraints of the flying car, and there are no unexecutable sharp turns or excessive gradients. The cumulative actual cost g_goal of the path is the total comprehensive cost from the starting point to the ending point, which includes the weighted sum of the no-fly zone safety cost, building wind farm energy consumption and stability cost, and take-off and landing point capacity and time cost. It is minimized among all feasible paths (guaranteed by A* optimality). The optimal path node sequence serves as the input control point sequence for subsequent B-spline smoothing and can also be directly used for route preview and feasibility analysis. Through a complete hybrid A* search process involving search initialization, priority queue node selection, kinematic constraint neighborhood generation and cost calculation, open set update and cyclic expansion, and endpoint backtracking, a kinematically feasible and comprehensively cost-optimal path node sequence was obtained on the multi-constraint cost graph, laying the foundation for subsequent smoothing and trajectory generation.

[0059] In some embodiments, the smoothing of the optimal path using B-spline curves includes: using nodes in the optimal path node sequence as control points of the B-spline curve, determining the order of the B-spline curve based on the number of control points; performing weighted interpolation on the control points according to the B-spline basis functions to generate a continuous and smooth B-spline curve; performing a curvature continuity check on the B-spline curve, and when the curvature of the B-spline curve exceeds the maximum turning curvature threshold of the flying car, adjusting the position of the corresponding control points and regenerating the B-spline curve until all curvatures of the B-spline curve satisfy the kinematic constraints of the flying car.

[0060] The nodes in the optimal path node sequence are used as control points for the B-spline curve, and the order of the B-spline curve is determined based on the number of control points. A B-spline curve is a piecewise polynomial parametric curve. Defined by control points and basis functions, it possesses excellent properties such as local controllability (moving a control point only affects a few adjacent curve segments, not the entire curve), high-order continuity (a p-th degree B-spline has a p-1 continuous derivative at the node, and a cubic B-spline is C2 continuous, meaning position, first derivative, and second derivative are all continuous), and convex hull property (the curve lies within the convex hull of the control points). These properties make it very suitable for smoothing discrete path points, generating continuous and smooth trajectories that flying cars can accurately track. Control point preparation: The 3D spatial coordinates (x_i, y_i, z_i) of each node in the optimal path node sequence obtained by hybrid A* search are used as the i-th control point P_i of the B-spline curve. The number of control points n+1 (i=0,1,...,n) is equal to the number of path nodes (city routes typically have tens to hundreds of nodes, so the number of control points is moderate). The control points are arranged in order from the start to the end of the path, with the first and last control points being the takeoff and landing points, respectively. If the number of path nodes is too large (e.g., more than 100), the Douglas-Puk algorithm can be used to thin out the path nodes (reducing the number of control points while keeping the path shape error less than the threshold) to avoid the increased computational load of B-spline and curve overfitting (the curve excessively following the jitter of discrete points). If the number of path nodes is too small (e.g., less than 4), linear interpolation can be used to add control points between nodes to ensure that the curve has sufficient degrees of freedom for smoothing. Determining the order: The order p (polynomial degree) of the B-spline curve is determined based on the number of control points and smoothing requirements. Cubic B-splines (p=3, i.e., 4th order polynomial, C2 continuous) are typically used because the continuity of the second derivative in a cubic B-spline implies continuity of curvature (curvature is related to the second derivative), which is crucial for trajectory tracking of flying cars: continuous curvature means that lateral acceleration changes continuously without abrupt changes during turns, allowing the flying car to turn smoothly without attitude abrupt changes. Quadratic B-splines (p=2, C1 continuous) have discontinuous curvature (the first derivative is continuous, but the second derivative may abruptly change at nodes), suitable for scenarios with lower requirements. Higher orders (p=4 or 5, C3 or C4 continuous) offer higher smoothness but increase computational complexity and reduce local controllability, and are generally not needed. When the number of control points n+1≤p+1 (insufficient control points to support a p-order B-spline), the order p needs to be reduced to n (a p-order B-spline requires at least p+1 control points), or the number of control points needs to be increased. The number of control points for typical urban flight routes is much greater than 4, therefore cubic B-splines are always feasible. Node vector settings: Use clamped node vectors to make the B-spline curve pass through the first and last control points (the start and end points are precisely matched with the takeoff and landing points). The clamped node vectors have p+1 repeated nodes at the beginning and end, and the intermediate nodes are evenly distributed (e.g., uniformly selected values ​​between 0 and 1) to ensure that the curve is interpolated to the control points at the endpoints.

[0061] Perform weighted interpolation on the control points according to the B-spline basis functions to generate a continuous and smooth B-spline curve. The parametric equation of the B-spline curve is: P(t)=Σ_{i=0}^{n}N_{i,p}(t)×P_i, where t is the curve parameter (the value range is [t_p,t_{n+1}], usually [0,1] under the clamped knot vector), P_i is the three-dimensional coordinate (vector) of the i-th control point, and N_{i,p}(t) is the i-th p-degree B-spline basis function (a scalar weight that determines the influence degree of the control point P_i on the curve point at the parameter t). The B-spline basis functions are calculated through the De Boor-Cox recurrence formula: for the p=0 degree basis function, N_{i,0}(t)=1 when t_i≤t<t_{i+1}, otherwise 0 (step function); for the p>0 degree basis function, N_{i,p}(t)=(t-t_i) / (t_{i+p}-t_i)×N_{i,p-1}(t)+(t_{i+p+1}-t) / (t_{i+p+1}-t_{i+1})×N_{i+1,p-1}(t) (recurrence formula, the term is zero when the denominator is zero), where t_i is the i-th knot value of the knot vector. Characteristics of the basis functions: for any t, the sum of all non-zero basis functions is 1 (partition of unity, which ensures that the curve is an affine combination of the control points); each basis function is non-zero only in the local parameter interval (local support property, moving the control point P_i only affects the curve segment in the interval t∈[t_i,t_{i+p+1}], that is, local controllability). Curve generation: perform uniform sampling on the parameter t within the range of [0,1] according to the set parameter interval Δt (for example, 0.001 to 0.01, the smaller the interval, the denser the curve sampling points, the finer the curve, and the larger the calculation amount). For each sampling parameter t_j=j×Δt (j=0,1,...,1 / Δt), calculate each basis function N_{i,p}(t_j) through the De Boor recurrence formula (or directly calculate all non-zero basis function values), and perform weighted summation according to the parametric equation to obtain the three-dimensional coordinates (x_j,y_j,z_j) of the point P(t_j) on the curve; all sampling points are arranged in the order of t to form a discrete representation (point sequence) of the continuous and smooth B-spline curve, and the density of the point sequence is sufficient to accurately describe the curve shape (the distance between adjacent points is much smaller than the positioning accuracy of the flying car). The first derivative of the curve P'(t) (tangent direction, corresponding to the velocity direction) and the second derivative P''(t) (corresponding to curvature and normal acceleration) can be calculated by the derivative recurrence of the basis functions or approximated by the finite difference of the sampling points, for subsequent curvature check. The generated B-spline curve is a C2 continuous three-dimensional space curve, which is smoother and more continuous than the original discrete path node sequence (polygonal line, with sudden direction change at the nodes), has no sudden direction change between adjacent segments, and is suitable for accurate tracking by flying cars.

[0062] A curvature continuity check is performed on the B-spline curve. When the curvature of the B-spline curve exceeds the maximum turning curvature threshold of the flying car, the positions of the corresponding control points are adjusted and the B-spline curve is regenerated until all curvatures of the B-spline curve satisfy the kinematic constraints of the flying car. Curvature is a geometric quantity describing the degree of curvature of a curve. The curvature formula for a three-dimensional curve is κ(t)=|P'(t)×P''(t)| / |P'(t)|³, where P'(t) is the first derivative (tangent vector), P''(t) is the second derivative, × is the vector cross product, and | | represents the vector magnitude; the unit of curvature κ is per meter (1 / meter). The larger the curvature, the sharper the curve turn. The radius of curvature R = 1 / κ (meters). The maximum turning curvature threshold κ_max of a flying car is determined by its minimum turning radius R_min: κ_max = 1 / R_min. R_min is the minimum turning radius of the flying car at the current speed (determined by the maximum lateral acceleration a_lat_max and the speed v: R_min = v² / a_lat_max. The higher the speed, the larger the minimum turning radius and the smaller the maximum curvature. Therefore, κ_max is a function of speed, with a smaller κ_max at high speeds and a larger κ_max at low speeds). A typical eVTOL flying car has a minimum turning radius of approximately 50 to 200 meters at cruising speeds (e.g., 100 to 200 km / h), corresponding to a κ_max of approximately 0.005 to 0.02 meters per meter. Curvature continuity check: For each sampling point P(t_j) along the B-spline curve, calculate the curvature κ(t_j) (using the cross product formula of the first and second derivatives, or approximating it using the circumradius of three adjacent sampling points). Compare κ(t_j) with the maximum curvature threshold κ_max(v_j) at the corresponding velocity at that point (v_j is the velocity of the curve at that point, which can be obtained by interpolation from the velocity information of the path nodes). If the curvature of all sampling points satisfies κ(t_j)≤κ_max(v_j), the curve passes the curvature check and is smoothed. If the curvature of a certain segment of the curve exceeds the threshold (κ(t_j)>κ_max), it means that the curve turns too sharply, exceeding the turning capability of the flying car. The flying car cannot accurately track this segment of the curve (forced tracking will lead to excessive lateral acceleration, sudden attitude changes, or even instability), and curvature correction is required.

[0063] Curvature Correction: Locate the parameter interval [t_a, t_b] corresponding to the curve segment with excessive curvature. Based on the local support of the B-spline, determine the set of control points affecting this curve segment (control points where the basis function is non-zero in [t_a, t_b]). Adjust the positions of these control points: move control points on the inside of the turn to the outside of the turn (increase the turning radius and reduce curvature), or increase the density of control points at sharp turns (insert new control points to make the curve smoother in that segment), or adjust the height of the control points (slow down the height change of sharp turns to reduce the three-dimensional curvature). The adjustment amount is adaptively determined according to the degree of curvature exceeding the limit (the magnitude of κ-κ_max), with larger adjustments for greater exceedance. After adjusting the control points, regenerate the B-spline curve (recalculate the basis function and sampling points), and perform curvature continuity checks again. Iterate in this way (control point position → regenerate curve → curvature check → readjust) until the curvature of all sampling points of the B-spline curve satisfies κ(t_j) ≤ κ_max(v_j), that is, the curvature of the entire curve is within the turning capability range of the flying car. The iterative process typically converges quickly (it only requires a few adjustments, as the original path generated by the hybrid A* already satisfies kinematic constraints, and curvature exceedances usually only occur at a few sharp turns). Curvature correction ensures that the smoothed B-spline curve is not only continuous and smooth but also accurately trackable by the flying car at all points (curvature is within the turning capability range), solving the problem of curvature abrupt changes in the flight path and the flying car's inability to track it in traditional simple polygonal line or circular arc transition methods. After smoothing, the sampling point sequence of the B-spline curve (including three-dimensional coordinates, tangent direction i.e., heading angle, curvature, and corresponding velocity) serves as the reference trajectory for subsequent local obstacle avoidance and spatiotemporal trajectory generation using the dynamic window method. Through a complete B-spline smoothing process involving control point preparation and order determination, curve generation via de Boer recursive basis function weighted interpolation, curvature continuity checking, and iterative correction, the discrete kinematically feasible path generated by the hybrid A* is transformed into a continuous, C2-smooth, three-dimensional spatial curve with curvature satisfying the flying car's turning constraints, providing the flying car with an accurately trackable reference trajectory.

[0064] In some embodiments, the step of combining the dynamic window method to perform local obstacle avoidance and generate a spatiotemporal trajectory point sequence includes: sampling along a smoothed B-spline curve at set time intervals to generate an initial spatiotemporal trajectory point sequence, where each spatiotemporal trajectory point contains spatial coordinates, a timestamp, and velocity information; during the execution of the flying car's trajectory, determining a dynamic velocity window based on the current velocity and acceleration constraints, and traversing all candidate velocity combinations within the dynamic velocity window; performing a forward simulation on each candidate velocity combination to predict the flying car's trajectory within a set prediction time domain, calculating the minimum distance between each predicted trajectory and the obstacle and the degree of deviation from the B-spline curve; selecting the candidate velocity combination with the maximum distance from the obstacle and the minimum degree of deviation from the B-spline curve as the current execution velocity, updating the spatiotemporal trajectory point sequence based on the current execution velocity; and using the updated spatiotemporal trajectory point sequence as the final output.

[0065] The B-spline curve is sampled at set time intervals to generate an initial spatiotemporal trajectory point sequence. Each spatiotemporal trajectory point contains spatial coordinates, a timestamp, and velocity information. The B-spline curve is a purely geometric curve (describing only spatial position and shape) and does not contain temporal dimension information. Therefore, it needs to be parameterized in time and sampled to generate the spatiotemporal trajectory point sequence (i.e., the velocity planning / time allocation step in trajectory planning), so that each spatial point corresponds to an arrival time and velocity, forming a four-dimensional spatiotemporal trajectory (x, y, z, t). Time interval setting: The B-spline curve is sampled in time parameterization according to a set time interval Δt_s (e.g., 0.5 seconds, 1 second; the smaller the interval, the denser the trajectory points, the higher the time resolution, and the more accurate the flying car tracking, but the larger the data volume; typical eVTOL routes use an interval of 0.5 to 1 second). Speed ​​planning: Based on the flying car's performance envelope (maximum speed v_max, minimum speed v_min, maximum acceleration a_max, maximum deceleration a_min) and mission time constraints (expected arrival time window), speed planning is performed along a B-spline curve to determine the desired speed at each position on the curve: During takeoff and climb, the speed gradually accelerates from hovering speed to cruise speed; during horizontal cruise, the speed is maintained (which can be fine-tuned according to wind field costs, with appropriate acceleration in tailwind and appropriate deceleration in headwind to optimize energy consumption); during descent and landing, the speed gradually decelerates to landing speed; and during turns, the speed automatically decelerates according to curvature (the greater the curvature, the lower the speed, to ensure that the lateral acceleration v² / R does not exceed the limit, i.e., v≤√(a_lat_max×R)). Speed ​​planning must satisfy that the speed change between adjacent points does not exceed the acceleration constraint (|v_{j+1}-v_j| / Δt_s≤a_max), ensuring that the speed curve is smooth and executable. Timestamp Calculation: Starting from the initial point t=0, the timestamp t_j for each sampling point is t_{j-1}+Δt_s (equal time interval sampling), or a non-uniform time interval is calculated based on the arc length and velocity (t_j=t_{j-1}+Δs_j / v_j, where Δs_j is the arc length of adjacent points and v_j is the average velocity of that segment, making the time allocation more consistent with the actual flight speed); the endpoint time t_end should fall within the expected arrival time window of the mission. If it exceeds this window, the cruise speed is adjusted (accelerated or decelerated) to ensure that the total flight time meets the time window constraint. Initial Spatiotemporal Trajectory Point Sequence Generation: Each trajectory point contains spatial coordinates (x_j, y_j, z_j) (B-spline curve sampling point), timestamp t_j, velocity information (velocity magnitude v_j and velocity direction, i.e., curve tangent direction / heading angle θ_j, climb angle φ_j), curvature κ_j, and cumulative flight distance s_j; the trajectory points are arranged in chronological order to form a complete four-dimensional spatiotemporal reference trajectory from takeoff to landing. The initial spatiotemporal trajectory point sequence serves as the reference execution trajectory for the flying car, and also as the tracking target for the local obstacle avoidance method of the dynamic window method (which avoids local dynamic obstacles while tracking the reference trajectory).

[0066] During the flight car's trajectory execution, a dynamic speed window is determined based on the current speed and acceleration constraints. All candidate speed combinations are traversed within this window. The Dynamic Window Approach (DWA) is a local reactive obstacle avoidance algorithm based on velocity space. Its core idea is to search for the optimal speed command within the achievable velocity space (dynamic window) under the current state, enabling the flight car to avoid obstacles while tracking the reference trajectory. It features low computational cost, high real-time performance, and suitability for dynamic environments. The DWA operates online at a high frequency (e.g., 10 Hz to 20 Hz, i.e., executed every 0.05 to 0.1 seconds), continuously sensing the local environment and adjusting speed commands as the flight car executes the reference trajectory. Current state acquisition: The current state of the flying car is acquired in each control cycle, including the current position (x_c, y_c, z_c) (from the navigation system such as GNSS+IMU integrated navigation), the current velocity vector (vx_c, vy_c, vz_c) or the velocity magnitude v_c and heading angle θ_c, the current climb angle φ_c, and the current acceleration (from the inertial measurement unit). Dynamic speed window determination: Based on the current speed and acceleration constraints, determine the achievable speed range within the next control cycle: the achievable speed range is [v_c-a_min×Δt_c,v_c+a_max×Δt_c] (a_min is the maximum deceleration, a_max is the maximum acceleration, and Δt_c is the control cycle, such as 0.1 seconds), while simultaneously clamping to the flying car's speed envelope [v_min,v_max] (hovering speed to maximum cruise speed); the achievable yaw rate (turning rate) range is [-ω_max,ω_max] (ω_max is the maximum turning rate, determined by the maximum lateral acceleration and the current speed; ω_max=a_lat_max / v_c, the higher the speed, the smaller the maximum turning rate); the achievable climb rate (vertical rate) range is [-descend_max,climb_max] (maximum descent rate to maximum climb rate). The dynamic speed window is a three-dimensional combination space of speed magnitude v, heading angular velocity ω, and climb angular velocity γ (or simplified to a two-dimensional space of v and ω, with the vertical direction controlled separately). Candidate speed combination traversal: The dynamic speed window is discretized according to a set resolution: speed magnitude step Δv (e.g., 0.5 m / s), heading angular velocity step Δω (e.g., 1 degree / s or 5 degrees / s), and climb angular velocity step Δγ (e.g., 0.5 m / s). All discrete (v, ω, γ) combinations within the window are traversed to form a candidate speed command set (typically tens to hundreds of combinations, such as v at level 5, ω at level 7, and γ at level 3, totaling 105 candidates). The purpose of this traversal is to search for the optimal speed command among all achievable speed commands. The design of the dynamic window ensures that only physically achievable speed commands under the current state are considered (unachievable speeds, such as instantaneous changes from low to high speed, are not included in the window), allowing the search result to be executed directly.

[0067] For each candidate speed combination, a forward simulation is performed to predict the trajectory of the flying car within a defined prediction time domain. The minimum distance between each predicted trajectory and obstacles, and the degree of deviation from the B-spline curve, are calculated. For each candidate speed combination (v_cmd, ω_cmd, γ_cmd), a forward simulation (trajectory extrapolation) is performed to predict the trajectory of the flying car over a future period under that speed command. Forward simulation: Using the current state of the flying car as the initial condition, and the candidate speed command (v_cmd, ω_cmd, γ_cmd) as a constant control input (assuming that the speed command is maintained in the prediction time domain, or updated in each simulation step), the flying car kinematic model is used to perform multi-step forward integration. The prediction time domain T_pred (e.g., 3 to 5 seconds; the longer the prediction time domain, the longer the trajectory planning, but the greater the computational load and the lower the prediction accuracy for dynamic environments. Typically, 3 to 5 seconds corresponds to 30 to 50 simulation steps, and the simulation step Δt_sim is 0.1 seconds). The position, velocity, heading, and climb angle are updated in each step of integration, resulting in a series of position points (e.g., 50 points) for the predicted trajectory, forming the predicted motion trajectory (a three-dimensional polyline or curve) corresponding to the candidate speed command. Obstacle information: Local obstacles include dynamic obstacles (other flying cars, drones, birds, etc., target detection and tracking from airborne perception systems such as radar, lidar, and vision, or traffic information broadcasts from air traffic control systems) and local details of static obstacles (building protrusions, temporary obstacles, sudden no-fly zone boundaries, from airborne perception or updated 3D maps); each obstacle has a location, size, speed of movement (dynamic obstacles), and safety radius (the minimum distance that flying cars must maintain, such as 30 to 50 meters, set according to the size of the flying car and positioning accuracy). Minimum distance calculation between predicted trajectory and obstacles: For each predicted trajectory, calculate the shortest distance from each point on the trajectory to each obstacle (for dynamic obstacles, the time dimension needs to be considered: the timestamp of the predicted trajectory point and the predicted position of the obstacle at that moment, calculate the spatiotemporal distance rather than the simple spatial distance, and avoid encountering moving obstacles in time). Take the minimum distance d_min among all points and all obstacles as the safe distance index of the predicted trajectory; if d_min is less than the safe radius (the predicted trajectory will collide with the obstacle or the distance is too close), the candidate speed command is directly marked as infeasible (discarded) and will not participate in the subsequent evaluation.

[0068] Calculation of deviation from B-spline curve: For each predicted trajectory, calculate its deviation from the reference B-spline curve (the geometric curve corresponding to the initial spatiotemporal trajectory point sequence). The deviation can be expressed as the shortest vertical distance from the predicted trajectory endpoint to the reference curve (lateral deviation), or the average vertical distance from all points on the predicted trajectory to the reference curve, or the time deviation between the predicted trajectory and the reference trajectory (the time difference between the predicted arrival at the corresponding reference point). Typically, the endpoint lateral deviation d_track = min|p_end - P(t)| (the shortest distance from the predicted trajectory endpoint to the B-spline curve) is used as the tracking accuracy indicator. The smaller d_track is, the better the speed command can track the reference route. At the same time, other auxiliary evaluation indicators are calculated: speed smoothness (the change between the candidate speed and the current speed |v_cmd - v_c| and the change in heading angular velocity; the smaller the change, the smoother the motion and the higher the passenger comfort), progress toward the target (the decrease in distance from the predicted trajectory endpoint to the landing point, or the distance traveled along the reference route; the larger the distance, the faster the progress toward the target), and energy consumption estimation (estimate the energy consumption of the predicted trajectory based on speed, acceleration, and climb angle; the lower the energy consumption, the more economical the trajectory). Evaluation function: The comprehensive evaluation score of the candidate speed command is obtained by weighted summation of the above indicators: Score = w_obs × d_min + w_track × (1 / d_track) + w_smooth × (1 / Δv) + w_progress × Δs_progress - w_energy × E_pred, where w_obs, w_track, w_smooth, w_progress, and w_energy are the weights of each indicator (configured according to task priority; w_obs is larger when safety is prioritized, and w_track is larger when tracking is prioritized). The higher the score, the better the speed command; the score of infeasible candidates (d_min < safety radius) is set to negative infinity. Through a complete process of forward simulation to predict trajectory, calculate spatiotemporal obstacle distance, calculate reference curve deviation, and perform multi-indicator weighted evaluation, each candidate speed command is quantitatively evaluated, providing a basis for optimal speed selection.

[0069] The candidate velocity combination that maximizes the distance from the obstacle and minimizes the deviation from the B-spline curve is selected as the current execution velocity. The spatiotemporal trajectory point sequence is updated based on the current execution velocity, and the updated spatiotemporal trajectory point sequence is used as the final output. Optimal velocity selection: Among all feasible candidate velocity commands (d_min≥safe radius), the candidate velocity command with the highest comprehensive evaluation score is selected as the execution velocity (v_exec, ω_exec, γ_exec) for the current control cycle. The highest score means that the velocity command has the highest tracking accuracy (minimum deviation) of the reference B-spline curve while ensuring safety (maximum distance from the obstacle or at least meeting the safe distance), and also has smooth speed, advances towards the target, and is energy-efficient. When safety and tracking conflict (e.g., strictly tracking the reference trajectory would lead to approaching the obstacle), the weight configuration ensures safety priority (when w_obs is large enough, a velocity that deviates from the reference trajectory but is far away from the obstacle is selected, i.e., local detour). Speed ​​Execution Application: The selected speed execution commands (v_exec, ω_exec, γ_exec) are sent to the flying car's flight control system. The flight control system controls the flying car to move at the specified speed within the current control cycle (adjusting thrust, control surfaces, and attitude to achieve the target speed, heading angular velocity, and climb angular velocity). Spatiotemporal Trajectory Sequence Update: After the flying car actually executes the speed command, its actual motion state (actual position, actual speed, actual heading, timestamp) is recorded. This records replace or supplement the reference points at the corresponding moments in the initial spatiotemporal trajectory point sequence, forming the actual executed spatiotemporal trajectory point sequence. Deviations between the actual trajectory and the reference trajectory (deviations caused by navigation errors, wind disturbances, and local obstacle avoidance) are continuously recorded. The dynamic window method recalculates the speed window and candidate speeds based on the actual state (not the reference state) in the next control cycle, forming a closed-loop tracking and obstacle avoidance system. The updated spatiotemporal trajectory point sequence contains the actual execution state (x, y, z, t, v, θ, φ) at each moment, which is a complete spatiotemporal record of the flying car's actual flight. It also serves as reference data for the route execution report after the mission is completed and for subsequent route optimization.

[0070] The dynamic window method operates at a high frequency of 10 to 20 Hz. Each control cycle acquires the current state, senses local obstacles, determines the dynamic window, traverses candidate speeds, performs forward simulation evaluation, selects the optimal speed, executes, and updates the trajectory. This cycle continues until the flying car reaches its destination and lands. The high-frequency operation ensures rapid response to dynamic environments (sudden appearance of other flying cars, sudden no-fly zones, and abrupt wind changes) (reaction time 0.05 to 0.1 seconds), achieving local dynamic obstacle avoidance and precise tracking under global reference trajectory constraints. Final output: The updated spatiotemporal trajectory point sequence (the actual executed four-dimensional spatiotemporal trajectory, including complete information such as spatial coordinates, timestamps, speed, heading, climb angle, and curvature) is used as the final output for execution by the flying car's flight control system, monitoring by the mission management system, and archiving by the post-event analysis system. Simultaneously, local obstacle avoidance event records (when, where, and why obstacles were bypassed, and the deviation in the bypass trajectory) are output for mission reporting and flight path optimization. By implementing a complete dynamic window method process, including initial spatiotemporal trajectory generation, dynamic velocity window determination, candidate velocity traversal and forward simulation evaluation, optimal safe tracking speed selection, execution and trajectory update, real-time avoidance of local dynamic obstacles and accurate tracking of the reference trajectory are achieved based on B-spline smoothed reference trajectory. This generates a spatiotemporal trajectory point sequence that is both globally optimal and locally safe and feasible, and can be accurately executed by flying cars. This solves the problems of traditional methods lacking local dynamic obstacle avoidance capabilities and being unable to cope with dynamic environmental changes.

[0071] In some embodiments, when a sudden no-fly zone or take-off / landing point capacity is detected to be lower than a preset threshold, reloading data and performing route replanning includes: receiving airspace control notices and take-off / landing point status updates in real time during route execution, detecting whether a sudden no-fly zone or take-off / landing point remaining capacity is lower than a preset capacity threshold; when a sudden no-fly zone or take-off / landing point remaining capacity is detected to be lower than the preset capacity threshold, freezing the current route execution, recording the current position, current speed, and remaining battery status of the flying car; reloading the updated static no-fly zone data, 3D building model data, real-time meteorological wind field data, and real-time take-off / landing point status data, using the current position of the flying car as the new search starting point, and the original landing point or alternative take-off / landing point as the new search endpoint, and re-executing the route planning; and sending the re-planned spatiotemporal trajectory point sequence to the flying car execution to restore the route execution status.

[0072] During flight path execution, the system receives real-time airspace control notices and landing point status updates, detecting sudden no-fly zones or landing point capacity falling below preset thresholds. Online monitoring is a prerequisite for dynamic replanning; the system continuously monitors environmental changes throughout flight path execution to promptly identify events affecting the feasibility of the current route. Airspace control notice reception: The system receives airspace control notices in real-time via aviation information service interfaces (such as flight information service systems, air traffic control data link communication, and push notifications from low-altitude traffic management platforms). Notices include newly issued temporary no-fly zones (activity control, emergency rescue airspace, VIP security airspace, construction and hoisting airspace, drone flight control areas, etc.), changes to the scope or time of existing no-fly zones, and changes to airspace usage restrictions (altitude restrictions, speed restrictions, access requirements). Notice data includes spatial boundaries (polygon vertex coordinates), altitude range, effective time (immediate or timed), expiration time, and control type and level. Landing Point Status Update Reception: The system continuously receives status updates from each landing point via real-time data interfaces of the landing point management system (such as MQTT message push, RESTAPI polling, and data link communication). Updates include changes in remaining capacity (slot occupancy or release), queue length changes (new aircraft joining the queue or aircraft completing takeoff / landing and departing), changes in landing point operational status (open / closed, equipment failure / maintenance, changes in weather airworthiness such as closure due to strong winds or thunderstorms), and changes in refueling facility status. Monitoring Frequency: Airspace control notices are received using an event-driven approach (new notices are received and processed immediately). Landing point status is detected through timed polling (e.g., polling every 10 seconds to 1 minute) or event-driven push, ensuring that the detection delay for status changes is within seconds to minutes.

[0073] Trigger condition detection: Parse and perform spatial intersection tests on the received airspace control notices: Calculate whether there is a spatiotemporal intersection between the newly issued or changed no-fly zone and the currently executed flight path (the future part of the spatiotemporal trajectory point sequence, i.e., the unflyed segment from the current time to the destination) (the spatial range of the no-fly zone intersects with the spatial range of the future segment of the flight path, and the effective time of the no-fly zone overlaps with the time when the flying car is expected to arrive at that location). If there is an intersection, it is determined that a sudden no-fly zone has occurred, affecting the current flight path (the flying car will enter the no-fly zone if it flies along the original flight path, violating airspace control); Perform capacity threshold judgment on the take-off and landing point status update: If the remaining capacity of the destination take-off and landing point (or the planned alternate landing point of the flight path) drops to 0 (no available slots) or the queuing waiting time exceeds the preset threshold T_max (capacity saturation, waiting time is unacceptable), it is determined that the take-off and landing point capacity is lower than the preset capacity threshold (the landing plan of the original flight path is not feasible). Simultaneously, other triggering conditions (expandable) are detected: real-time wind field changes (such as forecasted strong winds or turbulence exceeding the flying car's airworthiness range), airborne system failures (requiring nearby emergency landing), and mission changes (passengers requesting a change of destination), all of which can trigger replanning. When any triggering condition is met, the route replanning process is immediately initiated; if no triggering condition is detected, the original route continues to be executed, with continuous monitoring. Through an online monitoring mechanism that includes real-time multi-source data reception, spatiotemporal intersection testing, and capacity threshold judgment, environmental changes affecting route feasibility are ensured to be detected in a timely manner, providing trigger signals for dynamic replanning.

[0074] When a sudden no-fly zone or the remaining capacity of a take-off and landing point is detected to be lower than a preset capacity threshold, the current flight path execution is frozen, and the current position, current speed, and remaining battery status of the flying car are recorded. After triggering replanning, a safety freeze operation is first performed to ensure that the flying car remains in a safe state during the replanning period. Flight path execution freeze: A hover or level flight hold command is sent to the flying car's flight control system. The flying car suspends its forward movement along the original flight path and maintains its current position while hovering (hovering is possible in eVTOL vertical take-off and landing mode) or waiting in level flight within a safe airspace (if hovering is not possible in fixed-wing mode, it will circle within a designated waiting airspace). During the freeze, the flying car maintains stable altitude, speed, and heading, does not perform large maneuvers, and waits for a new flight path to be generated. The freeze operation must ensure that the current position of the flying car is safe (not in a no-fly zone, not in conflict with other aircraft, and not hovering above densely built-up areas). If the current position is not safe, it will first maneuver to the nearest safe waiting point before freezing. State Recording: Accurately records the complete current state information of the flying car, serving as the initial conditions for replanning: current position (x_now, y_now, z_now) (3D coordinates, from integrated navigation, meter-level accuracy), current speed (v_now, θ_now, φ_now) (speed magnitude, heading angle, climb angle, from atmospheric data computer and inertial navigation), current acceleration and attitude (roll angle, pitch angle, yaw angle, from inertial measurement unit), remaining battery status (battery remaining charge percentage SOC, remaining flight time or remaining range, estimated by the battery management system, considering current power consumption and wind conditions), current timestamp t_now, distance flown and flight time flown, and mission completion progress (e.g., percentage of completed range). This state information ensures that replanning starts from the flying car's true current state (rather than from a reference point on the original route), and the generated new route smoothly connects with the current state, avoiding trajectory jumps. It also records the reasons for triggering replanning and event details (the number and scope of sudden no-fly zones, specific data on changes in takeoff and landing point capacity) for updating the environmental model during replanning and for post-event analysis. The freeze and status recording operations are executed immediately upon triggering (with a delay of less than 1 second) to ensure that the flying car quickly enters a safe waiting state and avoids continuing to fly along the original route and entering dangerous areas after environmental changes.

[0075] The updated static no-fly zone data, 3D building model data, real-time weather and wind field data, and real-time take-off and landing point status data are reloaded. The current location of the flying car is used as the new search starting point, and the original landing point or alternative take-off and landing point is used as the new search endpoint. The route planning is then re-executed. Data Reloading and Environmental Model Updates: Based on triggering events and the latest data, update the 3D dynamic map and multi-constraint cost map. For airspace, add sudden no-fly zones to the no-fly zone database (spatial boundaries, altitude range, effective time), update the no-fly zone attributes of affected grids (grids within sudden no-fly zones have their costs set to positive infinity and impassable; grids within buffer zones have their costs updated according to the exponential potential field), and refresh the effective status of all no-fly zones (removing expired temporary no-fly zones). For take-off and landing points, update the latest remaining capacity, queue length, average service time, and operational status of each take-off and landing point, and recalculate the temporal capacity cost of take-off and landing points (grids related to saturated take-off and landing points have their costs set to positive infinity or high cost; waiting time costs can be updated using take-off and landing points). For wind fields, load the latest real-time meteorological wind field data (if the wind field has been updated), and refresh the grid wind field attributes and building wind field costs. For building models, if there are changes such as building construction or demolition, update the building model and related wind field costs (building models typically change infrequently, and existing models can be reused during replanning). Data loading adopts an incremental update method (updating only the affected grid areas and data items) to avoid the performance overhead of full map recalculation and ensure the speed of replanning.

[0076] New search start and end point settings: The search start point is set to the current position (x_now, y_now, z_now) and current state (speed v_now, heading θ_now, climb angle φ_now, time t_now) of the flying car, that is, the search starts from the actual state of the flying car, rather than from the original route start point (already flown segments do not need to be replanned); the search end point is determined according to the triggering reason and mission requirements: if the triggering reason is a sudden no-fly zone but the original landing point is still available (sufficient capacity, open status), then the end point is still the original landing point (the destination remains the same, only the unflown segments need to be replanned to bypass the sudden no-fly zone); if the triggering reason is that the take-off and landing point capacity is saturated (the original landing point is unavailable), then the end point is changed to an alternative take-off and landing point (searching for the optimal alternative take-off and landing point with sufficient capacity and acceptable waiting time nearby, the selection criteria are the same as the aforementioned alternative take-off and landing point search), and at the same time the mission scheduling system notifies passengers of the change in ground connection arrangements; if the original landing point and nearby alternative take-off and landing points are both unavailable, then the end point is changed to the nearest available emergency landing point (such as open ground, helipad, road emergency landing zone), prioritizing flight safety. Incremental replanning strategy: To accelerate the replanning process, incremental replanning is adopted instead of a full search from scratch. This involves reusing already flown segments (which do not require planning) and unaffected future segments from the original route that were not affected by the triggered event. (If a future segment of the original route does not intersect with the sudden no-fly zone and the take-off and landing points remain unchanged, that segment can be directly reused, and only the affected segment, i.e., the detour from the current position to the original route recovery point, needs to be replanned.) The search space for incremental replanning is much smaller than that for full replanning, and the computation time is significantly reduced (from tens of seconds to several seconds). If the impact of the triggered event is large (such as a large no-fly zone covering most of the future segments of the original route) or the take-off and landing points change, the process degenerates into full replanning (a complete search from the current position to the new destination). Re-execute route planning: Starting from the current state and targeting the new destination, re-execute the complete route planning process on the updated multi-constraint cost map: hybrid A* search for the optimal path, B-spline curve curvature constraint smoothing, dynamic window method for local obstacle avoidance parameter configuration, and generation of spatiotemporal trajectory point sequences; the planning of the new route must consider the remaining power constraints (estimate the maximum remaining range based on the remaining power, the total range of the new route must be less than the remaining range, and a safety margin such as 20% must be retained), time constraints (if the mission has a latest arrival time requirement, the new route must meet it), and flying car performance constraints (consistent with the initial planning). After the new route is generated, perform feasibility verification: check that the new route does not intersect with any no-fly zones (including sudden no-fly zones), the take-off and landing point capacity meets the requirements, the total range is within the remaining power range, the curvature and speed meet the performance constraints, and there is no conflict with other aircraft routes (if traffic information is integrated). If the verification passes, the new route is usable; otherwise, adjust the parameters and re-plan or expand the search range.Through a complete process of incremental data updates, setting the current starting point and new endpoint, incremental / full replanning, and feasibility verification, a new feasible route from the current state of the flying car was quickly generated after environmental changes, ensuring the continuity of the mission.

[0077] The newly generated spatiotemporal trajectory point sequence will be sent to the flying car for execution, restoring the flight path execution status. New flight path issuance: The newly generated spatiotemporal trajectory point sequence (containing the complete four-dimensional trajectory from the current position to the new destination, including the spatial coordinates, timestamp, speed, heading, climb angle, and curvature of each point) will be sent to the flying car's flight control system and mission management system via data link communication (such as low-altitude traffic management data link or dedicated command and control link). At the same time, a flight path change notification will be issued (reason for replanning, new destination (if changed), change in estimated arrival time, and remaining battery power assessment). Track Connection: The flying car smoothly transitions from its current frozen state (hovering or level flight) to the starting segment of a new flight path. The starting state (position, speed, heading) of the new flight path is consistent with the current state of the flying car (because replanning starts from the current state), so the trajectory connection is continuous (no jumps in position, speed, or heading). The flight control system only needs to start tracking the new flight path from the current state, without the need for special transition maneuvers. If the starting speed of the new flight path differs from the current speed (e.g., the current is hovering while the new flight path starts at cruise speed), a smooth acceleration phase is used for the transition (within acceleration constraints). Flight Execution Resumption: The flight control system unfreezes and begins executing the flight mission along the new flight path. The dynamic window method for local obstacle avoidance and the online monitoring mechanism simultaneously resume operation (continuing to monitor environmental changes; if triggering conditions reappear, replanning is performed again, forming a continuous dynamic adaptive closed loop). Execution Status Monitoring: The mission management system continuously monitors the flying car's execution status along the new route (position tracking error, speed tracking error, remaining battery power, and deviation from the estimated arrival time at each point on the route). If significant deviations or new environmental changes occur, they are addressed promptly. Simultaneously, the system monitors the capacity status of the take-off and landing points along the new route (if the capacity of the take-off and landing points changes again before arrival, the route is replanned or the arrival time is adjusted). Mission Completion: Once the flying car reaches the new destination (original landing point or alternative take-off and landing point) along the new route and completes landing, the mission ends, and a complete route execution report is output (including the original route, replanning events, new route, actual execution trajectory, obstacle avoidance events, energy consumption, and time statistics) for post-mission analysis and route optimization.

[0078] Replanning performance assurance: The entire replanning process (trigger detection → freeze → data update → search → smoothing → verification → deployment) must be completed within an acceptable time, typically a few seconds to tens of seconds (eVTOL hovering waits can be longer, such as 30 seconds to 1 minute, while fixed-wing hovering waits need to be shorter, such as 10 to 30 seconds); real-time replanning is ensured through incremental replanning, efficient search algorithms (hybrid A* heuristic guidance and state pruning), parallel computing (multi-threaded parallel computation of cost graph updates and searches), and pre-computation (pre-calculating backup route plans and directly calling fine-tuning upon triggering); if the replanning time exceeds the flying car's waiting time (e.g., insufficient remaining battery power to support waiting), an emergency strategy is executed first (moving to the nearest safe waiting point or emergency landing point), and the mission continues only after replanning is completed at a safe location. Through a complete process of issuing new routes, smoothly connecting trajectories, execution recovery and monitoring, and performance assurance, the newly replanned routes were ensured to be executed safely and smoothly by the flying car. This enabled online adaptive adjustment of routes under dynamic environmental changes, ensuring mission continuity and safety, and solving the problems of traditional methods lacking online dynamic replanning mechanisms and mission failure or safety risks after environmental changes.

[0079] In some embodiments, the method further includes a refined trajectory generation step based on the switching between vertical takeoff and landing (VTOL) and cruise modes of the flying car: dividing the flight mission into a vertical takeoff phase, a climb transition phase, a horizontal cruise phase, a descent transition phase, and a vertical landing phase, and configuring a corresponding flying car dynamics model and energy consumption model for each phase; generating trajectories using the VTOL dynamics model in the VTOL and vertical landing phases, generating trajectories using the mode-switching dynamics model in the climb and descent transition phases, and generating trajectories using the fixed-wing cruise dynamics model in the horizontal cruise phase; performing trajectory continuity verification at the mode-switching points, verifying the positional continuity, velocity continuity, and acceleration continuity of trajectories in adjacent phases, and adjusting the position of the mode-switching point and the trajectory parameters of the corresponding phase if the continuity requirements are not met; calculating the total energy consumption for the entire flight based on the energy consumption model of each phase, and incorporating the total energy consumption and the temporal capacity cost at the takeoff and landing points into the cost calculation of the multi-constraint cost map.

[0080] The flight mission is divided into five phases: vertical takeoff, climb transition, horizontal cruise, descent transition, and vertical landing. A corresponding flying car dynamics model and energy consumption model are configured for each phase. Flying cars (electric vertical takeoff and landing aircraft, eVTOL) typically have two flight modes: vertical takeoff and landing (VTOL) and fixed-wing cruise. The dynamic characteristics, control methods, and energy consumption characteristics of these two modes differ significantly: VTOL mode (multi-rotor / lift fan mode) relies on the thrust of the rotors or fans to directly overcome gravity, enabling vertical takeoff, landing, and hovering, but it has low cruise efficiency, high energy consumption, and low speed; fixed-wing cruise mode (wing-lift mode) relies on the aerodynamic lift generated by the wings to overcome gravity, resulting in high cruise efficiency, low energy consumption, and high speed, but it requires a runway or sufficient forward speed and cannot hover. A transition is required between the two modes (climb transition from vertical takeoff / landing mode to cruise speed and then to fixed-wing mode, descent transition from fixed-wing mode to hover speed and then to vertical takeoff / landing mode). The trajectory generation and continuity during the transition phases are crucial for the feasibility of the entire flight trajectory. Mission phase division: Based on the flying car's modal characteristics and mission profile, the complete flight mission is divided into five consecutive phases: Vertical takeoff phase (vertical takeoff from the ground at the takeoff / landing point to a safe altitude, such as 50 to 100 meters, in vertical takeoff / landing mode throughout, with the rotor / fan providing all lift, and the speed approaching zero or climbing at low speed); Climb transition phase (starting from a safe altitude, the flying car gradually accelerates forward while maintaining lift, increasing speed from hover speed to cruise speed, while lift is gradually transferred from the rotor / fan to the wings, completing the transition from vertical takeoff / landing mode to fixed-wing mode, and simultaneously climbing to cruise altitude, such as 100 meters). The flight path is divided into five phases: the horizontal cruise phase (fully converting to fixed-wing mode, with the wings providing all lift, flying horizontally at cruise speed; this is the main part of the flight, with the lowest energy consumption and highest efficiency), the descent transition phase (decelerating from cruise, decreasing altitude, with lift gradually transferring from the wings to the rotor / fan, completing the transition from fixed-wing mode to vertical takeoff and landing mode, with speed decreasing from cruise speed to hovering speed, and altitude decreasing to a safe landing altitude, such as 50 to 100 meters), and the vertical landing phase (fully converting to vertical takeoff and landing mode, descending vertically from a safe altitude to the ground at the takeoff and landing point, hovering at low speed throughout the descent for a precise landing). These five phases are sequentially connected, covering the complete flight profile from takeoff to landing.

[0081] Configuration of dynamic and energy consumption models for each stage: A corresponding flying car dynamic model (describing the relationship between force and motion in that mode) and an energy consumption model (describing the relationship between power consumption and flight state in that mode) are configured for each stage. The vertical takeoff and vertical landing stages use a vertical takeoff and landing dynamic model: This model uses the thrust of the multi-rotor / lift fan as the main control force. The dynamic equations are: vertical direction m×a_z=T_total-m×g-D_z (T_total is the total thrust of all rotors / fans, m is the flight mass, g is the gravitational acceleration, and D_z is the vertical drag), horizontal direction m×a_xy=T_horizontal-D_xy (horizontal component force generated by rotor / fan tilt or differential speed); attitude dynamics include torque balance in three channels: roll, pitch, and yaw; energy consumption of the vertical takeoff and landing mode. The model is P_vtol = P_hover + P_climb × |v_z| + P_maneuver × |a| (P_hover is the hovering power, which is related to mass and rotor efficiency, such as P_hover = m × g × v_induced, where v_induced is the induced velocity; P_climb is the additional power for climbing, which is the extra work done during vertical climb; P_maneuver is the additional power for maneuvering, and the greater the acceleration, the higher the energy consumption). The energy consumption of the vertical take-off mode is much higher than that of the cruise mode (hovering efficiency is low, and the typical eVTOL hovering power is 2 to 3 times the cruise power).

[0082] The climb and descent transition phases employ a mode-switching dynamics model: this model is a hybrid of a vertical takeoff model and a fixed-wing model. Lift is provided by the vertical component of the rotor / fan thrust and the wing lift. The dynamic equations are: m×a_z=T_vertical+L_wing-m×g-D_z (T_vertical is the vertical component of the rotor / fan thrust, L_wing is the wing lift, which increases with speed, L_wing=0.5×ρ×v²×S×C_L), m×a_xy=T_horizontal-D_wing-D_fuselage (T_horizontal is the horizontal component of the thrust, D_wing is the wing drag, and D_fuselage is the fuselage drag); during mode switching, the thrust distribution coefficient changes from 1.0 (all...). The lift gradually decreases from the rotor to zero (all lift is provided by the wing), and the speed increases from hovering speed to cruise speed (climb transition) or vice versa (descent transition). The energy consumption model for the transition phase is P_transition=(1-α)×P_vtol+α×P_cruise+P_transition_extra, where α is the mode transition progress coefficient (0 for pure vertical takeoff, 1 for pure cruise), and P_transition_extra is the additional energy consumption during the transition phase (aerodynamic and control losses during mode transition; the transition phase is usually one of the most energy-intensive phases of the entire flight). The time and distance of the transition phase need to be calculated precisely (a typical eVTOL climb transition takes 30 to 60 seconds and a distance of 500 to 1500 meters, and the descent transition is similar), which is an indispensable part of the flight path planning.

[0083] The horizontal cruise phase employs a fixed-wing cruise dynamics model: this model uses wing lift to overcome gravity and thrust to overcome drag, with the dynamic equations being m×a_z=L_wing-m×g (during cruise, a_z≈0, lift equals gravity), and m×a_xy=T_propulsion-D_total (T_propulsion is the thrust of the propeller / fan, D_total is the total drag including wing-induced drag and parasitic drag, fuselage drag, and interference drag; D_total=0.5×ρ×v²×S×C). _D, C_D includes zero lift drag and induced drag C_D=C_D0+k×C_L²); the energy consumption model for fixed-wing mode is P_cruise=T_propulsion×v / η_propulsion (thrust multiplied by speed divided by propulsion efficiency; during cruise, thrust equals drag, P_cruise=D_total×v / η). Cruise mode has the lowest energy consumption (high lift-to-drag ratio; typical eVTOL cruise power is 1 / 3 to 1 / 2 of hover power), and cruise speed is usually 150 to 300 km / h. The model parameters for each stage (mass, aerodynamic coefficient, rotor efficiency, propulsion efficiency) are loaded from the flying car's performance database. Changes in mass and center of gravity under different payloads (different passenger / cargo weights) will affect the model parameters and need to be adjusted in real time according to the mission payload.

[0084] The vertical takeoff and vertical landing phases utilize a vertical takeoff and landing dynamics model to generate the trajectory; the climb and descent transition phases utilize a mode-switching dynamics model to generate the trajectory; and the horizontal cruise phase utilizes a fixed-wing cruise dynamics model to generate the trajectory. The trajectory generation for each phase uses the corresponding phase's dynamics model as a constraint to ensure that the generated trajectory is physically executable within that mode. Vertical takeoff trajectory generation: Starting from the ground position of the takeoff and landing point (z = takeoff and landing point elevation, v = 0, hovering state), the vertical takeoff dynamics model is used to vertically climb to a safe altitude (e.g., 50 to 100 meters, higher than surrounding obstacles and buildings, ensuring sufficient safe altitude during the transition phase) at the maximum safe climb rate (e.g., 2 to 5 meters / second, to avoid excessively rapid climb leading to loss of ground effect and control instability). During the climb, the position is kept horizontally stable (hovering accuracy, horizontal deviation less than 1 to 2 meters). The trajectory is a vertically upward line segment (or a slight translation to align with the transition direction). The time t_tof = h_safe / v_climb (e.g., 100 meters / 3 meters / second ≈ 33 seconds). Energy consumption is calculated by integral of the vertical takeoff energy consumption model. Vertical landing trajectory generation: Starting from a hovering state at a safe landing altitude (e.g., 50 to 100 meters), a vertical takeoff dynamics model is used to descend vertically to the takeoff and landing point at a safe descent rate (e.g., 1 to 3 meters / second; descent that is too fast may lead to approach instability and abrupt changes in ground effect). During the descent, the horizontal position is maintained and aligned with the takeoff and landing point (precise landing, horizontal deviation less than 0.5 to 1 meter). In the final stage (5 to 10 meters above the ground), the vehicle decelerates and descends gently for a soft landing. The trajectory is a vertically downward line segment. The time t_landing≈h_safe / v_descend. Energy consumption is integrated according to the vertical takeoff energy consumption model (the rotor thrust during descent is slightly less than gravity, and the energy consumption is lower than during hovering but still significant).

[0085] Climb transition trajectory generation: Starting from the end of the vertical takeoff phase (safe altitude, hovering state), a mode-switching dynamics model is used to generate a composite trajectory that accelerates from hovering to cruise speed, simultaneously transitions from vertical takeoff mode to fixed-wing mode, and climbs to cruise altitude. Key parameters of the transition trajectory include: acceleration curve (speed increases from 0 to v_cruise, using a smooth acceleration curve such as trapezoidal or S-shaped acceleration, with the maximum acceleration limited to 0.1 to 0.3g to ensure passenger comfort and structural safety), climb profile (altitude increases from h_safe to h_cruise, with a larger climb rate in the early transition phase, which decreases as wing lift increases with speed), mode transition progress α(t) (gradually changes from 0 to 1 linearly or S-shaped; at α=0.5, the rotor and wing each provide half of the lift, which is the midpoint of the transition, typically corresponding to when the speed reaches 1.2 to 1.3 times the stall speed, the wing can provide sufficient lift), and thrust distribution (rotor / fan thrust gradually decreases from hovering thrust). The transition trajectory must satisfy force balance (lift = gravity, thrust = drag or acceleration) and attitude constraints (pitch angle changes during the transition, from hovering attitude (large pitch or horizontal) to cruise attitude (small pitch), with a limited rate of change). The generated trajectory is a three-dimensional spatial curve that simultaneously includes altitude increase, speed increase, and mode transition (usually slightly upward, with climb and acceleration synchronized). The transition time t_climb_transition ≈ 30 to 60 seconds, and the horizontal transition distance ≈ 500 to 1500 meters (depending on the rate of acceleration and cruise speed). Trajectory generation for the descent transition phase: Starting from the end of the horizontal cruise phase (cruise altitude, cruise speed, fixed-wing mode), a mode-switching dynamics model is used to generate a composite trajectory that decelerates from cruise speed to hover speed, simultaneously transitions from fixed-wing mode to vertical takeoff mode, and descends to a safe landing altitude. This is the reverse process of the climb transition (deceleration, descent, rotor restart and gradual takeover of lift, wing lift gradually decreasing to zero). The deceleration curve uses smooth deceleration (maximum deceleration 0.1 to 0.2g), and the mode transition rate α gradually changes from 1 to 0. The transition time and distance are similar to the climb transition. Finally, upon reaching the safe landing altitude, the speed is zero (hovering state), and the aircraft fully transitions to vertical takeoff mode. The transition phase trajectory is the most dynamically complex part of the entire flight and must be accurately generated using a mode-switching model; otherwise, problems such as trajectory discontinuities and abrupt attitude changes at the mode-switching points will occur.

[0086] Horizontal cruise phase trajectory generation: Starting from the end point of the climb transition phase (cruise altitude, cruise speed, fixed-wing mode) and ending at the start point of the descent transition phase (also with cruise altitude, cruise speed, fixed-wing mode), a three-dimensional trajectory for the horizontal cruise segment is generated using a fixed-wing cruise dynamics model. The cruise segment trajectory constitutes the main body of the flight (typically accounting for 80% to 95% of the total flight distance). The flight proceeds horizontally at a cruise altitude (e.g., 100 to 300 meters, selected based on urban low-altitude airspace regulations and wind field optimization; wind fields vary at different altitudes, so a tailwind and low-turbulence altitude can be selected) and a cruise speed (selected based on mission time requirements and energy consumption optimization; typically, 70% to 90% of the maximum cruise speed is the optimal lift-to-drag ratio speed with the lowest energy consumption). The spatial path of the cruise segment trajectory is determined by a hybrid A* search and B-spline smoothing (the optimal path searched on a multi-constraint cost map, passing through no-fly zone avoidance, building wind field optimization, and takeoff and landing point capacity constraints). The velocity profile of the cruise segment is automatically adjusted according to the path curvature (deceleration at turns to ensure lateral acceleration does not exceed limits, and cruise speed is restored on straight sections). The altitude profile maintains a safe obstacle clearance altitude based on terrain and buildings (50 to 100 meters above the highest building, or within the preset low-altitude corridor altitude layer). The cruise segment trajectory must meet fixed-wing dynamics constraints (speed between stall speed and maximum speed, turning radius greater than minimum turning radius, climb / descent rate within the allowable range of fixed-wing modes). Energy consumption is calculated by integral of the fixed-wing cruise energy consumption model (P_cruise=D×v / η, integral along the flight path to obtain the total energy consumption of the cruise segment). The trajectories of the five stages are sequentially spliced ​​to form a complete full-fledged trajectory from takeoff to landing. Each stage is generated using the corresponding mode's dynamics model, ensuring the physical executability of the trajectory in each mode.

[0087] At the mode switching points, trajectory continuity is verified, checking the positional, velocity, and acceleration continuity of the trajectories in adjacent stages. If the continuity requirements are not met, the position of the mode switching point and the trajectory parameters of the corresponding stage are adjusted. The mode switching point is the connection point between two adjacent stages (the switching point between vertical takeoff and climb transition, climb transition and level cruise, level cruise and descent transition, and descent transition and vertical landing). At these points, the flight mode of the flying car changes (or is about to change), and the trajectory must remain continuous. Otherwise, position jumps, sudden velocity changes, or sudden acceleration changes will occur, leading to drastic changes in the flying car's attitude, difficulty in control, or even instability. These are the key verification points for the feasibility of the trajectory throughout the entire flight. Continuity verification is performed at each mode switching point, verifying continuity at three levels: Positional continuity: The 3D coordinates (x_end, y_end, z_end) of the endpoint of the previous stage trajectory must be exactly the same as the 3D coordinates (x_start, y_start, z_start) of the starting point of the next stage trajectory (x_start, y_start, z_start) (or the deviation must be within the allowable range, such as less than 0.5 meters; in practice, they should be exactly equal because the splicing point is the same point), i.e., there is no positional jump; Velocity continuity: The velocity vector (vx_end, vy_end, vz_end) of the endpoint of the previous stage must be the same as the velocity vector (vx_start, vy_start, vz_start) of the starting point of the next stage (vx_start, vy_start, vz_start) (or the deviation must be less than 0.5 meters). 0.5 m / s², meaning no sudden changes in speed (the speed at the end of the climb transition should be the cruising speed, consistent with the speed at the start of the horizontal cruise; the speed at the end of the descent transition should be zero, i.e., hovering, consistent with the speed at the start of the vertical descent); acceleration continuity: the acceleration vector (ax_end, ay_end, az_end) at the end of the previous stage and the acceleration vector (ax_start, ay_start, az_start) at the start of the next stage should be the same or smoothly transition (the deviation is less than 0.5 m / s²), that is, no sudden changes in acceleration (a sudden change in acceleration means a sudden change in thrust or pull, which will lead to a sudden change in attitude and passenger discomfort. C2 continuity, i.e., acceleration continuity, is the ideal goal, and at least the rate of change of acceleration jerk should be limited).

[0088] Verification method: At each switching point, calculate the position, velocity (first derivative), and acceleration (second derivative) of the previous stage trajectory at the endpoint and the corresponding values ​​of the next stage trajectory at the starting point. Calculate the difference and compare it with the allowable threshold. The position difference should be strictly zero (because the same point should be used as the common endpoint of the two stages when stitching the trajectories), the velocity difference should be less than the threshold (e.g., 0.5 m / s), and the acceleration difference should be less than the threshold (e.g., 0.5 m / s²) or a smooth transition should be ensured through jerk constraints. If the continuity requirement is not met, adjustments are made: Adjust the position (spatial or temporal) of the mode switching point: for example, move the switching point from climb transition to cruise along the transition trajectory to make the velocity and acceleration at the switching point better match the initial conditions of the cruise segment; Adjust the trajectory parameters of the corresponding stage: for example, adjust the acceleration curve of the climb transition (increase or decrease the maximum acceleration, change the shape of the acceleration curve) to make the final velocity / acceleration consistent with the starting point of the cruise segment, or adjust the initial velocity / altitude of the cruise segment to match the end point of the transition segment; Increase the length of the transition segment: if the mode transition cannot be completed and continuity cannot be achieved within the existing transition distance, extend the time and distance of the transition segment (reduce the acceleration rate, increase the transition time) to make the mode transition smoother and the continuity easier to meet; Iterate and adjust until the position, velocity, and acceleration of all switching points meet the continuity requirement. Continuity verification and adjustment ensure that the trajectory is smooth and continuous at the mode switching points throughout the entire flight. The flying car can seamlessly transition from one mode to the next without attitude changes or control jumps. This solves the problems of traditional methods that use a single dynamic model throughout the entire flight, have discontinuous trajectories at mode switching points, and have the risk of attitude changes during actual execution.

[0089] The total energy consumption for the entire flight is calculated based on the energy consumption model for each stage, and the total energy consumption and the timing capacity cost at takeoff and landing points are included in the cost calculation of the multi-constraint cost map. Total energy consumption calculation for the entire flight: The energy consumption of each of the five stages is calculated by integrating along the trajectory according to the energy consumption model of the corresponding stage, and then summed to obtain the total energy consumption for the entire flight: E_total = E_vto + E_climb_trans + E_cruise + E_descend_trans + E_vland. Where E_vto is the energy consumption of the vertical takeoff stage (integration of the vertical takeoff model, climb work + hovering loss), E_climb_trans is the energy consumption of the climb transition stage (integration of the mode switching model, acceleration + climb + conversion loss), E_cruise is the energy consumption of the horizontal cruise stage (integration of the fixed-wing model, thrust × distance / efficiency), E_descend_trans is the energy consumption of the descent transition stage (integration of the mode switching model, deceleration + descent + conversion loss, there may be energy recovery during deceleration such as regenerative braking, and the energy of the eVTOL descent transition part can be recovered), and E_vland is the energy consumption of the vertical landing stage (integration of the vertical takeoff model, the thrust is slightly less than the gravity during descent, and the energy consumption is lower than that of takeoff). Total energy consumption calculations need to consider changes in the flying car's mass (fuel consumption or battery discharge; eVTOL, powered by batteries, maintains a relatively constant mass, but the payload may change after unloading), environmental conditions (air density ρ varies with altitude, affecting aerodynamic lift and drag and rotor efficiency; wind field affects ground speed and energy consumption; tailwinds save energy, headwinds consume energy), and flight parameters (the impact of speed, altitude, rate of climb, turning, etc., on energy consumption). The total energy consumption E_total is compared with the available battery energy E_battery (remaining charge × total battery energy × depth of discharge coefficient). If E_total ≤ E_battery × (1 - safety margin coefficient) (safety margin such as 20%, reserving emergency charge), then the flight is feasible within the endurance limit; if it exceeds this limit, the flight path needs to be adjusted (shortening the flight, selecting a more energy-efficient altitude / speed, adding alternate landing points) or an indication of insufficient endurance should be provided.

[0090] Incorporating total energy consumption into the multi-constraint cost graph: In the comprehensive cost calculation of the multi-constraint cost graph, in addition to the original no-fly zone obstacle avoidance cost, building wind field cost, and take-off and landing point time-series capacity cost, an energy consumption cost item C_energy is added. The energy consumption cost is obtained by integrating along the path based on the local energy consumption rate (energy consumption per unit distance, calculated from the wind field conditions (tailwind / headwind / turbulence), altitude layer, and flying car mode and speed) when the flight path passes through the grid, i.e., C_energy=w_energy×∫e_local(s)ds, where e_local(s) is the energy consumption per unit distance at position s on the path (joules / meter or kilowatt-hours / kilometer), and w_energy is the energy consumption cost weight; the energy consumption cost overlaps with the energy cost in the building wind field cost but is more comprehensive (the building wind field cost mainly considers the local wind field influence of building flow, while the energy consumption cost considers the modal energy consumption and global wind field of the entire flight path), and the weights can be combined or configured separately. Incorporating total energy consumption and takeoff / landing point time-series capacity costs into the cost calculation of the multi-constraint cost graph means that route search considers not only takeoff / landing point waiting time costs but also total flight energy consumption costs. This ensures that the searched routes are optimal across multiple objectives, including safety (no-fly zones), economy (energy consumption), timeliness (takeoff / landing point waiting time), and stability (building wind fields). The introduction of energy consumption costs enables hybrid A* search to proactively select more energy-efficient paths (such as tailwind altitudes, optimal cruise speeds, and reduced unnecessary climbs, descents, and turns), reducing operating costs and carbon emissions. At the same time, it ensures that the flight range remains within battery endurance (paths with high energy costs, i.e., high energy consumption, are suppressed to avoid searching for paths exceeding the flight range). By implementing a complete mechanism that includes five-stage dynamic model configuration, generation of corresponding model trajectories for each stage, continuity verification and adjustment of mode switching points, full-range energy consumption calculation, and multi-cost fusion, a refined trajectory generation for the entire flight range considering the mode switching characteristics of vertical take-off and landing and fixed-wing cruise is achieved. The trajectory is physically executable in each mode, the mode switching points are smooth and continuous, and the energy consumption is accurately calculable. This solves the problems of traditional methods, such as a single dynamic model throughout the entire flight, discontinuous trajectory at mode switching points, and the risk of attitude change.

[0091] In some embodiments, please refer to Figure 2The overall methodology is as follows: After starting, input flight mission parameters (takeoff point, landing point, waypoints, estimated time window); load multi-source data for the urban low-altitude corridor (static no-fly zone, 3D building model, real-time wind field, takeoff and landing point capacity, multi-source data fusion); create a 3D dynamic map and divide it into basic grid units; calculate in parallel the no-fly zone obstacle avoidance cost, building wind field cost, and takeoff and landing point temporal capacity cost for each grid unit; perform multi-constraint cost map fusion; use a hybrid A* search algorithm to search for the optimal path; use a B-spline curve smoothing algorithm to smooth the optimal path, combined with a dynamic window method to perform local obstacle avoidance; determine if replanning is needed (if a sudden no-fly zone is detected or the takeoff and landing point capacity is too low); if so, return to reload multi-source data and perform route replanning; if not, output the spatiotemporal trajectory point sequence; end. Please refer to [reference needed]. Figure 3 The building wind field cost calculation process is as follows: Input the city's 3D building outline model and real-time meteorological wind field (wind direction, wind direction, turbulence intensity); mark and classify each cell in the low-altitude corridor grid map; determine whether the current grid is located within the building's influence area, and if so, determine the building's relative location type (canyon area, roof area, leeward area); for canyon areas, calculate the wind speed amplification factor based on the funnel effect model, determine the wind direction / headwind / crosswind, and calculate the energy cost; for roof areas, introduce the downwash model, calculate the vertical wind shear intensity, and calculate the stability cost; for leeward areas, mark them as high turbulence areas, calculate the turbulent kinetic energy, and calculate the stability cost; calculate the comprehensive wind field cost of this grid (energy cost = k × wind direction / headwind coefficient, stability cost = k × turbulence / wind shear coefficient); if not located within the building's influence area (open area), then the wind field cost = basic atmospheric wind field cost (no building disturbance); determine whether all grids have been processed, if not, return to process the next grid, and if completed, output a 3D wind field cost map (each grid includes an energy gain coefficient + turbulence risk coefficient).

[0092] Please refer to Figure 4 The no-fly zone obstacle avoidance cost calculation process is as follows: Input static no-fly zone database, temporary no-fly zone notice, and grid map; traverse every grid cell in the current low-altitude corridor; determine whether the center point of the current grid cell is located inside any no-fly zone. If so, perform hard constraint processing (set the grid cost = +∞, absolutely impassable); otherwise, calculate the Euclidean distance d from the grid cell to the nearest no-fly zone boundary; set the safety buffer distance D_safe (configurable parameter); determine whether d is less than D_safe. If so, perform soft constraint processing (cost = A × exp(-σ × d), exponential potential field, the closer to the boundary, the higher the cost, A is the maximum cost coefficient, and σ is the attenuation factor); otherwise, perform safe zone processing (cost = 0); determine whether all grid cells have been processed. If not, return to process the next grid cell. If all are processed, output the 3D no-fly zone cost map.

[0093] Please refer to Figure 5 The calculation process for takeoff and landing point time-series capacity cost is as follows: Input the destination takeoff and landing points, warning arrival time window, and a list of nearby available takeoff and landing points from the flight plan; obtain the remaining capacity C_rem, queue length L_q, and average service time T_service of the destination takeoff and landing points in real time via API; calculate the estimated waiting time T_wait = L_q × T_service based on the M / M / 1 queuing theory model; determine whether T_wait is ≤ preset threshold T_max, if so, the capacity cost = w4 × T_wait (linearly mapping the waiting time to time cost); otherwise, search for nearby alternative takeoff and landing points, and filter those that meet the requirements of sufficient capacity and distance increase not exceeding d_max; determine whether there are available alternative takeoff and landing points, if so, change the destination to an alternative takeoff and landing point and recalculate the capacity cost; if not, execute forced waiting (capacity cost = +∞, takeoff and hover waiting strategy); output the takeoff and landing point capacity cost (time cost) and possible alternate landing point information.

[0094] Please see Figure 6 As shown, Figure 6 This is a schematic diagram of the structure of a refined flight path generation system 200 for urban low-altitude corridors provided in this application embodiment. This refined flight path generation system 200 for urban low-altitude corridors is used to execute the steps of the refined flight path generation method for urban low-altitude corridors shown in the above embodiments. The refined flight path generation system 200 for urban low-altitude corridors can be a single server or a server cluster, or it can be a terminal, such as a handheld terminal, laptop computer, wearable device, or robot. Figure 6 As shown, the refined flight path generation system 200 for flying cars oriented towards urban low-altitude corridors includes: Data loading unit 201 is used to load static no-fly zones, 3D building models, real-time meteorological wind fields and take-off and landing point capacity multi-source data of urban low-altitude corridors, to establish a 3D dynamic map and divide it into basic grid units. Cost calculation unit 202 is used to calculate the no-fly zone obstacle avoidance cost, building wind field cost and take-off and landing point temporal capacity cost for each grid cell, and merge them to generate a multi-constraint cost map. The no-fly zone obstacle avoidance cost sets the grid inside the no-fly zone as impassable and calculates the cost of the grid within the safe buffer distance according to the exponential potential field. The building wind field cost calculates the energy cost and stability cost according to the grid relative to the building position. The take-off and landing point temporal capacity cost calculates the waiting time cost according to the queuing theory model. The route planning unit 203 is used to search for the optimal path on the multi-constraint cost map using a hybrid A* search algorithm, smooth the optimal path using B-spline curves, perform local obstacle avoidance using a dynamic window method, and generate a spatiotemporal trajectory point sequence. When a sudden no-fly zone or take-off and landing point capacity is detected to be lower than a preset threshold, the data is reloaded and route replanning is performed.

[0095] It should be noted that those skilled in the art will understand that, for the sake of convenience and brevity, the specific working processes of the refined flight path generation system for urban low-altitude corridors and its modules described above can be found in the corresponding embodiments of the refined flight path generation method for urban low-altitude corridors, and will not be repeated here. The aforementioned refined flight path generation method for urban low-altitude corridors can be implemented as a computer program, which can be used in... Figure 6 It runs on the system shown.

[0096] Please see Figure 7 , Figure 7 This is a schematic block diagram of the computer device provided in an embodiment of this application. The computer device includes a processor, a memory, and a network interface connected via a device bus. The memory may include a storage medium and internal memory. The storage medium stores operating devices and computer programs. The computer program includes program instructions, which, when executed, cause the processor to execute any method for generating refined flight routes for flying cars oriented towards urban low-altitude corridors. The processor provides computational and control capabilities to support the operation of the entire computer device. The internal memory provides an environment for the execution of the computer program in the non-volatile storage medium; when executed by the processor, this computer program causes the processor to execute any method for generating refined flight routes for flying cars oriented towards urban low-altitude corridors.

[0097] This network interface is used for network communication, such as sending assigned tasks. Those skilled in the art will understand that... Figure 7The structures shown are merely block diagrams of a portion of the structure related to the present application and do not constitute a limitation on the terminal to which the present application is applied. Specific computer devices may include more or fewer components than shown in the figures, or combine certain components, or have different component arrangements. It should be understood that the processor may be a Central Processing Unit (CPU), or it may be other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The general-purpose processor may be a microprocessor or any conventional processor.

[0098] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any person skilled in the art can easily conceive of various equivalent modifications or substitutions within the technical scope disclosed in this application, and these modifications or substitutions should all be covered within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for generating refined flight routes for flying cars in urban low-altitude corridors, characterized in that, include: Load static no-fly zones, 3D building models, real-time meteorological wind fields, and multi-source data on take-off and landing point capacity of urban low-altitude corridors, establish a 3D dynamic map, and divide it into basic grid units; The no-fly zone obstacle avoidance cost, building wind field cost, and take-off and landing point temporal capacity cost are calculated for each grid cell. These costs are then fused to generate a multi-constraint cost map. The no-fly zone obstacle avoidance cost is calculated by setting the grid inside the no-fly zone as impassable and calculating the cost of the grid within the safe buffer distance according to the exponential potential field. The building wind field cost is calculated by classifying the grid relative to the building position and calculating the energy cost and stability cost. The take-off and landing point temporal capacity cost is calculated by queuing theory model to calculate the waiting time cost. The optimal path is searched on the multi-constraint cost map using a hybrid A* search algorithm, and the optimal path is smoothed using B-spline curves. Local obstacle avoidance is performed in combination with the dynamic window method to generate a spatiotemporal trajectory point sequence. When a sudden no-fly zone or take-off and landing point capacity is detected to be lower than the preset threshold, the data is reloaded and the route replanning is performed.

2. The method according to claim 1, characterized in that, The loading of static no-fly zones, 3D building models, real-time meteorological wind fields, and multi-source data on takeoff and landing point capacity for urban low-altitude corridors is used to establish a 3D dynamic map and divide it into basic grid units, including: Load the static no-fly zone database and temporary no-fly zone notices, and extract the spatial boundaries, altitude range, and effective time of the no-fly zones; Load the city's 3D building outline model and extract the building's geographical location, height, and shape parameters; Load real-time meteorological wind field data and extract wind direction, wind speed and turbulence parameters at different altitude levels; Load real-time status data of take-off and landing points, and extract parameters such as geographical location, remaining capacity, queue length, and average service time of the take-off and landing points; Static no-fly zone data, 3D building model data, real-time meteorological wind field data, and real-time status data of take-off and landing points are integrated into a unified coordinate system to create a 3D dynamic map. The 3D dynamic map is then divided into basic grid units according to a set spatial resolution.

3. The method according to claim 1, characterized in that, The calculation of the no-fly zone obstacle avoidance cost, building wind field cost, and takeoff and landing point time-series capacity cost for each grid cell is fused to generate a multi-constraint cost map, including: Traverse each basic grid cell in the 3D dynamic map and calculate the no-fly zone obstacle avoidance cost, building wind field cost, and take-off and landing point time sequence capacity cost for each basic grid cell. Configure corresponding fusion weight coefficients for obstacle avoidance cost in no-fly zones, building wind field cost, and take-off and landing point time-series capacity cost, respectively; The no-fly zone obstacle avoidance cost, building wind field cost, and take-off and landing point time-series capacity cost of each basic grid cell are weighted and summed according to the fusion weight coefficient to generate the comprehensive cost of each basic grid cell. The combined cost of all basic grid cells is combined according to their spatial location to generate a multi-constraint cost map.

4. The method according to claim 3, characterized in that, The obstacle avoidance cost in the no-fly zone is calculated by setting the grid inside the no-fly zone as impassable, calculating the cost of the grid within the safe buffer distance using an exponential potential field, calculating the energy cost and stability cost of the building wind field according to the grid's relative building position, and calculating the waiting time cost of the take-off and landing point time-series capacity cost using a queuing theory model, including: Determine whether the center point of the basic grid cell is located inside the no-fly zone. If it is located inside the no-fly zone, set the obstacle avoidance cost of the no-fly zone to be impassable. If it is not located inside the no-fly zone, calculate the distance from the basic grid cell to the nearest boundary of the no-fly zone. If the distance is less than the safe buffer distance, calculate the obstacle avoidance cost of the no-fly zone according to the exponential potential field. If the distance is greater than or equal to the safe buffer distance, set the obstacle avoidance cost of the no-fly zone to zero. Determine the location type of the basic grid cell relative to the building. When the location type is an open area, calculate the building wind field cost according to the basic atmospheric wind field. When the location type is a building canyon area, calculate the wind speed amplification factor and energy cost according to the septum effect model. When the location type is a building roof area, calculate the vertical shear intensity and stability cost according to the downwash flow model. When the location type is a building leeward area, calculate the turbulent kinetic energy and stability cost according to the high turbulence model. Based on the queue length and average service time of the take-off and landing points, the waiting time of the take-off and landing points is calculated according to the queuing theory model. The waiting time is linearly mapped to the time capacity cost of the take-off and landing points. When the waiting time exceeds a preset threshold, nearby alternative take-off and landing points are searched and the time capacity cost of the take-off and landing points is recalculated.

5. The method according to claim 1, characterized in that, The method of using a hybrid A* search algorithm to search for the optimal path on a multi-constraint cost graph includes: Use the takeoff point in the flight mission parameters as the search start point and the landing point in the flight mission parameters as the search end point to initialize the open node set and the closed node set; Select the node with the lowest overall cost value from the open node set as the current expansion node, and move the current expansion node from the open node set to the closed node set. The neighboring nodes of the current extended node are generated according to the kinematic constraints of the flying car. The actual cost and heuristic estimated cost of each neighboring node are calculated. The actual cost and heuristic estimated cost are added together to obtain the comprehensive cost of the neighboring node. Add neighboring nodes that are not in the closed node set to the open node set, and repeatedly select the node with the smallest comprehensive cost value for expansion until the search endpoint is reached. The search originates from the end point and traces back to the beginning point along the parent node pointer to generate the optimal path node sequence.

6. The method according to claim 5, characterized in that, The smoothing of the optimal path using B-spline curves includes: The nodes in the optimal path node sequence are used as control points of the B-spline curve, and the order of the B-spline curve is determined based on the number of control points. Weighted interpolation is performed on the control points according to the B-spline basis functions to generate a continuous and smooth B-spline curve; A curvature continuity check is performed on the B-spline curve. When the curvature of the B-spline curve exceeds the maximum turning curvature threshold of the flying car, the position of the corresponding control point is adjusted and the B-spline curve is regenerated until all curvatures of the B-spline curve satisfy the kinematic constraints of the flying car.

7. The method according to claim 6, characterized in that, The process of combining the dynamic window method to perform local obstacle avoidance and generate a spatiotemporal trajectory point sequence includes: Sample along the smoothed B-spline curve at set time intervals to generate an initial spatiotemporal trajectory point sequence. Each spatiotemporal trajectory point contains spatial coordinates, timestamps, and velocity information. During the flight car's trajectory execution, a dynamic speed window is determined based on the current speed and acceleration constraints, and all candidate speed combinations are traversed within the dynamic speed window; For each candidate velocity combination, a forward simulation is performed to predict the trajectory of the flying car within a set prediction time domain, and the minimum distance between each predicted trajectory and the obstacle and the degree of deviation from the B-spline curve are calculated. Select the candidate velocity combination that has the largest distance from the obstacle and the smallest deviation from the B-spline curve as the current execution velocity, and update the spatiotemporal trajectory point sequence according to the current execution velocity; The updated spatiotemporal trajectory point sequence is used as the final output.

8. The method according to claim 1, characterized in that, When a sudden no-fly zone or take-off / landing point capacity is detected to be below a preset threshold, the data is reloaded and route replanning is performed, including: During the flight route operation, receive airspace control notices and take-off and landing point status updates in real time, and detect whether there are sudden no-fly zones or whether the remaining capacity of the take-off and landing points is lower than the preset capacity threshold. When a sudden no-fly zone or take-off and landing point is detected to have a remaining capacity lower than the preset capacity threshold, the current flight path is frozen, and the current position, current speed, and remaining battery status of the flying car are recorded. Reload the updated static no-fly zone data, 3D building model data, real-time meteorological wind field data, and real-time take-off and landing point status data. Use the current location of the flying car as the new search starting point and the original landing point or alternative take-off and landing point as the new search endpoint, and re-execute the route planning. The newly generated spatiotemporal trajectory point sequence will be sent to the flying car for execution, restoring the flight path execution status.

9. The method according to claim 1, characterized in that, The method also includes a refined trajectory generation step based on the switching between vertical takeoff and landing mode and cruise mode of the flying car: The flight mission is divided into the vertical takeoff phase, the climb transition phase, the horizontal cruise phase, the descent transition phase, and the vertical landing phase, and a corresponding flying car dynamics model and energy consumption model are configured for each phase. The vertical takeoff and vertical landing phases use a vertical takeoff and landing dynamics model to generate the trajectory, the climb transition phase and descent transition phase use a mode switching dynamics model to generate the trajectory, and the horizontal cruise phase uses a fixed-wing cruise dynamics model to generate the trajectory. At the mode switching point, a trajectory continuity check is performed to verify the positional continuity, velocity continuity, and acceleration continuity of the trajectory in adjacent stages. If the continuity requirements are not met, the position of the mode switching point and the trajectory parameters of the corresponding stage are adjusted. The total energy consumption for the entire flight is calculated based on the energy consumption model for each stage, and the total energy consumption and the timing capacity cost at takeoff and landing points are included in the cost calculation of the multi-constraint cost map.

10. A refined flight path generation system for flying cars oriented towards urban low-altitude corridors, used to implement the method as described in any one of claims 1-9, characterized in that, include: The data loading unit is used to load static no-fly zones, 3D building models, real-time meteorological wind fields, and multi-source data on take-off and landing point capacity for urban low-altitude corridors, to create a 3D dynamic map and divide it into basic grid units. The cost calculation unit is used to calculate the no-fly zone obstacle avoidance cost, building wind field cost, and take-off and landing point temporal capacity cost for each grid cell, and merge them to generate a multi-constraint cost map. The no-fly zone obstacle avoidance cost sets the grid inside the no-fly zone as impassable and calculates the cost of the grid within the safe buffer distance according to the exponential potential field. The building wind field cost calculates the energy cost and stability cost according to the grid relative to the building position. The take-off and landing point temporal capacity cost calculates the waiting time cost according to the queuing theory model. The route planning unit is used to search for the optimal path on a multi-constraint cost map using a hybrid A* search algorithm, smooth the optimal path using B-spline curves, perform local obstacle avoidance using a dynamic window method, and generate a spatiotemporal trajectory point sequence. When a sudden no-fly zone or take-off and landing point capacity is detected to be lower than a preset threshold, the data is reloaded and route replanning is performed.

Images