Abstract
This paper proposes a new evolutionary planner to determine the trajectories of several Unmanned Aerial Vehicles (UAVs) and the scan direction of their cameras for minimizing the expected detection time of a nondeterministically moving target of uncertain initial location. To achieve this, the planner can reorient the UAVs cameras and modify the UAVs heading, speed, and height with the purpose of making the UAV reach and the camera observe faster the areas with high probability of target presence. Besides, the planner uses a digital elevation model of the search region to capture its influence on the camera likelihood (changing the footprint dimensions and the probability of detection) and to help the operator to construct the initial belief of target presence and target motion model. The planner also lets the operator include intelligence information in the initial target belief and motion model, in order to let him/her model real-world scenarios systematically. All these characteristics let the planner adapt the UAV trajectories and sensor poses to the requirements of minimum time search operations over real-world scenarios, as the results of the paper, obtained over 3 scenarios built with the modeling aid-tools of the planner, show.
1. Introduction
The strong research interest in UAVs trajectory planning takes advantage of the UAVs capabilities to perform different types of military and civil missions, such as georeference [1], wildlife monitoring [2], or target tracking [3]. Moreover, the current developments of their onboard sensorial systems make them also ideal for performing risky long-endurance reconnaissance and surveillance operations. This work focuses on a type of Probabilistic Target Search Problem (PTSP), named minimum time search (MTS) [4], since reducing the time required by the UAVs and their onboard sensors to detect the target is a critical objective of the mission. The selected problem has several applications that include looking for survivors after natural disasters (e.g., after fires or earthquakes), search and rescue operations, or searching for military targets.
Approaches that tackle PTSPs determine the trajectories of the UAVs in spite of the uncertainty associated with the target location and sensor capabilities. To do this, they probabilistically model the different uncertainty sources. In more detail, the information about the target position is modeled with an initial probability map (with the belief of target presence within the search area) and a probabilistic motion model, while the sensors uncertainty is modeled with detection likelihood functions. As an example, Figure 1 schematizes a PTSP where two UAVs equipped with electrooptic sensors (cameras) look for a static lost off-road vehicle in the mountains. The colored map at the bottom shows the altitude of the search area (with green being associated with valleys and brown with mountains), the grey shadowed map in the middle displays the target initial belief (darker areas, associated with lower altitudes, indicate a higher probability of target presence), and the blue/red polygons, respectively, show the current observed areas (within the cameras footprints) by the onboard sensor of each UAV. When the PTSP is also a MTS operation, the regions with higher target probabilities (darker grey areas) should be observed (falling within the cameras footprints) as soon as possible in order to minimize the target detection time.
One of the main goals of MTS planners is to reduce the target detection time, which can be achieved by optimizing the expected time of target detection [5–10]. Other PTSP approaches optimize alternative criteria, such as maximizing the probability of target detection [11–14] or minimizing its counterpart probability of nondetection [15, 16], maximizing the information gain [17], minimizing the system entropy [18], minimizing its uncertainty (areas with intermediate belief of target presence) [19], or optimizing normalized or discounted versions of the previous criteria [4, 20–22]. A common characteristic of the different approaches is that, in scenarios with bounded resources (e.g., limited flying time or fuel), they often obtain better results than predefined search patterns (e.g., spiral, lawnmower), as they adapt the UAV trajectories to the scenario specific target initial belief and motion [6, 20]. Besides, although the approaches that optimize the previous PTSP criteria can share the same elements and probabilistic models, MTS distinctiveness is the extreme influence of the visiting order of the high probability regions in the expected time of target detection, as prioritizing flying over high probability areas first increases the chances of finding the target earlier [4].
The NP-hard complexity of PTSP [23] is tackled with suboptimal algorithms and heuristics, such as gradient-based approaches [13, 15–17, 19], greedy methods [8, 12, 20], cross-entropy optimization [4, 7], Bayesian optimization algorithms [5], ant colony optimization [6, 9], or genetic algorithms [10]. Besides, streamlined formulations of the problem are typically accepted in order to further simplify the problem complexity. They range from considering static targets [8, 10, 13, 15, 16, 19, 20] instead of dynamic ones [4–7, 9, 11, 12, 14, 17, 21, 22] to modeling the sensors ideally [4–6, 8] instead of realistically (e.g., as radars [9–12] or downward-looking cameras [13, 14, 17, 19]) or to assuming that the UAVs fly following straight-lines according to the eight cardinal directions [4–8] or optimized waypoints [17, 21, 22] instead of considering the physical maneuverability constraints induced by the dynamic models of the UAVs [9–16, 19, 20]. Additionally, in some cases (e.g., in [10, 13, 16, 17, 19, 20]) the approach uses a receding horizon controller to divide the UAVs trajectory into sequentially optimized sections, narrowing the optimization search space at the expense of constructing suboptimal myopic solutions. Finally, it is worth noting that the approaches are often tested over synthetic scenarios, built by the authors, without a clear relation with a real-world problem.
The previous simplifications and the lack of analysis over real-world scenarios reduce the applicability of the approaches to real-world problems. This is especially relevant within the subset of MTS approaches [4–10], as the majority of them have been developed for UAVs flying accordingly to the eight cardinal directions [4–8] or only tested over ideal sensors and synthetic unreal scenarios [4–6]. To pave the path of using MTS methods for real-world scenarios, this work extends the capabilities of the planner introduced in [10] by contributing to the following fields:(i)MTS mission definition, by combining intelligence and geographical information to construct the target probability models that will be used in real-world missions. This feature is inspired by software tools that use geographical information to build the target initial belief or motion model [24–26], monitor search missions [26, 27], or evaluate predefined search patterns [24, 28]. Its main benefit, as the results of this paper will show, will be the substitution of the synthetic scenarios used in previous MTS planners by real-world inspired ones.(ii)MTS planning, by optimizing simultaneously the UAV trajectory (by means of changing the UAV heading, speed, and height) and camera pose (azimuth and elevation). This new competence combines the UAV trajectory optimization capabilities of previous MTS planners (which usually manipulate only the UAV heading) with the sensor moving capabilities of only a few PTSP approaches [21, 22]. As this innovation supports a higher moving capability of the camera footprint and a quicker coverage of the high probability areas, it is especially relevant for MTS where the target has to be detected as soon as possible.(iii)Sensor characterization, by incrementing the realism of the camera detection behavior, modeling the effects of the terrain elevation (occlusions and target-camera distance variation), camera orientation, and target and sensor size in the footprint dimensions and sensor likelihood. Not only is the realism of the likelihood model crucial to shorten the differences between the simulated and real behavior of the camera, but it will also be regarded during the mission definition presented in the Results section in order to set up the scenarios correctly.
In short, this paper presents a new planner that optimizes the UAV trajectories and onboard cameras orientations in real-world dynamic MTS scenarios, where the behavior of a dynamic target is probabilistically modeled from a novel perspective by combining intelligence and geographic information, and whose new camera likelihood model takes into account the effects of the terrain elevation in the camera detection capabilities. Besides, we extend the capabilities of the MTS planner presented in [10], which only manipulated the trajectory heading of a set of UAVs equipped with fixed radars in order to optimize the detection time of static targets in synthetically built scenarios, with the simultaneous optimization of UAVs trajectories and sensor poses over realistically built scenarios with dynamic targets. Moreover, the new planner also considers the extension to heterogeneous UAVs, which can have different parametrizations and start and leave the search mission at different times. Finally, this paper describes in detail, and in an algorithmic form, the functionality of the new planner with the purpose of clarifying the interaction between the different elements during the optimization process.
The remaining of this paper is organized as follows. The second section introduces the probabilistic formulation of the MTS problem, describes the novel approaches used to model the target initial state and motion behavior, and presents the new models used for the UAVs and their cameras. The third section details the new multistepped MTS planner presented in this paper, introducing its steps in an algorithmic form and analyzing its computational complexity. The fourth section compares the state-of-the-art of the closest related work with our new planner and highlights its differences with the previous planner in [10]. The fifth section analyzes its performance over three different real search scenarios and shows the benefits of letting the planner reorient the camera. And in the ending section, the conclusions are drawn and some open research questions discussed.
2. Minimum Time Search Definition
This section presents the probabilistic formulation of the MTS problem and describes the approaches used in this work to model the target, UAV, and sensor behavior.
2.1. General Problem Formulation
In the MTS problem presented in this paper, there is a set of UAVs overflying a search area (discretized into a grid of square cells) in order to detect the presence of a single target located within it. Besides, due to the uncertainty associated with the position and dynamics of the target and to the measurements of the onboard sensor of each UAV, the MTS problem is formulated within the probabilistic framework introduced in [5].
In more detail, the information about the target initial position (represented by random variable ) is modeled with the initial target belief , which is a probability mass function that represents the chances of finding the target at each . Besides, when the target is moving, the information about its dynamics is described with the target Markovian motion model , which states the probability that the target moves from any cell to any cell during a time lapse . Moreover, the detection capabilities of the onboard sensor of each UAV are described with the likelihood function , which returns the probability of detecting () or not detecting () a target placed at from the UAV and sensor pose stored in the state variable of the -th UAV.
For those readers familiarized with the Recursive Bayesian Filter (RBF, [12]), the previous probability models , , and are used to obtain the target belief , which represents the chances of finding the target at time step and at each cell , given the trajectory of the UAVs and sensors poses and the sensor measurements . The RBF process, stated in (1), calculates the current belief from the previous time step belief by (1) incorporating the current UAVs location and sensor poses and the last sensor measurements with the likelihood functions , and by (2) displacing and redistributing the target location probability with the target motion model . Besides, (1) assumes that the sensors provide measurements only at time steps that are multiple of and its normalization factor is used to ensure that .
A MTS mission can be formulated as an optimization problem, where some criteria related to the previous probabilistic models and other mission objectives are evaluated in order to determine the best UAV trajectories and sensor poses during the duration of the mission . A useful criterion for MTS is the Expected Time of Detection (ETD, [5–7, 9, 10]), which measures the average time of detection of the target when the UAVs and the sensors follow the trajectories and poses defined in a given . It is calculable with (2)-(4), setting to for the initial case (). The first equation is similar to RBF equation (1), but as it lacks the normalization term and all measurements are nondetection, it obtains the “unnormalized belief” assuming that the onboard sensors do not detect the target from the UAVs trajectories and sensor poses in . The second one obtains , which is the probability of not detecting the target from and whose value decays as the sensors make nondetection observations in regions with probability of target presence. The third expression obtains the ETD from when becomes for some and underestimates its value when . Finally, it is worth noting that minimizing is a better option for MTS than maximizing the probability of detection along the whole trajectory , as the MTS objective is to collect as much probability as possible sooner [4–6].
Besides, in order ensure that the UAV trajectories and sensor poses calculated by our MTS planner are feasible from the maneuverability point of view, we exploit the UAV and sensor deterministic dynamic model of each UAV-sensor pair , where stands for the set of control actions (e.g., commanded UAV heading, height & speed, and sensor elevation & azimuth) and for the time derivative of the UAV location and sensor poses. In particular, our MTS planner uses this dynamic model to obtain the solutions (best UAV trajectories and sensor poses ) from the initial UAVs locations and the sequence of sets of control actions proposed, manipulated, and evaluated by our approach in order to optimize MTS missions. Finally, it is worth highlighting that due to the deterministic nature of the UAV and sensor dynamic model, there is an unambiguous relation between the best sequence of control actions and the best UAVs trajectories and sensor poses.
The realism of the four models (i.e., of the initial target belief , of the target motion model , of the sensor likelihood , and of the UAV and sensor motion model ) is crucial to avoid discordances between the real ETD of a given and the ETD calculated by our approach. Hence, in the rest of this section we present the new models proposed in this paper to bring realism to the definition and solution of MTS missions performed by fixed-wing UAVs equipped with orientable cameras.
2.2. Initial Target Belief Definition
To construct the initial target belief , we merge knowledge coming from different sources, using a different probability layer for each information source and performing with (5) an addition of the probability layers weighted with their reliability/importance coefficients . In other words, considering the first term in (5) is a normalization coefficient that ensures that is a probability function (i.e., ), our initial target belief is calculated as the mixture of the beliefs associated with the different -th information sources.
The layers can be associated with geographical information (e.g., terrain altitude, road maps) or intelligence/user-defined clues (e.g., last/habitual location areas of the target). For the examples of this work, we consider the following two layers:(1)The terrain elevation probability layer , obtained with the following steps. First the digital elevation model (DEM, [29]) of is automatically resampled to the size of the cells in in order to obtain the average height of each cell in . Next the user/operator is required to divide the existing elevations within the cells in into consecutive ranges and to assign a chance of target presence to each range. Finally, the method automatically determines the cells in each elevation range and the probability associated with all , distributing the chances of target presence assigned by the operator. As this way of proceeding generates a geographical probability layer where areas with similar altitude share the same initial belief, it automatically spreads uniformly the belief over different regions of the search area.(2)The intelligence probability layer . For this layer the operator has to graphically define a mixture (weighted addition) of Gaussians (centered in eligible locations of the search area and spread according to selectable standard deviations) and of polygonal areas (defined by their external points placed in the desired locations of the search area ) with uniform probabilities (assigned by the operator). The weights of each element (each Gaussian and/or each polygonal area) within the intelligence probability layer are also selected by the operator according to the information gathered about the last known location of the target.
As an illustrative example, we show the initial belief defined by an operator when looking for a drifting boat next to the cost in Areia Branca, Brazil. To obtain the probability elevation layer shown at the bottom of Figure 2(b), where darker/lighter greys are associated with higher/lower probability cells, the operator analyzes the elevation map represented in Figure 2(a) and assigns high chances of target presence to the elevations by the coast (represented in dark green in the elevation map) to model a strong belief that the boat may have arrived a ground when the search mission starts, lower chances to the sea-level elevation (in blue) to model a moderate belief that the target is still navigating, and zero chances to higher inland elevations. In order to define the intelligence probability layer , represented at the top of Figure 2(b), the operator centers a moderate spread Gaussian in the last known position of the boat. Finally, the operator assigns a weight for each layer ( and ) in order to produce the initial target belief shown in Figure 2(c).
(a) Elevation map (m)
(b) Probability layers ( at bottom & at top)
(c) Initial belief
2.3. Target Motion Model Definition
To define that target motion model, we can distinguish two cases:(1)Scenarios with static targets. In this case, and due to the immobility of the target, if and only if and refer to the same cell in , and otherwise. Besides, in this case, the equations that contain the term can be simplified as .(2)Scenarios with dynamic targets. To construct this type of target motion model, we combine different types of information (e.g., elevation data and sea currents). In this work, the operator has to provide the following:(i)The elevation range where the target is allowed to move within . This allows the operator to indirectly assign the probabilities corresponding to the static target behavior to those cells that do not belong to the . In other words, , the target motion definition process automatically makes if and otherwise.(ii) vectors , each with 9 values that represent, for a few selected cells , the chances of moving to their 8 neighbor cells and of staying at the same cell . Besides, in order to force the target to stay within the search area , some of the 9 values of the vectors that correspond to cells in the borders of are automatically set to zero. With this information, the target motion definition process makes for all and computes applying the potential field method (according to the cardinal directions that make the target move from one cell to its neighbors) for all and . In more detail, if with stands for the adjacent cell to in the -th cardinal direction, represents stay in the same cell, is the distance between cell and cell , and is the -th element of , then and . This way of proceeding allows the operator to define the distribution of the target probability around the neighborhood cells in only a few cells and extend this definition to the remaining moving cells taking into account the distance between the selected cells and the others.(iii)The time period that has to be used to apply to the belief . As the speed of the target is related to the quotient of the size of the cell and , this value is used by the planner to relate the simulation of the motion of the target to the simulation of the motion of the UAVs and sensors.
As an illustrative example, we also obtain the target motion model of the boat drifting example next to the cost in Areia Branca, Brazil. In this case, the operator defines for three cells in the sea, making , for one of the cardinal directions towards the coast, and for the remaining directions . Figure 3 represents with the starting point of the red arrows the location of the selected and with the arrows themselves the direction where for each . Besides, the operator also makes the moving to only let the target (i.e., the drifting boat) move in the sea. To summarize the result of the target motion definition process, the green arrows in Figure 3 represent the direction obtained by weighting, for each cell , the arrows in the 8-cardinal directions with . Besides, the lack of inland green arrows indicates the static behavior of the target in this area. In short, the green arrows show how the defined target motion makes the belief over the sea move towards the coast and remain unchanged inland. Finally, it is worth highlighting that, thanks to our target motion model definition approach, this complex model is automatically built just considering the elevation range () and the target movement probabilities in three cells of .
2.4. Camera Likelihood with Terrain Occlusion
The likelihood of detecting the target with the onboard camera of each UAV in a given cell can be calculated scaling the Target Task Performance (TTP) metric [30] with the percentage of the selected cell within the camera footprint. This scaling, used to reduce the detection probability in the cells of the footprint border, can be modeled with (6), where represents the area (size of the surface) of the common region of the cell and of the footprint of the camera (oriented and placed according to ), is the total area of cell , and is the target task performance function.
To determine the footprint of the camera from the UAV location and camera pose in , we follow the process schematized in Figure 4(a) and consisting of the following two steps. First, we calculate geometrically the camera footprint at sea level taking into account the following pieces of information within : UAV location () and camera azimuth, elevation, and FOV (). To do this, we consider that the sensor location is the same as the UAV location (since its deviation from the location of the UAV mass center is negligible for the mission) and that its orientations are measured with respect to the vehicle coordinate frame (since the camera gimbal compensates the UAV attitude). Second, we approximate the corners of the real camera footprint by the intersections (obtained using the 3D Bresenham algorithm [31]) of the terrain elevation with the 4 lines that join the camera with the 4 corners of the sea-level camera footprint.
(a) Footprint determination
(b) TTPF curve vs.
(c) TTPF curves vs. relative height (for )
The value of the target task performance function TTPF for each cell within the sensor footprint from the UAV location and sensor pose in is calculated with (7), where is the cycle size of the target and (according to [32, 33]) is the critical cycle size for detecting a small target with a likelihood of 50% (since when , ). Moreover, the target cycle size is calculated with (8), where is the real size of the target, is the angle between the along-scan (vertical) and cross-scan (horizontal) directions of the footprint, and (with ) is the ground sample distance (corresponding size of a pixel of the camera over terrain, obtainable as the ratio of the length of the footprint in each direction to its corresponding number of pixels). Taking into account the trigonometric relations between the UAV location and camera footprint, the angle can be obtained solving (9) and the ground sample distance in each scan direction can be calculated with (10) and (11), where is the height of the terrain at the target location and is the number of pixels of the camera in the vertical or horizontal directions ( ).
As an illustrative example of the likelihood model, Figure 4(b) shows how the probability of detection, proportional to TTPF, is incremented as the cycle size of the target grows. Hence, and due to (8)-(11), the probability of detection grows as the target size is increased, as the UAV flying altitude with respect to the terrain elevation gets smaller (through ), as the FOV of the camera is reduced (through and ), or as the camera elevation gets closer to the vertical one (through and ). In more detail, Figure 4(c) shows several curves, corresponding to different and values (indicated in the legend), versus the UAV flying altitude with respect to the terrain elevation (i.e., ) for a target of (e.g., the boat in Areia Branca).
Finally, note that the camera likelihood function defined by (6)-(10) is able to observe partially the borders of the footprints and provide different values for targets of different sizes; cells observed from different UAV heights, camera azimuth, and elevation angles; etc. Hence, our model does not nullify the likelihood that those cells are partially within the footprint but their center falls outside the footprint as in [17, 19], and allows the planner to move the camera from the downward-looking pose assumed in other camera models tested in probabilistic target search problems [13, 14, 17].
2.5. UAV and Sensor Dynamic Models
The UAV dynamics corresponds to a fixed-wing UAV modeled with the upper part of the nonlinear parametrizable Simulink model represented in Figure 5 or the differential equations of the appendix. The UAV motion variables within , highlighted in light green on the right of the figure, are the UAV 3D location , 3D speeds , heading (), course angle (), air and ground velocity ( and ), and fuel consumption (). The UAV motion variables within the command , highlighted in cyan on the left, are the commanded UAV velocity (), heading (), and height (). Additional inputs to the UAV dynamical model are the wind speed () and direction (), highlighted on the left in pink. Besides, in order to let the reader identify the UAV dynamics within the Simulink model, the blocks associated with its height dynamics are colored in blue, with its velocity in grey, with its lateral displacement in white, with the wind in green, and with the fuel in magenta.
The sensor dynamic model is represented at the bottom of Figure 5 and corresponds to a gimbaled camera whose scan direction can be changed commanding its elevation () and azimuth () angles, with the model inputs highlighted in red on the left of the figure. Besides, the camera motion variables within are the outputs of the model highlighted in yellow on the right (camera elevation and azimuth ) and the blocks of the camera dynamics are colored in orange.
The whole model includes the usual limitations related to the UAV speed, height, and heading and to the camera elevation and azimuth. Its different parameters (provided in an external input file, eligible by the operator) allow adapting the movement of the UAV and of the camera to the behavior of different real-world aircraft and camera gimbals. Finally, the model integration, performed from each UAV and camera initial state using the values of a sequence of commands , allows obtaining the with a given simulation resolution time step .
3. MTS Planner
The MTS planner presented in this work is an extension of the Genetic Algorithm (GA) planner introduced in [10]. This section details the main elements of the new planner, introduces its algorithmic description following a top-down strategy, and discusses its computational cost.
3.1. Multistepped GA-Based Planner
The MTS planner obtains the UAVs trajectories for a given mission time following a receding horizon controller approach [10, 13, 16, 17, 19, 20]. More concretely, the planner divides into different sections of duration and optimizes each section ( with ) sequentially, using a GA whose inputs are the final state of the last section and its “unnormalized belief” , and whose outputs are the new and . Besides, it incorporates the possibility of letting each UAV engage and leave the MTS mission at different time instants and (without loss of generality, we assume that for at least one UAV, ) and of using different time steps for the measures of each UAV and for the target update . Hence, to obtain sequentially the sections of the UAV trajectories (and of the cameras orientations), the new planner obtains the overall mission ending time, calculates the number of sections required for the mission, and considers the UAVs and camera immobile and disabled for those time steps where they are not engaged in the MTS mission (i.e., and ). Besides, it also obtains the minimum time step required to update the “unnormalized belief” during the evaluation of the ETD taking into account the time steps of target motion model and of the camera of each UAV .
The algorithmic description of the new multistepped planner is sketched in the pseudocode of Algorithm 1. It starts computing the mission duration , the number of sections into which will be divided, and the time lapse required to update the “unnormalized belief” . Next, it initializes the “unnormalized belief” and the fitness criteria vector , which stores the different values of the criteria used by the GA to decide if a solution (proposal of UAV trajectories and sensor poses in ) is better than another. Finally, it sequentially optimizes each of the sections of the trajectory (i.e., each with and ) using the optimizer in , which will be described in detail in the following section.
| Require: ⊳Structure with scenario models and times: , , , , , | |
| Require: ⊳Structure with optimizer parameters | |
| Require: ⊳Initial UAVs locations and sensor poses | |
| Require: ⊳ Duration of each subsequence | |
| 1: ⊳ Determine the duration of the mission | |
| 2: ⊳ Determine the number of sections required for the multi-stepped planner | |
| 3: ⊳ Round the duration of the mission | |
| 4: ⊳ Determine minimum time step to update . | |
| 5: ⊳ Set starting time | |
| 6: ⊳ Initialize | |
| 7: ⊳ Initialize fitness criteria vector | |
| 8: for r=1:R do⊳ For each subsequence | |
| 9: ⊳ Obtain best subsequence | |
| 10: ⊳ Update time for the following section | |
| 11: end for | |
| 12: return ⊳ Complete UAV trajectory and sensor poses |
3.2. GA Optimizer
This section details the most relevant aspects of the GA used to determine the best UAV trajectory and sensor pose for each section ( with ) of the whole solution . We start presenting the GA genotype, continue with the main GA steps and operators, and finish with the evaluation of each possible solution .
3.2.1. GA Decision Variables
The decision variables directly manipulated by the genetic operators (i.e., the GA genotype) of the MTS planner are the subsequence of control actions used to obtain, by means of the UAV motion model, the -th subsection of the UAV trajectories (i.e., the GA phenotype or solution). Besides, the actions in are the commanded UAV heading , speed , and height and the camera elevation and azimuth . Additionally, the number of decision variables is reduced by applying each action constantly during a time lapse and by adapting the actions required for each UAV within each -th subsection to the time steps where the UAV is engaged in the MTS mission. In other words, the GA algorithm only obtains the values of sets of actions (UAV heading, speed, and height; camera elevation and azimuth) for each UAV engaged in the mission during the -th trajectory subsection.
3.2.2. GA General Description
The GA that sequentially optimizes each of the trajectory sections based on the set of command actions () and fitness criteria (which will be explained in the following section) maintains the operators, parameters, and structure of a standard GA and includes a new step to precalculate the number of actions required by each UAV.
To present it in an algorithmic form, we extend the notation of the paper with left superindexes ( or ) and left subindexes ( and ). The first superindex is used to indicate that there is a population or group of action sequences and fitness criteria vectors . The second one is used to identify the best action sequence , UAV trajectories and sensor poses , and fitness criteria vector within the population of solutions. And the left subindexes are incorporated to the previous notation when the algorithm needs to distinguish between the generic population elements (, ) and those associated with the parents selected by the GA () and with their children (, ).
The GA in Algorithm 2 performs the following steps. First, step 3 computes the number of actions of duration required by each UAV in the current subsection of the trajectory, taking into account the current time step and each UAV starting and ending engaging time ( and ). Step 4 initializes a random population of the J action sequences required for the UAVs engaged in the mission, using a parametrizable uniform distribution. Steps 5-7 simulate the UAV trajectories and sensor poses from the sequence of actions and evaluate them according to the criteria presented in the following section. In step 9, the GA enters an optimization loop which ends when a fixed computation time has passed. Next, step 10 selects parents of action sequences from the current population using binary tournament selection; step 11 creates the children by combining pairs of parents using a single-point crossover; step 12 slightly modifies the children with an incremental Gaussian noise mutation; steps 13-15 simulate and evaluate them; and step 16 selects the survival population from the previous and the mutated children using NSGA-II recombination [34]. Once the computation time has finished, step 19 identifies the best solution in the population; step 20 obtains the , trajectories and sensor poses, and “unnormalized belief” associated with it; and step 21 returns them.
| 1: procedure OPTIMGA(, , ) | |
| 2: Get population size | |
| 3: Obtain number of actions for each UAV | |
| 4: Generate J sequences with the required actions | |
| 5: for j=1:J do | |
| 6: Simulate and evaluate | |
| 7: end for | |
| 8: Initialize counter of the algorithm iterations (generations) | |
| 9: while do Main optimization loop | |
| 10: Select pair of parents (tournament selection) | |
| 11: Perform crossover of parents (single point crossover) | |
| 12: Mutate children (incremental gaussian mutation) | |
| 13: for j=1:J do | |
| 14: Simulate and evaluate | |
| 15: end for | |
| 16: | |
| 17: Increment the GA generation number | |
| 18: end while | |
| 19: Select best solution | |
| 20: Sim. and evaluate it | |
| 21: return Return best solution | |
| 22: end procedure |
It is worth noting that there are several steps whose behavior can be configured using the information in the configuration variable. The initialization parameters in step 4 allow, for each UAV and type of action, (1) selecting different bounds for the uniform distribution and (2) enabling/disabling its optimization. These options let our MTS planner generate trajectories with some fixed behavior (e.g., fixed height or fixed camera orientations). The stop condition parameter in step 9 allows fixing the computation time of the GA algorithm. The crossover parameter in step 11 allows selecting the percentage of parents that undergo the single-point crossover or are directly copied as children. Finally, the mutation parameter in step 12 allows selecting the Gaussian noise that is added to all variables in and the Gaussian noise that is added to a few variables uniformly selected in with a probability , where is the number of decision variables in .
Besides, there are three steps of the GA that use the fitness criteria variable to select solutions of the population to become pairs of parents (step 10), survive in the next generation (step 16), and become the best solution (step 19). Its behavior will be further detailed in the following section, after explaining the optimization criteria used in our MTS planner.
3.2.3. Simulation and Evaluation Process
In order to evaluate one of the solutions, subsequence of action commands , manipulated and proposed by the GA, first, the motion model in Figure 5 or in the appendix has to be simulated for each UAV, and, second, the fitness criteria have to be evaluated. The planner has also to take into account the UAVs mission engaging time and the possible differences among the time lapses of the problem (, , , ). In the remainder of this section we describe the simulation process and the fitness criteria evaluation, and the steps of the function, presented in Algorithm 3 and used by the GA in Algorithm 2.
| 1: procedure SIM&EVAL | |
| 2: | |
| 3: | |
| 4: Simulation of the proposed sequence of actions to obtain UAV trajectory and fuel consumption | |
| 5: for u=1:U do | |
| 6: [, Simulate each UAV and sensor movement | |
| 7: | |
| 8: end for | |
| 9: Evaluation of the trajectory obtained from the proposed sequence of actions | |
| 10: Update collision criterion | |
| 11: Update fuel consumption | |
| 12: Update smooth criterion | |
| 13: Update fuel consumption | |
| 14: Simulation of the target and evaluation of the probabilistic criteria | |
| 15: | |
| 16: while do | |
| 17: if then | |
| 18: Rename variable for the target transition | |
| 19: Target motion part of Equation (2) | |
| 20: end if | |
| 21: Measurement update part of Equation (2) | |
| 22: Equation (3) | |
| 23: Update ET, with equation (4) in sequential form | |
| 24: | |
| 25: end while | |
| 26: Calculate myopia reduction criterion | |
| 27: return , , | |
| 28: end procedure |
The simulation of each UAV is performed in step 6 of Algorithm 3 with the function, taking into account the Simulink or appendix model , the initial UAV location , and the timely-equispaced set of actions in . Moreover, during the simulation each set of actions proposed by the GA is considered constant during the time lapse and the resulting UAV trajectory and sensor poses are created with sampling simulation time of . We use two different time lapses in the simulator (for and ) to reduce the number of decision variables manipulated by the planner, while permitting a high time resolution for evaluating certain aspects of the UAVs trajectories . Finally, it is worth noting that also obtains, from the UAV motion model, the fuel consumption of the UAV during the corresponding trajectory subsection and checks if the UAV is engaged in the mission in order to only simulate its displacement in the simulation time steps that simultaneously belong to and .
The fitness criteria of the MTS planner consists in fulfilling two constraints (UAV collision and nonflying zones avoidance) and in optimizing four indexes (command smoothness, fuel consumption, expected detection time, and myopia reduction criterion), whose values are stored in the different elements of the Fitness Criterium Vector . Besides, in order to distinguish among the problem time periods and within the expressions of the evaluation criteria, we make stand for the time vector with elements equi-spaced by the selected and stand for the time steps within the set that fall within the limits . Moreover, the used in Algorithm 3 within expressions permits the identification of the element of the vector associated with each criterion. Finally, the details of each criterion are presented below:(i)UAV collision avoidance. This criterion counts the number of trajectory time instants at which the distance between the locations of all pairs of UAVs is smaller than the security collision distance . The criterion, evaluated in line 10 of Algorithm and detailed in (12), consists in checking the inter-UAV distances at time steps where both vehicles are engaged in the mission.(ii)Nonflying zones (NFZ) avoidance. This criterion counts the number of trajectory time instants at which the locations of any UAV overfly a selected set of forbidden overflying cells of . Again, the criterion, evaluated in line 11 of Algorithm 3 and detailed in (13), considers locations at the time steps where the UAVs are in the mission.(iii)Actions smoothness. This criterion, evaluated in line 12 of Algorithm 3, measures the smoothness of the commands for each UAV within . It accounts, with (14) and (15), for the consecutive changes of the five types of commands only during the time interval where the UAVs are engaged in the mission. The magnitude of the accounted change, calculated with (16), makes the smooth criterion only account for those changes in action that imply a modification of sign in the signal with a value proportional to the square of the difference. Therefore, the smoothness criterion penalizes big changes of values and sign in the sequence of each type of command.(iv)Fuel consumption. Its value is directly returned by the UAV motion model and affected by the UAV speed.(v)Expected time of detection. The main objective of the planner, described already in the general problem formulation section, is implemented in an iterative form within lines 16 and 25 of Algorithm 3. The process implements (2)-(4) after the following modifications. On the one hand, the expression in line 19 includes the target motion model transition operation for the time steps where the motion model has to be applied (i.e., when is true, because the target is dynamic and the last target motion update was performed seconds ago). On the other hand, the “unnormalized belief” update expression in line 21 only includes the likelihood models of the UAVs at the time steps where the UAV is engaged in the mission and its sensor has to take a measure (according to its ).(vi)Myopia reduction criterion. The planner can create myopic solutions due to its multistepped behavior, which only considers the probability regions achievable from the starting position of each section of the trajectory. In other words, by optimizing only the previous criteria, the planner becomes greedy, ignoring during the optimization of each trajectory section the benefits of getting closer to a high probability area which is currently unreachable. To reduce this myopic behavior, the planner evaluates, in line 26 of Algorithm 3, the suitability of the final UAVs positions and sensor poses of a trajectory section for the optimization of future trajectory sections by weighting the remaining unnormalized belief with a function that takes into account the of each cell in to the center of the camera footprint associated with the final location and camera pose of each UAV . Besides, this evaluation is only performed for the UAVs that are still engaged in the mission. In other words, the planner evaluates the suitability of the final camera footprints of the searching UAVs for optimizing the remaining sections of the trajectories.
As we have previously mentioned, the values obtained after the simulation of each trajectory section are later used in 3 different steps of the GA. The way of proceeding in each step is the following:(i)The binary tournament selection operator chooses randomly, with an uniform probability, two solutions of the population; among them it selects as a parent the one with lower constraint violation (i.e., with lower ), and when tied, the one that Pareto-dominates the other four optimization objectives (i.e., , , , ).(ii)The survival/recombination operator organizes the solutions in a similar way: first according to its constraint violation level and second according to the Pareto-dominance of the remaining four optimization objectives. Hence, as a whole, the GA tries to identify the solution that violates the constraints less times and, within them, those solutions that Pareto-dominate the other objective criteria.(iii)The selection of a unique best solution is required at the final step of the GA to let the multistepped planner use its final position as the starting point of next trajectory section. To do this, the GA selects, among the solutions in the best Pareto-front of the last population obtained by the recombination operator, the one with smallest rounded value of myopia reduction criterion (first), smallest ETD (second, if tied in the first), smallest command smoothness (third, if tied in the previous), and small fuel consumption (fourth). The priority levels used for the selections process reflect the importance order of the objectives. Moreover, the myopia criterion is preferred to ETD, because accounts both for the probability gathered by the trajectory (through ) and for the distance to all the cells with remaining “unnormalized belief”, while only considers the reachable probability within each trajectory section .
3.3. Computational Cost of the MTS Planner
Finally, we calculate the computational complexity of our planner. To do this, we will consider that all UAVs start and end the mission at the same time steps (, ) and that the values of quotients , , , , are natural numbers.
We start studying Algorithm 3 and estimate the computation complexity of the UAV and sensor motion simulation () and of the evaluation of each fitness criterion (). The computational cost of the UAV and sensor motion simulation is proportional to the number of UAVs and to the steps in which the simulator has to return UAV and sensor poses, . The computation cost of the collision avoidance criterion is proportional to the number of possible UAV combinations and to the steps in which the UAV trajectory and sensor poses have been discretized. Hence, . For the NFZ criterion, we can follow a similar process, but accounting only for the number of UAVs. Hence . In the case of the control smoothness, we have to consider the number of UAVs , the number of steps in which the action sequence is discretized, and the number of types of actions that are being optimized . Hence . As the fuel consumption is calculated during the simulation of the UAV and sensor motion, its computational complexity is already included in . To estimate the complexity associated with the ETD computation, we can observe that the cost associated with the target motion is proportional to and , and cost associated with the sensor update is proportional to . The first cost is associated to the operations in the target transition , which can be calculated, a number of times, as the multiplication of two matrixes of and size. The second cost is associated to the update of each sensor , which is calculated a number of times and is proportional to . Hence, . Finally, as the computational cost of the myopia reduction criterion is proportional to and to , its complexity . The addition of the previous complexities, after simplifying the final expression, makes the expression of the computational complexity of the simulation and evaluation function become (17).The expression shows that the complexity of the simulation and evaluation function grows with the increment of the duration of each evaluation, the number of UAVs , the number of decision variables manipulated by the planner , and the number of cells of the grid of the search space . Besides, also grows as the time steps of the simulation , of the actions , of the target , and of the camera measurements are decremented. Hence, in order to have a reasonable computational cost, it is necessary to achieve a compromise among all those parameters. Besides, the quadratic dependency of with (i.e., with the size of the grid of the search space) is worth noting. This quadratic cost is especially relevant, because a moderate size search space (e.g., of 30 x 30 cells in the Areia Branca example) has already significantly bigger values of (e.g., 900) than those used in other parameters (e.g., , , , , ).
To account for the computational complexity of the planner, we still have to consider the computational cost of the GA in Algorithm 2 and of the multistepped approach in Algorithm 1. For the GA, we compute the cost of the main loop of Algorithm 2. Considering a population size of solutions, decision variables, and iterations of the loop, the cost of the , , , and . This happens because computational cost of the tournament selection of each pair of parents is proportional to a constant, the cost of the single-point crossover and of the incremental mutation of a solution is proportional to the number of decision variables within it, and the cost of the recombination step of the population, which includes the sorting of the solutions, is proportional to . Hence, the overall GA complexity can be calculated as . Finally, as the multistepped GA-based planner in Algorithm 1 calls the GA times, its overall computation complexity becomes , calculable with (18), obtained by substituting within it by (17).
The overall computational complexity shows the same dependencies observed for the computational complexity of the simulation and evaluation and a (multiplying) growing behavior proportional to the population size , the number of generations in GA , and the number of subsections in which the trajectory is divided . Besides, as appears in the denominator of that multiplies the numerator of different quotients with the time steps, there is a part of the computational cost that can be considered independent of and of the number of sections in in which the trajectories are divided, and dependent on the quotient of the duration of the mission and the different time steps (, , , and ).
Finally, it is worth noting that the complexity analysis, which is useful to consider the effects of selecting certain variables of the system in the computation cost, does not provide information about (1) the time required for evaluating the solutions or (2) the convergence speed of the multistepped planner and of each GA run, whose behavior is highly dependent on the probability models (initial belief, target motion model, and sensor detection function) of each scenario. To analyze the convergence speed, in the Results section we will analyze the temporal performance of our new planner (implemented in Matlab and using its Parallel Toolbox in the evaluation loops, lines 5-7 and 13-15, of the GA algorithm) over different scenarios to be able to determine, empirically, the time needed to optimize them (in a 2.5 GHz Intel Core i7 with 6 GB RAM PC with Windows 7). To provide some insight related to the time required for evaluating the solutions within the planner, we have measured the values of the different steps of the GA algorithm in the scenarios of the following section and observed that the mean computation time of each population step of the GA (loop in line 9 of Algorithm 2) is 0.46 seconds for the first scenario, 0.49 seconds for the second one, and 1.05 seconds for the third one; that the evaluation of the solutions ( function) takes 99% of the time; and that the UAV motion simulation takes 43% of the time of the function, while the rest of it takes 53%.
4. Related Work
This section compares the new planner presented in this paper with existing approaches and formulations of PTSP. Instead of performing an exhaustive state-of-the-art, we present a compilation of the works that have motivated this paper. In detail, we focus on two types of approaches: methods that optimize the ETD or other criteria where the visiting order of the cells is important [4–10], or more generic PTSP approaches that obtain either smooth trajectories of the UAVs [11–17, 20] or that also optimize the sensor orientation [21, 22]. Hence, we leave out of this comparison other works such as [18, 35–42] that do not share those properties. The comparison is summarized in Table 1, where the first column identifies each work and the remaining ones show the properties highlighted in the comparison, and where the works are grouped row-wise by the type of criteria they optimize (which is shown in the tenth column). In more detail: Target column shows which approaches deal with multiple and/or moving targets. Although including both properties permits handling a wider range of scenarios, it is worth noting that none of the ETD strategies [5–7, 9, 10], at the bottom of the table, can still handle multiple targets. UAVs column shows which planners are developed for multiple UAVs and their type of dynamic model, grouped under the following labels: cell-to-cell indicates that the UAVs move from one cell of the belief to another cell according to a finite set of cardinal directions, waypoint indicates that the UAVs move according to segments defined by the waypoints selected by the planner, simplified indicates that the UAV motion model is linearized or consists in a constant speed steering model, and complex indicates that the planner works with a Simulink or differential motion model as ours. The benefits of multi-UAV approaches with complex UAV motion models are the possibilities of, on the one hand, observing different regions of the space simultaneously and, on the other hand, adapting the UAVs trajectories to the actual maneuvering capabilities of the UAVs. Sensor column shows the type of sensor, its observation range, and its efficiency, as well as the planner capability of moving it. The algorithms under comparison use mainly downward-looking (DL) cameras and radars as well as generic sensors, whose behavior is identifiable by the labels in its observation range (a single cell, a cell of the quadtree used to store the belief, and a limited or wide group of cells) and efficiency (ideal when the likelihood is always either 1 or 0, and FN and FP if it, respectively, considers false negatives and false positives). It is also worth noting that wide range sensors are often modeled with decaying likelihood functions with the distance [9–12, 16], while the remaining sensors are often modeled with constant likelihood functions [4–6, 8, 13, 19, 21]. Besides, although camera sensors are used in several works [13, 14, 17, 21, 22], only a few approaches (including the one presented in this paper) reorient them to observe different regions of the search space easier. Optimization column shows the criteria and the approach used to improve them within each planner. On the one hand, the optimization criteria are used to organize the works row-wise in four groups: optimizing the probability of detection (maximizing or minimizing its counterpart ), maximizing the information gain (IG), considering the trajectory distance or the time (by including a time discounting factor) in previous TPSP criteria, and minimizing the expected time of detection (ETD or the Monte Carlo Simulation Time of Detection MCSTD). Besides, in the table, the notation is extended with a right superindex in and in order to represent the corresponding probability of each () target (each with an initial belief and target motion model independent of the others). And to evaluate MCST for a given , several Monte Carlo simulations are run (each one consisting in sampling the target initial position from the initial belief and in measuring the time step of the trajectory where overflies the sample ) and the average value of the target detection times of all the simulation is calculated. In this regard it is worth noting that optimizing either , , IG, or the uncertainty by itself is not the best option for MTS, where it is necessary to collect the maximum probability of detection as soon as possible. Besides, although including the distance (and hence making the approaches prefer solutions that optimize the PTSP criterion with shorter trajectories) or the time (by including a disconting factor that gives higher importance to the optimization of the criterion in the initial steps) in the timeless objectives alleviates the problem, those solutions [4, 20–22] are biased by the strategies and the parameters used to include them. Hence, they do not optimize ETD as well as the approaches in the last group [5] do. On the other hand, the optimization approach column shows the variability of methods used to handle the high complexity of PTSP, ranging from greedy or local approaches to gradient or graph based techniques, max-sum algorithms, Sequential Quadratic Programming (QSP), Multiarmed Bandit algorithms (MAB), Cross-Entropy Optimization (CEO), Ant Colony Optimization (ACO), or Genetic Algorithms (GA). In this regard it is worth noting that the selected approach is often associated with the underlying type of decision variables (discrete when the UAV dynamic model is cell-to-cell, continuous otherwise), the type and number of criteria functions, and the authors expertise. The main benefit of the selected GA optimizer, inherited from the planner in [10], comes from its capabilities, tested by many researchers over different real-world problems, to handle many-decision multicriteria optimization problems. In this regard, it is worth noting that the selection of a GA has been especially useful to add several constraints and optimization criteria to ETD and to incorporate the five types of decision variables (actions) in the new planner.
If we analyze the table in detail, we can see that this work is the most general within the ETD group, as it allows the planner to obtain the best camera orientations for the best UAV trajectory, deal with moving target scenarios, and optimize multiple criteria. Moreover, the optimization approach in [8], which evaluates the expected time of detection via MCSTD, is directly targeted at initial Gaussian beliefs, while ours is targeted and tested against generic initial beliefs. Regarding the UAV motion model and the quantity of optimization criteria, this work is also more general than the ones classified in the other three groups. Its main drawback, inherited from existing ETD planners, is its single-target support. However, as we will show in one of the examples, this limitation can be shortened by using as initial belief the joint belief of the targets and making the UAVs continue the search mission until all the targets are found.
Besides, a more detailed comparison with its closest MTS state-of-the-art planners ([9, 10]) shows the following facts. First, the new planner optimizes five action types (UAV orientation, speed, and height; sensor azimuth and elevation) while the planners in [9, 10] only optimize the UAV orientation. Moreover, as the good performance of the ACO-based planner in [9] is supported by the use of a heuristic for generating promising UAV headings, its extension is not straightforward as it will require building heuristics for generating promising UAV speeds and heights, and promising camera orientations. Second, the camera introduced in this paper challenges the planner as its limited range likelihood enhances the stepped and multimodal behavior of , while the soft decaying wide range likelihood of the radars in [9, 10] smooths , helping those planners. Finally, although the new planner and the one in [10] share the multistepped GA-optimization support, UAV motion model, constraints, and objective functions, the new planner (and evaluation functions) extends the capabilities of the planner in [10] by allowing obtaining solutions (UAV trajectories plus sensor poses) for scenarios with dynamic targets, using five different types of actions, with different UAV motion model parameterizations, and with different UAV mission engaging and camera-capture times.
An additional contribution of this work to MTS search, not included in the table due to space limitation, is the methodology followed to build realistic scenarios for testing the planner. Moreover, the next section details how to set up the scenario correctly, as there is a strong influence of some parameters (e.g., the size of the target, the camera FOV, and the height and size of the search region) on others (such as the recommended UAV flying altitude, camera likelihood, and search cell size) and, hence, on the planner performance.
Finally, for comparison purposes with existing approaches, we can set up the new planner with the capabilities of the one presented in [10], as it is a direct extension of it and the new planner allows selecting which decision variables have to be optimized. We do not compare against the other works of the ETD group, as the planner in [9] is not capable of optimizing all the criteria used in this work, and the remaining ones [5–8] are designed for UAV motion models that move from cell to cell and for optimizing a single criterion.
5. Results
In this section the new planner is tested over three real-world inspired scenarios. We first explain their common properties and the analysis setup, and afterwards we detail the individual characteristics of each scenario and the solutions obtained by our planner under two different configurations: optimizing the UAVs trajectories and sensor poses, and only optimizing the UAV orientation (as in [10]). Therefore, this second configuration is used as a baseline to compare against and observe the benefits of this work proposal of optimizing multiple types of actions (including the sensor poses).
5.1. Common Scenario and Planner Setup
On the one hand, to study the benefits of optimizing multiple decision variables, we make the GA of our planner optimize: (1) only the commanded UAV heading (as in the previous planner in [10]) or (2) the commanded UAV heading , UAV speed , and camera azimuth . For shortening, hereafter the first configuration will be labeled as and the second as . Besides, in this work, we prefix the values of the remaining decision variables (commanded UAV height and commanded camera elevation ) to facilitate the analysis of the results. Finally, for comparison purposes, in we prefix the UAV speed and camera azimuth (i.e., the two extra decision variables of ) to the mean value of the permitted ranges in .
On the other hand, and due to the stochastic nature of a GA, we run each planner configuration 20 times and represent the evolution (against the computation time) of the mean and variance of the obtained at different runs. Note that we only show the ETD as it is in fact the most relevant objective for this planner. Besides, as is not available until the last section of the trajectory is being optimized, for the computation times where the multistepped planner has not yet started the final section optimization, we approximate the value of with the following expression, which accounts for the ETD up to the section that is being optimized plus the worst case (associated with the case in which the UAVs do not collect more beliefs after ) for the remaining time ().It is worth noting that the approximation in (19) induces the stepped behavior of the ETD curves presented in the paper in Figures 7(d), 9(d), and 10(d). This happens because the contribution of the worst case term in the , which is associated with the belief that can/has not yet been collected by the UAVs, is usually decremented after a new section of the trajectory is being optimized, because as soon as the new probability of nondetection is lower than the one of the previous section , there is at least a decrement of in . Moreover, in the initial sections, where is lower, for a similar reduction in the nondetection probability, the decrement in the represented value of the is often bigger.
The results of the planner for each scenario are, respectively, schematized in Figures 7, 9, and 10. In all of them, the graphic at the left-top (a) represents the initial and UAVs initial locations (colored airplanes), the graphic underneath (d) shows the evolution of the mean and standard deviation, and the second and third columns display additional information of a representative solution obtained by each planner configuration (only heading in the second column, and heading + speed + azimuth in the third one). In more detail, the top row graphics of the second and third column show the , , and the probability of detection of a representative solutions of the corresponding planner configuration, while the graphics at the bottom row also represent the UAVs trajectories (with the red and blue lines), the camera footprints (with the shadowed red and blue polygons), and the camera aiming point without terrain elevation (by the end of the yellow line that starts at the UAV location of each measurement).
Additionally, the elevation information used to construct the initial beliefs for each scenario is available in [43]. Besides, the wind speed in all the scenarios and all the cameras has a resolution of 2160x3840 pixels and a fixed elevation angle of 45°. Moreover, the camera FOV and UAV altitude are precalculated and maintained constant for each mission to ensure a high detection likelihood during the whole search. To achieve this, the FOV and are adjusted to make acceptable the value of for the worst possible case, which occurs at the highest distance between the target and camera and, hence, at the lowest terrain elevation. Additionally, the size of the cells in is calculated for each mission to ensure that the footprint covers more than one cell in the worst case (at the highest terrain elevation). Finally, all missions share the frequency of the camera captures and command changes (s), the simulation step (s), the duration of the trajectory sections (s), and the collision security distance (m).
Lastly, we inherit several values of the planner parameters from [10]. The population size , the probability of crossover is 0.8, the probability of mutation , and the high-lower additive mutation standard deviation is 1.0-0.1. Besides, the myopia parameter .
5.2. Looking for a Static Off-Road Vehicle in the Mountains
The target of the first scenario (proposed by Airbus) is an off-road vehicle broken down in Gador mountains (Spain), within a search area of 10x10km2 whose elevation varies, according to the elevation map in Figure 6(a), from 626 to 2037m. The elevation probability layer, represented at the bottom of Figure 6(b), is obtained by dividing the elevation range into four fractions and assigning lower probabilities to higher elevations. The intelligence layer, displayed at the top of Figure 6(b), consists in two Gaussians centered at the valleys oftener visited by the vehicle. The initial belief, displayed in Figure 6(c), is obtained merging both layers with and , since this weights combination makes the two Gaussians be only 1/16 (i.e., 6.25%) of . Besides, the only UAV engaged in this MTS operation, from to s, starts the mission at the right top corner of the scenario and can fly at (with respect to the sea level) and m/s, with camera azimuth . We select a camera FOV=25° and UAV m to ensure that the sensor likelihood varies, according to the elevations in and target size (m), from 0.90 to 0.99. Besides, as the estimated smallest complete footprint, according to the UAV constant height and maximum terrain elevation, is 590x482m2, we divide into a grid of 30x30 cells, each of 333x333m2.
(a) Elevation (m)
(b) Probability layers
(c) Scenario (, )
(a) Scenario (, )
(b) of
(c) of
(d)
(e) Solution of
(f) Solution of
The results of running both configurations of the planner for this scenario are presented in Figure 7. The graphic of the ETD evolution in Figure 7(d), which represents the approximated up to the optimization of the last section of the trajectory, clearly shows the stepped behavior associated with the inclusion of new trajectory sections in the optimization of the planner. It is worth noting that the duration of the ETD steps can be easily clipped reducing the computational time of the GA stop condition, as each GA run within the multistepped planner converges to its ETD plateau long before its designated time has finished. Besides, the ETD evolution graphic also shows that the produces better results (with lower ETD) during all the optimization processes than .
Furthermore, if we analyze the additional information of the representative solutions of both configurations, presented in Figures 7(b), 7(c), 7(e), and 7(f), we can observe that permits the UAV + camera to collect more than (51% vs. 43%), because the UAV adapts its speed to collect slowly the two Gaussians and arrive quicker to the southeast region of . Besides, the footprints obtained with are more scattered than those obtained by , due to the possibility of changing the camera azimuth. Finally, it is worth noting how the shape of the footprints changes and gets closer to the camera aiming point at zero elevation in the valleys (at the north and at the south west) than in the mountains.
5.3. Looking for Survivors after an Earthquake
In the second scenario, inspired by the earthquake of Haiti in 2010, we search survivors at Carrefour (Port-au-Prince) in a rectangular area of 10x7km2, whose elevation varies, according to the elevation map in Figure 8(a), from sea level to 875m. To use our single-target search planner for a multitarget scenario, we only have to make the initial belief of our scenario the joint-probability of all the survivors and make the UAVs perform the search until all the targets are found [16]. The elevation probability layer, represented at the bottom of Figure 8(b), is defined by dividing the elevation range into three fractions and assigning 0 probability of survivors presence at the sea, high probability at the coast, and lower probability to elevations bigger than 400m. The intelligence layer, represented at the top of Figure 8(b), indicates the probability of survivors based on building damage assessment. This information, available in [44] and used in [45] for another search problem, was obtained by several public and private entities, comparing previous and a posteriori aerial imagery of the earthquake area. The initial belief, represented in Figure 8(c), is obtained merging both layers with and , in order to make the intelligence information represent 1/2 (i.e., 50%) of . The first UAV (in red in Figure 8(c)) arrives at the search area from the East at s, while the second UAV (in blue) does it from the south at s. Both UAVs remain in the mission up to s and can fly at and m/s, with camera azimuth . However, due to the survivors size (m), both UAVs are forced to fly at lower altitudes m with a smaller FOV=20° than in the first scenario, in order to ensure an acceptable camera likelihood that varies from 0.82 to 0.99 (respectively, at the lowest and highest elevations of the search region). Besides, as the estimated smallest complete footprint in this scenario is 383x313m2, the scenario is divided into a grid of 40x28 cells, each of them of 250x250m2.
(a) Elevation (m)
(b) Probability layers
(c) Scenario (, )
(a) Scenario (, )
(b) of
(c) of
(d)
(e) Solution of
(f) Solution of
(a) Scenario (, , NFZ)
(b) of
(c) of
(d)
(e) Solution of
(f) Solution of
The results of running both configurations of the planner for this scenario are presented in Figure 9. The ETD evolution graph in Figure 9(d) shows again that the proposed outperforms regarding the ETD. This happens because both UAVs of can arrive quicker to the high probability areas (in the damaged buildings). It is worth noting also that the probabilities of detection at the end of the mission of two representative solutions of each GA are similar (54% in vs. 55% in ), because both solutions manage to get the same belief during the development of the mission. In more detail, the analysis of and of the trajectories in the second and third column of Figure 9 shows that the solution by collects quicker the concentrated intelligence belief than the elevation belief, while the solution by gets quicker the elevation belief than part of the intelligence belief. However, the similarity in the values of the final probability of detection does not imply that the solutions are equivalent for our planner, as its purpose is to minimize the expected time of detection (which is better for than for ) and not the probability of detection.
Another interesting result of this scenario is the higher variance in the ETD of . This is caused by the uniformity of elevation probability layer, which hardens the scenario for multistepped planners based on generic heuristics as ours. The effect is bigger in the ETD of than in the ETD of , since its capability to manipulate more decision variables allows the planner to obtain more solutions that are similar to the myopic criterium and different from the ETD perspective.
5.4. Looking for a Drifting Boat in the Sea
The target of the last scenario is a drifting boat (m) close to the coast (Areia Branca, Brazil) within a square area of 15x15km2, whose elevation varies from sea level to 147m. The construction of its initial belief and target motion models has already been described in the Minimum Time Search Definition section.
Besides in this scenario, there is a NFZ (e.g., restricted military area) over the coast cells marked in red in Figure 10(a), where the green arrows summarize the tendency of the target motion model. Both UAVs enter the search area from the west at and stay in the mission up to s. Their flying capabilities are different: the first UAV (in red) can do it at and m/s, while the second (in blue) can do it at and m/s. The camera azimuth of both . Besides, both UAVs are required to fly at their highest altitude with a FOV=30° to ensure that the sensor likelihood of the first UAV varies from 0.95 to 0.97 and of the second one from 0.89 to 0.91. Additionally, as the estimated smallest complete footprint is 990x821m2, we divide the search area into a grid of 30x30 cells of 500x500m2. Finally, as the drifting boat is a dynamic target, we also have to select its motion lapse. Taking into account the size of the cells, we make s with the purpose of letting the target advance at approximately 10m/s.
The results of running both configurations of the planner for this scenario are presented in Figure 10. Again, the ETD evolution graph in Figure 10(d) shows that produces overall better ETD results than , at the expense of a bigger variability of its ETD. The initial closer mean values of the ETD results in and are due to the fact that the of this scenario and the initial locations of the UAVs let their cameras initially observe areas with the same probability independently of the UAVs speed and cameras azimuth. The benefit on ETD comes later, when the UAVs with changing velocity can arrive quicker to the high probability areas.
For this case, the of the representative solution of is also higher than its counterpart in (74% vs. 71%). That is, the planner with lets the UAVs collect quicker more probability. Several benefits of changing and are clearly observed in the solutions displayed in the second and third columns of Figure 10. For configuration, not only does the red UAV arrive sooner to the moving Gaussian, but also the blue UAV orients the camera to cover the coast at its right while flying straightforward, while in the red UAV follows the track left by the Gaussian and the blue UAV has to change its orientation several times to be able to collect the belief on the coast. Hence, obtains safer straighter flights, where the UAV headings changes proposed by are substituted by camera reorientations. It is also worth noting that in both solutions the red footprint is smaller than the blue, because the first UAV flies at lower altitude than the second.
6. Conclusions
This paper presents a new planner for real-world scenarios that optimizes the UAV trajectories and poses of its onboard camera for minimum time search missions with uncertain UAV target locations and movements, and noisy camera-based target detection. This new planner extends the capabilities of the one presented in [10] by (1) permitting the change of UAV heading, speed, and height and of camera azimuth and elevation angles; (2) supporting the use of a DEM and of intelligence information of the search area to define the target behavior (initial belief and motion); (3) modeling realistically (considering the target size, sensor specifications, and DEM) the cameras used by the UAVs to search the target; and (4) handling scenarios with moving targets. Besides, we also present an algorithmic description of the planner and analyze the effects of its different elements in its computational complexity. Finally, we (1) show how to model in detail 3 real-world inspired scenarios (looking for an off-road vehicle in the mountains, survivors in an earthquake, and a drifting boat by the coast) taking into account the influence of different variables (e.g., target size, camera FOV, search area and grid size, and UAV height) and (2) analyze the performance and benefits of the new planner capabilities to tackle them.
As future work, we plan to further extend the capabilities of the planner. On the one hand, we would like to let it handle the definition of the probabilistic behavior of multiple targets and the optimization of the expected detection time of all of them. On the other hand, we are considering treating each action type differently, by including weights in the smoothness functions and spending different computational resources in each type (e.g., using a different action step for each type or permitting the multistepped planner to modify each of them in different runs of the GA with a different time computation limit for each type). Besides, we are also exploring the possibility of substituting the GA that supports the optimization of our planner by other metaheuristics, to include new heuristics in the planner for speeding up the optimization, or to generate solutions supported partially by predefined patterns (e.g., spirals, lawnmower) often used in general probabilistic target search problems.
Appendix
In this appendix we sketch the differential equations that model the UAV motion ((A.2)-(A.8)) and sensor motion ((A.2)-(A.8)). In more detail, we consider that the velocity of fuel consumption is proportional to the UAV air speed ; that the horizontal speeds and are related to the UAV air speed and orientation , and to the wind speed and orientation ; that the derivatives of the UAV air acceleration , the UAV angular acceleration , and the UAV height acceleration are related to their time constants () and commanded speed , orientation , and altitude ; and that the camera elevation and azimuth speeds () are related to the commanded camera poses (). Besides, although not stated in the expressions, the UAV air speed , angular speed , and vertical speed and the camera elevation and azimuth speeds () are bounded, and the rate of change of the last three () is limited (see the Simulink model in Figure 5 for further details).
The model parameters ( and ) and the limits in the speeds and rate changes allow adapting the movement of the UAV and of the camera to the behavior of different real-world aircraft and camera gimbals. Finally, this differential model can be simulated, after rewriting the expression in a state-space fashion, using different tools (such as the ODE solvers in Matlab or a Runge-Kutta algorithm).
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by Airbus under SAVIER AER-30459 project.