Optimization Dubins Path of Multiple UAVs for Post-Earthquake Rapid-Assessment

In the last decade, with the wide application of UAVs in post-earthquake relief operations, the images and videos of affected areas obtained by UAVs immediately after a seismic event have become an important source of information for post-earthquake rapid assessment, which is crucial for initiating effective emergency response operations. In this study, we first consider the kinematic constraints of UAV and the Dubins curve is introduced to fit the shortest flyable path for each UAV that meets the maximum curvature constraint. Second, based on the actual requirements of post-earthquake rapid assessment, heterogeneous UAVs, multi-depot launching, and targets allowed access to multiple times, the paper proposes a multi-UAV rapid-assessment routing problem (MURARP). The MURARP is modeled as the multi-depot revisit-allowed Dubins TOP with variable profit (MD-RDTOP-VP) which is a variant of the team orienteering problem (TOP). Third, a hybrid genetic simulated annealing (HGSA) algorithm is developed to solve the problem. The result of numerical experiments shows that the HGSA algorithm can quickly plan flyable paths for heterogeneous UAVs to maximize the expected profit. Finally, a case study based on real data of the 2017 Jiuzhaigou earthquake in China shows how the method can be applied in a post-earthquake scenario.


Introduction
Rapid assessment after a catastrophic event, such as earthquake, is crucial to initiating effective emergency operations [1]. Rapid assessment after an earthquake in particular is a serious and difficult problem involving quickly and accurately obtaining disaster area information, updating the post-earthquake information database, and shortening the post-earthquake black box period. These efforts can provide decision support for theory-based and effective development of relief operations to reduce the losses associated with earthquake disasters. The primary purpose of rapid assessment is to obtain a general understanding of the effects of the earthquake in the shortest time and to determine the degree of damage in different directions, particularly the location of the heavy disaster area. It is an effective way to obtain the damage degree of buildings in different directions around the epicenter.
we design a hybrid genetic simulated annealing (HGSA) algorithm to solve it and the flyable path of each UAV can be obtained quickly. The main innovations of this paper are as follows: (1) The MURARP proposed in this paper can be regarded as a variant of the TOP in the real application of multi-UAV for post-earthquake rapid assessment task, in which UAVs equipped with sensors access a group of potential targets, such that the total expected profit can be maximized. In the optimization process of the MURARP, various factors, such as the UAV endurance, dynamic constraints, multi-depot, the weight of potential targets, the heterogeneity of the UAV swarm and sensor errors, are considered. In this paper, the MURARP is modeled as a multi-depot revisit-allowed Dubins TOP with variable profit (MD-RDTOP-VP). The similarity between this model and the TOP model is that it is necessary to select the appropriate potential targets and determine the sequence of access. The differences are reflected as follows: (i) Multiple depots-In the classic TOP, all UAVs depart from the same start depot and return to the same end depot after completing the task. However, a UAV carrying out the rapid-assessment task in MD-RDTOP-VP needs to start from one of multiple depots and return to its departure depot after completing the task. (ii) The objective function-The optimization objective of the classic TOP is to maximize the total profit, while the optimization objective of MD-RDTOP-VP is to maximize the expected profit, which is calculated by multiplying the potential target's weight by the probability that valid information about that target was successfully collected. (iii) Target revisit-allowed strategy-The classic TOP model stipulates that each potential target can only be accessed at most once. However, in MD-RDTOP-VP, to increase the probability that valid information about a target was successfully collected, a UAV is allowed to visit a potential target more than once. (2) In this paper, the influence of the minimum turning radius of the UAV and heading angle discretization on the completion quality of the rapid assessment is studied. The instances of MD-RDTOP-VP are generated by modifying the classic benchmark instances of the TOP model. The numerical experiments are used to prove the influence of the minimum turning radius of the UAV on the completion quality. Then, the influence of the heading angle discretization is analyzed through simulation experiments. The results show that the high-quality MD-RDTOP-VP solution can be accelerated by selecting the appropriate heading angle discretization. In previous reports on the UAV path planning using the Dubins curve, most studies have assumed that the initial heading angle of the UAV is specified, but we also optimize that parameter. (3) This paper designs an HGSA algorithm for solving the MD-RDTOP-VP model. Since the MURARP proposed in this paper is a novel problem, the existing algorithms cannot be directly applied here. Therefore, a hybrid algorithm is developed to solve this problem. The performance of the algorithm is shown through numerical and simulation experiments. Finally, a case study based on real-world data of the 2017 Jiuzhaigou earthquake in China shows how to use our method in a post-earthquake scenario. The results show that the HGSA algorithm can obtain a high-quality feasible solution of the MD-RDTOP-VP model in a short time, which can meet the actual demand of quick path planning for each UAV in a post-earthquake scenario.
The remainder of this paper is structured as follows. In Section 2, the mathematical model of the MURARP is given. The HGSA algorithm for solving the MURARP mathematical model is discussed in Section 3. In Section 4, the performance of the HGSA algorithm is evaluated through numerical and simulation experiments, and a case study is also presented. Finally, a concluding summary of this paper and the potential future work are discussed in Section 5.

Problem Modeling
In this section a formulation of the MURARP is presented. A directed graph G = (T ∪ D ∪ D , A) is given, where T is the potential targets set, D is the set of UAV launching depots, and D is the set of UAV landing depots. The arc set A consists of all possible connections between any two vertices in T Appl. Sci. 2020, 10, 1388 5 of 24 and all directed arcs connecting all vertices in D ∪ D and all vertices in T. A group of heterogeneous UAVs, represented by U, launched from their associated launching depots and return to their associated landing depots after visiting a subset of T. Each UAV is constrained by its endurance. When a UAV accesses a potential target, the profit associated with the weight of the target is successfully obtained with a certain probability. The optimization objective of the MURARP is to maximize the expected profit. Table 1 defines the major notation of relevant sets, indices, parameters, and variables of MURARP. Table 1. List of major notations of sets, indices, parameters, and variables of MURARP.
R k feasible route of UAV k, a configuration sequence of vertices visited by UAV k, Scheme feasible routing scheme for the MURARP, ρ k turning radius of UAV k, k ∈ U.
T k max endurance of UAV k, k ∈ U.
w i weight of a potential target i, i ∈ T.
p k error probability of the sensor carried by UAV k, p k ∈ (0, 1). P i probability of successfully obtaining a profit for visiting potential target i, i ∈ T. d k ij distance between vertices i and j under the turning radius of UAV k, i, j ∈ V,k ∈ U.
t k ij travel time between vertices i and j for UAV k, i, j ∈ V,k ∈ U.
y k i times that potential target i is visited by UAV k, i ∈ T,k ∈ U.
x k ij binary variable representing whether UAV k travels from point i and j, i, j ∈ A, k ∈ U.

Heterogeneous UAVs
Let U ={1, · · · , k, · · · , K} be a set of K UAVs which are heterogeneous. The UAV heterogeneity is manifested in terms of different airspeeds v k , different turning radii ρ k , different endurances T k max and different detection errors of the carried sensors p k . Since UAVs are relatively scarce disaster-relief resources, we assume that each rescue team has only one UAV and that each UAV is launched from a different depot. In reality, most UAVs are equipped with automatic obstacle avoidance devices and stabilizing devices. Although the wind gusts and obstacles will affect the UAVs' motion trajectory, making the UAVs' actual flight trajectory slightly different from the pre-planned path, it will not cause the UAVs to fail to perform the rapid assessment task. Therefore, we also assume that the UAVs have the ability to avoid obstacles automatically and can perform rapid assessment tasks safely, and the resulting path deviation from the total flight path length is very small and can be ignored.
In this paper, the UAV motion state is described by the Dubins car model [9]. The state of a Dubins car is represented by q = (x, y, θ),which consists of its position in the plane (x, y) and its heading angle θ. One of the specifics of the Dubins car model is that the minimal turning radius ρ influences the length of the shortest path between two spatial configurations. The main purpose of modeling UAV with the Dubins car model in this paper is to qualitatively analyze whether the turning radius and the heading angle have an impact on the UAV path planning. The kinematic model of a Dubins car with a constant velocity v and a control input u ∈ [−1, 1] can be described as

Potential Targets
Typical potential target of a post-earthquake rapid-assessment task is a building whose size is generally smaller than the UAV's vision field. The sensor of a UAV can completely cover the target when it flies over the building, so the potential target is defined as the point-target. Let T ={1, · · · , i, · · · , N} be the set of N potential targets with known spatial configuration. Since the Dubins curve is used to fit the UAV path in this paper, the spatial configuration of a potential target can be defined as (x i , y i , θ i ).
Let w i ∈ {1, 2, . . . , 10} be the weights of potential targets. A larger weight value means that the target is more important, and the higher the weight of the target visited by the UAV, the greater is the profit. Due to the detection error of the sensor, when the UAV accesses the potential target i ∈ T, the valid information about the target can be successfully obtained with a certain probability, which is recorded as P i . The probability P i is related to the following three factors: i) the weight of the potential target i, w i ; ii) the error probability of the sensor carried by UAV k, p k ; and iii) the number of times that the potential target i is visited by UAV k, y k i . The formula of P i is

Feasible UAV Path
Dubins curves are widely used to generate the feasible paths of minimal length for a UAV between any two vertices. Since the length of shortest feasible path between any two vertices is related to the UAV's minimum turning radius, the path length of the heterogeneous UAVs is usually different between the same two points. Reference [20] demonstrated that the sensors on UAVs can capture an area with a width of at least 2ρ, so it is a reasonable assumption that any two targets are at least 2ρ apart. There are two kinds of Dubins curves, the first one is with a terminal heading constraint and the other is without any terminal heading constraint [35]. The potential targets studied in this paper have no terminal heading constraints. Therefore, by optimizing the UAV heading angle when accessing each target, the UAV path length can be shortened, resulting in a savings of UAV energy to access more targets and in turn improve the expected profit. Figure 1 shows the heading angle discretization when the UAV accesses a potential target, with the heading angle discretization set defined as where N θ > 0 is an integer defining the heading angle dispersion rate. For example, for N θ = 8 of Figure 1, the discretization set is where 0 N θ > is an integer defining the heading angle dispersion rate. For example, for 8 N θ = of A feasible path of UAV k is an ordered set, Rk, of at least . In the scheme shown in Figure 2, UAV1 starts from D1 at a 90° heading angle, first accesses the #1 target at a 90° heading angle, then accesses the #4 target at a 225° heading angle, and finally returns to D1 at a 315° heading angle; UAV2 starts from D2 at a 315° heading angle, first accesses the #2 target at a 315° heading angle, then accesses the #5 target at a 225° heading angle, accesses the #3 target at a 135° heading angle, and finally returns to D2 at a 45° heading angle. The scheme can be expressed as ) ) ) ) ) ) ) , (3,135 ,45 } ) ,( target appears twice in Scheme, which means that the #5 target is revisited. A feasible path of UAV k is an ordered set, R k , of at least r + 2(r ≥ 1) vertices of the form . . , r, corresponds to the index of r potential targets being accessed in that sequence by UAV k. A feasible multi-UAV rapid-assessment routing scheme is composed of the flight paths of all UAVs, Scheme = {R 1 ; . . . ; R k ; . . . ; R K }, k ∈ U. In the scheme shown in Figure 2, UAV1 starts from D 1 at a 90 • heading angle, first accesses the #1 target at a 90 • heading angle, then accesses the #4 target at a 225 • heading angle, and finally returns to D 1 at a 315 • heading angle; UAV2 starts from D 2 at a 315 • heading angle, first accesses the #2 target at a 315 • heading angle, then accesses the #5 target at a 225 • heading angle, accesses the #3 target at a 135 • heading angle, and finally returns to The #5 target appears twice in Scheme, which means that the #5 target is revisited.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 25 According to the research results of Pham, T.A. [36], the revisiting modes are divided into restricted node revisiting and node revisiting without any restrictions. Considering the actual needs of the task completion time of rapid assessment, this paper adopts the first revisiting method, i.e., to revisit target i, UAV k has to visit another target before it can return to target i. In this paper, we take the constraint of an overnight stay being forbidden, which refers to the naming of the first revisiting mode in [36], as an important constraint of the MURARP model.  According to the research results of Pham, T.A. [36], the revisiting modes are divided into restricted node revisiting and node revisiting without any restrictions. Considering the actual needs of the task completion time of rapid assessment, this paper adopts the first revisiting method, i.e., to revisit target i, UAV k has to visit another target before it can return to target i. In this paper, we take the constraint of an overnight stay being forbidden, which refers to the naming of the first revisiting mode in [36], as an important constraint of the MURARP model.
The feasible path of UAV k consists of the following three parts: (i) the path from UAV's configuration at launching depot q D k to the entry configuration at the first target q k target V k l+1 , l = 1, 2, . . . , r − 1; and (iii) the path from the exit configuration q V k r of the last target V k r to the UAV's configuration at landing depot q D k . The total path length, d R k , associated with the feasible path R k of UAV k can be calculated as In this paper, we assume that each UAV airspeed v k is constant, so the task completion time of UAV k, can be calculated as

Mathematical Model
To solve the MURARP, the objective function is designed to maximize the total expected profit of the routing scheme. The mathematical model of MURARP is formulated as Equation (6) is the objective function, corresponding to maximizing the expected profit of the multi-UAV post-earthquake rapid-assessment routing scheme, where w i is the weight of target i, p k is the error probability of the sensor carried by UAV k, y k i is times that target i is visited by UAV k. Equation (7) ensures that all UAV must departure from their own depot and ultimately return to the same depot. Equation (8) ensures the connectivity of each UAV's route. Equation (9) is the endurance constraint, that is, the task execution time of the UAV k cannot exceed its maximum safe flight duration. Equations (10) and (11) are the constraints of an overnight stay being forbidden, thus ensuring that the UAV cannot stay at the same potential target, and a target can be revisited only after the UAV visits another target. Equation (12) involves the values of the binary decision variable whether UAV k flies to target j after leaving target i; if so, x k ij = 1, and x k ij = 0 otherwise.

HGSA Algorithm
A hybrid heuristic algorithm for efficiently solving the MURARP will be presented in this section. The GA has excellent ability to find new solutions quickly, while its hill climbing ability is poor and easily falls into local optimum in later iterations. It has been proved that simulated annealing (SA) can converge to the global optimum with a probability of approximately 1 by controlling the cooling process [37]. However, the convergence rate of SA is too slow to meet actual needs. In this paper, we specially develop a two-stage hybrid heuristic algorithm for the characteristics of the MURARP. In the first stage, a better solution is obtained quickly by GA, which is used as the initial solution of the SA algorithm. In the second stage, the optimal solution of the MURARP is obtained through SA operation. To balance the convergence rate and global optimization ability of the HGSA algorithm, a switching criterion from the GA to SA is adaptively designed. The pseudocode of the algorithm is shown in Algorithm 1.
Algorithm 1: HGSA Algorithm 1: Input: P, a population of N P chromosomes; 2: Output: C. best, a chromosome with the best fitness found; 3: begin 4: initialize the parameters of the HGSA algorithm; 5: initialize the population P of the GA (see Section 3.1.2); 6: constraint checking and adjustment for each chromosome in the population P (see Section 3.1.3); 7: evaluate the fitness of the population, and update C. best (see Section 3.1.4); 8: while the switch point between the GA and SA is not reached (see Section 3.5); 9: crossover (see Section 3.2); 10: mutation (see Section 3.3); 11: update population (see Section 3.4); 12: endwhile 13: T iter =T 0 ; 14: while T iter > T end 15: improve C. best through SA iteration (see Section 3.6); 16: T iter = T iter * R T ; 17: endwhile 18: the chromosome with the best fitness C. best is selected as the final scheme; 19: end

Initialization
The key parameters of the HGSA algorithm are listed in Table 2. We do not set the maximum iteration of the GA because an automatic switching mechanism from the GA to SA is designed. When the optimization ability of the GA shows a decay trend, the algorithm automatically switches to SA to avoid falling into a local optimum.

Integer-Encoded Chromosome with Variable Length
A chromosome represents a feasible routing scheme for the MURARP. Since the MURARP allows the target to be revisited, each scheme may have a different number of access vertices, i.e., the length of each chromosome may be different. In addition, the MURARP has two other features: (i) The UAVs start from different stations, and (ii) the UAVs have multiple optional heading angles when accessing each target. To facilitate the operation of the HGSA algorithm, this paper designs integer-encoded chromosomes with variable lengths consisting of three lines, in which the first line is the index of visited vertices accessed by the UAVs, the second line is the UAV indices of the heading angle at each vertex, and the third line is the index of UAVs corresponding to the first line. The chromosome shown in Figure 3 corresponds to the UAV routing scheme shown in Figure 2. 3.1.1. Integer-Encoded Chromosome with Variable Length A chromosome represents a feasible routing scheme for the MURARP. Since the MURARP allows the target to be revisited, each scheme may have a different number of access vertices, i.e., the length of each chromosome may be different. In addition, the MURARP has two other features: (i) The UAVs start from different stations, and (ii) the UAVs have multiple optional heading angles when accessing each target. To facilitate the operation of the HGSA algorithm, this paper designs integer-encoded chromosomes with variable lengths consisting of three lines, in which the first line is the index of visited vertices accessed by the UAVs, the second line is the UAV indices of the heading angle at each vertex, and the third line is the index of UAVs corresponding to the first line. The chromosome shown in Figure 3 corresponds to the UAV routing scheme shown in Figure 2. The generation method of integer-encoded chromosomes is introduced in the next section.

Population Initialization
Complete the population initialization by the following five steps: Step 1: Draw a " max k T circle" with the depot of UAV k, Dk, as the center and the endurance of UAV k, max k T , as the radius. Then delete the indices of potential targets outside the max k T circle from the target set T and obtain set Tk; Step 2: Randomly arranged the indices in set Tk to obtain a random array. Then add '0' to the front and back of the random array to obtain a temporary routing scheme of UAV k, denoted as Ak; Step 3: Randomly generate the same number of heading angle indices according to the number of elements in Ak, place them in the second row of Ak, and place UAV index, k, in the third row of Ak, then the routing scheme of UAV k, denoted as Rk is obtained; Step 4: Repeat steps 2 and 3 according to the number of UAVs, K, until all UAVs obtain a routing scheme. An initial chromosome is obtained by combining all Rk by row.
Step 5: Repeat steps 2-4 according to the population size to obtain the initial population.
It is not necessary that all chromosomes of the initial population satisfy all constraints of the MURARP model. Therefore, constraint checks should be performed on the population and the chromosomes that do not satisfy the constraints must be adjusted.

Constraint Check and Adjustment
There are two constraints in the MURARP model: the UAV endurance constraint represented by Equation (9) and the overnight stay forbidden constraint represented by Equations (10) and (11). In this paper, chromosomes that do not satisfy Equation (9) are called class A violators, chromosome that do not satisfy Equations (10) and (11) are called class B violators, and chromosomes that do not satisfy Equations (9)-(11) are called class C violators.
In view of the above three classes of violator chromosomes, three chromosome adjustment strategies are designed as follows:  The generation method of integer-encoded chromosomes is introduced in the next section.

Population Initialization
Complete the population initialization by the following five steps: Step 1: Draw a "T k max circle" with the depot of UAV k, D k , as the center and the endurance of UAV k, T k max , as the radius. Then delete the indices of potential targets outside the T k max circle from the target set T and obtain set T k ; Step 2: Randomly arranged the indices in set T k to obtain a random array. Then add '0' to the front and back of the random array to obtain a temporary routing scheme of UAV k, denoted as A k ; Step 3: Randomly generate the same number of heading angle indices according to the number of elements in A k , place them in the second row of A k , and place UAV index, k, in the third row of A k , then the routing scheme of UAV k, denoted as R k is obtained; Step 4: Repeat steps 2 and 3 according to the number of UAVs, K, until all UAVs obtain a routing scheme. An initial chromosome is obtained by combining all R k by row.
Step 5: Repeat steps 2-4 according to the population size to obtain the initial population.
It is not necessary that all chromosomes of the initial population satisfy all constraints of the MURARP model. Therefore, constraint checks should be performed on the population and the chromosomes that do not satisfy the constraints must be adjusted.

Constraint Check and Adjustment
There are two constraints in the MURARP model: the UAV endurance constraint represented by Equation (9) and the overnight stay forbidden constraint represented by Equations (10) and (11). In this paper, chromosomes that do not satisfy Equation (9) are called class A violators, chromosome that do not satisfy Equations (10) and (11) are called class B violators, and chromosomes that do not satisfy Equations (9)-(11) are called class C violators.
In view of the above three classes of violator chromosomes, three chromosome adjustment strategies are designed as follows: After the constraint checking and adjustment on the population, all chromosomes are feasible routing schemes for MURARP.

Fitness Evaluation
The fitness of a chromosome represents the expected profit of a routing scheme. According to the objective function of the MURARP, Equation (6), the higher the expected profit, the better the routing scheme. However, UAVs have multiple optional heading angles when they access a target, and chromosomes with the same fitness may have different task completion times. The routing scheme shown in Figure 4a and the routing scheme shown in Figure 4b are exactly the same in both the targets to be accessed and the number of times that targets are accessed, so the expected profits of the two schemes are also the same. However, the total flight time of the routing scheme shown in Figure 4b is shorter, which helps make room in UAV routing, thus allowing UAVs to visit more potential targets [39]. According to the above analysis, the quality of the routing scheme is related to both the total expected profit and the total flight time. Therefore, the paper designs the chromosome evaluation mechanism of a double fitness function, that is, the total expected profit fitness function Fit1 and the total flight time fitness function Fit2. The calculation formulas are The value of Fit1 is the expected profit of the routing scheme represented by the chromosome. The bigger value of Fit1, the higher is the fitness of the chromosome. The value of Fit2 is the total flight time of the routing scheme represented by the chromosome. The smaller the value of Fit2, the higher is the fitness of the chromosome. In this paper, Fit1 is used as the main fitness function for crossover and SA operations, and Fit2 is used to evaluate two chromosomes with the same expected profit, i.e., the values of Fit1 are the same. Fit2 is adopted in the mutation operation.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 12 of 25 Strategy 3: For class C violators, keep only one gene with the same target indices and delete others first, and then randomly delete a gene and repeat this deletion operation until the chromosome is no longer a violator.
After the constraint checking and adjustment on the population, all chromosomes are feasible routing schemes for MURARP.

Fitness Evaluation
The fitness of a chromosome represents the expected profit of a routing scheme. According to the objective function of the MURARP, Equation (6), the higher the expected profit, the better the routing scheme. However, UAVs have multiple optional heading angles when they access a target, and chromosomes with the same fitness may have different task completion times. The routing scheme shown in Figure 4a and the routing scheme shown in Figure 4b are exactly the same in both the targets to be accessed and the number of times that targets are accessed, so the expected profits of the two schemes are also the same. However, the total flight time of the routing scheme shown in Figure 4b is shorter, which helps make room in UAV routing, thus allowing UAVs to visit more potential targets [39]. According to the above analysis, the quality of the routing scheme is related to both the total expected profit and the total flight time. Therefore, the paper designs the chromosome evaluation mechanism of a double fitness function, that is, the total expected profit fitness function Fit1 and the total flight time fitness function Fit2. The calculation formulas are  The value of Fit1 is the expected profit of the routing scheme represented by the chromosome. The bigger value of Fit1, the higher is the fitness of the chromosome. The value of Fit2 is the total flight time of the routing scheme represented by the chromosome. The smaller the value of Fit2, the higher is the fitness of the chromosome. In this paper, Fit1 is used as the main fitness function for crossover and SA operations, and Fit2 is used to evaluate two chromosomes with the same expected profit, i.e., the values of Fit1 are the same. Fit2 is adopted in the mutation operation.

Crossover
The crossover operation is performed through the classic roulette wheel method, which is based on the selection of the best fit solutions. Thus, the higher the Fit1 value of the chromosome, the higher the probability is that its superior gene will be inherited by the next generation. Because this paper adopts an integer-coded chromosome with variable length, traditional GA single point crossover and double point crossover may lead to the failure of crossover operation. Therefore, a segmented single point stitching crossover method is specifically designed with the following specific operation steps: Step 1: Select two chromosomes to be crossed from the parent population: parent A and parent B.
Divide parent A and parent B into K segments according to the third row of the chromosome, with each segment represents a routing scheme for a UAV. For example, the kth segment of parent A, denoted by , represents the routing scheme of UAV k in parent A.
Step 2: Stitch the chromosomes' segments with the same index in parent A and parent B, denoted as Ak and Bk, respectively. The stitching rules are as follows: (i) if there is only one shared target index, g, in Ak and Bk, then the genes behind g are exchanged to obtain the new segments Ck and Dk; (ii) if there is more than one index in Ak and Bk, randomly select one index, g, and then exchange the genes behind g to obtain the new segments Ck and Dk; (iii) if there is no shared index g in Ak and Bk, then stitch Bk to the end of Ak to obtain a new segment Ck, and stitch Ak to the end of Bk to obtain a new segment Dk. A stitching example is illustrated in Figure 5.

Crossover
The crossover operation is performed through the classic roulette wheel method, which is based on the selection of the best fit solutions. Thus, the higher the Fit1 value of the chromosome, the higher the probability is that its superior gene will be inherited by the next generation. Because this paper adopts an integer-coded chromosome with variable length, traditional GA single point crossover and double point crossover may lead to the failure of crossover operation. Therefore, a segmented single point stitching crossover method is specifically designed with the following specific operation steps: Step 1: Select two chromosomes to be crossed from the parent population: parent A and parent B.
Divide parent A and parent B into K segments according to the third row of the chromosome, with each segment represents a routing scheme for a UAV. For example, the kth segment of parent A, denoted by A k , represents the routing scheme of UAV k in parent A.
Step 2: Stitch the chromosomes' segments with the same index in parent A and parent B, denoted as A k and B k , respectively. The stitching rules are as follows: (i) if there is only one shared target index, g, in A k and B k , then the genes behind g are exchanged to obtain the new segments C k and D k ; (ii) if there is more than one index in A k and B k , randomly select one index, g, and then exchange the genes behind g to obtain the new segments C k and D k ; (iii) if there is no shared index g in A k and B k , then stitch B k to the end of A k to obtain a new segment C k , and stitch A k to the end of B k to obtain a new segment D k . A stitching example is illustrated in Figure 5. The single point stitching crossover method can preserve the superior genes of the chromosome as much as possible and increase the diversity of the population, which is helpful to find a new solution. It is not necessary that all chromosomes of the child population after crossover operation satisfy all constraints of the MURARP model. Therefore, constraint checks should be performed on the child population and the chromosomes that do not satisfy the constraints must be adjusted. See Section 3.1.3 for details. After constraint checking and adjustment, the chromosomes are a feasible routing scheme of the MURARP.

Mutation
For chromosomes with the same expected income and different total flight times as shown in Figure 4, we designed a heading angle index mutation operator for the second line of the chromosome. The pseudocode of the mutation operator is shown in Algorithm 2. The single point stitching crossover method can preserve the superior genes of the chromosome as much as possible and increase the diversity of the population, which is helpful to find a new solution. It is not necessary that all chromosomes of the child population after crossover operation satisfy all constraints of the MURARP model. Therefore, constraint checks should be performed on the child population and the chromosomes that do not satisfy the constraints must be adjusted. See Section 3.1.3 for details. After constraint checking and adjustment, the chromosomes are a feasible routing scheme of the MURARP.

Mutation
For chromosomes with the same expected income and different total flight times as shown in Figure 4, we designed a heading angle index mutation operator for the second line of the chromosome. The pseudocode of the mutation operator is shown in Algorithm 2. Input: a chromosome C; 2: Output: a new chromosome C m with no worse fitness than C; 3: begin 4: calculate the total flight duration fitness of C using Equation (14), denoted as Fit2 C ; 5: calculate the number of columns in C, denoted as N C ; 6: for i=1 to 0.1 × N C 7: randomly select a gene in C, denoted as C g ; 8: obtain the index of the heading angle of C g , denoted as h Cg (h Cg = the second line of C g ); 9: randomly generate an index of the heading angle, denoted as h r ; 10: while h r = = h Cg 11: randomly generate an index of the heading angle, denoted as h r ; 12: endwhile 13: replace the index of the heading angle at the location of h Cg in C with h r ; 14: endfor 15: obtain a mutated chromosome, denoted as C m ; 16: calculate the total flight duration fitness of C m using Equation (14), denoted as Fit2 Cm ; 17: if Fit2 Cm > Fit2 C 18: Since the mutation operation is used only for the heading angle, the targets accessed by the UAV and the access sequence are not changed, so there is no influence on the value of Fit1. Furthermore, only when Fit2 is worth raising, that is, only the variation is small, is the variation is accepted so that the fitness of the new chromosome obtained by the mutation operation is not worse than that of the original chromosome.

Population Update
The population update operation is accomplished by replacing some of the chromosomes in the newly generated offspring with the same number of elite chromosomes in the parent population. The formula for calculating the number of replacing chromosomes is N F = N P × (1 − Gap), where N p is the population size and Gap is the generation gap, which are both preset. A new population for the next iteration of the GA is obtained by using the chromosomes of the pre-N F position in the parent population to replace the chromosomes of the post-N F position in the child population. In addition, perform a fitness evaluation of the optimal chromosome C. Nbest in the new population and the optimal chromosome C. best in the current population. C. best is replaced by C. Nbest only if C. Nbest has a higher fitness.

Adaptive Switching from the GA to SA
For a two-phase hybrid GA and SA algorithm, the switching point between the GA and SA is crucial and is affected by the number of iterations. Setting the switch point too early will reduce the quality of the solution. However, if the switching point is too late, the efficiency of the algorithm will be affected [32]. In general, the GA must run a fixed number of generations which depends mainly on the users' experience or preference rather than adaptive self-computation. Because of the strong optimization ability of the GA in the early stage, the growth rate of the fitness value is faster, but in the later stage, the GA exhibits a decaying rate of increase. Therefore, in this paper, the switching point between the GA and SA is set when the growth rate of the Fit1 value of C. best decreases for 10 generations continuously. When the adaptive switching condition from the GA to SA is satisfied, the GA is terminated, and C. best is used as the initial solution of SA.

Simulated Annealing
The core idea of SA is to generate multiple neighborhood solutions by disturbing the current solutions. The algorithm needs to maintain a suitable balance between disturbing the existing solutions and keeping the good parts of the existing solutions. After the iterative optimization of the GA, we consider that the initial solution C. best has found a better UAV routing scheme that selects the appropriate targets to access and that determines the UAV access order and heading angle when accessing these targets. The main goal of SA is to fine-tune C. best to increase the expected profit for the routing scheme by adding or replacing a certain access target. In this paper, a disturbance strategy based on the profit distance ratio R rd is designed with the following calculation where ∆w i is the increase in the profit and ∆d i represents the additional distance after adding target i to the routing scheme. In this paper, two disturbance strategies are designed for the MURARP: Disturbance strategy 1: Insert a target. Randomly select a gene site, and then find a target to maximize R rd after insertion, which means inserting a target with the highest possible weight at the lowest possible cost.
Disturbance strategy 2: Replace a target. Randomly select a gene site, and then find a target to maximize R rd after replacing the current target.
The main purpose of these disturbance strategies is to improve the expected profit of the routing scheme by maximizing the additional profit per distance. One of the above two strategies will be randomly used to generate a new solution in each annealing operation, and then the operations described in Section 3.1.3 should be performed on the new solution. Finally, the metropolis criterion is used to accept the new solution with a certain probability.

Experiments and Discussions
Since the MD-RDTOP-VP model is a generalization of the classic TOP, the former can be simplified to the traditional TOP model by special processing, and the model and algorithm of this paper can be tested by using the TOP's benchmark instances. In Section 4.1, we use the classic TOP's benchmark instances to carry out two numerical experiments, in which the first experiment shows the impact of the UAV turning radius on the collected profit and the second experiment shows the impact of heading angle discretization. In Section 4.2, according to the characteristics of the post-earthquake rapid assessment, we generate 42 simulation instances by transforming six sets of TOP benchmark instances. Then, the MD-RDTOP-VP model is used to solve these instances by the HGSA algorithm, and the relevant analysis is carried out. Finally, we conducted a case study based on real data as described in Section 4.3 to illustrate the application of our method in real-life scenarios.
All experiments were carried out using a desktop computer with a single core of i5-6500 3.2 GHz CPU. Since the performance of the algorithm can be improved by properly setting the values of parameters [38], we refer to the parameter setting methods of related literatures [32][33][34], and consider the actual requirements of the MURARP problem to set the key parameters of the HGSA algorithm.

Numerical Experiments
The TOP can be regarded as a special case of MD-RDTOP-VP while meeting the following five conditions at the same time: (i) all UAVs are launched from the same depot and return to the same destination, which may be different from the launching depot; (ii) each potential target can be accessed at most once; (iii) the sensor detection error probability is 0; (iv) all UAVs are homogeneous; and (v) the UAV kinematic constraint is omitted, i.e., the turning radius of UAV is 0. Therefore, we can use the TOP benchmark instances which is downloaded from https://www.mech.kuleuven.be/en/cib/op to test the model and algorithm proposed in this paper. We select six sets of TOP benchmark instances [40] for two numerical experiments. In the first numerical experiment, we selected Set 6 with 64 vertices to show the impact of the UAV turning radius on the collected profit. To further demonstrate the impact of heading angle dispersion on the collected profit, numerical experiment 2 was performed with six sets of TOP benchmark instances.

Numerical Experiment 1
To evaluate the influence of the turning radius of UAV on the routing scheme, we select 10 benchmark instances with a routing number of k = 2 for Set 6. The number of vertices of these benchmark instances is N = 64, and the vertex positions and weights are shown in Figure 6. collected profit, numerical experiment 2 was performed with six sets of TOP benchmark instances.

Numerical Experiment 1
To evaluate the influence of the turning radius of UAV on the routing scheme, we select 10 benchmark instances with a routing number of k = 2 for Set 6. The number of vertices of these benchmark instances in Figure 6. Since the distance between any two nodes is assumed to be greater than 2ρ in Section 2.3, where ρ is the minimum turning radius of the UAV. Since the minimum distance between any two vertices in the above figure is 2 , we select the typical turning radii  Table 3; BEST R represents the known traditional TOP optimal solutions of these 10 instances. Table 3 records the best results for 10 instances of Set 6, which was run 10 times under the same experimental conditions. Since the distance between any two nodes is assumed to be greater than 2ρ in Section 2.3, where ρ is the minimum turning radius of the UAV. Since the minimum distance between any two vertices in the above figure is √ 2, we select the typical turning radii ρ ∈ {0, 0.1, 0.3, 0.5, 0.7} to perform numerical experiments on the above data set, where ρ = 0 is a solution of the traditional TOP. We set the heading angle discretization N θ = 8 in other cases of turning radii. The results for the proposed method are presented in Table 3; R BEST represents the known traditional TOP optimal solutions of these 10 instances. Table 3 records the best results for 10 instances of Set 6, which was run 10 times under the same experimental conditions. The experimental results in Table 3 show that the HGSA algorithm can find the known optimal solution of the 10 TOP benchmark instances when ρ = 0, which demonstrates that the HGSA algorithm which is specially designed for the MD-RDTOP-VP model in this paper can also obtain the optimal solution of the classic TOP benchmark instances. We also find that for the turning radius ρ = 0.1, the optimal solution of all 14 instances is as same as that for ρ = 0. For other turning radius the collected profit decreases gradually with increasing ρ. Figure 7 shows the routing schemes of maximum profit under two different turning radii of the p6.2.e instance. solution of the 10 TOP benchmark instances when , which demonstrates that the HGSA algorithm which is specially designed for the MD-RDTOP-VP model in this paper can also obtain the optimal solution of the classic TOP benchmark instances. We also find that for the turning radius 0 .1 ρ = , the optimal solution of all 14 instances is as same as that for 0 ρ = . For other turning radius the collected profit decreases gradually with increasing ρ . Figure 7 shows the routing schemes of maximum profit under two different turning radii of the p6.2.e instance. In Figure 7a, UAV1 and UAV2 each visit 10 targets, while in Figure 7b, UAV2 visits 10 targets and UAV1 only visit nine targets. The results indicate that increasing the turning radius will lead to a longer path; therefore, the optimal routing scheme for the classic TOP violates the UAV endurance constraint. For example, in Figure 7b, if the 61# target is visited by UAV1, the path length of UAV1 will reach 17.6784, which exceeds TMAX = 17.5 of UAV1. The results of Numerical Experiment 1 show that if the impact of the turning radius on the path is not considered in the route planning, the path planned is likely to be impossible for the UAV to execute, which leads to task incompletion or even failure.

Numerical Experiment 2
To present the influence of the heading angle dispersion on the routing scheme, we select six benchmark instances with routing number K = 2 and maximum endurance Tmax = 20 in the TOP benchmark instances for numerical experiments. In each instance, the UAV turning radius is set to In Figure 7a, UAV1 and UAV2 each visit 10 targets, while in Figure 7b, UAV2 visits 10 targets and UAV1 only visit nine targets. The results indicate that increasing the turning radius will lead to a longer path; therefore, the optimal routing scheme for the classic TOP violates the UAV endurance constraint. For example, in Figure 7b, if the 61# target is visited by UAV1, the path length of UAV1 will reach 17.6784, which exceeds T MAX = 17.5 of UAV1. The results of Numerical Experiment 1 show that if the impact of the turning radius on the path is not considered in the route planning, the path planned is likely to be impossible for the UAV to execute, which leads to task incompletion or even failure.

Numerical Experiment 2
To present the influence of the heading angle dispersion on the routing scheme, we select six benchmark instances with routing number K = 2 and maximum endurance T max = 20 in the TOP benchmark instances for numerical experiments. In each instance, the UAV turning radius is set to ρ = 0.1, and the set of heading angle dispersions is N θ ∈ {3, 6,8,15,24,36,60}; thus, we obtain 42 benchmark instances. Table 4 records the best results for the 42 benchmark instances, which was run 10 times under the same experimental conditions. The optimal results of the 10 runs are listed in the R MAX column and the average computational time of the 10 runs in seconds are listed in the CPU AVG column. For comparison, we also list the known traditional TOP optimal solutions in the R BEST column. For each instance, except p2.2.j, R MAX tends to increase until N θ = 8; with increasing N θ , the computational time increases continuously. The main reason is that greater discretization will lead to more accurate heading angle for UAV, which will effectively shorten the path length of the UAV, in turn it will help the UAV to access more targets to obtain greater profit. However, a larger heading angle dispersion increases the computational complexity of the problem and increases the computational time. The results of Numerical Experiment 2 show that under the premise of ensuring a high-quality solution, setting appropriate heading angle discretization can effectively improve the time efficiency of the algorithm.

Simulation Experiments
In this paper, we construct two types of datasets based on the TOP benchmark instances for simulation experiments. The first simulation scene is reconstructed on Set 6 of the TOP benchmark, which is called the even distribution dataset, denoted as E. The second simulation scene, which is reconstructed on Set 2 of the TOP benchmark, is called the clustered distribution dataset, denoted as C. The vertex positions and weights are shown in Figure 8.
in turn it will help the UAV to access more targets to obtain greater profit. However, a larger heading angle dispersion increases the computational complexity of the problem and increases the computational time. The results of Numerical Experiment 2 show that under the premise of ensuring a high-quality solution, setting appropriate heading angle discretization can effectively improve the time efficiency of the algorithm.

Simulation Experiments
In this paper, we construct two types of datasets based on the TOP benchmark instances for simulation experiments. The first simulation scene is reconstructed on Set 6 of the TOP benchmark, which is called the even distribution dataset, denoted as E. The second simulation scene, which is reconstructed on Set 2 of the TOP benchmark, is called the clustered distribution dataset, denoted as C. The vertex positions and weights are shown in Figure 8. For the above two scenarios, K heterogeneous UAVs visit potential targets after launching from different depots whose indices corresponds to the UAV indices, and finally return to their respective departure depots. The number, K, of UAVs used increases from 2 to 4. The relevant parameter values of the UAVs used in the simulation experiment are shown in Table 5. For the above two scenarios, K heterogeneous UAVs visit potential targets after launching from different depots whose indices corresponds to the UAV indices, and finally return to their respective departure depots. The number, K, of UAVs used increases from 2 to 4. The relevant parameter values of the UAVs used in the simulation experiment are shown in Table 5. To demonstrate the influence of the number of UAVs and their depot position on the collected profit, we designed 22 simulation experiments. Each case is named according to the type of dataset, the number of UAVs, and the depot indices. For example, the E2-12 case indicates that the UAVs launch from D1 and D2 and perform a rapid assessment in the E-type dataset. We modeled this problem as a MURARP model and HGSA is used to obtain the routing scheme for rapid assessment. Table 6 records the results for the 22 simulation instances, which was run 10 times under the same experimental conditions. The maximum profit of the 10 runs are listed in the R MAX column, the average profit of the 10 runs are listed in the R AVG column, the minimum profit of the 10 runs are listed in the R MIN column, and the average computational time of the 10 runs in seconds are listed in the CPU AVG column.  Table 6 shows that with additional UAVs, the collected profit increases continuously. Moreover, the profits collected by the same number of UAVs launched from different depots may also be different, which indicates that the selection of the launching depot has a certain impact on the quality of mission completion. The results inspired us to improve the completion quality of the task by optimizing the location of the depot that launches the UAV in the case of limited UAV resources. Table 6 also shows that the stability of the HGSA algorithm is strong. The average gap between R MAX and R AVG is only 5.28% and between R MAX and R MIN is only 11.11%. From the experimental results shown in Table 6, it can also be seen that the average value of the R MAX column is only 5.28% different from the average value of the R AVG , and the average value of the R MAX column is only 11.11% different from the average value of the R MIN column. This demonstrates that the stability of the HGSA algorithm is strong.

Case Study
To illustrate how our method can be applied in a post-earthquake scenario, a case study based on the actual data of the 2017 Jiuzhaigou earthquake in China is conducted. An earthquake with a magnitude of 7.0 struck in Jiuzhaigou County, China at 21:19:46 on 8 August 2017. The epicenter was located at 33.20 • N, 103.82 • E, and the focal depth was 20 km. Soon after, the rescue organization used a UAV to quickly survey the buildings in the disaster area. Data analysis and field surveys based on UAVs were subsequently conducted for the disaster area to assess the geological hazards generated by the earthquake [41]. National Earthquake Response Support Service (NERSS) participated in the post-earthquake survey and rescue operations of the Jiuzhaigou earthquake. Experts of NERSS determined several potential targets, such as schools, hospitals, and other important buildings, that needed UAVs to assess quickly and to capture images and videos. Then, combined with the distance from the epicenter and the population size of these buildings, the possible damage degree of these potential targets was estimated in the form of weights. Then, four rescue teams were sent to carry out rapid assessment of the disaster area from different directions by UAVs. This paper makes a case study from the actual situation of relief operation. Thus, 50 potential targets and the UAV depots are selected for case studies, as shown in Figure 9.  We modeled the above rapid assessment by using the MD-RDTOP-VP model and solving the problem with the HGSA algorithm. The results are shown in Table 8. A multi-UAV routing scheme is found within 10 s using the same parameters as in Section 4.2. In the routing scheme, the UAVs accessed 42 targets, of which 15 were revisited. The expected profit of the routing scheme is 253.540. The relevant parameters values of the UAVs used by the rescue teams are shown in Table 7. We modeled the above rapid assessment by using the MD-RDTOP-VP model and solving the problem with the HGSA algorithm. The results are shown in Table 8. A multi-UAV routing scheme is found within 10 s using the same parameters as in Section 4.2. In the routing scheme, the UAVs accessed 42 targets, of which 15 were revisited. The expected profit of the routing scheme is 253.540. Specifically, Table 8 records the number of visit targets, N VT , the number of revisit targets, N RT , the task execution time, T E , and the utilization rate, R U , of the UAV endurance. From the experimental results shown in Table 8, it can also be seen that the routing scheme generated by our model and algorithm maximizes the endurance utilization rate of each UAV under the premise of target revisit-allowed strategy.

Conclusions and Potential Future Work
In this paper, a novel routing problem for multiple-UAV post-earthquake rapid-assessment is presented. Considering endurance and kinematic constraints, a multi-depot launching heterogeneous UAVs rapid-assessment routing scheme is optimized with the objective of maximize the expected profit. We also specifically design a two-stage hybrid heuristic algorithm, called HGSA, to solve this problem and three sets of experiments are carried out to verity the performance of the HGSA. The experimental results indicate that the HGSA algorithm can obtain a high-quality solution of MURARP stably. For future research, the problem introduced in this paper can be further studied from two aspects. On the one hand, we can study the cooperation between vehicles and UAVs while each rescue team carries different numbers of UAVs and vehicles; on the other hand, the MURARP for heterogeneous targets-such as point targets, line targets, and area targets-can be studied.