Mountainous area electric power emergency multi-machine and ground station cooperative task allocation processing method and system

Abstract

The invention provides a mountainous area electric power emergency multi-machine and ground station cooperative task allocation processing method and system, and relates to the technical field of electric power emergency, and the method comprises the steps: 1, analyzing an original task request, generating a structured task description containing the coordinates of a first-aid repair target point, the quality of a to-be-transported material, the envelope size of the to-be-transported material, and a task time window, based on the first-aid repair target point coordinates, elevation and meteorological data of the corresponding area are obtained, and a three-dimensional environment model of the task area is generated; 2, performing spatial rasterization sampling on the three-dimensional environment model to obtain a spatial discrete sampling field, and identifying a target first-aid repair structure sampling cluster and an airspace channel sampling cluster to obtain a reference direction vector and an offset direction vector; through three-dimensional environment modeling, risk quantification airspace network construction, adaptive formation screening and precise control and cooperative correction, mountainous area complex terrain and airflow interference are effectively avoided, and the safety, precision and execution efficiency of mountainous area power emergency material cooperative transportation are improved.

Description

Technical Field

[0001] This invention relates to the field of power emergency technology, and in particular to a method and system for the collaborative task allocation and processing of multiple generators and ground stations in mountainous power emergency situations. Background Technology

[0002] In emergency power repair in mountainous areas, when it is necessary to transport insulators and other components to a damaged transmission tower located on a steep mountain ridge, it is feasible to use multiple drones to coordinate the hoisting through the load-bearing network. Currently, when planning such coordinated hoisting missions, a three-dimensional environmental model of the mission area is usually performed, and a flight path from the starting point to the target point is planned. However, existing methods may have insufficient consideration when constructing the final approach airspace near the target point (such as the tower). Specifically, most of these methods treat the target point as a simple three-dimensional coordinate point without fully considering the disturbance effects of the structure of the repair target itself (such as the truss of the tower) on the surrounding wind field, such as shielding, acceleration, or flow around it.

[0003] For example, the leeward side of a steel tower may form a specific low-speed zone or vortex zone due to structural obstruction. If existing planning methods rely solely on large-scale average wind field data for airspace safety assessment, they may fail to identify such local and specific airflow patterns caused by the target structure. This could result in the planned final approach path and hovering point being located in an airspace with actual unstable disturbances. When the drone formation carrying supplies arrives at this point to prepare for precise delivery or hovering docking, it may encounter unexpected local airflow impacts, increasing the difficulty and risk of flight control. Summary of the Invention

[0004] The technical problem to be solved by this invention is to provide a method and system for the collaborative task allocation and processing of multiple emergency power units and ground stations in mountainous areas, so as to improve the safety and reliability of collaborative transportation.

[0005] To solve the above-mentioned technical problems, the technical solution of the present invention is as follows: Firstly, a method for coordinating task allocation between multiple emergency power units and ground stations in mountainous power supply emergencies, the method comprising: Step 1: Parse the original task request and generate a structured task description that includes the coordinates of the target point to be repaired, the quality and envelope size of the materials to be transported, and the task time window. Based on the coordinates of the target point to be repaired, obtain the elevation and meteorological data of the corresponding area and generate a three-dimensional environment model of the task area. Step 2: Spatial rasterization sampling is performed on the three-dimensional environment model to obtain a spatial discrete sampling field, and the target emergency repair structure sampling cluster and the spatial channel sampling cluster are identified to obtain the reference direction vector and the offset direction vector. Step 3: Based on the centroid of the target emergency repair structure sampling cluster, the reference direction vector and the offset direction vector, generate a three-dimensional spatial cone, and based on the envelope size of the material to be transported and the three-dimensional spatial cone, generate a three-dimensional feasible spatial network with cone risk weights. Step 4: Based on the three-dimensional feasible airspace network and the quality of the materials to be transported, candidate UAV formation configurations are selected, and based on the candidate UAV formation configurations and the three-dimensional feasible airspace network, a set of suitable formation configurations and safe operation boundaries corresponding to each path are generated. Step 5: Based on the three-dimensional feasible airspace network, mission time window and adaptive formation configuration set, generate mission execution plan and continuous motion path and formation configuration transition manifold of each UAV in the execution UAV set; Step 6: Based on the continuous motion path and formation configuration transition manifold of each UAV, generate the control command sequence of each UAV and send it to the execution UAV set to perform the collaborative transportation task.

[0006] Secondly, the mountainous area power emergency multi-unit and ground station collaborative task allocation and processing system includes: The parsing module is used to parse the original task request, generate a structured task description that includes the coordinates of the repair target point, the quality and envelope size of the materials to be transported, and the task time window. Based on the coordinates of the repair target point, it obtains the elevation and meteorological data of the corresponding area and generates a three-dimensional environment model of the task area. The identification module is used to perform spatial rasterization sampling on the three-dimensional environment model to obtain a spatial discrete sampling field, and to identify the target emergency repair structure sampling cluster and the spatial channel sampling cluster to obtain the reference direction vector and the offset direction vector. The module is used to generate a three-dimensional spatial cone based on the centroid of the sampling cluster of the target emergency repair structure, the reference direction vector and the offset direction vector, and to generate a three-dimensional feasible spatial network with cone risk weights based on the envelope size of the material to be transported and the three-dimensional spatial cone. The adaptation module is used to screen candidate UAV formation configurations based on the three-dimensional feasible airspace network and the quality of the goods to be transported, and to generate a set of adapted formation configurations and safe operation boundaries for each path based on the candidate UAV formation configurations and the three-dimensional feasible airspace network. The execution module is used to generate a mission execution plan and the continuous motion path and formation transition manifold of each UAV in the execution UAV set based on a three-dimensional feasible airspace network, mission time window and adaptive formation configuration set; The control module is used to generate a sequence of control commands for each UAV based on the continuous motion path and formation configuration transition manifold of each UAV, and then send them to the UAV set to perform the collaborative transportation task.

[0007] Thirdly, a computer-readable storage medium storing a program that, when executed by a processor, implements the method.

[0008] The above-described solution of the present invention has at least the following beneficial effects: By constructing a 3D environmental model incorporating terrain elevation and wind field vector attributes, and combining this with specialized identification of the target repair structure and airspace channels, the system captures local airflow characteristics around the target caused by structural obstruction and bypassing, preventing UAVs from encountering sudden airflow impacts when approaching the target or hovering for docking, thus reducing flight control risks. Furthermore, by generating a 3D feasible airspace network with cone-shaped risk weights, the system precisely matches material transportation needs with airspace safety status, providing quantifiable flight space data for UAV formation planning, effectively avoiding obstacles in complex mountainous terrain and unstable airflow areas, and improving the safety and rationality of flight paths. By selecting suitable UAV formation configurations based on mission requirements and airspace characteristics, and clearly defining the safe operation boundaries under each path, the system ensures that the formation's payload capacity, anti-interference capability, and mission scenario are highly compatible. This avoids low transportation efficiency or safety hazards caused by improper formation selection and improves the reliability of multi-UAV collaborative transportation. By generating continuous motion paths and smooth formation configuration transition manifolds, the system enables smooth switching of UAV formations at different flight stages, reducing the control difficulty caused by attitude changes. At the same time, combined with the issuance of control command sequences, the system ensures the consistency and synchronization of multi-UAV collaborative actions, improving the accuracy of material transportation and mission execution efficiency. Attached Figure Description

[0009] Figure 1 This is a flowchart illustrating the method for collaborative task allocation and processing between multiple generators and ground stations in mountainous power emergency response, provided by an embodiment of the present invention. Figure 2 This is a schematic diagram of a multi-machine and ground station collaborative task allocation and processing system for emergency power supply in mountainous areas, provided by an embodiment of the present invention. Detailed Implementation

[0010] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.

[0011] like Figure 1 As shown, embodiments of the present invention propose a method for collaborative task allocation and processing between multiple generators and ground stations in mountainous power emergency response, the method comprising the following steps: Step 1: Parse the original task request and generate a structured task description that includes the coordinates of the target point to be repaired, the quality and envelope size of the materials to be transported, and the task time window. Based on the coordinates of the target point to be repaired, obtain the elevation and meteorological data of the corresponding area and generate a three-dimensional environment model of the task area. Step 2: Spatial rasterization sampling is performed on the three-dimensional environment model to obtain a spatial discrete sampling field, and the target emergency repair structure sampling cluster and the spatial channel sampling cluster are identified to obtain the reference direction vector and the offset direction vector. Step 3: Based on the centroid of the target emergency repair structure sampling cluster, the reference direction vector and the offset direction vector, generate a three-dimensional spatial cone, and based on the envelope size of the material to be transported and the three-dimensional spatial cone, generate a three-dimensional feasible spatial network with cone risk weights. Step 4: Based on the three-dimensional feasible airspace network and the quality of the materials to be transported, candidate UAV formation configurations are selected, and based on the candidate UAV formation configurations and the three-dimensional feasible airspace network, a set of suitable formation configurations and safe operation boundaries corresponding to each path are generated. Step 5: Based on the three-dimensional feasible airspace network, mission time window and adaptive formation configuration set, generate mission execution plan and continuous motion path and formation configuration transition manifold of each UAV in the execution UAV set; Step 6: Based on the continuous motion path and formation configuration transition manifold of each UAV, generate the control command sequence of each UAV and send it to the execution UAV set to perform the collaborative transportation task.

[0012] In this embodiment of the invention, specialized identification of airspace channels captures local airflow characteristics around the target caused by structural obstruction and bypass flow, preventing UAVs from encountering sudden airflow impacts when approaching the target or hovering for docking, thus reducing flight control risks. By generating a three-dimensional feasible airspace network with cone-shaped risk weights, the material transportation needs are precisely matched with the airspace safety status, providing a quantifiable flight space basis for UAV formation planning, effectively avoiding obstacles in complex mountainous terrain and unstable airflow areas, and improving the safety and rationality of flight paths. Based on mission requirements and airspace characteristics, suitable UAV formation configurations are selected, and the safe operation boundaries under each path are clearly defined, ensuring that the formation's load capacity, anti-interference capability, and mission scenario are highly compatible, avoiding low transportation efficiency or safety hazards caused by improper formation selection, and improving the reliability of multi-aircraft collaborative transportation. By generating continuous motion paths and smooth formation configuration transition manifolds, the smooth switching of UAV formations in different flight stages is achieved, reducing the control difficulty caused by attitude changes. At the same time, combined with the issuance of control command sequences, the consistency and synchronization of multi-aircraft collaborative actions are ensured, improving the accuracy of material transportation and mission execution efficiency.

[0013] In a preferred embodiment of the present invention, step 1 includes: Step 1.1: Extract text and data fields from the original task request, and parse them to obtain the initial repair target point coordinates, initial material quality and initial envelope size information, and initial task time requirements. Specifically, the original task request may take various forms, such as electronic forms submitted by the ground station, text instructions issued by the emergency command center, or speech-to-text files reported by on-site personnel. To unify the subsequent parsing standards, these different types of requests first need to be formatted. For unstructured data such as voice and images, speech recognition technology is used to convert speech into text, and image recognition technology is used to extract... The text and numerical information in the images were ultimately converted into directly parsable text and data field formats. Then, a joint extraction algorithm combining keyword matching and entity recognition was used to accurately locate three types of core information from the standardized content: first, information related to the coordinates of the repair target point, focusing on extracting numerical combinations including longitude, latitude, and altitude; second, if the original information did not explicitly provide specific values, but only vague descriptions such as "near XX tower" or "XX line tower area," a pre-built power facility database was used to preliminarily determine the standard coordinates of the power facility by matching information such as facility name and location. This serves as the initial coordinates of the target repair point; secondly, it provides initial material mass and initial envelope size information. First, it identifies specific material names such as insulators, conductor joints, and surge arresters, then extracts the corresponding mass descriptions (e.g., 50 kg, 3 tons, two boxes, etc.). For quantity-based mass information like two boxes or three bundles, it needs to be matched against a pre-set material type-single box (single bundle) mass correspondence database to determine the mass of a single unit of that type of material. The initial material mass is then obtained by multiplying the quantity by the single unit mass. Simultaneously, it extracts dimensional descriptions such as length × width × height = XX × XX × XX, diameter XX, and height XX to form the initial envelope size information; thirdly, it provides the initial... The initial task time requirements are specifically extracted, including time descriptions such as immediate execution, delivery within 2 hours, completion before 6 PM today, and execution between 9 AM and 11 AM tomorrow. These clearly define the initial constraints (such as the earliest executable time) and the deadline constraints (such as the latest completion time), thus obtaining the initial task time requirements. After information extraction, the above three types of fields need to be validated. Information that is obviously invalid or contradictory, such as coordinate values ​​exceeding the geographical range of the target mountainous area or the same material being labeled as both 10 kg and 50 kg, needs to be removed. Missing key fields need to be supplemented (through database matching) to ensure the completeness and rationality of the preliminary analysis results.

[0014] Step 1.2: Convert the initial target point coordinates from the geodetic coordinate system to the local engineering coordinate system to obtain the target point coordinates; standardize and unify the units of the initial material mass and initial envelope size information to obtain the material mass and envelope size to be transported; perform time windowing processing on the initial task time requirement to obtain the task time window; integrate the target point coordinates, the material mass and envelope size to be transported, and the task time window to generate a structured task description, specifically including: The Gauss-Kruger projection algorithm is used to transform the geodetic coordinate system to the local engineering coordinate system, ensuring that the coordinates are adapted to the on-site UAV operation planning requirements. The specific operation is as follows: First, determine the origin of the local engineering coordinate system by selecting the UAV take-off and landing point or a nearby power facility with known precise coordinates (such as an intact tower) as the origin. Simultaneously, record the longitude (denoted as L0) and latitude (denoted as B0) data of this origin in the geodetic coordinate system. Next, calculate the positional difference between the initial repair target point and the origin, and use the initial repair... Subtracting the origin's longitude L0 from the longitude of the target point (L1) yields the longitude difference ΔL = L1 - L0; subtracting the origin's latitude B0 from the initial target point's latitude (B1) yields the latitude difference ΔB = B1 - B0; then substituting the longitude difference ΔL and latitude difference ΔB into the Gauss-Kruger projection formula, we calculate the plane rectangular coordinate increments (ΔX is the increment in the X direction, ΔY is the increment in the Y direction). The Earth's radius is taken as a fixed value R = 6378137 meters. The specific formula is as follows. Finally, the X coordinate of the emergency repair target point is obtained by adding the X coordinate (denoted as X0) of the origin in the local engineering coordinate system to the X-direction increment ΔX (X=X0+ΔX); the Y coordinate of the emergency repair target point is obtained by adding the Y-direction increment ΔY to the Y coordinate (denoted as Y0) of the origin in the local engineering coordinate system (Y=Y0+ΔY); the elevation value (denoted as H1) in the initial emergency repair target point coordinates is directly used as the H coordinate in the local engineering coordinate system (H=H1), and the three together constitute the emergency repair target point coordinates (X, Y, H) in the local engineering coordinate system.

[0015] For the initial material mass, first identify the unit type (e.g., ton, kilogram, jin), and uniformly convert it to kilograms. The conversion formula is: when the unit is ton, the mass of the material to be transported = initial material mass × 1000; when the unit is jin, the mass of the material to be transported = initial material mass × 0.5. If there are multiple types of materials, the total mass is calculated by summing, i.e., the mass of the material to be transported = the sum of the masses of all types of materials. For the initial envelope size, the unit is also standardized to meters. The cylindrical dimensions, such as diameter × height, are converted into the envelope size of the circumscribed cuboid (length = diameter, width = diameter, height = height), ultimately forming a standardized envelope size for the material to be transported. The initial task time requirements are converted into a clearly defined start and end time range. If the initial requirement is immediate execution, the start time is directly set to the current system time, and the end time needs to be calculated based on a pre-set mountain emergency task duration benchmark table. This table uses the weight of the materials to be transported and the distance between the repair target point and the take-off and landing point as core classification dimensions. Combining the complexity of mountainous terrain, the characteristics of drone swarm operations, and actual repair case data, it clarifies the standard operation duration values ​​for different scenarios. The standard operation duration covers drone swarm inspection and debugging, material loading and fixing, route planning, round-trip flight, hovering docking and delivery, and emergency redundancy time. The standard operating time is as follows: For materials weighing ≤10 kg, the standard operating time is 30 minutes if the distance is ≤3 km; 50 minutes if the distance is 3 to 8 km; and 90 minutes if the distance is >8 km. For materials weighing 10 to 50 kg, the standard operating time is 45 minutes if the distance is ≤3 km; 70 minutes if the distance is 3 to 8 km; and 120 minutes if the distance is >8 km. For materials weighing >50 kg, the standard operating time is 60 minutes if the distance is ≤3 km; 100 minutes if the distance is 3 to 8 km; and 150 minutes if the distance is >8 km. The current system time is used as the starting point for calculating the cutoff time. The standard operation time for the corresponding scenario is superimposed to obtain the latest completion time of the task; if the delivery time is within XX hours, the deadline = current system time + required time; if there is a specific time node requirement (such as completion before 12 o'clock tomorrow), the node is directly used as the deadline, and the start time is determined according to the drone formation preparation time (including equipment inspection and material loading), usually with a 30 to 60 minute allowance; after completing the above processing, the coordinates of the repair target point, the quality of the materials to be transported, the envelope size of the materials to be transported, and the task time window are integrated according to the preset template, the field names and data types of each parameter are clarified, and a structured task description is generated.

[0016] Step 1.3: Based on the coordinates of the target repair point in the structured task description, retrieve the digital elevation data and historical meteorological data set of the corresponding area from the pre-set geographic information database. Convert the digital elevation data into a regular grid terrain elevation matrix. Extract the prevailing wind direction and average wind speed data from the historical meteorological data set, and interpolate to generate a wind field vector matrix spatially aligned with the regular grid terrain elevation matrix. Specifically, extract the hourly prevailing wind direction and average wind speed data of the target area for the past three months from the historical meteorological data set. During this period, the historical prevailing wind direction is assumed to be mainly concentrated between 35 and 55 degrees northeast, with the most frequent wind direction being 45 degrees northeast. The historical average wind speed shows obvious time differences. The average wind speed during the day (6:00 AM to 6:00 PM) is 3.2 m / s, while the average wind speed during the night (6:00 PM to 6:00 AM the following day) is 1.8 m / s. The hourly data completeness throughout the statistical period exceeds 98%. If a real-time weather monitoring station is located within 5 kilometers of the target area, its real-time data is collected simultaneously as a supplement. For example, real-time monitoring might show the prevailing wind direction as 55 degrees northeast, with a real-time average wind speed of 4.5 m / s. To preserve the statistical representativeness of historical data while incorporating the dynamic timeliness of real-time data, a weighted average method is used to calculate the regional benchmark wind field parameters. The weighting rule is 70% for historical data and 30% for real-time data. Specifically, the two types of data are first matched according to the same hourly interval, and then calculated step by step... In the initial calculation, the baseline prevailing wind direction is calculated by multiplying the historical wind direction for each hour by 0.7 and then adding the real-time wind direction multiplied by 0.3 to obtain the weighted wind direction value for that hour. Then, the weighted wind direction values ​​for all hours over the past three months are averaged to obtain the baseline prevailing wind direction, assumed to be 48 degrees northeast. Special attention must be paid to the circular continuity of wind direction during the calculation; that is, the wind direction cycles 360 degrees, with both 0 degrees and 360 degrees representing true north. It is necessary to correct for the shortest arc length before substituting it into the weighted formula to avoid logical errors. The following uses a historical wind direction of 355 degrees and a real-time wind direction of 5 degrees as an example to explain the complete process in detail. First, the numerical difference between the two wind directions is calculated: |355 degrees - 5 degrees| = 350 degrees. Since the shortest arc length of the wind direction is... The angle between the two wind directions is 180 degrees (semicircle), and 350 degrees is far beyond this range. This means that direct calculation will misjudge the actual wind field relationship and must be corrected. Based on this judgment, the correction must follow the principle of making the difference between the two wind directions fall within 0 to 180 degrees. It is observed that 355 degrees is close to 360 degrees (essentially equivalent to 0 degrees), and 5 degrees is also close to 0 degrees. Therefore, 355 degrees is corrected to 360 degrees, which does not change the actual wind direction and can intuitively show the shortest arc length relationship. After the correction, the two wind directions become 360 ​​degrees (corrected historical wind direction) and 5 degrees (real-time wind direction). At this time, the difference is |360 degrees - 5 degrees| = 10 degrees, which meets the shortest arc length requirement. Next, we substitute it into the weighted formula for calculation. The weighted wind direction value = corrected historical wind direction × 0.7 + real-time wind direction × 0.3. After calculation, result validity processing is required. If the result exceeds 360 degrees, subtract 360 degrees; if it is negative, add 360 degrees to ensure the result is within the valid range of 0 to 360 degrees. This coherent process—first determining if correction is needed, then precisely adjusting the angle according to principles, substituting into the weighting formula for calculation, and finally verifying the validity of the result—avoids deviations caused by numerical misjudgments. When calculating the baseline average wind speed, the logic is consistent with the wind direction calculation. The weighted wind speed value for each hour is the historical wind speed of that period multiplied by 0.7, plus the real-time wind speed multiplied by 0.3. The arithmetic mean of the weighted wind speed values ​​for all periods is then taken to obtain the baseline average wind speed.

[0017] Through this calculation, the fused baseline wind field parameters retain both the regional wind field patterns reflected in historical data and the current wind field changes reflected in real-time data. However, the spatial resolution of meteorological data is typically only 1 to 5 kilometers per monitoring point, far lower than the 1 to 3-meter grid resolution of topographic data. Direct use of this data can lead to inaccurate matching between wind field information and topographic details. For example, the meteorological data for narrow canyons captured by the topographic grid is still a large-scale average wind field, which cannot reflect the wind speed acceleration effect within the canyon. Therefore, it is necessary to use bilinear interpolation to spatially refine the baseline wind field parameters and generate a wind field vector matrix that is perfectly aligned with the regular grid topographic elevation matrix. The specific interpolation operation is performed sequentially as follows: for each grid cell in the regular grid topographic elevation matrix (e.g., a 3-meter × 3-meter grid), the intersection of the two diagonals of the cell is taken as the center coordinate. This coordinate is completely consistent with the planar coordinate system of the topographic matrix, ensuring that the interpolation result accurately corresponds to the topographic details. The terrain grid, within a pre-set set of meteorological data points, uses the previously determined center coordinates as the core to filter out the four nearest known meteorological data points. The filtering process involves first considering the horizontal direction, finding the horizontal coordinates of the two nearest known points on either side of the center coordinates; then considering the vertical direction, finding the vertical coordinates of the two nearest known points on either side of the center coordinates. This results in a rectangular area formed by these four coordinates, with the four vertices representing the known meteorological data points used for interpolation. The wind field data for these four points is well-defined; for example, the wind direction at the top left vertex is 45 degrees, and the wind speed is 3.2 m / s. Following the principle that closer points have higher weights, the distance between the center coordinates and these four known points is first calculated using the planar distance formula. The distances of the four known points are then added together to obtain the total distance. Next, the weight of each point is calculated by subtracting the distance of a single point from the total distance and then dividing by the total distance. The resulting value is the interpolation weight. Finally, these weights are normalized (ensuring that the sum of the four weights equals 1), ultimately determining the weights of the four points.

[0018] Wind direction and wind speed are calculated separately. The interpolated wind direction is calculated by multiplying the wind direction of each of the four known points by its respective weight, and then summing the four products. The interpolated wind speed is calculated by multiplying the wind speed of each of the four known points by its respective weight and then summing the products. Combining the calculated interpolated wind direction and interpolated wind speed forms the wind field vector of the current terrain grid cell. This vector contains two core parameters: wind direction angle and wind speed magnitude, which can fully reflect the wind field characteristics of the grid location. By repeating the above operations of determining the target point, filtering known points, calculating weights, and interpolating for all grid cells in the regular grid terrain elevation matrix, a wind field vector matrix is ​​finally formed. The grid division method and spatial location of this matrix are completely consistent with the regular grid terrain elevation matrix.

[0019] Step 1.4 involves spatially overlaying and fusing the regular grid terrain elevation matrix and the wind field vector matrix to construct a 3D environment model of the task area containing both terrain elevation and wind field vector attributes. Specifically, this includes: using the unified row and column numbers of the regular grid as a spatial association benchmark, ensuring each grid cell is precisely matched with two sets of core data: terrain height information from the terrain elevation matrix and wind force and direction data from the wind field vector matrix, achieving a one-to-one correspondence between terrain and wind field information; in the 3D modeling stage, first, the horizontal position of each grid cell in space is locked using its planar coordinates (x-axis, y-axis), and then the terrain elevation corresponding to that grid cell is... Using the z-axis data as the data, a three-dimensional terrain structure is constructed through layer-by-layer rendering and 3D stitching to fully reproduce the real landforms such as mountain valleys, steep slopes, and gentle slopes. Then, the wind field vector data bound to each grid cell is accurately superimposed on the corresponding terrain location in the form of visual arrows. The direction of the arrows directly matches the wind direction, and the length or thickness of the arrows is proportional to the wind force, intuitively presenting the wind force differences in different areas. The final three-dimensional environment model of the task area not only fully preserves the undulating features of the mountain terrain, but also clearly reflects the wind force strength and wind direction distribution patterns at each spatial location, realizing the integrated visualization of terrain and wind field information.

[0020] This embodiment solves problems such as ambiguous coordinates, inconsistent units, and unclear time requirements in the original task request by extracting and verifying multi-source information and standardizing data, thus avoiding deviations in subsequent modeling and planning due to data errors. It specifically retrieves terrain and meteorological data and generates an alignment matrix to construct a multi-attribute 3D environment model, breaking the limitation of separating terrain and wind field data, and particularly highlighting the correlation between complex mountain terrain and wind fields. Coordinate transformation adopts a local engineering coordinate system suitable for mountain operations, material processing considers the diverse forms of power repair materials, and time windowing combines the urgency of emergency tasks, making the entire data processing process more in line with the actual needs of mountain power emergency scenarios, improving the practicality and operability of the method. The 3D environment model can present key information such as mountain terrain obstacles and wind field changes in advance, avoiding unpredictable terrain obstructions or airflow impacts encountered by drone formations during flight, thus improving the safety of collaborative transportation tasks from the source.

[0021] In a preferred embodiment of the present invention, step 2 includes: Step 2.1: Based on the preset spatial sampling resolution, the 3D space containing the 3D environment model is divided into a uniform 3D regular grid, with the geometric center of each grid cell serving as a sampling point. For each sampling point, the terrain elevation attribute value corresponding to the sampling point is extracted from the 3D environment model, and the mean magnitude of the wind field vector attribute at all locations within the grid cell containing the sampling point is calculated. The terrain elevation attribute value and the mean magnitude of the wind field vector are combined to form the sampling feature vector of the sampling point. The set of sampling feature vectors of all sampling points constitutes a spatial discrete sampling field, specifically including: First, the preset spatial sampling resolution is defined as 1 meter, which directly corresponds to the side length of the 3D grid cell. Each grid cell is a 1m x 1m x 1m cube. Based on this resolution, the complete 3D space containing the 3D environment model is uniformly divided into a regularly arranged 3D grid. The geometric center of each grid cell is the sampling point. The coordinates are calculated as follows: extract the maximum and minimum values ​​of the grid cell in the x, y, and z axes respectively, add them together, and divide by 2 to obtain the center coordinates of each axis. The combination of the three axis coordinates is the complete spatial coordinates of the sampling point (for example, for a grid with an x-axis range of 0 to 1 meter, the center x-coordinate is 0.5 meters). For each sampling point, the terrain elevation data corresponding to its location is extracted from the 3D environment model. This type of data reflects the undulations of the Earth's surface. The core parameters, measured in meters, directly correspond to the actual altitude or relative height of the sampling point. It's important to clarify that the core data supporting the 3D terrain structure in the 3D environment model is a pre-established regular grid terrain elevation matrix. The model creates the visualized terrain undulation effect by converting the elevation data of each grid in the matrix into z-axis coordinates in 3D space. Therefore, the model's terrain and elevation matrix are completely identical in spatial location and correspond one-to-one. During extraction, the planar coordinates (x-axis and y-axis coordinates) of the sampling points are used as the core retrieval basis. That is, the row number of the terrain elevation matrix is ​​directly associated with the grid range of the x-axis coordinate, and the column number is precisely matched with the grid range of the y-axis coordinate. Through this set of... Row and column indexing allows locating the grid cell in the matrix that perfectly overlaps with the spatial location of the sampling point. The value recorded in this cell is the terrain elevation data corresponding to the current sampling point. At the same time, the mean wind field vector magnitude of the grid cell containing the sampling point is calculated. First, the wind field vectors of all sub-sampling points in the grid are obtained (the magnitude of the wind field vector is equivalent to the wind speed). After summing all the magnitude values, the result is divided by the total number of sub-sampling points in the grid to obtain the mean wind field vector magnitude of the grid. The terrain elevation data of each sampling point is bound to the mean wind field vector magnitude of the corresponding grid to form a sampling feature vector containing two types of core information. The feature vectors of all sampling points are collected to finally form a spatial discrete sampling field.

[0022] Step 2.2: Threshold segmentation is performed based on the terrain elevation attribute values ​​of each sampling point in the spatial discrete sampling field. A first-class candidate point set with elevations higher than a preset elevation threshold is selected. Euclidean distance clustering is then performed on this first-class candidate point set to obtain multiple independent clusters. From these independent clusters, the cluster with the smallest spatial distance to the target repair point is selected and defined as the target repair structure sampling cluster. Specifically, this includes: performing threshold segmentation based on the terrain elevation data of the spatial discrete sampling field. The preset elevation threshold is 10 meters, which is set with reference to the typical height difference between common repair structures such as power towers and communication towers in mountainous areas and the surrounding terrain (mountainous repair structures are typically 50 to 80 meters high, forming a significant elevation difference with the surrounding terrain). Each point is compared... The elevation data and threshold of each sampling point are used to retain sampling points with elevations higher than 10 meters, forming the first type of candidate point set. Euclidean distance clustering is performed on the first type of candidate point set, with two preset core parameters: a distance threshold of 2 meters (twice the spatial sampling resolution to ensure that the clustering accuracy matches the sampling accuracy) and a minimum number of cluster points of 10 (to avoid a small number of isolated points forming invalid clusters). To accelerate the efficiency of neighborhood point retrieval, a spatial index of the first type of candidate point set is constructed using a KD-tree. The KD-tree is an efficient nearest neighbor query structure adapted to high-dimensional data. Compared with the traversal method (where the retrieval efficiency is linearly proportional to the total amount of data), the KD-tree can improve the neighborhood search efficiency to the logarithmic level, which is proportional to the logarithm of the total amount of data, thus shortening the retrieval time.

[0023] Traverse the unlabeled sampling points in the first candidate point set. Starting from the current point, search for all neighboring points with a spatial Euclidean distance of less than 2 meters using a KD-tree. If the number of neighboring points is ≥10, the current point and its neighboring points are grouped into an initial cluster. Then, recursively search for neighboring points of each point within the cluster and continuously add them until no new points can be included. If the number of neighboring points is less than 10, they are marked as isolated points and do not participate in clustering. Repeat this process until all points have been processed, resulting in multiple independent clusters that do not overlap. Calculate the cluster center coordinates of each independent cluster by summing the x, y, and z coordinates of all sampling points within the cluster and dividing by the total number of sampling points within the cluster to obtain the three-dimensional coordinates of the cluster center. Calculate the spatial Euclidean distance between each cluster center and the target repair point, and select the cluster with the smallest distance as the target repair structure sampling cluster.

[0024] Step 2.3: Based on the wind field vector magnitude attribute values ​​of each sampling point in the spatial discrete sampling field, threshold segmentation is performed to filter out a second type of candidate point set whose magnitude is lower than a preset wind resistance threshold. Based on the second type of candidate point set, points located in the spatial corridor between the target repair structure sampling cluster and the preset flight start point are selected to form a third type of candidate point set. Specifically, this includes: performing threshold segmentation based on the wind field vector magnitude data of the spatial discrete sampling field, with a preset wind resistance threshold of 6 m / s. This value is set with reference to the upper limit of safe flight wind speed for small and medium-sized operational UAVs. The magnitude data of each sampling point is compared with the threshold, and sampling points with a magnitude lower than 6 m / s are retained to form the second type of candidate point set; constructing the spatial corridor between the target repair structure sampling cluster and the preset flight start point. The spatial corridor is defined by taking the coordinates of the flight start point and the cluster center coordinates of the target repair structure sampling cluster as endpoints and determining a connecting line segment. The preset corridor width threshold is 20 meters. The spatial corridor is a cylindrical three-dimensional region with a radius of 10 meters and the line segment as the central axis (total width 20 meters), and its axial range is between the two endpoints. Each sampling point in the second type of candidate point set is traversed, and it is determined whether it is located in the spatial corridor by vector projection method. First, the vertical distance from the sampling point to the connecting line segment is calculated (by solving the coordinates of the perpendicular foot, and then calculating the straight distance between the sampling point and the perpendicular foot). If the distance is ≤20 meters and the sampling point is located between the two endpoints, it is determined to meet the conditions. All sampling points that meet the conditions are collected to form the third type of candidate point set.

[0025] Step 2.4: Perform density clustering on the third type of candidate point set, extract the connected regions with the highest data point density, and define the corresponding regions as spatial channel sampling clusters; perform principal component analysis on the spatial coordinates of all sampling points in the target repair structure sampling cluster to obtain the first principal component direction, and use the corresponding direction as the reference direction vector; perform principal component analysis on the spatial coordinates of all sampling points in the spatial channel sampling cluster to obtain the first principal component direction, and use the corresponding direction as the offset direction vector. Specifically, this includes: performing density clustering on the third type of candidate point set, with preset core parameters, namely, a neighborhood radius of 2 meters (consistent with the Euclidean distance clustering threshold to ensure parameter consistency), and a minimum number of core points of 6 (i.e., 2 × 10⁻⁶). In terms of data dimensions (3D data corresponds to a minimum of 6 core points), each sampling point is traversed. A sphere is drawn with that point as the center and a radius of 2 meters. The number of other sampling points within the sphere is counted. If the number is ≥6, it is a core point; if the number is <6 but it is located in the neighborhood of a core point, it is a boundary point; the rest are noise points. Starting from the core point, all density-reachable core points and boundary points are grouped into a connected region (density-reachable means that there is a series of core points between two points, and the distance between adjacent core points is <2 meters). This process is repeated until all core points are classified, resulting in multiple connected regions. The point density of each region is calculated, and the total number of sampling points in the region is divided by the spatial volume of the region. The connected region with the highest point density is selected and defined as the spatial channel sampling cluster.

[0026] Principal component analysis (PCA) is performed based on the coordinates of sampling points in the target emergency repair structure sampling cluster. The process is as follows: First, data centralization is performed by calculating the mean of the three-axis coordinates of all sampling points. The mean of the corresponding axis is subtracted from the coordinates of each point to eliminate positional bias. Second, a 3×3 covariance matrix is ​​calculated. The covariance of each pair of dimensions is calculated by multiplying the corresponding centralized data, summing the results, and then dividing by the total number of sampling points minus one, reflecting the correlation of data in each dimension. Third, eigenvalue decomposition is performed on the covariance matrix. The covariance matrix is ​​a 3×3 symmetric matrix (the element in the i-th row and j-th column is equal to the element in the j-th row and i-th column). The core of eigenvalue decomposition is to find a set of eigenvalues ​​that satisfy covariance matrix × eigenvector = eigenvalue × eigenvector (denoted as Cv = λv, where C is the 3×3 covariance matrix). The solution to the covariance matrix (where v is the eigenvector and λ is the corresponding eigenvalue) is calculated as follows: First, the core equation of the eigenvalue decomposition is defined. For the covariance matrix C, we need to find a non-zero vector v and a constant λ such that the formula Cv = λv holds. To solve this equation, it can be transformed into (C - λI)v = 0 (where I is a 3×3 identity matrix, i.e., a matrix with 1s on the main diagonal and 0s elsewhere). The necessary and sufficient condition for this homogeneous linear system of equations to have a non-zero solution v is that the determinant of the coefficient matrix is ​​0, i.e., det(C - λI) = 0. This equation is called the characteristic equation of the covariance matrix. Next, the eigenvalues ​​are calculated. By expanding the determinant det(C - λI), a cubic polynomial equation in λ can be obtained (the general form of a cubic equation is...). ), where a, b, c, and d are constants calculated from the original elements of the covariance matrix. Solving this cubic equation yields three roots, which are the three eigenvalues ​​of the covariance matrix, denoted as λ1, λ2, and λ3. Then, the corresponding eigenvectors are solved. Substituting each eigenvalue into the transformed equation (C-λI)v=0, three different homogeneous linear equation systems are obtained. Taking λ1 as an example, solving the equation system (C-λ1I)v=0 yields the non-zero solution vector, which is the eigenvector v1 corresponding to λ1. Similarly, v2 corresponding to λ2 and v3 corresponding to λ3 can be obtained. During the solution process, the obtained eigenvectors need to be normalized (i.e., by calculating the magnitude of the vector, dividing each element of the vector by the magnitude to make the vector magnitude 1). Ensure the eigenvectors are unit vectors; finally, through the above calculations, obtain the three eigenvalues ​​(λ1, λ2, λ3) and three independent unit eigenvectors (v1, v2, v3) corresponding to the covariance matrix. The magnitude of the eigenvalues ​​directly reflects the dispersion of the data (i.e., the variance) in the direction of the corresponding eigenvector, while the eigenvectors represent the main directions of the data distribution in three-dimensional space. The three together cover the core dimensions of the data distribution. In the fourth step, select the eigenvector corresponding to the largest eigenvalue, whose direction is the first principal component direction (the direction in which the data distribution is most concentrated), as the reference direction vector; use the same process to process the coordinate data of the spatial channel sampling cluster to obtain its first principal component direction, as the offset direction vector.

[0027] This embodiment, through uniform sampling and data fusion, constructs a spatial discrete sampling field that can capture the spatial distribution characteristics of terrain elevation and wind field modulus within the task area; through threshold segmentation and Euclidean distance clustering, it achieves precise positioning of the target repair structure area, effectively eliminating interference from irrelevant terrain areas; combined with wind resistance threshold and spatial corridor screening, it selects candidate points with low wind resistance and located within a reasonable flight path range, providing a safe and feasible basic range for UAV flight path planning; through density clustering, it determines the optimal airspace channel, and combined with principal component analysis, it obtains the key direction vector, providing a core reference for the accurate planning of UAV flight paths, ensuring wind resistance safety during flight, and efficiently reaching the target repair structure, thereby improving operational efficiency and mission success rate.

[0028] In a preferred embodiment of the present invention, step 3 includes: Step 3.1: Calculate the arithmetic mean of the spatial coordinates of all sampling points in the target emergency repair structure sampling cluster to obtain the centroid of the target emergency repair structure sampling cluster. Using the centroid of the target emergency repair structure sampling cluster as the vertex and the reference direction vector as the central axis direction, calculate the projection vector of the offset direction vector onto the plane perpendicular to the central axis. Calculate the angle between the projection vector and the central axis to obtain the spatial azimuth angle. Using the central axis as the rotation axis and the spatial azimuth angle as the cone half-angle, generate a three-dimensional spatial cone. Specifically, this includes: first, calculating the centroid of the target emergency repair structure sampling cluster. The centroid is the core representative of the spatial position of the sampling cluster. The calculation method is to sum the x-axis, y-axis, and z-axis coordinates of all sampling points in the target emergency repair structure sampling cluster, and then divide the sum of coordinates on each axis by the total number of sampling points in the sampling cluster to obtain the arithmetic mean of the x-axis, y-axis, and z-axis. The combination of these three averages is the target emergency repair structure. The centroid coordinates of the sampling cluster are determined. After centroid localization, the next step is to determine the opening characteristics of the spatial cone based on the previously obtained reference direction vector and offset direction vector, i.e., to calculate the projection vector. The specific operation consists of three steps: First, calculate the dot product of the offset direction vector and the reference direction vector. The dot product result is equal to the magnitude of the offset direction vector multiplied by the magnitude of the reference direction vector, and then multiplied by the cosine of the angle between the two vectors. Second, calculate the parallel component. Divide the above dot product result by the square of the magnitude of the reference direction vector to obtain a scaling factor. Then multiply this scaling factor by the reference direction vector to obtain the parallel component of the offset direction vector on the reference direction vector. Third, subtract this parallel component from the offset direction vector. The resulting difference vector is the projection vector of the offset direction vector on the plane perpendicular to the central axis. This vector is completely located in the vertical plane and only reflects the direction in the plane, directly determining the opening direction of the spatial cone.

[0029] With the projection vector, the opening angle of the cone needs to be quantified using the spatial azimuth. The calculation uses the projection of the central axis onto the vertical plane as a reference, and is solved using the vector angle formula. First, the dot product of the projection vector and the reference vector is calculated, then the magnitudes of the two vectors are calculated separately. The dot product result is divided by the product of the two magnitudes to obtain the cosine of the angle. The angle calculated using the inverse cosine function is the spatial azimuth, which is the half-angle of the cone and determines the size of the opening of the spatial cone. Finally, when constructing the three-dimensional spatial cone, the centroid of the sampling cluster is used as the vertex, the reference direction vector as the central axis, and the spatial azimuth as the half-angle. The cone extends along the axis away from the centroid, forming a cone with its opening facing the direction of the projection vector. Its range precisely covers the potential flight path from the starting point to the target area.

[0030] Step 3.2: The internal space of the 3D spatial cone is uniformly divided radially, circumferentially, and axially, discretized into a series of continuous voxel elements. For each voxel element, based on the envelope size of the transported material, an envelope space cuboid is constructed, aligned with the center point of the voxel element and with its edges aligned with the envelope size of the transported material. The minimum spatial distance between the outer surface of the envelope space cuboid and the cone surface of the 3D spatial cone is calculated and denoted as the first distance. The minimum spatial distance between the outer surface of the envelope space cuboid and the surfaces of all obstacles in the 3D environment model is calculated and denoted as the second distance. The minimum value between the first distance and the second distance is taken as the minimum directed distance of the voxel element, specifically including: First, the three-dimensional spatial cone is discretized. The core step is to determine the preset voxel resolution. Combining the typical envelope size of the transported material (0.8m × 0.5m × 0.5m) with the path planning accuracy requirements, the voxel resolution is set to 0.2m. This value accurately captures the material's placement details while controlling the computational load. After determining the resolution, the spatial cone needs to be uniformly divided along the three-dimensional direction. Radially, starting from the cone apex, it is divided towards the cone base at fixed intervals of 0.2m, ensuring the radial division density matches the voxel size. Circumferentially, it is divided around the central axis at fixed angles of 15 degrees, ensuring uniform voxel distribution and avoiding blind spots. Axially, it is divided along the central axis at fixed lengths of 0.2m, consistent with the voxel resolution. Through this three-dimensional division, the continuous spatial cone is broken down into a series of 0.2m × 0.2m × 0.2m cubic voxel units. The units are continuous and non-overlapping, and the spatial position of each unit can be uniquely determined by the three-dimensional coordinates of its center point. After completing the spatial cone discretization, in order to accurately assess the flight safety of the materials in each voxel unit, it is necessary to construct an envelope space cuboid for each voxel unit to simulate the space occupied by the materials. The envelope size of the materials to be transported is known to be 0.8m × 0.5m × 0.5m (length × width × height). During the operation, the coordinates of the center point of the voxel unit are first obtained through the spatial index of the voxel unit (e.g., the center point of a certain unit is (x0, y0, z0)). Then, a cuboid is built with this point as the geometric center. The center of the cuboid coincides with the center point of the voxel unit. The length, width, and height strictly match the envelope size of the materials, and the edge direction is also aligned with the actual placement direction of the materials. This ensures that the cuboid can restore the actual space occupation of the materials in the unit position in a 1:1 ratio, thus avoiding the distortion of the safety distance calculation due to size deviation from the source.

[0031] After constructing the envelope space cuboid of each voxel unit, two critical safety distances need to be calculated. Both use the spatial two-point distance formula (the sum of the squared differences of the coordinates of two points on the x, y, and z axes, and then the square root of the sum). The only difference is the object being compared. The first distance calculation focuses on the feature points on the outer surface of the envelope cuboid that are most likely to touch the cone surface. It traverses the 8 vertices of its outer surface and the midpoints of its 12 edges (a total of 20 feature points). At the same time, discrete points on the cone surface of the three-dimensional spatial domain are collected at 0.1-meter intervals. After calculating the distances between all cuboid feature points and cone surface sampling points, the minimum value is selected as the first distance. For example, if the first distance calculated for a certain unit is 0.3 meters, this means that there is a 0.3-meter safety redundancy between the materials in that unit and the boundary of the spatial cone. The second distance calculation uses the surface of obstacles (mountains, transmission towers, buildings, etc.) in the three-dimensional environment model as the reference. Similarly, for the object, traverse the 20 feature points on the outer surface of the envelope cuboid, and collect discrete points on the obstacle surface at 0.1-meter intervals. Calculate the distance between all point pairs using the same distance formula, and take the minimum value as the second distance. This value directly reflects the safety margin between the material and the obstacle. For example, if the calculated second distance is 0.5 meters, it means that the material has at least 0.5 meters of safe space between itself and the obstacle at this position. After calculating the two safe distances, compare the first distance and the second distance for each voxel unit, and take the smaller value as the minimum directed distance of that unit. Taking the above unit as an example, its first distance of 0.3 meters is less than its second distance of 0.5 meters, so the minimum directed distance is 0.3 meters. The magnitude of this value is directly related to the flight safety level. The smaller the value, the tighter the safety buffer space when the material flies at this position; the larger the value, the higher the safety redundancy.

[0032] Step 3.3: Based on the preset piecewise linear distance and risk mapping function, the minimum directed distance of each voxel unit is mapped to a risk weight. All voxel units and their corresponding risk weights together constitute a three-dimensional feasible spatial network with cone-shaped risk weights. Specifically, this includes: First, a piecewise linear distance and risk mapping function needs to be preset. This function is the core bridge for converting the minimum directed distance (safety distance index) of voxel units into quantifiable risk weights. Its rules use a safety threshold and a danger threshold as the core boundaries, dividing the minimum directed distance into three continuous intervals. Each interval corresponds to a clear risk judgment logic, and the danger threshold is set to 3 meters (materials and obstacles, spatial cone). The critical safety distance at the boundary (the risk of collision increases dramatically when the distance is less than this value) is defined as follows: the safety threshold is 8 meters (a sufficient safety buffer distance for cargo flights; there is no risk of collision when the distance is greater than this value); the risk weight ranges from 0 to 10 (0 represents absolute safety, and 10 represents extremely high collision risk). The specific mapping formulas and rules for the three intervals are as follows: When the minimum directed distance > 8 meters (safety threshold), the risk weight = 0 (fixed value), indicating that the area has sufficient safety redundancy; when 3 meters (danger threshold) ≤ minimum directed distance ≤ 8 meters (safety threshold), the linear formula risk weight = 16 - proportional coefficient 2 × minimum directed distance is used for calculation, where 2 is the core proportional coefficient of the linear function, and its value is determined by the interval... The weighting requirements at the endpoints are derived from the formula, specifying that 8 meters corresponds to a weight of 0 and 3 meters to a weight of 10. The slope of the linear function is (10-0) ÷ (3-8) = -2, where 2 is the absolute value of the slope. This directly reflects the correspondence between distance changes and weight changes. Within the 3-meter to 8-meter range, for every 1-meter decrease in the minimum directed distance, the risk weight increases linearly by 2. This ratio ensures weight matching at both endpoints and makes the gradual change in risk easier to quantify. The formula is derived based on the endpoints of the range, ensuring that the weight is 0 at 8 meters and 10 at 3 meters, achieving a gradual increase in risk weight with decreasing distance. When the minimum directed distance is < 3 meters (danger threshold), the risk weight... A value of 10 (fixed value) indicates that there is a direct collision risk in the area. This segmented design not only fits the actual safety logic of drone flight, but also simplifies risk calculation through linear relationships. After clarifying the mapping rules, a risk weight mapping operation is performed on each voxel unit. First, the interval to which the minimum directed distance of the unit belongs is determined by numerical comparison. Then, the calculation rules of the corresponding interval are applied. For example, if the minimum directed distance of a voxel unit is 5 meters, it is determined to be in the middle interval of 3 meters to 8 meters. Substituting 5 meters into the linear formula, we can get the risk weight = 16 - 2 × 5 = 6. This value intuitively reflects the risk level of the location. A score of 6 represents a medium to high risk, and a better path should be found to bypass it.After calculating the risk weights for all voxel units, a three-dimensional feasible airspace network is constructed. This involves arranging all voxel units in an ordered manner according to their three-dimensional spatial coordinates (x, y, z axis positions), with each unit firmly bound to its corresponding risk weight. This forms a structured dataset where spatial location and risk level correspond one-to-one. This network not only fully preserves the spatial distribution characteristics of the three-dimensional airspace cone, ensuring that the flight path remains within its bounds, but also clearly marks the safety level of each location through quantified risk weights, transforming the originally abstract concept of airspace safety into a concrete and verifiable numerical indicator.

[0033] This embodiment uses centroid calculation to pinpoint the core spatial location of the target repair structure. Combined with a three-dimensional airspace cone generated by direction vector projection and azimuth calculation, it achieves precise airspace delineation from the starting point to the target area, avoiding subsequent computational redundancy caused by an excessively large airspace range and ensuring the analysis focuses on the core feasible region. The airspace cone is discretized into voxel units, and the minimum directed distance is calculated based on the material envelope size. This considers the impact of the actual space occupied by the materials on the path and comprehensively quantifies safety redundancy through dual distance assessment, providing precise quantitative basis for risk assessment. By converting distance into risk weights through a piecewise linear mapping function, a weighted three-dimensional feasible airspace network is constructed, making flight risks visible and structured, effectively improving the safety and efficiency of path planning, and providing a clear target direction for path optimization.

[0034] In a preferred embodiment of the present invention, step 4 includes: Step 4.1: Based on the weight of the goods to be transported, select all formation configurations with a nominal total payload capacity greater than the weight of the goods from the pre-set formation configuration library to form an initial candidate formation configuration set. Specifically, the pre-set formation configuration library is constructed based on the performance parameters of typical multi-rotor UAVs and formation flight test data. It stores a variety of preset multi-rotor UAV formation configurations. Each configuration clearly indicates the nominal total payload capacity parameter calculated by the number of UAVs in the formation, the nominal payload capacity of a single UAV, and the load distribution coefficient based on safety redundancy. The load distribution coefficient is preset according to the formation size and anti-interference requirements. Generally, the larger the formation size, the higher the safety redundancy reserved to cope with the complexity of cooperative control and airflow disturbances, and the lower the corresponding load distribution coefficient value. For example, the library includes, but is not limited to, a 2-aircraft diamond configuration, consisting of two UAVs with a single payload of 20 kg. The configurations are as follows: A 4-unit rectangular configuration, consisting of four UAVs with a single payload of 20 kg, employs a redundant payload distribution scheme (distribution coefficient 0.8, with 20% reserved for anti-interference redundancy), with a nominal total payload capacity of 2 × 20 × 1.0 = 40 kg; a 6-unit hexagonal configuration, consisting of six UAVs with a single payload of 20 kg, employs a high-redundancy payload distribution scheme (distribution coefficient 0.7, with 30% reserved for anti-interference redundancy), with a nominal total payload capacity of 6 × 20 × 0.7 = 84 kg; and an 8-unit ring configuration, consisting of eight UAVs with a single payload of 20 kg, employs a super-strong redundant payload distribution scheme (distribution coefficient 0.6, with 40% reserved for anti-interference redundancy), with a nominal total payload capacity of 8 × 20 × 0.6 = 96 kg.

[0035] After clarifying the configuration library parameters, the initial screening stage begins. First, the actual weight of the materials to be transported is obtained. Combining this with the weight range of common emergency repair materials, it is assumed that the weight of the materials to be transported is 50 kg, and this value is used as the core screening benchmark. Then, all four formation configurations in the pre-set formation configuration library are traversed, and the nominal total load capacity of each configuration is compared with 50 kg one by one. The nominal total load capacity of the 2-aircraft diamond configuration is 40 kg, which is less than 50 kg and cannot meet the basic load requirement, so it is excluded. The nominal total load capacity of the 4-aircraft rectangular configuration (64 kg), the 6-aircraft hexagonal configuration (84 kg), and the 8-aircraft ring configuration (96 kg) are all greater than 50 kg, indicating that these configurations have the basic conditions to safely carry the materials, and they are included in the initial candidate formation configuration set. The final initial candidate formation configuration set only includes the three configurations that meet the load requirements: 4-aircraft rectangular, 6-aircraft hexagonal, and 8-aircraft ring.

[0036] Step 4.2: For each pre-planned candidate flight path in the three-dimensional feasible airspace network, perform path adaptability analysis to screen suitable formation configurations for the candidate flight paths and determine the safe operational boundaries of the suitable formation configurations on the candidate flight paths. The path adaptability analysis includes: calculating the interference force spindles at each waypoint on the candidate flight path based on the cone risk weights corresponding to all voxel units traversed by the candidate flight path and the wind field vector attributes in the three-dimensional environment model; for each formation configuration in the initial candidate formation configuration set... The system calculates the maximum set of anti-disturbance dynamic spindles for the formation configuration at each waypoint of the candidate flight path. It then performs an inclusion test between the interference force spindles at each waypoint and the maximum set of anti-disturbance dynamic spindles for the formation configuration at the corresponding waypoint to screen suitable formation configurations for the candidate flight path. For each screened suitable formation configuration, based on the minimum Euclidean distance from the boundary of the maximum anti-disturbance dynamic spindle set of the formation configuration at each waypoint of the candidate flight path to the interference force spindle, it determines the safe operating boundaries of the suitable formation configuration on the candidate flight path, specifically including: The comprehensive environmental interference at each waypoint is quantified, with interference force spinor being the core quantification indicator. It comprehensively reflects the intensity and form of interference, consisting of interference force and interference torque. The calculation process requires integrating the cone risk weights of voxel elements, the wind field vector attributes of the three-dimensional environmental model, and the inherent parameters of different formation configurations (including formation windward area, lever arm, etc.). Based on the initial candidate formation configurations determined in step 4.1, the key inherent parameters of each configuration have been derived through the characteristics of individual UAVs and the formation arrangement logic, namely, a single UAV fuselage length of 0.5 meters and a width of 0.4 meters, corresponding to a single fuselage windward area of ​​0.2 square meters and a maximum thrust of 50 Newtons. The formation consists of four aircraft arranged in a rectangular 2x2 configuration with a frontal area of ​​0.8 square meters and a lever arm (vertical distance from the formation's center of mass to the point of application of the interference force) of 0.3 meters; six aircraft arranged in a hexagonal configuration with a centrally symmetrical configuration with a frontal area of ​​1.0 square meter and a lever arm of 0.4 meters; and eight aircraft arranged in a ring configuration with a frontal area of ​​1.2 square meters and a lever arm of 0.5 meters. Basic data on the association between waypoints and formations are collected. Each candidate flight path consists of multiple consecutive waypoints, and each waypoint corresponds precisely to a voxel unit in the three-dimensional feasible airspace network. The cone risk weight of this unit (with a weight range of 0 to 10, and a maximum risk weight of 10, directly reflecting environmental interference) needs to be obtained first. The uncertainty of the disturbance); at the same time, the wind field vector attributes of the waypoint are extracted, including wind speed (e.g., a wind speed of 8 m / s at a waypoint determines the intensity of the wind force) and wind direction (determines the direction of the wind); finally, based on the current formation configuration being analyzed, its corresponding preset parameters are called (e.g., when analyzing a 4-aircraft rectangular configuration, 0.8 square meters of windward area and 0.3 meters of lever arm are called simultaneously) to form a complete set of calculation data for risk weights, wind field parameters, and formation parameters; with the basic data, the calculation of the corrected disturbance force can be carried out. Wind is the most important environmental disturbance source on the path. The first step is to calculate the direct force of the wind on the formation (basic disturbance force). The calculation formula is as follows: The basic disturbance force is calculated as follows: wind speed × air density × formation frontal area. The air density is taken as 1.225 kg / m³ under standard atmospheric conditions to ensure consistent calculation standards across different scenarios. The formation frontal area is determined based on the current analysis configuration. Considering that the risk weight of the cone reflects the uncertainty of environmental disturbances (the higher the weight, the stronger the additional disturbances from turbulence and obstacles), the basic disturbance force needs to be corrected using the weight. The correction logic is that the higher the risk weight, the greater the superposition ratio of disturbance intensity. The specific formula is: corrected disturbance force = basic disturbance force × (1 + the ratio of risk weight to the maximum risk weight), thereby quantifying the additional disturbances caused by environmental uncertainty.

[0037] After the corrected interference force is determined, the interference moment is further calculated. When the interference force acts on the formation, it will cause the formation attitude to deflect. The interference moment is the quantitative indicator of this deflection effect. Its calculation requires combining two core parameters: the corrected interference force and the lever arm. The lever arm has been determined by the formation structure design (0.3 meters for a 4-aircraft rectangular configuration, 0.4 meters for a 6-aircraft hexagonal configuration, and 0.5 meters for an 8-aircraft ring configuration). The calculation formula is: Interference moment = Corrected interference force × Lever arm. The direction of the moment is determined by the wind direction (the direction of the interference force) and the direction of the lever arm, which fully reflects the impact of interference on the formation attitude. Finally, the corrected interference force and interference moment of each waypoint are used as two core components to form the interference force spin of that waypoint. This indicator includes both the straight-line thrust interference of the wind on the formation and the deflection moment interference on the attitude, which can comprehensively and accurately describe the comprehensive effect of the environment on the current formation configuration being analyzed.

[0038] Step 4.3: Summarize all suitable formation configurations and their safe operating boundaries for each candidate flight path, generating a set of suitable formation configurations and safe operating boundaries for each candidate flight path. Specifically, this includes: After clarifying the quantitative standards for environmental interference, the next core task is to screen suitable formation configurations that can withstand this interference. The suitability judgment is essentially a matching verification between the formation's maximum anti-interference capability and the path interference; that is, the formation's anti-interference limit must completely cover the comprehensive interference of all waypoints on the path. This is achieved through two steps: calculating the formation's anti-interference limit and performing a full-path interference coverage verification. To ensure the reliability of the screening results, the first step is to calculate the maximum interference force and upper limit of the interference torque that the initial candidate 4-drone rectangular, 6-drone hexagonal, and 8-drone ring configurations can withstand, based on their power system and structural design parameters. This will form the maximum anti-interference power spinor set, i.e., the anti-interference capability boundary. When calculating the maximum anti-interference force, the anti-interference force of the formation originates from the thrust coordination of each UAV. The formula is: Maximum anti-interference force = Maximum thrust of a single UAV × Number of UAVs in the formation × Power redundancy coefficient. For example, the maximum thrust of a single UAV is fixed at 50 Newtons, and the power availability coefficient (e.g., 0.8 for the 4-drone configuration, 0.7 for the 6-drone configuration, and 0.6 for the 8-drone configuration) is a factor less than 0.8. The coefficient of 1 is obtained through hovering and anti-disturbance tests on the formation configuration. Specifically, the formation is allowed to hover stably in a windless environment, and the average percentage of each UAV's thrust relative to its maximum thrust is recorded. Then, a known horizontal disturbance force is applied, gradually increasing until the formation begins to drift uncontrollably. The percentage of additional thrust used to resist the disturbance at this critical state is recorded. The available power coefficient is the sum of the hovering thrust percentage and the available additional anti-disturbance thrust percentage. It comprehensively reflects the thrust share necessary to maintain basic flight and attitude stability; the remaining portion can be used to counter sudden disturbances along the flight path. When calculating the maximum anti-disturbance moment, it is considered in conjunction with the formation... The anti-interference torque is related to the moment of inertia and attitude control capability. The formula is: Maximum anti-interference torque = Formation moment of inertia × Maximum allowable angular acceleration. The formation moment of inertia is the moment of inertia of the formation about its center of mass axis. It is calculated based on the mass of each UAV constituting the formation and its relative position in the configuration. During the calculation, the mass of each UAV (including its own mass and the shared load mass) is regarded as concentrated at its own center of mass, and its position coordinates in the formation coordinate system are obtained. The overall moment of inertia of the formation (e.g., the yaw moment of inertia about the vertical axis) is obtained by summing the contributions of each UAV to the moment of inertia of that axis. The specific calculation formula is: Formation moment of inertia = In the formula, This represents the summation over all participating UAVs in the formation; the horizontal distance to the formation's center of mass is calculated based on the UAV's coordinates and the overall formation's center of mass coordinates. Using this formula, characteristic values ​​for different configurations can be calculated, for example, approximately 0.8 kg·m² for a 4-UAV rectangular configuration, approximately 1.5 kg·m² for a 6-UAV hexagonal configuration, and approximately 2.2 kg·m² for an 8-UAV ring configuration; the maximum permissible angular acceleration is a preset safety upper limit, obtained through standard anti-interference flight tests on UAV formations. In a simulated mountainous wind field environment, a series of gusts of known intensity are applied to the flying formation, and its attitude response is monitored; dividing this interference torque value by the formation's moment of inertia yields the critical angular acceleration that the formation can overcome without instability; averaging the critical values ​​obtained from multiple tests and multiplying them by an engineering safety factor less than 1 (e.g., 0.8) finally determines this upper limit value (e.g., 5 radians / second²), ensuring a rapid and smooth attitude recovery process.

[0039] Using the maximum anti-interference force and maximum anti-interference moment of each configuration as the horizontal and vertical boundaries, a closed region is constructed in two-dimensional space, namely the maximum anti-interference dynamic spinor set of that configuration, clearly marking the upper limit range of its anti-interference capability. The second step involves verifying, for each formation configuration, whether the interference force spinor of all waypoints on the candidate flight path falls within its maximum anti-interference dynamic spinor set. The verification strictly follows the principle of dual component compliance to ensure that all interference is covered throughout the entire flight path. Specifically, during single waypoint verification, for a given waypoint, if the two components of its interference force spinor—that is, the corrected interference force ≤ the formation's maximum anti-interference force and the interference moment ≤ the formation's maximum anti-interference moment—then the interference at that waypoint is covered by the formation's anti-interference capability. For example... For example, if the interference force spin of a waypoint to a 4-aircraft rectangular configuration is (11.76 N, 3.528 N·m), and compared with its maximum anti-interference spin set (240 N, 4 N·m), both components are within the limits, indicating that the interference at that waypoint can be resisted. During the full-path test, single waypoint coverage does not guarantee full-path safety. Only when the interference force spin of all waypoints meets both criteria can the formation configuration resist interference throughout the entire path and be judged as a suitable formation configuration. If the interference force or interference moment of any waypoint exceeds the upper limit, such as a waypoint interference moment of 4.2 N·m, which exceeds the upper limit of 4 N·m for a 4-aircraft rectangular configuration, then the configuration cannot guarantee flight safety throughout the entire path and is directly excluded.

[0040] After selecting a suitable formation configuration, it is necessary to further clarify its safe operating boundaries, i.e., the limits of power output and attitude adjustment when the formation flies along this path. The core purpose is to ensure that the formation's anti-interference capability always has sufficient redundancy to avoid instability under extreme conditions. The boundary determination is based on the safety redundancy of anti-interference capability and actual interference, and is achieved by quantifying the redundancy of single waypoints and integrating the boundaries of the entire path to ensure that the boundary rules are specific and executable. The first step is to regard the interference force spinor of each waypoint as a point in two-dimensional space (the horizontal axis is the interference force, and the vertical axis is the interference torque). The boundary of the maximum anti-interference power spinor set of the formation is regarded as a closed rectangle in two-dimensional space (the length of the rectangle is the maximum anti-interference force, and the width is the maximum anti-interference torque). The core of the safety redundancy is the distance from the interference point to the anti-interference boundary. The minimum value of this distance is calculated using the Euclidean distance formula, which is: Minimum Euclidean distance = The interference force difference is calculated as follows: maximum anti-interference force - actual interference force. The interference torque difference is calculated as: maximum anti-interference torque - actual interference torque. Taking a 4-aircraft rectangular configuration at a certain waypoint as an example, the interference force difference is 240 - 11.76 = 228.24 N·m, and the torque difference is 4 - 3.528 = 0.472 N·m. Substituting these values ​​into the formula, the minimum Euclidean distance is assumed to be 228.24 N·m. The larger this value, the more sufficient the safety redundancy of the formation's anti-interference capability relative to the actual interference. The second step is to determine the safe operating boundary based on the redundancy distribution along the entire path. The redundancy of a single waypoint cannot reflect the needs of the entire journey. It is necessary to comprehensively consider the minimum Euclidean distances of all waypoints along the entire path and formulate differentiated operating boundary rules to ensure that the boundary and redundancy distribution are accurately matched. In areas with sufficient redundancy (minimum Euclidean distance > 100 N·m), the formation has sufficient anti-interference margin and can operate safely. The boundaries can be appropriately relaxed, allowing the formation to adjust its attitude within a relatively large range (e.g., attitude angle fluctuation ±3 degrees) and flexibly adjust power output (e.g., thrust fluctuation ±10%), balancing flight efficiency and safety. In the redundancy-tight region (minimum Euclidean distance < 50 N·m), the formation's anti-disturbance margin is limited, and the operational boundaries must be strictly tightened, limiting the attitude adjustment range (e.g., attitude angle fluctuation ≤ ±1 degree) and strictly controlling power output fluctuation (≤ ±5%) to avoid exhausting redundancy due to excessive operational range. Finally, boundary integration is completed, and the operational restrictions of each waypoint are connected in sequence according to the path, forming a continuous correspondence between waypoint position and allowable operational range. This integrated safe operational boundary ensures that when the formation flies at any position on the path, the power output and attitude adjustment will not push the anti-disturbance capability close to the upper limit, always maintaining safe redundancy, and providing clear operational specifications for actual flight control.

[0041] This embodiment, through load capacity screening, quickly eliminates formation configurations that cannot meet basic transportation requirements, reducing the scope of subsequent analysis and improving the efficiency of the overall process; by quantifying the matching relationship between environmental interference and formation anti-interference capability, it achieves accurate screening of suitable formation configurations, and at the same time, based on the safety redundancy-determined safety operation boundary, it provides clear safety constraints for formation flight, effectively ensuring the stability and safety of the flight process; by summarizing and forming structured path, formation, and boundary correspondence data, it avoids blind decision-making and ensures that the final selected flight path and formation configuration can adapt to complex environmental conditions and transportation needs.

[0042] In a preferred embodiment of the present invention, step 5 includes: Step 5.1: Based on the three-dimensional feasible airspace network, the mission time window, and the set of suitable formation configurations and safe operation boundaries corresponding to each of the candidate flight paths, multi-objective programming is performed to generate a mission execution plan with the goal of shortening the total flight time and improving the average safe operation boundary. The mission execution plan includes the selected set of execution UAVs, the selected flight trajectory from the preset starting point to the repair target point and traversing the three-dimensional feasible airspace network, and the formation configuration specified for each waypoint on the selected flight trajectory. Specifically, it includes: first, clarifying the objective function and quantification method, transforming the abstract efficiency and safety objectives into calculable mathematical indicators. Among them, the quantification of the total flight time is based on the path length and formation speed. The calculation formula is: total flight time = total spatial length of candidate flight paths ÷ average flight speed of the corresponding formation configuration. The specific logic and parameters of the formation flight speed are as follows: the maximum level flight speed of a single UAV is determined by its power system (e.g., the maximum level flight speed of a single UAV is 20 m / s, which is pre-stored in the UAV performance parameter library). However, when flying in formation, the speed needs to be adjusted in combination with power redundancy retention and anti-interference capability matching. The core calculation logic is as follows. The formula is as follows: Formation flight speed = (maximum thrust of a single UAV × number of UAVs in the formation × power redundancy coefficient × thrust efficiency coefficient) ÷ (formation air drag coefficient × air density × formation frontal area). In this formula, the maximum thrust of a single UAV and the thrust efficiency coefficient (reflecting the energy conversion efficiency of the power system, with values ​​ranging from 0.85 to 0.95, provided by the UAV manufacturer) are fixed performance parameters. The number of UAVs in the formation, the power redundancy coefficient (previously stated as 1.2 for 4 UAVs, 1.3 for 6 UAVs, and 1.4 for 8 UAVs), and the formation frontal area (0.8 m² for 4 UAVs, 1.0 m² for 6 UAVs) are also considered. The values ​​of 1.2 m² (8 aircraft 1.2 m²) are derived from the inherent parameters of the adapted formation configuration; the formation air drag coefficient (related to the configuration shape, 0.9 for rectangle, 0.7 for hexagon, and 0.6 for ring) is calibrated using wind tunnel test data; the air density is taken as 1.225 kg / m³ under standard atmospheric conditions; the quantification of the average safe operating boundary is achieved by statistical path redundancy, and the formula is: average safe operating boundary = sum of safe operating boundary values ​​of all waypoints on the selected path ÷ total number of waypoints. The larger this value is, the more sufficient the formation operating redundancy is, and the higher the flight safety is.

[0043] Next, three core constraints were determined to ensure that the planning results met the actual execution requirements. The airspace constraint required that the selected flight trajectory must be entirely within the three-dimensional feasible airspace network and must not exceed the safe range formed by voxel units. The time constraint stipulated that the total flight time must not exceed the mission time window; for example, if supplies need to be delivered within 4 hours, the total planned time must strictly be ≤4 hours. The configuration constraint stipulated that the formation configuration specified at each waypoint must belong to the set of suitable formation configurations corresponding to that path, and the total number of configuration changes throughout the journey must be ≤3 to avoid frequent changes leading to formation attitude instability. Finally, the dual-objective problem was transformed into a single-objective optimization problem using a weighted summation method, and solved using a genetic algorithm, with weights set according to task priority. Since different tasks have different priorities regarding time and safety, the weighting directly reflects these priorities. Taking an emergency repair mission as an example, materials need to be delivered as quickly as possible to minimize losses. Therefore, shortening the total flight time takes precedence over increasing the average safe operating boundary. The weight of the total flight time is set at 0.6, and the weight of the average safe operating boundary is set at 0.4. This weighting ensures both the priority of efficiency and the safety redundancy is not neglected through the 0.4 weight, avoiding sacrificing flight safety for speed. The average safe operating boundary is normalized. The core of normalization is to convert the actual value of the average safe operating boundary into a dimensionless value between 0 and 1 using the formula (actual target value - minimum target value) ÷ (maximum target value - minimum target value). The minimum and maximum target values ​​are derived from the statistical results of the average safe operating boundaries of all candidate paths (e.g., the minimum safety boundary is 50 N·m and the maximum is 200 N·m across all paths). (The actual value of a certain path is 150 N·m). After normalization, the value of the average safe operating boundary is unified to the same scale as the normalized value of the total flight time, ensuring the fairness of the weighted calculation. Then, by processing the normalized value of the average safe operating boundary (1 - average safe operating boundary), the target direction is unified. Since the average safe operating boundary needs to be maximized, the closer its normalized value is to 1, the more sufficient the safety redundancy. The normalized value of 1 will decrease as the safety boundary increases, which just transforms the goal of maximizing the safety boundary into the goal of minimizing this part of the value, which is consistent with the goal of minimizing the total flight time. Based on this, the final single-objective optimization function is: total optimization goal = (0.6 × total flight time) + (0.4 × (1 - normalized value of average safe operating boundary)). When the total optimization goal of a certain combination is minimized, it means that it has reached the optimal state of mission requirements in terms of efficiency and safety, and can be directly adapted to the minimization solution logic of the genetic algorithm.

[0044] The objective function was clearly defined, and a genetic algorithm was chosen for solving it. The core reason for this choice is that the mission planning of UAV formations is essentially a dual selection of the optimal path and the optimal configuration combination. The genetic algorithm, by simulating the selection, crossover, and mutation processes of biological evolution, can efficiently traverse the combination space and lock in the optimal solution under multiple constraints such as airspace, time, and configuration. The candidate flight path and the suitable formation configuration combination are encapsulated into the basic computational unit of the algorithm, namely chromosomes. Each chromosome corresponds to a complete preliminary planning scheme, and its encoding logic is completely aligned with the actual needs of the mission. That is, the first half of the chromosome is the unique identifier of the candidate flight path (for example, path 3 corresponds to three-dimensional feasibility). The first part is a continuous sequence of voxel units in the airspace from the starting point to the target point. The second part is the adaptive formation configuration code for each waypoint on the path (pre-defined: 01 represents a 4-aircraft rectangle, 02 represents a 6-aircraft hexagon, and 03 represents an 8-aircraft ring). All configuration codes are strictly selected from the adaptive configuration set of the path to ensure that the configuration switching number is ≤3 times. For example, a chromosome can be interpreted as path 3 + waypoints 1-5 using configuration 01, waypoints 6-10 using configuration 02. From path selection to configuration allocation, it is clear at a glance, forming a directly verifiable planning scheme. The entire solution process has an iteration limit of 50 generations (to ensure solution accuracy and avoid computational redundancy). The evolutionary process (EC) progresses through three core operations: selection, crossover, and mutation. Each step builds upon the results of the previous step to generate better combinations, forming a clear evolutionary logic. Selection is the fundamental screening operation, its core being the preservation of high-quality solutions and the transmission of advantages. Before the selection operation, the objective function value of all chromosomes in the current generation (i.e., all preliminary planned solutions) is calculated. The smaller the value, the better the solution balances efficiency and safety. Then, the chromosomes are sorted from smallest to largest by objective function value. A combination strategy of elite preservation and roulette wheel selection is used to directly send the top 20% of the best chromosomes (elite individuals) to the next generation, ensuring that core advantages are not lost. The remaining 80% of chromosomes are selected using the roulette wheel selection method. The specific process is as follows: Since the objective function is to minimize, the objective function value is first converted into fitness (to avoid negative selection probability). The formula is: Fitness of a chromosome = Maximum objective function value of the remaining 80% of chromosomes in the current generation - Objective function value of the chromosome. The better (smaller) the objective function value, the greater the fitness. Selection probability of a chromosome = Fitness of the chromosome ÷ Sum of fitness of the remaining 80% of chromosomes in the current generation. The sum of selection probabilities of all chromosomes is 1. The roulette wheel is divided into sector areas according to the selection probability of each chromosome. The higher the probability, the larger the sector area. The roulette wheel is randomly rotated. The chromosome corresponding to the sector area pointed to by the pointer is the individual selected to enter the next generation.

[0045] After selecting high-quality chromosomes through selection, the next step is to achieve synergistic effects through crossover operations, generating a more competitive new scheme. During this operation, the selected high-quality chromosomes are first randomly paired. Then, a crossover point is determined among the continuous waypoints in the path (e.g., between waypoints 5 and 6). The paired chromosomes exchange the configurational coding segments after the crossover point. Taking paired chromosomes A and B as an example, chromosome A uses a 01 configuration for path 3 + waypoints 1-5 and a 02 configuration for 6-10 (its advantage is a shorter path distance in the first half). Chromosome B uses a path 4 + waypoints... Points 1-5 use a 0-2 configuration, and points 6-10 use a 0-1 configuration (the advantage is high safety redundancy in the latter half). If the intersection point is set as waypoint 5, two new chromosomes will be generated after the intersection: path 3+1-5 uses a 0-1 configuration, path 6-10 uses a 0-1 configuration, path 4+1-5 uses a 0-2 configuration, and path 6-10 uses a 0-2 configuration. This retains the path advantages of the original chromosomes while incorporating the configuration advantages of the other, thus improving combinatorial performance. After the crossover operation generates a new scheme, to avoid the algorithm getting trapped in a local optimum (i.e., a certain scheme is optimal in the current generation but not globally optimal), further optimization is needed. Furthermore, a mutation operation is required to introduce moderate random changes. During the operation, a portion of the new chromosome is randomly selected with a low probability of 1%. Then, 1-2 mutation bits (i.e., the configuration code of a certain waypoint) are randomly determined in the configuration coding segment of the chromosome and replaced with other codes in the set of configurations that are suitable for that path. For example, the 01 configuration of waypoint 5 is replaced with the 03 configuration that is suitable for the same path. The key is that after mutation, it is necessary to verify whether the total number of configuration switching times is ≤3 times to ensure that the new scheme still meets the constraints. This operation not only retains the main advantages of the scheme, but also opens up space for exploring better configuration combinations. After 50 iterations, the algorithm will output the chromosome with the smallest objective function value. The planning scheme corresponding to it is the global optimal solution, that is, the final task execution scheme. This scheme clearly includes the set of execution drones (determined by the number of drones with the selected configuration, such as 4 numbered drones corresponding to the 01 configuration), the selected flight trajectory from the starting point to the repair target point (i.e., the specific path corresponding to the first half of the chromosome), and the specified formation configuration of each waypoint on the trajectory (the actual configuration corresponding to the second half of the chromosome code).

[0046] Step 5.2: Obtain consecutive waypoints A and B on the selected flight trajectory. Obtain the 3D position coordinate set A for each UAV in the execution UAV set under the specified formation configuration at waypoint A, and the 3D position coordinate set B for each UAV in the execution UAV set under the specified formation configuration at waypoint B. Specifically, this includes: first, locking consecutive waypoint pairs from the mission execution plan, determining adjacent waypoints A and B according to the flight sequence. For example, the third waypoint in the flight trajectory is A, and the following fourth waypoint is B. Simultaneously, record the specified formation configuration for each point, such as a 4-UAV rectangular configuration at point A and a 6-UAV hexagonal configuration at point B. The core of this step is to clarify the position-configuration correspondence, establishing an index for subsequent coordinate retrieval. Then, retrieve the position coordinates based on the execution UAV set, which is determined by the mission plan. For example, for the four UAVs numbered U1 to U4 corresponding to a 4-aircraft rectangular configuration, two sets of data are retrieved from the preset formation configuration-position coordinate library. One set is set A, which is the three-dimensional position coordinates of U1 to U4 under the 4-aircraft rectangular configuration at waypoint A. The coordinates are based on a unified coordinate system of the three-dimensional feasible airspace. For example, U1 is (100, 200, 50), and the unit is meters. The other set is set B, which is the three-dimensional position coordinates of U1 to U4 under the 6-aircraft hexagonal configuration at waypoint B. If the configuration is expanded (e.g., from 4 aircraft to 6 aircraft), the coordinates of the newly added UAVs (e.g., U5, U6) are added. For example, the coordinates of U1 at point B are (105, 203, 52). The data in the coordinate library are all pre-calculated by the formation arrangement logic. For example, in the 4-aircraft rectangular configuration, U1 is located in the upper left corner. The coordinates of all UAVs are symmetrical with respect to the configuration center to ensure that the positional relationship meets the requirements of formation flight.

[0047] Step 5.3: For each drone in the drone set, project the corresponding 3D position coordinates of the drone in set A and the 3D position coordinates in set B onto a virtual unit sphere to obtain the spherical projection point pair corresponding to the drone. Specifically, this involves transforming the 3D positional relationship into a spherical geometry problem by projecting the drone position coordinates onto the unit sphere. The specific implementation process is as follows: First, define the unit sphere. Construct a virtual sphere with the origin of the 3D coordinate system as the center and a fixed radius of 1. The projection points of all drone coordinates will fall on this sphere, thereby eliminating distance scale differences. Next, perform projection calculations on the coordinates of each drone. The processing objects are the 3D coordinates of the corresponding drones in sets A and B obtained in step 5.2. The projection calculation is divided into two steps. The first step is to calculate the coordinate modulus. Let the 3D coordinates of the drone be (x, y, z), and the modulus r be the distance from the coordinate point to the origin. The calculation formula is r = The second step is to calculate the coordinates of the projection point. Divide the three components of the original coordinates by the modulus r, so that the coordinates of the unit spherical projection point are (x / r, y / r, z / r). Taking U1 as an example, its coordinates in set A are (100, 200, 50). First calculate the modulus, and then calculate the projection point. Similarly, calculate the projection point of U1 in set B to form the spherical projection point pair of the UAV. Process the coordinates of all UAVs in the UAV set in this way to obtain the spherical projection point pair of each UAV.

[0048] Step 5.4: On the unit sphere, calculate the shortest great circle path connecting the corresponding UAV spherical projection point pairs. Based on the shortest great circle path, backproject it back into three-dimensional space to construct a smooth spatial transition curve for each UAV in the execution UAV set, from the configuration position of waypoint A to the configuration position of waypoint B. The set of smooth spatial transition curves of all UAVs in the execution UAV set collectively represents the formation configuration transition manifold. Specifically, this includes: First, calculating the shortest great circle path on the unit sphere. The shortest path between two points on the unit sphere is the plane passing through the center of the sphere and the sphere surface. For the great circle arc formed by the intersection, in specific calculations, first assume the spherical projection points of a certain UAV are PA and PB. The angle θ between the two points is calculated using the dot product formula: θ = arccos(PA·PB). Since the projection points are unit vectors, the dot product result is directly equal to cosθ. Then, spherical linear interpolation (SLERP) is used to obtain continuous points on the great circle arc. The interpolation formula is P(t) = [sin((1-t)θ) / sinθ]×PA + [sin(tθ) / sinθ]×PB, where t is the interpolation value representing the transition progress from projection point PA to PB. The parameter t ranges from [0, 1]. t=0 corresponds to PA, and t=1 corresponds to PB. To ensure a smooth path, t is set to intervals of 0.1, 0.2…1.0, resulting in 10 interpolation points. The second step involves back-projecting the spherical interpolation points into three-dimensional space. The back-projection logic is to maintain the direction and restore the distance scale. First, the distance scale at different times is determined. Let the modulus of the UAV at point A be rA, and the modulus at point B be rB. The modulus r(t) at time t is calculated using a linear transition, with the formula r(t) = rA + t × (rB - r A) Ensure a smooth change in distance from A to B; then multiply the spherical interpolation point P(t) at time t by the corresponding modulus r(t) to obtain the three-dimensional coordinates of the UAV at that time; the third step is to construct the formation configuration transition manifold, connecting the back-projected coordinates of each UAV in t-value order to form a smooth spatial transition curve; the set of all UAV transition curves together constitutes the transition manifold of the formation configuration from A to B. This manifold ensures both the stability of the movement of a single UAV and the absence of collision risk when the overall formation configuration changes, clearly presenting the switching rules of the formation configuration.

[0049] Step 5.5: Based on each smooth spatial transition curve in the formation configuration transition manifold, generate a continuous motion path from waypoint A configuration to waypoint B configuration for the corresponding UAV in the execution UAV set. Specifically, this includes: first, determining the kinematic constraint parameters by retrieving core physical constraint indicators from the UAV parameter library, including maximum flight speed vmax (e.g., 15 m / s), maximum acceleration amax (e.g., 2 m / s²), and maximum angular velocity ωmax (e.g., 0.5 radians / s). These parameters are the core criteria for determining whether the path is executable; then, discretizing the transition curve in time and calculating a reasonable time step Δt based on the maximum acceleration, using the formula Δt = Where Δs is the straight-line distance between adjacent interpolation points of the transition curve, for example, when Δs = 1 meter, Δt = =1 second, according to this Δt, the interpolation process of t∈[0,1] is discretized into 10 time nodes (t0=0, t1=0.1, ..., t10=1.0), so that the theoretical curve is transformed into discrete time-position nodes; finally, a continuous motion path is generated. For each UAV, the corresponding three-dimensional coordinates are extracted sequentially according to the time nodes to form a preliminary time-position correspondence; the speed of each segment of the path is further verified. The speed v=Δs / Δt. If v exceeds vmax, then Δt is reduced until v≤vmax. For example, when Δs=1 meter and vmax=1 meter / second, Δt is adjusted to 1 second. Finally, the continuous motion path of each UAV output has its position, speed and acceleration at each time node strictly in accordance with kinematic constraints, ensuring that the formation can smoothly and safely complete the configuration switch and flight mission from waypoint A to B.

[0050] This embodiment, by quantifying objectives and defining constraints, avoids the problem of simply pursuing speed while neglecting safety. The generated mission execution plan not only meets the time window requirements but also ensures that the formation remains in a suitable and safe configuration throughout the entire process. It converts the complex positions in three-dimensional space into a standardized projection of a unit sphere, providing a simplified model for shortest path calculation. By inferring the three-dimensional space curve through the spherical shortest path, the UAVs transition smoothly from one configuration to another without abrupt changes in trajectory, effectively avoiding collisions within the formation and reducing the difficulty of attitude control. Through discretization and constraint verification, the theoretical curves are transformed into motion commands that conform to the physical limits of the UAVs, avoiding flight instability caused by exceeding speed and acceleration limits, and improving the safety and reliability of mission execution. From the macro-level mission plan to the specific path of a single UAV, the data at each stage is interconnected, preserving the safety characteristics of the three-dimensional airspace while ensuring the stability of formation operation through refined calculations, providing strong support for the efficient completion of material transportation tasks by UAV formations.

[0051] In a preferred embodiment of the present invention, step 6 includes: Step 6.1 involves discretizing the continuous motion path of each drone in the drone ensemble according to a preset control period, generating a discrete-time position and attitude point sequence for each drone. Specifically, this includes: first, determining the preset control period. The control period must match the hardware response frequency of the drone flight control system and be an integer multiple of the minimum sampling period of the flight control system to avoid data conflicts. Based on the performance parameters of common industrial-grade drones, the preset control period is set to 0.1 seconds, meaning the flight control system receives the control target every 0.1 seconds. This period ensures a control accuracy with a discrete error ≤ 0.1 meters without exceeding the flight control system's limits. The data processing capability is then used to generate a discrete time node sequence. First, the total duration T of the continuous motion path of a single UAV is extracted (for example, the total duration from waypoint A to B is 2 seconds). Then, the total duration is divided according to the control period Δt = 0.1 seconds to obtain a discrete time node sequence of t0 = 0 seconds, t1 = 0.1 seconds, t2 = 0.2 seconds, ..., t20 = 2 seconds, a total of 21 time nodes, which completely cover the entire motion process. The calculation of discrete position and attitude both adopt the linear interpolation method to ensure smooth transition. For any discrete time node tI on the continuous path, its position is obtained by interpolating the positions of two adjacent known times. Let... ( <tI< The position at that time is ( , , ), The position of the moment is ( , , ),but The formula for calculating the position PI at time 1 is: , , Calculated using the same logic; attitude is based on roll angle. - Pitch angle -Yaw angle This indicates that discrete attitudes are also obtained through linear interpolation, let... At any given moment, the attitude is ( , , ), Time for ( , , ),but Roll angle of time Pitch and yaw angles are calculated using the same logic. After completing the position and attitude calculations for all discrete time nodes, the discrete time nodes, corresponding position coordinates, and corresponding attitude angles of each UAV are combined sequentially in chronological order to form a discrete time position and attitude point sequence for that UAV.

[0052] Step 6.2: Based on the predefined UAV dynamics model, the discrete-time position and attitude point sequence of each UAV is converted into a low-level flight control command sequence containing throttle control, pitch control, roll control, and yaw control. Specifically, this includes: First, constructing a six-degree-of-freedom UAV dynamics model (a predefined UAV dynamics model) to provide a theoretical basis for control calculation. The first step is to define the coordinate system, establishing an inertial coordinate system (with the ground station location as the origin, horizontal forward as the forward-backward axis, horizontal right as the left-right axis, and vertical upward as the altitude axis) and a body coordinate system (with the UAV's center of mass as the origin, the fuselage forward as the longitudinal axis, right as the transverse axis, and vertical upward as the longitudinal axis). The two coordinate systems are transformed through attitude angles (roll angle, pitch angle, and yaw angle). The first step is to ensure that force and torque calculations are performed under a unified benchmark. The second step is to construct position dynamics equations to describe the translational motion of the UAV in the inertial coordinate system. The translational motion of the UAV in the forward, backward, left, and right directions all follow the rule that mass multiplied by acceleration equals the resultant force of all forces acting in that direction. Among them, the forces acting in the altitude direction include the downward gravity (equal to the UAV's mass multiplied by the gravitational acceleration of 9.8 m / s²) and the upward lift (related to throttle control and determined by the single motor thrust and lift coefficient). The forces acting in the forward and backward and left and right directions are mainly air resistance (equal to the air resistance coefficient multiplied by the square of the velocity in the corresponding direction, in the opposite direction to the direction of motion). The relationship between position acceleration and forces is established through this equation. The specific equation is as follows: For the forward direction, UAV mass × forward acceleration = Left and right direction: Drone mass × left and right acceleration = In the altitude direction, the drone's mass multiplied by its altitude acceleration equals... The third step is to construct attitude dynamics equations to describe the rotational motion of the UAV in the body coordinate system. The rotational motion of the UAV around the three axes of roll, pitch, and yaw follows the rule that the moment of inertia multiplied by the angular acceleration equals the resultant torque of all torques on that axis. The roll torque is generated by the thrust difference between the left and right motors, the pitch torque by the thrust difference between the front and rear motors, and the yaw torque by the thrust difference between the diagonal motors. At the same time, the air damping torque (proportional to the angular velocity and opposite to the rotation direction) is considered. The relationship between attitude angular acceleration and torque is established through this equation. The specific equations are as follows: Roll direction: Roll inertia × Roll angular acceleration = Roll torque - Air damping coefficient × Roll angular velocity; Pitch direction: Pitch inertia × Pitch angular acceleration = Pitch torque - Air damping coefficient × Pitch angular velocity; Yaw direction: Yaw inertia × Yaw angular acceleration = Yaw torque - Air damping coefficient × Yaw angular velocity.

[0053] After completing the construction of the six-degree-of-freedom dynamic characteristic model of the UAV, the core parameters of the model are retrieved. These parameters must be consistent with the previously established formation configuration parameters, including UAV mass (the mass of a single UAV is fixed at 4 kg), moments of inertia in each direction (roll moment of inertia 0.3 kg·m², pitch moment of inertia 0.4 kg·m², yaw moment of inertia 0.5 kg·m²), lift coefficient 0.8, drag coefficient 0.6, and maximum thrust of a single motor 50 N. Simultaneously, the PID controller parameters are defined (proportional coefficient 0.5, integral coefficient 0.1, derivative coefficient 0.2). These parameters are the basis for calculating the control quantities and ensure matching with the actual physical characteristics of the UAV. Next, the discrete time parameters of each UAV are extracted. The system generates a sequence of position and attitude points. From this sequence, the target position (including forward / backward position, left / right position, and altitude) and target attitude (including roll angle, pitch angle, and yaw angle) for each discrete time node are obtained. Simultaneously, the actual position and attitude at the current time node are obtained through real-time feedback data from the UAV. The deviation data for each dimension are calculated one by one. That is, the forward / backward position deviation is equal to the target forward / backward position minus the actual forward / backward position; the left / right position deviation is equal to the target left / right position minus the actual left / right position; the altitude deviation is equal to the target altitude position minus the actual altitude position; the roll angle deviation is equal to the target roll angle minus the actual roll angle; the pitch angle deviation is equal to the target pitch angle minus the actual pitch angle; and the yaw angle deviation is equal to the target yaw angle minus the actual yaw angle.

[0054] Then, the throttle control value is calculated. This throttle control value is used to adjust the UAV's lift, ensuring altitude stability and vertical movement. First, the integral value of the altitude deviation is calculated, which is the sum of all altitude deviations from the initial time point to the current time point multiplied by the control period (0.1 seconds). Next, the rate of change of the altitude deviation is calculated, which is the altitude deviation at the current time point minus the altitude deviation at the previous time point, divided by the control period. Finally, the target's vertical acceleration is calculated, which equals the proportional coefficient multiplied by the altitude deviation, plus the integral coefficient multiplied by the integral value of the altitude deviation, plus the derivative. The coefficient is multiplied by the rate of change of altitude deviation; then, the required lift is calculated based on the vertical force balance relationship. The lift is equal to the mass of the UAV multiplied by the gravitational acceleration (valued at 9.8 m / s²) plus the mass of the UAV multiplied by the vertical acceleration of the target, i.e., lift = UAV mass × 9.8 + UAV mass × target vertical acceleration; finally, the throttle control amount is calculated based on the correspondence between lift and throttle control amount. Throttle control amount = lift ÷ (maximum thrust of a single motor × lift coefficient) × 100%. The result must be limited to between 0% and 100% to avoid exceeding the range of motor thrust.

[0055] Next, calculate the pitch control input. The pitch control input is used to adjust the UAV's forward and backward attitude and directional movement. It is related to the forward and backward position deviation and the pitch angle deviation. First, calculate the integral value of the pitch angle deviation, which is the sum of all pitch angle deviations from the initial time node to the current time node multiplied by the control period. Then, calculate the rate of change of the pitch angle deviation, which is the pitch angle deviation at the current time node minus the pitch angle deviation at the previous time node, divided by the control period. Next, calculate the target pitch angle acceleration. The target pitch angle acceleration is equal to the proportional coefficient multiplied by the pitch angle deviation, plus the integral coefficient multiplied by the integral value of the pitch angle deviation, plus the derivative coefficient multiplied by the rate of change of the pitch angle deviation, while also considering the forward and backward position deviations. Deviation correction is performed by multiplying the proportional coefficient by the forward and backward position deviation. The final target pitch angle acceleration is calculated as: proportional coefficient × pitch angle deviation + integral coefficient × integral value of pitch angle deviation + differential coefficient × rate of change of pitch angle deviation + proportional coefficient × forward and backward position deviation. Then, the required pitch torque is calculated based on the relationship between torque and angular acceleration: pitch torque = pitch direction rotational inertia × target pitch angle acceleration. Finally, the pitch control quantity is calculated based on the correspondence between pitch torque and pitch control quantity: pitch control quantity = pitch torque ÷ (maximum thrust of a single motor × lever arm coefficient). The lever arm coefficient is set to 0.3, and the result is limited to between -30° and 30° (corresponding to the pitch servo control range).

[0056] Next, the roll control variable is calculated. The roll control variable is used to adjust the UAV's left and right attitude and left and right directional movement. It is related to the left and right position deviation and roll angle deviation. The calculation logic is the same as that of the pitch control variable. First, the integral value of the roll angle deviation (the sum of all roll angle deviations multiplied by the control period) and the rate of change of the roll angle deviation (the current roll angle deviation minus the previous roll angle deviation divided by the control period) are calculated. Then, the target roll angle acceleration is calculated: target roll angle acceleration = proportional coefficient × roll angle deviation + integral coefficient × integral value of roll angle deviation + derivative coefficient × rate of change of roll angle deviation + proportional coefficient × left and right position deviation. Next, the required roll torque is calculated: roll torque = moment of inertia in the roll direction × target roll angle acceleration. Finally, the roll control variable is calculated: roll control variable = roll torque ÷ (maximum thrust of a single motor × lever arm coefficient). The lever arm coefficient is taken as 0.3, and the result is limited to between -30° and 30° (corresponding to the roll servo control range).

[0057] Finally, the yaw control variable is calculated. This variable adjusts the UAV's heading and is related to the yaw angle deviation. First, the integral value of the yaw angle deviation (the sum of all yaw angle deviations multiplied by the control period) and the rate of change of the yaw angle deviation (the current yaw angle deviation minus the previous yaw angle deviation divided by the control period) are calculated. Then, the target yaw angle acceleration is calculated: Target yaw angle acceleration = proportional coefficient × yaw angle deviation + integral coefficient × integral value of yaw angle deviation + derivative coefficient × rate of change of yaw angle deviation. Next, the required yaw moment is calculated: Yaw moment = yaw angle... The rotational inertia is multiplied by the target yaw acceleration. Finally, the yaw control quantity is calculated: yaw control quantity = yaw torque ÷ (maximum thrust of a single motor × lever arm coefficient). The lever arm coefficient is 0.4, and the result is limited to between -45° and 45° (corresponding to the yaw servo control range). The throttle control quantity, pitch control quantity, roll control quantity, and yaw control quantity calculated at each discrete time node are arranged in chronological order to form the underlying flight control command sequence for each UAV, ensuring that each control command corresponds one-to-one with the discrete time position and attitude point sequence.

[0058] Step 6.3: Through the wireless data link between the ground station and the UAVs, the underlying flight control command sequence of each UAV in the UAV group is sent to the corresponding UAV to drive all UAVs to complete the collaborative transportation task according to the mission execution plan. Specifically, this includes: First, configuring the wireless data link. For a formation of 4 to 8 UAVs, a master-slave link architecture is preferred. The ground station is equipped with a master link module (operating frequency 2.4GHz, transmission rate 500kbps, latency ≤50ms), and each UAV is equipped with a slave link module. The master link establishes point-to-point communication with all slave links to ensure that commands can be sent independently to each UAV. At the same time, the link adopts a CRC cyclic redundancy check mechanism. The check formula is: check code = bitwise XOR operation of command data bytes. The receiving end verifies the integrity of the command through the check code. If the verification fails, it automatically requests a retransmission. To ensure reliable command transmission, the ground station processes commands in three steps: First, it encodes the control command sequence for each UAV according to the flight control protocol, for example, encoding 30% throttle as hexadecimal byte 0x1E to ensure correct decoding and recognition by the UAV. Second, to ensure coordinated formation flight, the ground station synchronously triggers the issuance of commands to all UAVs at a control cycle of 0.1 seconds. At each discrete time point, it first verifies the command execution feedback of all UAVs from the previous cycle (e.g., command signals have been received), and only after confirming that there are no errors does it simultaneously issue the command at time tᵢ, avoiding formation position deviations due to command issuance time differences. After the command is issued, the ground station receives real-time status data such as position, attitude, and battery voltage from the UAVs. If the actual state of a UAV deviates from the target by more than a preset threshold (e.g., position deviation > 1 meter), an emergency command adjustment is immediately triggered to ensure flight safety.

[0059] The process of a UAV receiving and executing commands is as follows: After receiving commands from the link module, the flight control chip decodes the commands, extracts the throttle, pitch, roll, and yaw control quantities, and converts them into drive signals for the actuators (e.g., throttle control quantities correspond to motor PWM signals, and attitude control quantities correspond to servo drive signals). The actuators (motors and servos) act according to the drive signals. At the same time, the flight control system collects its own position and attitude data through GPS and IMU (Inertial Measurement Unit), calculates the deviation between the actual state and the target, and uses this data for its own closed-loop control correction. On the other hand, the state data is encoded and fed back to the ground station. If a UAV deviates slightly due to external interference, the flight control system will first obtain the state data of 2 to 3 adjacent UAVs in real time through the short-range communication link within the formation, including the actual position, attitude, and control command execution feedback of the adjacent UAVs. Then, it compares its own deviation with the state of the adjacent UAVs to determine whether the deviation is unique to itself (rather than a shift in the overall formation). If the deviation is confirmed to be self-inflicted, the average position and attitude of adjacent UAVs are used as a reference benchmark. Combined with the preset formation spacing requirements (e.g., horizontal spacing error between adjacent UAVs does not exceed 0.5 meters, attitude angle error does not exceed 1°), the required additional position and attitude corrections are calculated. The position correction is calculated separately in three directions: forward / backward, left / right, and altitude. First, the forward / backward positions of all adjacent UAVs are added together and divided by the number of adjacent UAVs to obtain the average forward / backward position. Then, the preset forward / backward spacing between the self and adjacent UAVs in the formation is subtracted from this average position to obtain the self's correct forward / backward target reference position. Finally, the forward / backward position correction equals the forward / backward target reference position. The reference position is subtracted from the actual forward / backward position. Position corrections in the left / right and altitude directions are calculated using the same logic, with a maximum single position correction of 0.2 meters to avoid excessive adjustments. Attitude corrections are also calculated separately for roll, pitch, and yaw angles. First, the roll angles of all adjacent UAVs are summed and divided by the number of adjacent UAVs to obtain the average roll angle. The roll direction correction equals the average roll angle minus the actual roll angle. If the absolute value of this correction exceeds 1°, ±1° is taken as the final roll correction (sign consistent with the calculation result); otherwise, the calculated result is used directly. Attitude corrections for pitch and yaw angles are calculated using the same logic. This correction is then superimposed on the next round of control commands, fine-tuning the output amplitude of throttle and attitude control to gradually reduce the state difference with adjacent UAVs. Throughout the collaborative correction process, the previously determined safe operating boundaries are strictly followed to avoid instability due to excessive adjustments, while ensuring that correction actions do not affect the flight status of adjacent UAVs, ultimately maintaining the stability of the overall formation configuration and driving all UAVs to complete the collaborative transport mission according to the mission execution plan.

[0060] This embodiment discretizes the control cycle, matching the flight control response frequency, to convert a continuous path into a target sequence executable by the flight control, avoiding the mismatch between continuous commands and hardware processing capabilities. Linear interpolation for discrete position and attitude calculation ensures smooth transitions between adjacent target points, reducing impact and jitter during UAV movement. Inverse kinematics calculation of control quantities based on the UAV's six-degree-of-freedom dynamics model ensures that throttle and attitude control commands strictly conform to the UAV's physical motion limits, avoiding risks such as stall and overload due to commands exceeding performance limits. The application of PID control algorithms enables real-time compensation for attitude deviations, improving the UAV's resistance to interference and ensuring that the actual flight state closely approximates the target. The combination of a master-slave wireless link and CRC check mechanism reduces command transmission latency and loss probability, guaranteeing the real-time performance and reliability of control commands. Synchronous distribution and status monitoring from the ground station, along with the UAV's collaborative correction logic, ensures that all UAVs execute tasks at a unified pace, avoiding the accumulation of individual deviations within the formation and providing core assurance for the safe completion of collaborative transportation missions.

[0061] like Figure 2 As shown, embodiments of the present invention also provide a system information processing system based on power data, comprising: The parsing module is used to parse the original task request, generate a structured task description that includes the coordinates of the repair target point, the quality and envelope size of the materials to be transported, and the task time window. Based on the coordinates of the repair target point, it obtains the elevation and meteorological data of the corresponding area and generates a three-dimensional environment model of the task area. The identification module is used to perform spatial rasterization sampling on the three-dimensional environment model to obtain a spatial discrete sampling field, and to identify the target emergency repair structure sampling cluster and the spatial channel sampling cluster to obtain the reference direction vector and the offset direction vector. The module is used to generate a three-dimensional spatial cone based on the centroid of the sampling cluster of the target emergency repair structure, the reference direction vector and the offset direction vector, and to generate a three-dimensional feasible spatial network with cone risk weights based on the envelope size of the material to be transported and the three-dimensional spatial cone. The adaptation module is used to screen candidate UAV formation configurations based on the three-dimensional feasible airspace network and the quality of the goods to be transported, and to generate a set of adapted formation configurations and safe operation boundaries for each path based on the candidate UAV formation configurations and the three-dimensional feasible airspace network. The execution module is used to generate a mission execution plan and the continuous motion path and formation transition manifold of each UAV in the execution UAV set based on a three-dimensional feasible airspace network, mission time window and adaptive formation configuration set; The control module is used to generate a sequence of control commands for each UAV based on the continuous motion path and formation configuration transition manifold of each UAV, and then send them to the UAV set to perform the collaborative transportation task.

[0062] It should be noted that this system is a system corresponding to the above method. All implementation methods in the above method embodiments are applicable to this embodiment and can achieve the same technical effect.

[0063] Embodiments of the present invention also provide a computer-readable storage medium storing instructions that, when executed on a computer, cause the computer to perform the method described above. All implementations in the above method embodiments are applicable to this embodiment and can achieve the same technical effects.

[0064] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for collaborative task allocation and processing between multiple emergency power units and ground stations in mountainous areas, characterized in that: The method includes: Step 1: Parse the original task request and generate a structured task description that includes the coordinates of the target point to be repaired, the quality and envelope size of the materials to be transported, and the task time window. Based on the coordinates of the target point to be repaired, obtain the elevation and meteorological data of the corresponding area and generate a three-dimensional environmental model of the task area. Step 2: Spatial rasterization sampling is performed on the three-dimensional environment model to obtain a spatial discrete sampling field, and the target emergency repair structure sampling cluster and the spatial channel sampling cluster are identified to obtain the reference direction vector and the offset direction vector. Step 3: Based on the centroid of the target emergency repair structure sampling cluster, the reference direction vector and the offset direction vector, generate a three-dimensional spatial cone, and based on the envelope size of the material to be transported and the three-dimensional spatial cone, generate a three-dimensional feasible spatial network with cone risk weights. Step 4: Based on the three-dimensional feasible airspace network and the quality of the materials to be transported, candidate UAV formation configurations are selected, and based on the candidate UAV formation configurations and the three-dimensional feasible airspace network, a set of suitable formation configurations and safe operation boundaries corresponding to each path are generated. Step 5: Based on the three-dimensional feasible airspace network, mission time window and adaptive formation configuration set, generate mission execution plan and continuous motion path and formation configuration transition manifold of each UAV in the execution UAV set; Step 6: Based on the continuous motion path and formation configuration transition manifold of each UAV, generate the control command sequence of each UAV and send it to the execution UAV set to perform the collaborative transportation task.

2. The method for collaborative task allocation and processing of multi-unit and ground station power emergency response in mountainous areas according to claim 1, characterized in that, Step 1 includes: Extract text and data fields from the original task request, and parse them to obtain the initial emergency repair target point coordinates, initial material quality and initial envelope size information, and initial task time requirements; The coordinates of the initial emergency repair target point are converted between the geodetic coordinate system and the local engineering coordinate system to obtain the coordinates of the emergency repair target point; the initial material mass and initial envelope size information are standardized and the units are unified to obtain the mass of the material to be transported and the envelope size of the material to be transported; the initial task time requirement is processed by time windowing to obtain the task time window; the coordinates of the emergency repair target point, the mass and envelope size of the material to be transported, and the task time window are integrated to generate a structured task description; Based on the coordinates of the target repair point in the structured task description, digital elevation data and historical meteorological statistics data of the corresponding area are retrieved from the pre-set geographic information database. The digital elevation data is converted into a regular grid terrain elevation matrix. The prevailing wind direction and average wind speed data are extracted from the historical meteorological statistics data set and interpolated to generate a wind field vector matrix that is spatially aligned with the regular grid terrain elevation matrix. By spatially overlaying and fusing the regular grid terrain elevation matrix and the wind field vector matrix, a three-dimensional environment model of the task area containing terrain elevation attributes and wind field vector attributes is constructed.

3. The method for collaborative task allocation and processing of multi-unit and ground station emergency power systems in mountainous areas according to claim 2, characterized in that, Step 2 includes: Based on a preset spatial sampling resolution, the 3D space containing the 3D environment model is divided into a uniform 3D regular grid, with the geometric center of each grid cell serving as a sampling point. For each sampling point, the terrain elevation attribute value corresponding to the sampling point is extracted from the 3D environment model, and the mean magnitude of the wind field vector attribute at all locations within the grid cell containing the sampling point is calculated. The terrain elevation attribute value and the mean magnitude of the wind field vector are combined to form the sampling feature vector of the sampling point. The set of sampling feature vectors of all sampling points constitutes a spatial discrete sampling field. Threshold segmentation is performed based on the terrain elevation attribute values ​​of each sampling point in the spatial discrete sampling field. The first type of candidate point set with elevations higher than the preset elevation threshold is selected. Euclidean distance clustering is then performed on the first type of candidate point set to obtain multiple independent clusters. The cluster with the smallest spatial distance to the target point coordinates is selected from the multiple independent clusters and the corresponding cluster is defined as the target repair structure sampling cluster. Threshold segmentation is performed based on the wind field vector magnitude attribute value of each sampling point in the spatial discrete sampling field. A second set of candidate points with a magnitude lower than the preset wind resistance threshold is selected. Based on the second set of candidate points, points located in the spatial corridor between the target emergency repair structure sampling cluster and the preset flight starting point are selected to form a third set of candidate points. Density clustering is performed on the third type of candidate point set to extract the connected regions with the highest data point density, and the corresponding regions are defined as spatial channel sampling clusters. Principal component analysis is performed on the spatial coordinates of all sampling points in the target repair structure sampling cluster to obtain the first principal component direction, and the corresponding direction is used as the reference direction vector. Principal component analysis is performed on the spatial coordinates of all sampling points in the spatial channel sampling cluster to obtain the first principal component direction, and the corresponding direction is used as the offset direction vector.

4. The method for collaborative task allocation and processing of multi-unit and ground station power emergency response in mountainous areas according to claim 3, characterized in that, Step 3 includes: Calculate the arithmetic mean of the spatial coordinates of all sampling points in the target emergency repair structure sampling cluster to obtain the centroid of the target emergency repair structure sampling cluster. Then, with the centroid of the target emergency repair structure sampling cluster as the vertex and the reference direction vector as the direction of the central axis, calculate the projection vector of the offset direction vector on the plane perpendicular to the central axis. Calculate the angle between the projection vector and the central axis to obtain the spatial azimuth angle. Generate a three-dimensional spatial cone with the central axis as the rotation axis and the spatial azimuth angle as the cone half-angle. The internal space of the three-dimensional spatial cone is uniformly divided into a series of continuous voxel units along the radial, circumferential and axial directions. For each voxel unit, an envelope space cuboid is constructed according to the envelope size of the material to be transported, with its edges aligned with the center point of the voxel unit and the envelope size of the material to be transported. The minimum spatial distance between the outer surface of the envelope space cuboid and the cone surface of the three-dimensional spatial cone is calculated and denoted as the first distance. The minimum spatial distance between the outer surface of the envelope space cuboid and the surfaces of all obstacles in the three-dimensional environment model is calculated and denoted as the second distance. The minimum value between the first distance and the second distance is taken as the minimum directed distance of the voxel unit. Based on the preset piecewise linear distance and risk mapping function, the minimum directed distance of each voxel unit is mapped to a risk weight. All voxel units and their corresponding risk weights together constitute a three-dimensional feasible spatial network with cone-shaped risk weights.

5. The method for collaborative task allocation and processing of multi-unit and ground station emergency power systems in mountainous areas according to claim 4, characterized in that, Step 4 includes: Based on the weight of the goods to be transported, all formation configurations with a nominal total load capacity greater than the weight of the goods to be transported are selected from the pre-set formation configuration library to form an initial candidate formation configuration set. For each pre-planned candidate flight path in the three-dimensional feasible airspace network, perform path adaptability analysis, screen the adaptable formation configurations of the candidate flight paths, and determine the safe operation boundaries of the adaptable formation configurations on the candidate flight paths. Summarize all suitable formation configurations and safe operating boundaries for each candidate flight path, and generate a set of suitable formation configurations and safe operating boundaries for each candidate flight path.

6. The method for collaborative task allocation and processing of multi-unit and ground station emergency power systems in mountainous areas according to claim 5, characterized in that, The path adaptability analysis includes: Based on the cone risk weights corresponding to all voxel units traversed by the candidate flight path and the wind field vector properties in the 3D environment model, the interference force spindle of each waypoint on the candidate flight path is calculated. For each formation configuration in the initial candidate formation configuration set, the maximum anti-disturbance dynamic spinor set of the formation configuration at each waypoint of the candidate flight path is calculated. By performing an inclusion test between the interference force spinor at each waypoint and the maximum anti-disturbance dynamic spinor set of the formation configuration at the corresponding waypoint, the suitable formation configuration for the candidate flight path is selected. For each selected suitable formation configuration, the safe operating boundary of the suitable formation configuration on the candidate flight path is determined based on the minimum Euclidean distance from the boundary of the maximum anti-disturbance dynamic spinor set of the formation configuration to the disturbance force spinor at each waypoint of the candidate flight path.

7. The method for collaborative task allocation and processing of multi-unit and ground station emergency power systems in mountainous areas according to claim 6, characterized in that, Step 5 includes: Based on the three-dimensional feasible airspace network, the mission time window, and the set of suitable formation configurations and safe operation boundaries corresponding to each of the candidate flight paths, a multi-objective programming solution is performed with the goal of shortening the total flight time and improving the average safe operation boundary to generate a mission execution plan. The mission execution plan includes the selected set of execution UAVs, the selected flight trajectory from the preset starting point to the repair target point and through the three-dimensional feasible airspace network, and the formation configuration specified for each waypoint on the selected flight trajectory. Obtain waypoints A and B that are consecutively flying on the selected flight path. Obtain the three-dimensional position coordinate set A of each UAV in the UAV set under the specified formation configuration of waypoint A, and the three-dimensional position coordinate set B of each UAV in the UAV set under the specified formation configuration of waypoint B. For each drone in the set of drones, project the three-dimensional position coordinates of the corresponding drone in set A and the three-dimensional position coordinates in set B onto a virtual unit sphere to obtain the spherical projection point pair corresponding to the drone. On a unit sphere, calculate the shortest great circle arc path connecting the corresponding UAV spherical projection point pairs. Based on the shortest great circle arc path, back-project back into three-dimensional space to construct a smooth spatial transition curve from the configuration position of waypoint A to the configuration position of waypoint B for each UAV in the execution UAV set. The set of smooth spatial transition curves of all UAVs in the execution UAV set collectively characterizes the formation configuration transition manifold. Based on each smooth spatial transition curve in the formation configuration transition manifold, a continuous motion path from waypoint A configuration to waypoint B configuration is generated for the corresponding UAV in the execution UAV set.

8. The method for collaborative task allocation and processing of multi-unit and ground station emergency power systems in mountainous areas according to claim 7, characterized in that, Step 6 includes: The continuous motion path of each drone in the drone ensemble is discretized in time according to a preset control cycle to generate a discrete time position and attitude point sequence for each drone. Based on the predefined UAV dynamics model, the discrete-time position and attitude point sequence of each UAV is converted into a low-level flight control command sequence containing throttle control, pitch control, roll control and yaw control. Through the wireless data link between the ground station and the drones, the underlying flight control command sequence of each drone in the drone ensemble is sent to the corresponding drone to drive all drones to complete the collaborative transportation task according to the mission execution plan.

9. A multi-unit and ground station collaborative task allocation and processing system for emergency power supply in mountainous areas, wherein the system implements the method as described in any one of claims 1 to 8, characterized in that, include: The parsing module is used to parse the original task request, generate a structured task description that includes the coordinates of the repair target point, the quality and envelope size of the materials to be transported, and the task time window. Based on the coordinates of the repair target point, it obtains the elevation and meteorological data of the corresponding area and generates a three-dimensional environment model of the task area. The identification module is used to perform spatial rasterization sampling on the three-dimensional environment model to obtain a spatial discrete sampling field, and to identify the target emergency repair structure sampling cluster and the spatial channel sampling cluster to obtain the reference direction vector and the offset direction vector. The module is used to generate a three-dimensional spatial cone based on the centroid of the sampling cluster of the target emergency repair structure, the reference direction vector and the offset direction vector, and to generate a three-dimensional feasible spatial network with cone risk weights based on the envelope size of the material to be transported and the three-dimensional spatial cone. The adaptation module is used to screen candidate UAV formation configurations based on the three-dimensional feasible airspace network and the quality of the goods to be transported, and to generate a set of adapted formation configurations and safe operation boundaries for each path based on the candidate UAV formation configurations and the three-dimensional feasible airspace network. The execution module is used to generate a mission execution plan and the continuous motion path and formation transition manifold of each UAV in the execution UAV set based on a three-dimensional feasible airspace network, mission time window and adaptive formation configuration set; The control module is used to generate a sequence of control commands for each UAV based on the continuous motion path and formation configuration transition manifold of each UAV, and then send them to the UAV set to perform the collaborative transportation task.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a program that, when executed by a processor, implements the method as described in any one of claims 1 to 8.

Images