Multi-UAV path planning methodology for postdisaster building damage surveying

For effective disaster relief decision-making, responders require extensive and rapid information on the damage situation in affected areas. Areas with unknown conditions pose a high risk of injury, and working on the ground limits the coverage and speed of information acquisition. An alternative is to exploit aerial observations and, in particular, unmanned aerial vehicles (UAVs). UAVs can be rapidly deployed to access remote areas without risking survey teams. Moreover, large-scale disasters impact wide areas, and multiple UAVs are needed to increase coverage without compromising resolution or speed. Of particular importance for evaluation are assets such as hospitals, shelters and essential infrastructures. UAVs can survey such structures to construct three-dimensional (3D) models for inspection. A structure-from-motion (SfM) survey generates 3D models from multiple images. However, most path planning algorithms for SfM focus on points of interest taken from an individual UAV and consider a single structure. Here, we propose a path design method for multi-UAV SfM surveys. By designing flight paths with sufficient overlap and sidelap ratios for all faces of the surveyed objects, more precise 3D models can be constructed than with conventional methods. The fuzzy C-means method is adopted to reduce the UAV flight loads to a uniform minimum to ensure full battery utilization.


Introduction
When a disaster such as an earthquake or tsunami occurs, obtaining a comprehensive picture of the damage situation is one of the most important tasks in disaster response. To assess the extensive damage produced by a disaster such as the Great East Japan Earthquake, which occurred on March 11, 2011, responders relied on aerial photographs and satellite remote sensing to avoid delays in information collection caused by destabilization and limitations on ground transportation 18 . Especially when the accident at the Fukushima Daiichi Nuclear Power Plant occurred, satellite remote sensing and unmanned aerial vehicles (UAVs) were crucial platforms to avoid human exposure to radiation on the premises of the power plant 5,20 . Significant progress has been made in the estimation of damage sustained from a disaster using satellite remote sensing 7,19 ; however, the details of damage to buildings typically remain unclear because of the scarcity of information about the condition of building walls. For this reason, increasing attention is being paid to structure-from-motion (SfM) and multiview stereo (MVS) photographic surveys conducted using UAVs 11,21,22 .
Indeed, UAVs are currently applied in several fields, such as aerial assessment of social infrastructure, aerial photography, aerial surveying, and pesticide spraying. In addition, UAVs are used for disaster site surveys because of their low cost and low measurement time requirements. Meanwhile, SfM technology has progressed by virtue of advances in computer vision. Currently, SfM techniques are able to extract many characteristic points from overlapping photographic images and simultaneously estimate the three-dimensional (3D) coordinates of elements and characteristic points of a camera image. Furthermore, MVS is a method for generating a high-precision 3D point group or mesh model via a matching function among multiple photographic images based on the 3D geometric information collected through SfM. SfM and MVS processing is, in general, conducted simultaneously in various commercial software programs; for simplicity, we call this overall process an SfM survey. Software tools such as Pix4D 16 and PhotoScan 1 , among others, have become popular among UAV users. In fact, private companies and individuals can now conduct SfM surveys quite easily. However, in most commercial software, the mission planning feature is limited to a few standard configurations, such as grid missions with a nadir-looking camera (overhead flight) or a fly-around of a point of interest. The limited flexibility of these path planning algorithms necessitates several flights when capturing multiple targets in a complex environment 4 . This issue is critical in the case of disaster response, in which biases and a lack of information in surveys are of particular concern. Moreover, planning multiple flight missions and paths is a time-consuming process that can delay disaster relief activities. A precise 3D model of a complicated structure can be produced using a method implemented in commercial software. This method combines images acquired by a UAV with images acquired from the ground 16 . Alternatively, one can find locations from which to capture images by first using conventional methods for initial model construction and then solving the photogrammetric network design problem (NDP) to calculate the optimal photography locations for building a 3D model 4 . Although these methods enable the construction of quite precise 3D models compared with those produced using conventional methods, they can be time consuming and may require multiple visits to the survey location. To address these shortcomings, previous research has proposed shortening the time needed for surveying by assigning jobs to multiple UAVs using the K-means method 9 . Nevertheless, some difficulties with this approach have been reported. One is a lack of examination of the differences between UAV tasks for a simple environmental structure and those for a complicated and nonsymmetric structure.
This paper proposes a path design method for an SfM survey executed using multiple UAVs in a complex environment. After the simultaneous generation of sufficient photography locations or waypoints for the survey, the degree of attribution of each waypoint with respect to each UAV is calculated using the fuzzy C-means method. Next, we solve a multiple traveling salesman problem (MTSP) to infer the paths that minimize the UAV cost while making the path lengths for all UAVs mutually uniform. This problem is formulated by solving the traveling salesman problem (TSP) for the flight distances and distributing the waypoints in accordance with the degree of attribution (the degree to which a certain waypoint belongs to a particular cluster of waypoints). Finally, we apply the A* algorithm to the obtained UAV paths to improve their obstacle avoidance. This methodology is applied offline since the objects to be avoided are static. For online methods of avoiding moving objects, see 17 . This paper is organized as follows. Section 2 briefly introduces the SfM and MVS methods. Section 3 describes the methodology used for multipath planning. Section 4 examines the differences between conventional methods and the one proposed here. Finally, Section 5 concludes this study.

Background
An SfM survey restores the geometrical shape of an object that has been photographed in multiple images from various camera locations 2, 10, 12-14 . A typical process of treatment by conventional SfM software is depicted in Figure 1. First, pixels containing contrasting density information are extracted from the input images as tie points. Then, matching tie points are associated among multiple images using the estimated camera poses. The point locations are triangulated, and through bundle adjustment, the 3D coordinates are calculated for the locations and angles in each image and for the tie points among images. The camera matrix P P P i i i for m camera locations in 3D space and n points in photographs is described as shown below.
Here, K K K i i i represents the internal matrix of the camera (i.e., the camera settings), and R R R i i i and t t t i i i represent the rotation matrix and translation vector corresponding to the camera orientation and location, respectively. When the 3D coordinate vector of point j ( j = 1, . . . , n) is described as q q q j j j = [X j ,Y j , Z j ] T , the coordinates of camera i (i = 1, . . . , m) in image z z z i j = [u u u i j , v v v i j ] T are given by the following equation: z z z i j 1 ∝ P P P i i i q q q j j j 1 Here, ∝ means that the vectors are proportional. Regarding the parameters of the camera matrix, P P P i i i and K K K i i i are known or partially unknown, but R R R i i i and t t t i i i are completely unknown. A vector p p p i i i can be calculated from the external matrices of 2/14 unknown parameters R R R and t t t by placing them side by side. Then, the unknown parameters p p p 1 1 1 , . . . , p p p m m m of the cameras and the 3D coordinates q q q 1 1 1 , . . . , q q q n n n of the tie points can be estimated based on Equation 2. These are the outlines for bundle adjustment 14 . Because the tie points may include many erroneous points, these outlines are used to delete points that do not satisfy the constraints. The specific treatment and estimation of the parameters differ between software programs. Then, MVS processing is applied to generate dense points. Although many different MVS methods have been proposed, they all fundamentally proceed as presented below.
1. Calculation of the matching costs (similarities) of the tie points.
2. Aggregation of the costs in surrounding domains.
3. Calculation and optimization of parallax (i.e., displacement in the apparent positions of objects).

Improvement of parallax (filtering).
Finally, the point cloud generated in the MVS step is used to create a triangular mesh, and textures are pasted onto the generated faces of the mesh. Examples of SfM-MVS results are shown in section 4.

Methodology
In conventional studies, a low-precision initial model of the objects in the target area is first generated using a traditional path planning method such as the overhead pattern approach. Based on this model, a detailed plan is then calculated to generate the waypoints or photography points 4, 6 that will result in a suitable camera constellation for sufficient coverage.
In this study, first, we built several virtual environments representative of urban areas using NetLogo3D 24 (see Figure 2). The lot sizes and building heights at the sites (expressed in meters) were set randomly, and the structures to be surveyed were visualized in darker colors. Based on these environments, photography points and flight paths were to be calculated. Note that the environment depicted on the left in this figure is later used as site A (see section 4).

Generation of waypoints
First, waypoints separated every 2.5 m horizontally and 5 m vertically were generated around the objects of interest. The camera was presumed to face the center of the nearest object. In addition, waypoints satisfying any of the following conditions were deleted from the candidate waypoint list.
• The distance to the object of interest is less than a certain user-defined value (7.5 m in our case).
• A structure exists between the object and the camera.
• The coordinates of the waypoint are inside of a structure.
We selected a camera visibility angle of 90 • horizontally; additionally, the camera aspect ratio was set to 4:3, and the distance between the camera and the object was set to 7.5 m. Moreover, the images taken at each waypoint were required to have a forward overlap ratio of 83% and a side overlap ratio of 56% to ensure sufficient multiplicity.

Sorting of waypoints
Sorting is among the most studied problems in machine learning. Various sorting methods are available with different objectives. Some typical methods include support vector machines (SVMs) or deep learning models, which require training data, and the K-means method, which uses no training data. In the examination undertaken for this study, to sort the waypoints for photographs based on distance in 3D space, the authors used the K-means method based on the distance from each cluster's center of gravity and the fuzzy C-means method, which is an application of the K-means method. Since the generated waypoints are shared among all UAVs, a sorting method for allocating specific waypoints to specific UAVs is necessary. Here, the fuzzy C-means method was applied to sort the waypoints among the available UAV units -note that all UAVs were assumed to start from the same initial point.
The basic concepts of the K-means and fuzzy C-means methods are described below.

K-means method
The K-means method is an algorithm for nonhierarchical clustering. Clustering is the process of grouping a series of data. As a kind of unsupervised learning method, nonhierarchical clustering is a method of directly classifying similar data into groups using an evaluation function. The K-means method functions on the basis of a simple algorithm; therefore, it is a typical method of clustering that is widely used. Generally, the algorithm is implemented as follows 15 : 1. Randomly allocate all data x x x i i i (i = 1, . . . , n) to clusters.
2. Find the center-of-gravity vector v v v j j j ( j = 1, . . . , k) for each cluster.
3. Find the distances between each x x x i i i and each v v v i i i and re-sort the data into their closest clusters.
4. When the change in v v v j j j is less than a predetermined threshold value, terminate the operation. Otherwise, return to step 2.
In mathematical terms, the process described above is applied to minimize the value given by Equation 3: The performance of this method is known to strongly depend on the positions of the initial clusters. Therefore, the first application of this method might not yield the best sorting results. In addition, it is assumed that the clusters are spherical and that the number of clusters into which the data should be classified is known, and classification is difficult if these assumptions are violated.

Fuzzy C-means method
When each datum x x x i i i is sorted into just one cluster, such clustering is referred to as hard clustering. By contrast, when data are partially sorted into multiple clusters, such ambiguous clustering is known as fuzzy clustering. The K-means method belongs to the former type of clustering; the fuzzy C-means method belongs to the latter. The fuzzy C-means method provides each x x x i i i with an attribution degree (from 0 to 1) with respect to each cluster based on the distance to the center-of-gravity vector v v v i i i . This method is implemented as described below 23 .
1. Set s = 0 and determine the initial value U U U (0) of U U U appropriately.

Determine the centroid vector
4. Determine the attribution degree matrix u u u i j of each x x x i i i with respect to each cluster j using Equation 5.
4/14 5. When the change in v v v j j j is less than a predetermined threshold value, terminate the operation. Otherwise, return to step 2.
Here, U U U is the set of attribution degree matrices u u u i j , d is the Euclidean norm of a datum and a cluster, and m is a constant that represents the extent of fuzziness of the attribution degrees. The value of m is usually set to a real number between 1.3 and 2. If m = 1, then this method reduces to the K-means method of hard clustering. In mathematical terms, this procedure minimizes the value given by Equation 6 below:

Route optimization
For the multiple traveling salesman problem (MTSP) solved in this study, each UAV visits all waypoints in its cluster only once and then returns to the starting point. For all UAVs, the starting point and the end point after the flight were assumed to be the same based on practical considerations. The MTSP for the UAV routes was solved to select one solution for which the cost was minimal and uniform across all devices. To this end, the procedure presented below was applied to ensure that the costs for the UAVs were minimal and uniform, referring to the attribution degree matrices. a brief explanation of other methods supporting our methodology are described in the Supplementary Material. 7. From among the waypoints of UAV max , select the one with the maximum attribution degree with respect to UAV min and re-sort it to UAV min .
Here, tick is the number of trials in a calculation run, f i is the cost function of UAV i , and θ is a constant (threshold value). Note that f i is either the total flight distance, as given by Equation 7, or the total flight time, as given by Equation 8.
Each UAV is presumed to accelerate to its maximum speed with uniform acceleration and to hover at the waypoint for a predetermined duration. Additionally, d k denotes the distance between waypoints k and k + 1, t k denotes the flight time between waypoints k and k + 1, t h denotes the hovering time at a photography point, n is the total number of waypoints, v max is the maximum speed of a UAV, and a represents the uniform acceleration of a UAV.

Obstacle avoidance
The routes found through the procedure presented in the preceding section may pass through obstacles; thus, the A * algorithm was implemented for obstacle avoidance. We selected the A * algorithm for finding a path between two points because it is more efficient in terms of memory than the Dijkstra method. The process of generating waypoints such that obstacles are avoided is as follows: 1. Let the current waypoint be W 0 .
3. i f there is an obstacle between W 0 and W 1 , proceed to the next step; else, go to step 5.
4. Apply the A * algorithm to find the shortest route to W 1 via neighboring nonobstacle points.
6. Repeat from the start until the last waypoint on the path is reached.
Through the flow described above, the SfM survey routes for multiple UAVs were generated. Figure 3 presents the output UAV flight routes for city site A. Each point represents a waypoint, and the color corresponds to the particular UAV to which that waypoint is assigned. The waypoints are linked by lines in accordance with the flight progression.

Examination of the TSP
To examine the effectiveness of the proposed method using the NN, 2-opt, and 1.5-opt algorithms, the precision of the path solutions is compared with that achieved using the genetic algorithm (GA) proposed by 8 . Site A, introduced in section 3, was used for the comparative examination. The GA parameters used here are presented in Table 1. Because both methods use random numbers, to ensure a fair comparison, 100 trials of finding solutions were conducted using both methods. The cost function given in Equation 7 was used here.
The results are presented in Table 2. Figure 4 illustrates the convergence of the cost function for both methods. Figure 4-left shows the case with the maximum converged value of the cost function among the 100 trials for each method, whereas Figure  4-right shows the case with the minimum converged value. Table 2 shows that the combination of the NN, 2-opt, and 1.5-opt methods, as explained in section 3.3, produced more precise results. Basically, the time for one tick is shorter in the 2-opt and 1.5-opt methods than in the GA method. Ultimately, the method proposed in this study is approximately five times faster than the GA method. Figure 4-right demonstrates that even the initial solution provided by the NN method is more precise than the final solution of the GA method. Moreover, with the 2-opt and 1.5-opt methods, convergence was attained in 20 ticks, whereas the GA method took 600 ticks to reach convergence. Therefore, the combination of the NN, 2-opt, and 1.5-opt methods is superior to the GA method in terms of both precision and convergence speed. The reason might be that the combination of the NN, 2-opt, and 1.5-opt methods may be particularly well suited to the current problem under examination, in which the waypoints are distributed regularly but with some deviations. A method in which the initial solution is obtained via the NN method but the subsequent improvement is achieved using the GA method requires more time for convergence than the current method. Improving the initial solution using the 2-opt and 1.5-opt methods might be more compatible with the NN method than improvement via the GA method.

Comparative examination of sorting methods
To confirm the benefits of the fuzzy C-means sorting method used in the proposed method, it was compared to the K-means method 8 . Site A, introduced in section 3, was used for the comparative experiment: five UAVs were used to survey 10 randomly selected structures at the site. After waypoints were generated for each structure, they were sorted using the K-means and fuzzy C-means methods to make the UAV costs minimal and uniform via the procedure explained in section 3.3. In K-means sorting, the waypoints did not have attribution degree matrices. Therefore, the cost of the TSP was instead minimized for each UAV. After 100 trials of the above-described procedures, the results were examined for comparison. The cost functions provided in both Equations 7 and 8 were considered. The parameters were as follows: the constant m in the fuzzy C-means method was 2, the UAV hovering time t h was 2 s, the maximum UAV speed v max was 5 m/s, and the uniform UAV acceleration a was 2 m/s 2 . The total cost F was given as shown in Equation 10. The subtour averaging coefficient ρ was given as shown in Equation 11 for another evaluation function of the sorting methods 23 . Figure 5 shows the cost functions for each UAV in the fuzzy C-means method. Figure 5-left shows the cost function calculated using Equation 7, whereas Figure 5-right shows the results with Equation 8. These figures suggest that the redistribution of the waypoints among the UAVs based on the attribution degree function of the fuzzy C-means method drives the cost functions for all UAVs to approach the average. Tables 3 and 4 present the average values of the evaluation functions for the cost functions given by Equations 7 and 8, respectively. These tables show that the fuzzy C-means method yielded a smaller subtour averaging coefficient and more uniform solutions, whereas the K-means method found a slightly lower total cost. One reason for this result is that with the hard clustering of the K-means method, long paths between structures might not be generated. Another is that f min increased more rapidly than f max decreased with redistribution based on the attribution degree, as shown in Figure 5. However, the max f i results in Table 4 show that the K-means method resulted in a cost of approximately 20 min on average, whereas redistribution with the fuzzy C-means method reduced the cost to less than 15 min. The current capacity limit of UAV batteries is typically 10-30 min, depending on the equipment. The termination of UAV operation is determined by the UAV with the maximum cost. Therefore, it is better to have a smaller maximum cost. From the discussions presented above, one can conclude that the redistribution of the waypoints in accordance with the fuzzy C-means method is quite effective in making the costs for multiple UAVs used to execute an SfM survey minimal and uniform when solving the MTSP.

Comparison of 3D model outputs between the proposed and conventional methods
The proposed method of waypoint generation was examined in terms of its ability to generate more precise 3D models than those produced by means of conventional overhead flight. As the target problem for the comparative examination, city site B with an area of 85 m 2 , as shown in Figure 6, was prepared. The same site was also generated using the Unity3D game engine 3 ,  figure, where the red oval indicates the object to be surveyed. A flight route for one UAV was designed using the proposed method, including the coordinates of the waypoints, the camera orientation, and the order of the waypoints to be visited. Then, this flight route was used in Unity3D to photograph the object. An additional survey of the same object using an overhead flight path was also conducted. In contrast to the flight path designed using the proposed method, the altitude of the overhead flight needed to be sufficiently high to achieve sufficient overlap and sidelap ratios. The UAV took pictures at every waypoint on the flight path. The camera parameters are presented in Table 5. Then, after the SfM survey process was conducted using Pix4d to generate 3D models for both flights, their respective precisions were compared. Table 6 presents the parameters applied in Pix4d. For comparative examination, graffiti and damage details are illustrated in Figures 7 and 8 to allow the precision of their reproduction to be examined. For graffiti drawn on the sunny side of the building with readily available tie points, both methods produced quite similar results. However, for wall objects placed in sunny but obscured areas, the qualities of the images were quite different between the two methods. Even on a shady side where only a few tie points could be extracted, the proposed method precisely reproduced cracks and graffiti. For damage features such as those shown in Figure 8, which are in obscured areas and cannot be adequately interpreted using conventional methods, the proposed method provided clear images of cracks. However, for the buildings shown in Figure 7-left, extra meshes were generated at the building edges. Possible reasons for this behavior are given below.
• Errors in estimating the camera position and orientation.
• Errors in MVS processing.
• False recognition of the boundary between a roof and the sky as a tie point in an image acquired while looking toward the border.

Discussion and Future Work
Despite the promising results, some room for improvement remains for the proposed methods of generating waypoints and redistributing waypoints using the attribution degree. Regarding the former, the authors have proposed a rule-based waypoint generation method with sufficient overlap and sidelap ratios to achieve speedy generation of waypoints. As a result of this method, erroneous extraction of tie points can occur. Additionally, tendencies toward very poor estimation of the camera position and orientation for short distances and multiple photographs are apparent. Therefore, the optimal photography points for the target objects must be generated after solving the NDP problem 4 . Regarding the redistribution of waypoints, the results show that the UAV costs are more uniform than those achieved by sorting using the K-means method. However, the K-means method is superior in terms of the total and average costs. In the examination performed above, waypoint redistribution was performed for five UAVs. If waypoints are to be sorted among more UAVs, groups of clusters that are not mutually close will tend to appear, and long paths might be generated through redistribution using the same cost function and degree of attribution settings used in the present study. To avoid this difficulty, one might add procedures for renewing the center-of-gravity vectors v v v of the clusters in the algorithm used for the MTSP. By a means such as providing a weight coefficient w w w i i i for the Euclidean norm of the K-means method depending on the cost calculated for the MTSP, one could implement such additional procedures rather easily. In this way, the method proposed herein could be more usefully extended to more than five UAVs, allowing it to yield better results in terms of the total and average costs compared to the method of solving the TSP for each UAV for K-means sorting. Another shortcoming of the present proposal is the practical problem that it cannot respond quickly to a cost greater than the prescribed cost. This shortcoming derives from the fact that a long computation time is necessary to achieve convergence of the UAV costs. Only in the final stage of processing can it be known whether the cost exceeds the prescribed cost. To overcome this difficulty, one might predetermine the number of UAVs corresponding to the generated number of waypoints. However, if the converged cost value exceeds the prescribed value, then the calculations must still be repeated for an increased number of UAVs. The next research goal of the authors in this context is to combine the algorithms for solving the MTSP and clustering the waypoints.
A method that can uniformly minimize the costs for several tens of UAVs would be useful for developing methods of surveying much wider areas. The ability to provide highly precise 3D information about disasters over a wide area is expected to contribute greatly to the preparation of countermeasures to mitigate the effects of future disasters.

Concluding Remarks
We have proposed a path planning method for multiple UAVs to aid in 3D building damage surveys or disaster situations. The proposed methodology combines the fuzzy C-means method for assigning waypoints to each UAV with a route optimization algorithm for calculating the visit order of the waypoints for each UAV by solving the multiple traveling salesman problem. Finally, the selected paths are corrected to avoid obstacles between waypoints using the A* algorithm. To assess the effectiveness of the proposed method, a comparative study of the precision of SfM surveys based on waypoints was conducted, and another comparative study was performed concerning solving the TSP to optimize the routes between the photographed locations. In addition, a comparative examination was conducted of the results obtained for the MTSP to minimize and unify the costs for multiple UAVs. In all respects, the proposed method produced more effective results than conventional methods.