Clustering-Based Multi-Region Coverage-Path Planning of Heterogeneous UAVs

: Unmanned aerial vehicles (UAVs) multi-area coverage-path planning has a broad range of applications in agricultural mapping and military reconnaissance. Compared to homogeneous UAVs, heterogeneous UAVs have higher application value due to their superior ﬂexibility and efﬁciency. Nevertheless, variations in performance parameters among heterogeneous UAVs can signiﬁcantly amplify computational complexity, posing challenges to solving the multi-region coverage path-planning problem. Consequently, this study studies a clustering-based method to tackle the multi-region coverage path-planning problem of heterogeneous UAVs. First, the constraints necessary during the planning process are analyzed, and a planning formula based on an integer linear programming model is established. Subsequently, this problem is decomposed into regional allocation and visiting order optimization subproblems. This study proposes a novel clustering algorithm that utilizes centroid iteration and spatiotemporal similarity to allocate regions and adopts the nearest-to-end policy to optimize the visiting order. Additionally, a distance-based bilateral shortest-selection strategy is proposed to generate region-scanning trajectories, which serve as trajectory references for real ﬂight. Simulation results in this study prove the effective performance of the proposed clustering algorithm and region-scanning strategy.


Introduction
With the rapid development of artificial intelligence and automated control [1][2][3], UAVs have been widely used in both military and civilian fields, including reconnaissance and strike operations [4,5], target tracking [6][7][8][9], forest fire prevention [10], regional surveillance [11][12][13], etc. Due to limitations in performance and payload capacity, it is often challenging for one UAV to complete complex missions [14].Therefore, the study of multi-UAV systems, which have good scalability and cooperation capabilities, has gained significant attention in current research.In order to maximize the overall effectiveness of the multi-UAV system, it is necessary to research task planning to obtain suitable task orders and flight paths of UAVs.
As an important branch of multi-UAV task planning, coverage-path planning includes region allocation and path planning.It has been studied by numerous scholars from various aspects, such as region shape [15,16], energy constraints [17][18][19][20], and obstacle avoidance [21,22].Nielsen [15] tackled the issue of region coverage for non-convex polygons by dividing the area into numerous separate convex sub-polygons and utilizing a scanning pattern to ensure complete coverage.Huang [17] determined the energy consumption of UAVs in various flight modes and presented a coverage path-planning algorithm Drones 2023, 7, 664 3 of 20 programming.The problem of multi-area coverage-path planning is divided into two subproblems: allocating the task regions and determining the execution order.This decomposition significantly reduces the complexity associated with solving the problem; 2.
Based on the iterative idea of the k-means algorithm and the requirement of multiregion coverage, a novel clustering algorithm based on spatiotemporal similarity is proposed, and a clustering center iteration method is designed to complete regions clustering; 3.
A novel method is proposed to minimize the flight distance of a single UAV when scanning multiple regions.The method involves selecting entry points and scanning patterns based on the shortest flight distances.
1.An optimization model is employed to tackle the problem of multi-region coveragepath planning of heterogeneous UAVs, and it is formulated using integer linear programming.The problem of multi-area coverage-path planning is divided into two subproblems: allocating the task regions and determining the execution order.This decomposition significantly reduces the complexity associated with solving the problem; 2. Based on the iterative idea of the k-means algorithm and the requirement of multiregion coverage, a novel clustering algorithm based on spatiotemporal similarity is proposed, and a clustering center iteration method is designed to complete regions clustering; 3. A novel method is proposed to minimize the flight distance of a single UAV when scanning multiple regions.The method involves selecting entry points and scanning patterns based on the shortest flight distances.
This study is organized as follows.Section 2 establishes the model and formulas for the multi-region coverage path-planning problem of heterogeneous UAVs.Section 3 proposes a novel region clustering algorithm and optimizes the visited order of regions.Section 4 introduces a trajectory-planning method based on the shortest flight distance.Section 5 presents simulation experiments and comparative analysis.Finally, Section 6 concludes this work.

System Model
The coverage path-planning problems of heterogeneous UAVs are often classified as non-deterministic polynomial hard (NP-hard); it is difficult to obtain an accurate solution directly [33].So, we first analyze the constraints that need to be satisfied when the UAV system performs tasks and obtain the exact formulation.Then, in order to improve the solution efficiency, we intend to solve this problem from two aspects: regional allocation and regional sequence optimization, which are achieved through regional clustering and subregion reordering, respectively.Finally, we present a multi-region coverage-path planning strategy specifically for the practical flight trajectory of UAVs.

Problem Description
This study deploys n fixed-wing heterogeneous UAVs ; the task scenario diagram is shown in Figure 2. A list of key symbols used in this study and their definitions are provided in Table 1.In this work, the main differences between the UAVs are their flight speed, endurance, and scanning width.So, each UAV is characterized as , , , , This study is organized as follows.Section 2 establishes the model and formulas for the multi-region coverage path-planning problem of heterogeneous UAVs.Section 3 proposes a novel region clustering algorithm and optimizes the visited order of regions.Section 4 introduces a trajectory-planning method based on the shortest flight distance.Section 5 presents simulation experiments and comparative analysis.Finally, Section 6 concludes this work.

System Model
The coverage path-planning problems of heterogeneous UAVs are often classified as non-deterministic polynomial hard (NP-hard); it is difficult to obtain an accurate solution directly [33].So, we first analyze the constraints that need to be satisfied when the UAV system performs tasks and obtain the exact formulation.Then, in order to improve the solution efficiency, we intend to solve this problem from two aspects: regional allocation and regional sequence optimization, which are achieved through regional clustering and subregion reordering, respectively.Finally, we present a multi-region coverage-path planning strategy specifically for the practical flight trajectory of UAVs.

Problem Description
This study deploys n fixed-wing heterogeneous UAVs U = {U 1 , U 2 , . . . ,U n } to carry out reconnaissance and scanning missions across m rectangular regions R = {R 1 , R 2 , . . . ,R m }; the task scenario diagram is shown in Figure 2. A list of key symbols used in this study and their definitions are provided in Table 1.In this work, the main differences between the UAVs are their flight speed, endurance, and scanning width.So, each UAV is characterized as U i = Uid i , U pos i , V i , T i , d i , where Uid i and U pos i represent the serial number and the takeoff coordinate of U i , respectively.V i denotes the cruising flight speed, and we assume that the cruising for each UAV maintains a constant throughout the mission.T i represents the maximum endurance for cruising flight, while d i represents the scan width of sensors installed in U i .
represent the serial number and the takeoff coordinate of i U , respectively.i V denotes the cruising flight speed, and we assume that the cruising for each UAV maintains a constant throughout the mission.i T represents the maximum endurance for cruising flight, while i d represents the scan width of sensors installed in i U .An m -row m -column matrix is adopted to represent the distance between regions j R and k R , which is calculated using the Euclidean distance between their center points.In this study, heterogeneous UAVs are required to travel from their takeoff positions i Upos to the designated mission area.Upon finishing their assigned area scanning tasks, the UAVs must return to their respective starting points.As a result, , , i j k TF is used to represent the time consumption of i U flying from one region j R to another region k R .Since we do not know the order in which the UAV scans regions dur- ing the allocation phase, the flying distance is roughly estimated using the distance between the centers of these two regions.
, i j

TS
is used to represent the time consumption of i U scanning the region j R ; each UAV will use the back-and-forth scanning method to scan a single area.
, , i j k TF and , i j TS can be calculated using , , , ,

Symbol Definitions U
The set of UAVs  The ith UAV in the set U n The number of UAVs in the set U U pos i The takeoff position of U i (also called the base of U i ) The cruising flight speed of U i T i The maximum flight endurance of U i d i The scan width of onboard sensors of U i R The set of regions R j The jth region in the set R m The number of regions in the set R Rpos j The center coordinate of R j Rbeta j The angle between the major axis and the X-axis of the map L j The length of the major axis of R j W j The length of the minor axis of R j D j,k The distance between R j and R k TF i,j,k The time consumption of U i in flying from R j to R k TS i,j The time consumption of U i in scanning R j C The set of clusters C i The ith cluster in the set C STS i,j The spatiotemporal similarity between C i and R j Que i The regions queue of C i Seq i The sequence of regions that U i is going to coverage {Lin, Win}  The optional entry points set of a region In this study, each region's features are shown as R j = Rpos j , Rbeta j , L j , W j .In this notation, Rpos j represents the coordinates of the center point in the rectangular region.Rbeta j signifies the angle formed between the major axis of the region and the X-axis of the map, with a counterclockwise direction identified as positive.L j represents the length of the major axis, and W j denotes the length of the minor axis.
An m-row m-column matrix D = D j,k is adopted to represent the distance between regions R j and R k , which is calculated using the Euclidean distance between their center points.In this study, heterogeneous UAVs are required to travel from their takeoff positions U pos i to the designated mission area.Upon finishing their assigned area scanning tasks, the UAVs must return to their respective starting points.As a result, TF i,j,k is used to represent the time consumption of U i flying from one region R j to another region R k .Since we do not know the order in which the UAV scans regions during the allocation phase, the flying distance is roughly estimated using the distance between the centers of these two regions.TS i,j is used to represent the time consumption of U i scanning the region R j ; each UAV will use the back-and-forth scanning method to scan a single area.TF i,j,k and TS i,j can be calculated using Drones 2023, 7, 664 5 of 20

Exact Formulation
This section analyzes the constraints of the multi-area coverage task-allocation problem of heterogeneous UAVs, and mathematical formulations for each constraint are established.Furthermore, an optimization model for task allocation is established, aiming to minimize the overall completion time of the system.
To obtain appropriate allocation results, the main constraints in this study are as follows: C1: The constraint on the takeoff and landing of UAV.Once a UAV is selected to perform the scanning task, it needs to depart from the starting point and return to the starting point after completing all tasks.
C2: The constraint on the scanning of regions.To avoid performing things twice, each task region should only be scanned by one UAV.It also means that only one UAV can fly to and from one mission area.This will prevent repeated scanning of the same area.
C3: The constraint on the maximum endurance of an UAV.During the execution of missions, the total flight time U i cannot exceed its maximum endurance requirement.
C4: The constraint on the number of regions covered by a single UAV.The number of scanning regions conducted by each UAV must not surpass the total number of regions.
C5: The constraint on the number of UAVs engaged in all tasks.This constraint ensures that the number of UAVs carrying out the task remains below the total quantity and is greater than one.
Constraint (C1) refers to the restriction that all UAVs can only perform takeoff and landing operations at most once.This also means that if the UAV departs from its base, it is obliged to return to the base upon completion of the coverage task.To describe constraint (C1) through mathematical expressions, this study adopts a Boolean array represent the decision variables of UAVs' planned paths, where the subscripts j, k represent U i flies from R j to R k .It is worth noting that only if U i is chosen to fly from R j to R k , Boolean variable x i,j,k = 1, otherwise, x i,j,k = 0.If j = 0 or k = 0, it represents that the U i departs from or returns to the starting point.Through the aforementioned description, constraint (C1) can be expressed as: Constraint (C2) pertains to the situation in which only one UAV is capable of traversing a given region.Therefore, constraint (C2) can be represented as: Constraint (C3) indicates that the overall flight duration of a chosen UAV, including departure from and return to the base and coverage of the specified task regions, should not exceed its maximum flight time, i.e., ∀i ∈ The variable y i,j in the above expression is a Boolean variable, which has a similar meaning with x i,j,k .If U i needs to perform a coverage scanning task in R j , y i,j = 1, otherwise y i,j = 0.The relationship between y i,j and x i,j,k can be expressed as: Drones 2023, 7, 664 6 of 20 Constraint (C4) implies that each UAV cannot perform more coverage tasks than the total number of regions, i.e., ∀i ∈ Constraint (C5) shows that the number of UAVs taking off from the bases to carry out the tasks must not exceed the total number of UAVs, which can be written as: This article aims to seek the sequence of UAV access to the task regions in the scenario where heterogeneous UAV systems take off from different bases, scan and cover corresponding areas, and return to the bases.The goal is to minimize the overall system's task completion time f , and the constraints are constraints from (C1) to (C5).Since each UAV takes off simultaneously and flies towards their respective task regions, the task completion time of the entire system can be equivalently represented as the time consumed by the UAV, which returns to the base the last.Therefore, the optimization objective function f can be expressed as: The exact formulation of this system can be written as: In the above programming, the unknown variables are a mixture of integer (e.g., elements of X) and real variables (e.g., the maximum time f ), and all constraints are linear.Therefore, this problem belongs to the class of MILP problems.Since the coveragepath planning of heterogeneous UAVs is NP-hard, although the precise solutions can be obtained using the proposed MILP formulation, it requires searching the entire solution space.Furthermore, the computational time will drastically increase with the growth of both the number of UAVs and the number of regions, resulting in significant consumption of computational time and cost.Inspired by the concept of clustering, we devised a clustering-based approach to tackle the task-allocation problem during the coverage-path planning of heterogeneous UAVs in the following sections.Implementing region allocation through clustering, and then optimizing the region access sequence, can greatly reduce the complexity of the problem.

Coverage Scanning Clustering Algorithm-Based Coverage-Path Planning
The target clustering methods normally consist of two steps: target area clustering and cluster target allocation.Target area clustering involves clustering the regions based on specific characteristics; cluster target allocation requires assigning each sub-cluster to a specific UAV.
The k-means clustering algorithm is one of the iterative and classical target clustering methods.It has the advantages of simplicity, wide applicability, and fast convergence speed.Therefore, it is frequently applied to solve multi-UAV task assignment problems [34][35][36].The basic process of the k-means clustering algorithm is commonly represented as: Step 1: Determine a value k, which represents the number of sub-clusters aiming to obtain by clustering the dataset.
Step 2: Randomly select k data points from the dataset as centroids (the centers of the sub-clusters).
Step 3: For each point in the dataset, calculate its distance to each centroid and assign it to the sub-cluster with the nearest centroid.
Step 4: Calculate the mean coordinates of all points within the sub-cluster and take it as the new centroid.
Step 5: If the distance between the newly calculated centroid and the previous centroid is less than a predefined threshold, it can consider the clustering to have achieved the desired result, and the algorithm terminates.Otherwise, we need to iterate Step 3~Step 5.
From the flow of the k-means clustering algorithm, it can be seen that it is highly sensitive to the initial selection of centroids, and different random seeds can yield completely different clustering results, significantly influencing the outcome.Moreover, the k-means clustering algorithm focuses on clustering data points, while the research background of this study involves multiple heterogeneous UAVs flying to multiple rectangular regions to perform scanning tasks.Therefore, the dataset in this study consists of planes rather than points, making the k-means algorithm unsuitable for this study.Inspired by the iterative clustering idea of the k-means algorithm, this study proposes an algorithm called the coverage-based scanning clustering algorithm (CSCA), which includes three major phases: the initial cluster centers selection phase, the multi-regional initial clustering phase, and the clustering regions update phase.Additionally, we investigated the method for optimizing the access order of regions following the completion of the clustering process.

Initial Cluster Centers Selection Phase
The purpose of clustering the various regions is to allocate each clustered subset area to its corresponding UAVs to complete the overall task in the shortest possible time.This means that each UAV will be assigned to a sub-clustered area.Afterward, the UAVs will perform scanning tasks of the regions within their respective sub-clusters in a certain order.
The task completion time includes both flight time and scanning time.The clustering centers have spatiotemporal similarities (which will be explained in the second phase) with their corresponding subsets of regions.Consequently, the duration of the UAV flight time is partially influenced by the distance between the UAV base and the clustering center.In other words, the shorter the distance between the clustering center and the UAV base, the less time the UAV will spend on flight.Therefore, we take the bases of each UAV as the initial clustering center points for the CSCA algorithm.A set C = {C 1 , C 2 , . . . ,C n } is used to store the clustering centers, where the serial number of cluster centers is equal to the serial number of UAVs involved in the given tasks.Each clustering center is defined as C i = x i , y i , where x i and y i represent the coordinates on the X-axis and Y-axis of the map, respectively.For each C i , a queue Que i = que 1 i , que 2 i , . . ., que e i is utilized to indicate the indexes of regions that have been grouped into C i and would be scanned by UAV U i ; additionally, e represents the total number of regions to be covered by U i .

Multi-Regional Initial Clustering Phase
After determining the initial cluster centers, the initial clustering methods for the regions will be proposed in this phase.A n-row and m-column matrix STS = STS i,j is employed to characterize the connection between the regions targeted for clustering and the clustering centers.STS i,j is defined as a spatiotemporal similarity between cluster C i and region R j .It consists of two parts: the time taken by U i to travel from the center R j to the cluster center C i , and the time taken by U i to coverage R j , i.e., The above equation d R j , C i refers to the Euclidean distance between the regional center R j and the cluster center C i .It is necessary to calculate the spatiotemporal similarity of each R j with all C i in the set C to find the minimum value index ind i , and add the index j corresponding to R j to the queue Que ind i , i.e.,

Clustering Regions Update Phase
This phase is the core of the proposed CSCA algorithm, where the cluster centers will continuously change with each iteration until the desired outcome is achieved.For the cluster C i , the mean coordinates of all regional centers in this cluster are defined as the updated coordinates for the cluster center.Additionally, TR(C i ) is used as an approximate substitute for the task-completion time requirements of C i , which can be written as, From the above formulation, it can be deduced that the more regions C i contains, the larger TR(C i ) will be.Consequently, this will lead to an overall increase in the total task completion time for UAV U i .Our goal in this study is to minimize the overall completion time of the entire heterogeneous UAV system, meaning that there is minimal difference in completion time between the UAVs.Thus, this can be equivalently represented as minimizing the difference between the maximum completion time max{TR(C i )} and the minimum completion time min{TR(C i )} of tasks.
In order to meet the maximum endurance constraint of UAVs, it is imperative to determine the remaining flight time RFT(U i ) for each individual UAV U i .RFT(U i ) can be calculated via, If there is any RFT(U i ) value less than 0, the TR(C i ) of the RTF(U i ) with the lowest negative value should be assigned as max{TR(C i )}.In order to reduce the regions of the cluster corresponding to max{TR(C i )} and increase the regions of the cluster corresponding to min{TR(C i )}, we proposed a region-transfer strategy based on the ranking of clustering center distance, which consists of mainly three steps.In the following steps, to enhance convenience and facilitate ease of understanding, the symbol maxind and minind are employed as representations of arg(max{TR(C i )}) and arg(min{TR(C i )}), respectively: Step 1: Calculate the distances between the cluster center C maxind and all the remaining cluster centers and sort them in ascending order to obtain the sequence Π = {π 1 , π 2 , . . . ,π n−1 }.Π(i) is equal to π i , and each π i represents the index number of a cluster center.
Step 2: Compute the distance between all regions que 1 maxind , que 2 maxind , . . ., que e maxind in the Que maxind and the cluster center C Π(1) , then assign the region index with the minimum distance to the Que Π(1) of cluster C Π(1) .
Step 3: If Π(1) = minind, it indicates the completion of the region transfer process, otherwise, define maxind as Π(1), and remove Π(1) from the sequence Π.Update Π by the distance between C maxind and other remaining cluster centers, then proceed to Step 2.
The schematic diagram of the algorithmic process described above is illustrated in Figure 3. Once the clustered regions have been updated, recalculate the centroids for each cluster and TR(C i ).The centroid coordinate of C i can be calculated by Subsequently, compute the disparity between the max{TR(C i )} and min{TR(C i )} values and ascertain whether it falls below the designated threshold.If the disparity is less than the threshold, the clustering assignment for all regions is considered complete.If not, proceed with the above region-transfer strategy.
Drones 2023, 7, 664 9 of 20 Subsequently, compute the disparity between the values and ascertain whether it falls below the designated threshold.If the disparity is less than the threshold, the clustering assignment for all regions is considered complete.If not, proceed with the above region-transfer strategy.Calculate STS i,j by Equation ( 11); 6: end for 7: end for 8: for j ← 1 to m do 9: ind i ← argmin STS i,j ; 10: Que ind i ← Que ind i ∪ j ; 11: end for 12: while the difference between max{TR(C i )} and min{TR(C i )} has not dropped below the specified threshold or the maximum number of iterations has not been reached, do 13: Calculate the TR(C i ) for each C i according to Equation ( 13) and determine max{TR(C i )} and min{TR(C i )}; 14: Update the RFT(U i ) for each U i according to Equation ( 14) and determine the minimum value min{RFT(U i )};

Regional Sorting
After executing the clustering assignment for all regions, if the number of regions in a sub-cluster is greater than one, it is necessary to sort the order of scanning regions.Therefore, the scanning order of regions assigned to UAVs will be addressed in this phase.There are many strategies that can be employed to solve this problem, such as the Genetic algorithm (GA), particle swarm optimization (PSO), and various other heuristic algorithms.However, these algorithms have the problems of non-uniqueness in calculation results and slow computation.Therefore, we adopt the nearest-to-end (NE) policy, which has significant efficiency and reliability.In this policy, all the unallocated regions in each Que i will be sorted.The policy involves two parts: initialization and interpolation.
(1) Initialization A sequence Seq i is used to denote the sequential order of tasks executed by U i , and the region in the unallocated regions that is closest to the base U pos i is chosen as the tail (seq e i ) of the sequence Seq i .In addition, U pos i is utilized as the head (seq 1 i ) of the sequence.Then, choose the unallocated region nearest to either the head or tail of the sequence, R best , which can be represented as, (2) Interpolation After initialization, R best will be stored in Seq i .If R best is closer to the tail than to the head, R best is placed in the end of Seq i and U pos i will be deleted from Seq i ; thus, R best will be the new tail; otherwise, put the R best in the first place of Seq i and delete U pos i , thus R best becomes the new head.As a result, the unallocated regions in Que i can be sorted using Equation ( 16) and placed in the Seq i by the head-tail update method described above until all regions in Que i are allocated.
Once all the regions have been sorted, the mission duration of all UAVs can be calculated effortlessly using Equations ( 1) and (2), and the maximum duration signifies the overall mission duration of the entire heterogeneous UAV system.The pseudo-codes of the coverage-based scanning clustering algorithm with a nearest-to-end policy are shown in Algorithm 2, and the algorithm flowchart is shown in Figure 4.

Algorithm 2: Pseudo-codes of coverage-based scanning clustering algorithm with Nearest-to-End Policy
Input: set of regions R, set of UAVs U Output: The regions scan sequence Seq i for each U i 1: get the set of clusters C and regional queue Que according to Algorithm 1; 2: for i ← Once all the regions have been sorted, the mission duration of all UAVs can be calculated effortlessly using Equations ( 1) and ( 2), and the maximum duration signifies the overall mission duration of the entire heterogeneous UAV system.The pseudo-codes of the coverage-based scanning clustering algorithm with a nearest-to-end policy are shown in Algorithm 2, and the algorithm flowchart is shown in Figure 4.

Bilateral Shortest-Selection Strategy-Based Trajectory Planning
In the previous section, we discussed the multi-region allocation method and obtained the regional scanning sequence for each UAV.Hence, this section will introduce a pathplanning method for the UAV flights.
The back-and-forth scanning method [21] is widely employed for regional scanning.When scanning a rectangular region R j using the back-and-forth scanning method, the scanning can be performed along either its long side or its short side.Currently, the majority of studies concentrate on scanning a single region.Consequently, scanning along the longer side indeed yields a shorter path compared to scanning along the shorter side.However, in this study, a single UAV must sequentially scan multiple regions.Its route encompasses the scanning distance within each region, as well as the flight distance between regions.As a result, scanning along the longer side does not necessarily minimize the overall distance, as it may increase the inter-region distance.Therefore, this study proposes a multi-region scanning path-planning method based on minimizing the total flight distance called the bilateral shortest-selection strategy (BSSS).
As defined earlier in this study, the scanning width of the sensors carried by U i is represented by d i .Additionally, R j has a total of 8 optional entry points {Lin, Win}.They are classified into categories {Lin 1 , . . . ,Lin 4 } and {Win 1 , . . . ,Win 4 }, depending on whether they are located on the long side or the short side, as shown in Figure 5.The red dots represent entry points on the long side, which means the UAV will scan along the shorter side of the region, while the blue dots represent entry points on the short side, which means the UAV will scan along the longer side of the region.
The current position of the UAV U i is defined as Outpos i , and the areas to be scanned are R j .The entry points InL and InW for the longer and shorter sides of region R j will be determined by calculating the closest points in {Lin 1 , . . . ,Lin 4 } and {Win 1 , .
After defining the entry points of j R , i U would fly from i Outpos to InL or InW , using the back-and-forth scanning method and scanning along either the short or long sides.Upon completion of the scanning process, the departure points OutL and OutW will be identified for each region.Then, it is necessary to separately calculate the shortest distances between OutL and OutW with the optional entry points   ' ' , Lin Win of the next to-be-scanned region.Therefore, the flying distance _ Len W and _ Len L of the aforementioned process can be represented individually as: Since the subsequent region can be scanned either along its longer or shorter edge, the last terms in Equations ( 19) and ( 20) are expressed as a weighted average.If _ _ Len L Len W  , choose InW as the entry point and scan along the longer side of the region.Otherwise, choose InL as the entry point and scan along the shorter side of the region.Finally, the current exit point of the region will be defined as new i Outpos , and we can employ the previously mentioned approach to determine both the entry point and the scanning method for the subsequent region.

Simulation Experiments
This section conducted numerical simulation experiments on the proposed CSCA algorithm and generated actual flight reference paths for UAVs based on the BSSS method.The simulation scenario was established, and subsequently, we validated the accuracy and effectiveness of our algorithm by conducting a comparative analysis with the spatial- After defining the entry points of R j , U i would fly from Outpos i to InL or InW, using the back-and-forth scanning method and scanning along either the short or long sides.Upon completion of the scanning process, the departure points OutL and OutW will be identified for each region.Then, it is necessary to separately calculate the shortest distances between OutL and OutW with the optional entry points {Lin , Win } of the next to-be- scanned region.Therefore, the flying distance Len_W and Len_L of the aforementioned process can be represented individually as: Since the subsequent region can be scanned either along its longer or shorter edge, the last terms in Equations ( 19) and ( 20) are expressed as a weighted average.If Len_L < Len_W, choose InW as the entry point and scan along the longer side of the region.Otherwise, choose InL as the entry point and scan along the shorter side of the region.Finally, the current exit point of the region will be defined as new Outpos i , and we can employ the previously mentioned approach to determine both the entry point and the scanning method for the subsequent region.

Simulation Experiments
This section conducted numerical simulation experiments on the proposed CSCA algorithm and generated actual flight reference paths for UAVs based on the BSSS method.The simulation scenario was established, and subsequently, we validated the accuracy and effectiveness of our algorithm by conducting a comparative analysis with the spatialtemporal clustering-based algorithm (STCA) cited in [31], which provided an effective method for the assignment of heterogeneous UAVs multi-area scanning tasks.Utilizing the regional allocation outcomes obtained via CSCA, the actual feasible flight routes can be obtained using the BSSS method while accounting for the constraints imposed by UAV maneuverability.

Parameter Setting
For this simulation, the mission area is defined as a square with a side length of 50 km.There are a total of eight UAVs available for the mission, with their specific parameter settings presented in Table 2.The setting of the rectangular regions needs to obey the following rules: 1.
Randomly assign coordinates to the center point of each rectangular region in order to achieve a uniform distribution across the entire map; 2.
Set the range of values for the major axis L ∈ [3, 3.5, 4] and the range of values for the minor axis W ∈ [2, 2.5, 3], with units in kilometers; 3.
Iterate through the rectangular region in order and assign values to the major and minor axes according to the index of the region.When the index is divisible evenly by 3, assign the values L(1) and W(1); when the index leaves a remainder of 1, assign the values L(2) and W(2); otherwise, assign the values L(3) and L(3); 4.
Set the rotation angle of the rectangle.Divide π into m equal parts (i.e., π/m), where m is the total number of regions.The rotation angle for each region is obtained by multiplying the index of the region by the equal parts angle (i.e., index × π/m).The unit is radians, and it is defined that the counterclockwise direction is positive; 5.
The total number of regions ranges from 5 to 40, with an increment of 5.

Task Completion Time Simulation
Under the premise of the same parameters and regional distribution of UAVs, we conducted a simulation comparison between the CSCA algorithm and the STCA algorithm in this section.The maximum task completion time of UAVs is used to represent the overall task completion time of the heterogeneous UAV system.Furthermore, we explored the reliability of the algorithm from the perspectives of different numbers of regions and different numbers of UAVs.
In Figure 6, we used two different region clustering and assignment algorithms, CSCA and STCA [31].Region sequencing methods based on the nearest-to-end policy and genetic algorithm are used to compare the applicability of the CSCA and STCA for different region sequencing methods.The number of UAVs was three, and the X-axis represents the number of regions, with a step size of 5, while the Y-axis represents the task completion time of the UAV system."STCA-NE", "CSCA-NE", "STCA-GA", and "CSCA-GA" correspondingly indicate the allocation of region-scanning sequences based on the STCA algorithm and the CSCA algorithm utilizing the nearest-to-end policy and allocation of region-scanning sequences based on the STCA algorithm and the CSCA algorithm employing the genetic algorithm.Under different region numbers, the statistical results of task completion time obtained using different algorithms are shown in Table 3.
From Figure 6, it can be observed that with the increase in the number of regions, the overall task completion time of the four methods shows an upward trend.Additionally, the CSCA-NE algorithm proposed in this study surpasses the STCA-NE algorithm.Furthermore, CSCA-GA consistently achieves the shortest task completion time, which demonstrates the outstanding performance of our clustering algorithm.Furthermore, it can be found that CSCA-GA and STCA-GA are superior to CSCA-NE and STCA-NE, respectively, due to the higher accuracy of genetic algorithms compared to the nearest-end strategy.From Table 3, it can further be observed that under different numbers of regions, when using the average task completion time as a measurement standard, CSCA-NE's result is 19.6% lower than STAC-NE, and CSCA-GA's result is 9.7% lower than STAC-GA.Moreover, the average difference between CSCA-NE and CSCA-GA is 60.6% lower than the average difference between STCA-NE and STCA-GA.This indicates that under different region quantities, the CSCA algorithm has a higher adaptability for the sequential allocation of regions in clustering.As a result, the above analysis demonstrates that under the premise of a fixed number of UAVs scanning multiple regions, the CSCA algorithm exhibits higher reliability and accuracy compared to the STCA algorithm.spectively, due to the higher accuracy of genetic algorithms compared to the nearest-end strategy.From Table 3, it can further be observed that under different numbers of regions, when using the average task completion time as a measurement standard, CSCA-NE's result is 19.6% lower than STAC-NE, and CSCA-GA's result is 9.7% lower than STAC-GA.Moreover, the average difference between CSCA-NE and CSCA-GA is 60.6% lower than the average difference between STCA-NE and STCA-GA.This indicates that under different region quantities, the CSCA algorithm has a higher adaptability for the sequential allocation of regions in clustering.As a result, the above analysis demonstrates that under the premise of a fixed number of UAVs scanning multiple regions, the CSCA algorithm exhibits higher reliability and accuracy compared to the STCA algorithm.In Figure 7, the number of regions is set to 20, and we attempt to compare the results of the four algorithms under different numbers of UAVs.To better understand the results, the statistical results of task completion time obtained using different algorithms are  In Figure 7, the number of regions is set to 20, and we attempt to compare the results of the four algorithms under different numbers of UAVs.To better understand the results, the statistical results of task completion time obtained using different algorithms are shown in Table 4.It can be seen that with the increase of UAVs, the overall task completion time of the four methods shows a downward trend.In addition, under different numbers of UAVs, the average task completion time calculated using CSCA-NE is 21.1% lower than STCA-NE, and CSCA-GA is 13.1% lower than STCA-GA.The results of CSCA-NE even outperform those of STCA-GA consistently throughout the entire process, with the exception of UAVs; the number is four.Furthermore, it is apparent that there is a minimal disparity in outcomes between CSCA-NE and CSCA-GA, with their average difference being 60.4% lower than the average difference between STCA-NE and STCA-GA.This reinforces the notion of the CSCA algorithm's remarkable adaptability in sequentially allocating regions in clustering.
outperform those of STCA-GA consistently throughout the entire process, with the exception of UAVs; the number is four.Furthermore, it is apparent that there is a minimal disparity in outcomes between CSCA-NE and CSCA-GA, with their average difference being 60.4% lower than the average difference between STCA-NE and STCA-GA.This reinforces the notion of the CSCA algorithm's remarkable adaptability in sequentially allocating regions in clustering.

Execution Time Simulation
This section further compares the differences in calculation time between algorithms based on the previous subsection.Figure 8 demonstrates the comparative results of simulation time for four algorithms as the number of regions increases from 10 to 40.Table 5 presents the statistical results of the time consumed by various algorithms for different numbers of regions.The computation time for CSCA-GA and STCA-GA is within the same order of magnitude, around 0.1 s, while the computation time for CSCA-NE and STCA-NE is in the millisecond range.In Table 5, despite STCA-NE having the shortest computation time, the maximum difference between STCA-NE and CSCA-NE is approximately 3 milliseconds.The computation time of CSCA-NE and STCA-NE is significantly lower than that of CSCA-GA and STCA-GA due to the higher complexity of GA in

Execution Time Simulation
This section further compares the differences in calculation time between algorithms based on the previous subsection.Figure 8 demonstrates the comparative results of simulation time for four algorithms as the number of regions increases from 10 to 40.Table 5 presents the statistical results of the time consumed by various algorithms for different numbers of regions.The computation time for CSCA-GA and STCA-GA is within the same order of magnitude, around 0.1 s, while the computation time for CSCA-NE and STCA-NE is in the millisecond range.In Table 5, despite STCA-NE having the shortest computation time, the maximum difference between STCA-NE and CSCA-NE is approximately 3 milliseconds.The computation time of CSCA-NE and STCA-NE is significantly lower than that of CSCA-GA and STCA-GA due to the higher complexity of GA in contrast to the nearest-end strategy.Moreover, the increase in the number of regions displays a gradual rise in computation time, suggesting a minimal impact of the number of regions on the algorithm's computation time.
Table 6 displays the time consumed by various algorithms as the number of UAVs increases from 3 to 8; Figure 9 is a graphical representation of Table 6.As can be seen in Figure 9, the execution times of CSCA-GA and STCA-GA remain highly similar, and this similarity also holds for CSCA-NE and STCA-NE.The calculation times for both CSCA-NE and STCA-NE scarcely increase with an increasing number of UAVs; however, CSCA-GA and STCA-GA exhibit a gradual increase.This occurs because the increase in UAVs results in a rise in the number of calls made to GA, which, in turn, results in an increase in time.The computation time of STCA-NE is the shortest, while the computation time of CSCA-NE is also within the same scale of magnitude as STCA-NE, with a maximum difference of 5 milliseconds between them.
Drones 2023, 7, x FOR PEER REVIEW 17 of 21 contrast to the nearest-end strategy.Moreover, the increase in the number of regions displays a gradual rise in computation time, suggesting a minimal impact of the number of regions on the algorithm's computation time.9 is a graphical representation of Table 6.As can be seen in Figure 9, the execution times of CSCA-GA and STCA-GA remain highly similar, and this similarity also holds for CSCA-NE and STCA-NE.The calculation times for both CSCA-NE and STCA-NE scarcely increase with an increasing number of UAVs; however, CSCA-GA and STCA-GA exhibit a gradual increase.This occurs because the increase in UAVs results in a rise in the number of calls made to GA, which, in turn, results in an increase in time.The computation time of STCA-NE is the shortest, while the computation time of CSCA-NE is also within the same scale of magnitude as STCA-NE, with a maximum difference of 5 milliseconds between them.
All of the above results demonstrate that, when compared to the STCA algorithm, the CSCA algorithm proposed in this work can achieve more accurate results while ensuring the computational speed at the same scale, and it has better adaptability for optimizing subsequent access sequences.All of the above results demonstrate that, when compared to the STCA algorithm, the CSCA algorithm proposed in this work can achieve more accurate results while ensuring the computational speed at the same scale, and it has better adaptability for optimizing subsequent access sequences.

Flighting-Path Simulation
In this section, the top three UAVs from Table 2 will be selected, and the number of regions will be set to 15.In order to enhance the visibility of simulation results, the range of values for the long axes of the regions is [6][7][8], while the range for the short axes is [4, 4.5, 5].Furthermore, the minimum turning radius constraint of fixed-wing UAVs is considered.In this simulation, the minimum turning radius of each UAV is set to 200 m, and the flight trajectory is optimized using the Dubins curve optimization method [37].Figure 10 shows the allocation sequence obtained using the CSCA-GA algorithm, and it displays UAVs represented by lines of different colors, with black arrows indicating their respective flight directions.

Flighting-Path Simulation
In this section, the top three UAVs from Table 2 will be selected, and the number of regions will be set to 15.In order to enhance the visibility of simulation results, the range of values for the long axes of the regions is [6][7][8], while the range for the short axes is [4, 4.5, 5].Furthermore, the minimum turning radius constraint of fixed-wing UAVs is considered.In this simulation, the minimum turning radius of each UAV is set to 200 m, and the flight trajectory is optimized using the Dubins curve optimization method [37].Figure 10 shows the allocation sequence obtained using the CSCA-GA algorithm, and it displays UAVs represented by lines of different colors, with black arrows indicating their respective flight directions.The flight trajectories of the UAVs obtained using the BSSS algorithm and long edge scanning strategy (in this study called LESS) are depicted in Figure 11a and Figure 11b, respectively.Upon comparing the two figures, it becomes evident that the UAV employs distinct scanning methods for the area enclosed by the red dashed circle in each figure.The flight trajectories of the UAVs obtained using the BSSS algorithm and long edge scanning strategy (in this study called LESS) are depicted in Figures 11a and 11b, respectively.Upon comparing the two figures, it becomes evident that the UAV employs distinct scanning methods for the area enclosed by the red dashed circle in each figure.Unlike LESA, BSSS does not carry out scanning along the longer side.In contrast to scanning along the longer side, BSSS achieves a reduction of approximately 8 km in total distance compared to LESA.As a result, this also substantiates the superiority of the proposed BSSS, as presented in this study.The flight trajectories of the UAVs obtained using the BSSS algorithm and long edge scanning strategy (in this study called LESS) are depicted in Figure 11a and Figure 11b, respectively.Upon comparing the two figures, it becomes evident that the UAV employs distinct scanning methods for the area enclosed by the red dashed circle in each figure.Unlike LESA, BSSS does not carry out scanning along the longer side.In contrast to scanning along the longer side, BSSS achieves a reduction of approximately 8 km in total distance compared to LESA.As a result, this also substantiates the superiority of the proposed BSSS, as presented in this study.

Conclusions
This study researched the multi-region coverage path-planning problem to identify appropriate region coverage schemes for UAVs with varying flight speeds and scanning widths.The main conclusions are as follows: 1.A novel clustering algorithm was proposed to assign regions to UAVs, which effectively reduces the computational complexity during the task allocation process; 2. A regional iterative strategy was designed to ensure the balancing of task completion time among UAVs;

Conclusions
This study researched the multi-region coverage path-planning problem to identify appropriate region coverage schemes for UAVs with varying flight speeds and scanning widths.The main conclusions are as follows: 1.
A novel clustering algorithm was proposed to assign regions to UAVs, which effectively reduces the computational complexity during the task allocation process; 2.
A regional iterative strategy was designed to ensure the balancing of task completion time among UAVs; 3.
A path-planning method was devised to select the entry points and scanning modes of the regions by considering the shortest flight distance, offering a reference trajectory for the actual flight of UAVs; 4.
The simulation results demonstrated that the clustering algorithm (CSCA) proposed in this study surpasses others in terms of both quality and flexibility.Compared to utilizing STCA for clustering, the average task completion time of CSCA has decreased by 9-21%.In addition, the discrepancy between CSCA-NE and CSCA-GA was approximately 60% lower than the discrepancy between STCA-NE and STCA-GA.Furthermore, the path-planning method (BSSS) proposed in this study can yield shorter flight paths.
Although the method proposed in this study has the capability of swiftly achieving results within a brief timeframe, we also want to enhance the algorithm's accuracy by increasing the coupling between subproblems, which will be the primary focus of our future research.

Figure 1 .
Figure 1.The structure of this study.

Figure 1 .
Figure 1.The structure of this study.

Figure 2 .j
Figure 2. Task scenario diagram.In this study, each region's features are shown as , , ,j j j j j R Rpos Rbeta L W . In this

Figure 3 .
Figure 3. Illustrative diagram of region transfer.Figure 3. Illustrative diagram of region transfer.

Figure 3 .
Figure 3. Illustrative diagram of region transfer.Figure 3. Illustrative diagram of region transfer.To better understand the overall description of the coverage-based scanning clustering algorithm (indicated in Sections 3.1-3.3),the entire pseudo-codes of this algorithm are shown in Algorithm 1.

Algorithm 1 :
Pseudo-Codes of Coverage-Based Scanning Clustering AlgorithmInput: set of regions R, set of UAVs U Output: the final set of clusters C and regional queue Que 1: Take the base of each U i as the initial cluster center C i ; 2: Build Que i for each C i to store the index of regions; 3: for i ← 1 to n do 4: for j ← 1 to m do 5:

15 :
if min{RFT(U i )} < 0 do 16: Define the TR(C i ) associated with min{RFT(U i )} < 0 as the max{TR(C i )}; 17: end if 18: Obtain the sequence Π = {π 1 , π 2 , . . . ,π n−1 } according to the Step 1 of region-transfer strategy; 19: while Que minind and Que maxind have not completed the update, do 20: Calculate the distance between all regions in the Que maxind and C Π(1) ; 21: Assign the regional index que l maxind with the minimum distance to the Que Π(1) end while 29: Recalculate the centroid coordinate of each C i according to Equation (15) 30: end while Drones 2023, 7, 664 10 of 20

Figure 5 .
Figure 5. Illustrative diagram of scanning entry points in a rectangular region.

Figure 6 .
Figure 6.Comparison of task completion time achieved using the four algorithms when the region number increases from 10 to 40.

Figure 6 .
Figure 6.Comparison of task completion time achieved using the four algorithms when the region number increases from 10 to 40.

Figure 7 .
Figure 7.Comparison of task completion time achieved using the four algorithms when the UAVs number increases from 3 to 8.

Figure 7 .
Figure 7.Comparison of task completion time achieved using the four algorithms when the UAVs number increases from 3 to 8.

Figure 8 .
Figure 8.Comparison of execution time among the four algorithms when the region number increases from 10 to 40.

Figure 8 .
Figure 8.Comparison of execution time among the four algorithms when the region number increases from 10 to 40.

Figure 9 .
Figure 9.Comparison of execution time among the four algorithms when the UAV number increases from 3 to 8.

Figure 9 .
Figure 9.Comparison of execution time among the four algorithms when the UAV number increases from 3 to 8.

21 Figure 10 .
Figure 10.Diagram of UAVs executing tasks in sequential order.

Figure 10 .
Figure 10.Diagram of UAVs executing tasks in sequential order.

Figure 10 .
Figure 10.Diagram of UAVs executing tasks in sequential order.

Figure 11 .
Figure 11.Illustrative diagram of UAVs flight path: (a) UAVs flight path based on BSSS; (b) UAVs flight path based on LESS.

Figure 11 .
Figure 11.Illustrative diagram of UAVs flight path: (a) UAVs flight path based on BSSS; (b) UAVs flight path based on LESS.
. . ,Win 4 } to Outpos i , respectively, i.e., InL = argmin(d(Outpos i , Lin 1 ), . . ., d(Outpos i , Lin 4 )) i , Win 1 ), . . ., d(Outpos i , Win 4 )) (18) Figure 5. Illustrative diagram of scanning entry points in a rectangular region.The current position of the UAV i U is defined as i Outpos , and the areas to be scanned are j R .The entry points InL and InW for the longer and shorter sides of re- gion j R will be determined by calculating the closest points in   i Outpos , respectively, i.e.,

Table 2 .
Settings for UAV parameters.

Table 3 .
The statistical results (in seconds) of task completion time under different region numbers.

Table 3 .
The statistical results (in seconds) of task completion time under different region numbers.

Table 4 .
The statistical results (in seconds) of task completion time under different UAV numbers.

Table 4 .
The statistical results (in seconds) of task completion time under different UAV numbers.

Table 5 .
The statistical results (in seconds) of the time consumed using various algorithms under different region numbers.

Table 5 .
The statistical results (in seconds) of the time consumed using various algorithms under different region numbers.

Table 6 .
The statistical results (in seconds) of the time consumed by various algorithms under different UAV numbers.

Table 6 .
The statistical results (in seconds) of the time consumed by various algorithms under different UAV numbers.